Statistics & Algorithms

jsmetrics contains metrics of three types:

1. jet statistics – methods for isolating individual quantities synonymous with the jet stream from upper-level wind speed (e.g. latitude, speed, width).

2. jet core algorithms – methods that return a mask of coordinates related to the jet location throughout the horizontal and/or vertical plane.

3. waviness metrics – methods for determining the “waviness” of upper-level mean flow.

For specification details of each metric’s DOI and original study area please see all metrics.

For progress on their completion see issues.

Jet statistics

Statistics for isolating individual quantities synonymous with the jet stream from upper-level wind speed within a given time window (e.g. latitude, speed, width).

Overview table:

Metric	Status		Metric	Status
Archer & Caldiera 2008	Complete		Woollings et al. 2010	Complete
Barnes & Polvani 2013	Complete		Screen & Simmonds 2013	In progress*
Grise & Polvani 2014	Complete		Barnes & Polvani 2015	Complete
Barnes & Simpson 2017	Complete		Chenoli et al. 2017	In progress
Molnos et al. 2017	In progress*		Bracegirdle et al. 2018	Complete
Ceppi et al. 2018	Complete		Zappa et al. 2018	Complete
Rikus 2018	In progress		Kerr et al. 2020	To verify

* == help needed

Documentation:

Statistics for isolating individual quantities synonymous with the jet stream from upper-level wind speed within a given time window (e.g. latitude, speed, width).

The following statistics each return a xarray.Dataset and are ordered by paper publish year.

jsmetrics.metrics.jet_statistics.sort_xarray_data_coords(coords, ascending=True)

Sort the coordinates of the input data (e.g. lat or lon). Sort in ascending order by default. Used in all jet statistics and jet core algorithms included in this package.

Parameters:

coordsarray-like

Coords to sort in xarray data

ascendingBoolean (default=True)

Whether the coords should be sorted in ascending order

jsmetrics.metrics.jet_statistics.archer_caldeira_2008(data)

This method calculates three monthly-averaged jet stream properties (windspeed, pressure and latitude) via integrated quantities from u-component wind.

This method returns three properties:

1. weighted-average wind speed – jet stream wind speed (\(WS\)), calculated by:

\[WS_{i, j} = \frac{\sum\limits_{k=400hPa}^{k=100hPa} m_{k} \times \sqrt{u^{2}_{i, j, k} + v^{2}_{i, j, k}}}{\sum\limits_{k=400hPa}^{k=100hPa} m_{k}}\]

where \(u_{i,j,k}\) and \(v_{i,j,k}\) are the monthly-average horizontal wind components at grid point (i,j,k), and \(m_{k}\) is the mass at level k.

2. mass flux weighted pressure – the average pressure of flows by the tropopause (\(P\)), calculated by:

\[P_{i, j} = \frac{\sum\limits_{k=400hPa}^{k=100hPa} \left(m_{k} \times \sqrt{u^{2}_{i, j, k} + v^{2}_{i, j, k}}\right) \times p_k}{\sum\limits_{k=400hPa}^{k=100hPa} m_{k} \times \sqrt{u^{2}_{i, j, k} + v^{2}_{i, j, k}}}\]

where \(p_k\) is the pressure at level \(k\).

3. mass flux weighted latitude – Latitude of the Northern Hemisphere jet (\(L^{NH}\)), calculated by:

\[L_{i}^{NH} = \frac{\sum\limits_{j=15°N}^{j=70°N} \left[\sum\limits_{k=400hPa}^{k=100hPa} \left(m_{k} \times \sqrt{u^{2}_{i, j, k} + v^{2}_{i, j, k}}\right) \right] \times \phi_{i,j}}{\sum\limits_{j=15N}^{j=70N} \sum\limits_{k=400hPa}^{k=100hPa} m_{k} \times \sqrt{u^{2}_{i, j, k} + v^{2}_{i, j, k}}}\]

where \(\phi_{i,j}\) is the grid cell latitude.

This method was originally introduce in Archer & Caldiera (2008) (https://doi.org/10.1029/2008GL033614) and is described in Section 3 of that study.

Note: this method does not explicitly limit inputted wind to 100-400 hPa, see ‘Notes’ for more information about the implementation of this method to this package.

Parameters:

dataxarray.Dataset

Data which should containing the variables: ‘ua’ and ‘va’, and the coordinates: ‘lon’, ‘lat’, ‘plev’ and ‘time’.

Returns:

outputxarray.Dataset

Data containing the three jet properties: ‘mass_weighted_average_ws’, ‘mass_flux_weighted_pressure’ and ‘mass_flux_weighted_latitude’

Notes

While the initial methodology provides limits for pressure level (100-400 hPa), here the mass weighted outputs will be calculated for any pressure levels passed into the method.

The latitude calculation is limited to 15-70N (as we only provide a way to extract the Northern Hemisphere jet), but you may find it easy enough to edit this method to calculate outputs for a different region.

Examples

import jsmetrics
import xarray as xr

# Load in dataset with u and v components:
uv_data = xr.open_dataset('path_to_uv_data')

# Subset dataset to range used in original methodology for the NH jet (100-400 hPa & 15-70 N)):
uv_sub = uv_data.sel(plev=slice(100, 400), lat=slice(15, 70))

# Run statistic:
archer_outputs = jsmetrics.jet_statistics.archer_caldiera_2008(uv_sub)

# Subset mass weighted wind by a windspeed threshold
windspeed_threshold = 15 # remember this is for monthly averages
strong_jet = archer_outputs['mass_weighted_average_ws'].where(lambda row: row > windspeed_threshold)

jsmetrics.metrics.jet_statistics.woollings_et_al_2010(data, filter_freq=10, window_size=61)

This method follows an in-text description of 4-steps describing the algorithm of jet-stream identification from Woollings et al. (2010).

This method returns four outputs:

1. jet_lat – latitude of maximum speed within low-pass filtered zonally averaged wind profile

2. jet_speed – speed at the ‘jet_lat’

3. ff_jet_lat – Fourier-filtered ‘jet_lat’ by season

4. ff_jet_speed – Fourier-filtered ‘jet_speed’ by season

This method was first introduce in Woollings et al (2010) (http://dx.doi.org/10.1002/qj.625) and is described in section 2 of that study.

Please see ‘Notes’ below for any additional information about the implementation of this method to this package including Step 5 of the methodology.

Parameters:

dataxarray.Dataset

Data which should containing the variables: ‘ua’, and the coordinates: ‘lon’, ‘lat’, ‘plev’ and ‘time’.

filter_freqint

number of days in filter (default=10 days)

window_sizeint

number of days in window for Lancoz filter (default=61 days)

Returns:

outputxarray.Dataset

Data containing the four output variables: ‘jet_lat’, ‘jet_speed’, ‘ff_jet_lat’ and ‘ff_jet_speed’

Notes

In the original paper, a further step (Step 5) is carried out to express the values of jet latitude and jet speed anomalies from the seasonal cycle, this is shown in the ‘Examples’

Examples

import jsmetrics
import xarray as xr

# Load in dataset with u component wind:
ua_data = xr.open_dataset('path_to_u_data')

# Subset dataset to range used in original methodology (700-925 hPa & 20-70 N, 300-360 W)):
ua_sub = ua.sel(plev=slice(700, 925), lat=slice(20, 70), lon=slice(300, 360))

# Run statistic with a filter frequency and window size used in the original methodology:
w10 = jsmetrics.jet_statistics.woollings_et_al_2010(ua_sub, filter_freq=10, window_size=61)

# Express jet latitude and speed as anomalies from smoothed seasonal cycle (Step 5 of methodology)
w10_seasonal_anomalies = w10.groupby('time.season').apply(lambda row: row['jet_lat'] - row['ff_jet_lat'])

jsmetrics.metrics.jet_statistics.barnes_polvani_2013(data, filter_freq=10, window_size=41)

This method constructs the ‘eddy-driven jet’ by performing a pressure-weighted average of zonal winds. The winds are then low-pass frequency filtered at each grid point using a 10-day Lanczos filter with 41 weights by default. Finally a 0.01 degree quadratic function is fitted to the peak of the subsequent wind speed profile for each time step.

This method returns three outputs:

1. jet_lat – latitude of maximum speed at an interval of 0.01 degree within the subseqeunt profile

2. jet_speed – speed at the ‘jet_lat’

3. jet_width – full width at half of the maximum ‘jet_speed’ within the 0.01 quadratic fitted to the peak of the wind speed profile

Note: ‘jet_width’ is undefined where the ‘jet_speed’ never drops below half maximum.

This method was originally introduce in Barnes & Polvani (2013) https://doi.org/10.1175/JCLI-D-12-00536.1 and is described in Section 2b and 2c of that study.

Please see ‘Notes’ below for any additional information about the implementation of this method to this package.

Parameters:

dataxarray.Dataset

Data which should containing the variables: ‘ua’, and the coordinates: ‘lon’, ‘lat’, ‘plev’ and ‘time’.

filter_freqint

number of days in filter (default=10 timeunits)

window_sizeint

number of days in window for low-pass Lancoz filter (default=41 timeunits)

Returns:

outputxarray.Dataset

Data containing the three output variables: ‘jet_lat’, ‘jet_speed’, ‘jet_width’

Notes

Whereas the original analysis using this method focuses on three distinct sections of the globe, the method here does not make any distinction. Instead we highlight how to do subset the input data and calculate this metric for those three sections in ‘Examples.

This method was based on the method from Woollings et al. (2010) (http://dx.doi.org/10.1002/qj.625)

Examples

import jsmetrics
import xarray as xr

# Load in dataset with u component wind:
ua_data = xr.open_dataset('path_to_u_data')

# Subset dataset to three sectors of the globe used in original methodology:
ua_sh = ua.sel(plev=slice(700, 850), lat=slice(-90, 0)) # the Southern Hemisphere
ua_na = ua.sel(plev=slice(700, 850), lat=slice(0, 90), lon=slice(300, 360) # the North Atlantic
ua_np = ua.sel(plev=slice(700, 850), lat=slice(0, 90), lon=slice(135, 235) # the North Pacific

# Run statistic with a filter frequency and window size used in the original methodology:
bp13_sh = jsmetrics.jet_statistics.barnes_polvani_2013(ua_sh, filter_freq=10, window_size=41)
bp13_na = jsmetrics.jet_statistics.barnes_polvani_2013(ua_na, filter_freq=10, window_size=41)
bp13_np = jsmetrics.jet_statistics.barnes_polvani_2013(ua_np, filter_freq=10, window_size=41)

jsmetrics.metrics.jet_statistics.grise_polvani_2014(data)

This method calculates the latitude of the midlatitude eddy-driven jet (‘jet_lat’) by finding the peak value of the input u-wind field. A polynomial fit is then applied to get an appropriate value of ‘jet_lat’ at a resolution 0.01 degrees. As opposed to the original method, this implementation also returns the speed at the ‘jet_lat’: the ‘jet_speed’

This method was originally introduce in Grise & Polvani (2014) https://doi.org/10.1002/2013GL058466 and is described in Section 2 of that study.

Please see ‘Notes’ below for any additional information about the implementation of this method to this package.

Parameters:

dataxarray.Dataset

Data which should containing the variables: ‘ua’, and the coordinates: ‘lon’, ‘lat’, ‘plev’ and ‘time’.

Returns:

outputxarray.Dataset

Data containing the two outputs: ‘jet_lat’ and ‘jet_speed’

Notes

This method was originally developed for the jet streams in the Southern Hemisphere.

The original paper also includes two other metrics for zonal mean atmospheric circulation.

Examples

import jsmetrics
import xarray as xr

# Load in dataset with u component wind:
ua_data = xr.open_dataset('path_to_u_data')

# Subset dataset to range used in original methodology (850 hPa &  -65--30N)):
ua_sub = ua.sel(plev=850, lat=slice(-65, -30))

# Run statistic:
gp16 = jsmetrics.jet_statistics.grise_polvani_2014(ua_sub)

jsmetrics.metrics.jet_statistics.barnes_polvani_2015(data)

This method calculates the jet positon and wind speed at that position by fitting a parabola around the maximum of zonally average u-component wind speed, using the magnitude at the maximum (‘jet_speed’) and latitude at that maximum (‘jet_lat’).

This method was originally introduce in Barnes & Polvani (2015) https://doi.org/10.1175/JCLI-D-14-00589.1 and is described in Section 2b of that study.

Please see ‘Notes’ below for any additional information about the implementation of this method to this package.

Parameters:

dataxarray.Dataset

Data which should containing the variables: ‘ua’, and the coordinates: ‘lon’, ‘lat’, ‘plev’ and ‘time’.

Returns:

outputxarray.Dataset

Data containing the two outputs: ‘jet_lat’ and ‘jet_speed’

Notes

This methodology make an assumption that the a parabola can be fit to windspeed profile, so it performs quite different from other jet latitude methods available in the package in cases where the windspeed profile is more complex (multiple jets), and on data with finer temporal resolution in general.

Examples

import jsmetrics
import xarray as xr

# Load in dataset with u component wind:
ua_data = xr.open_dataset('path_to_u_data')

# Subset dataset to range used in original methodology (700-925 hPa &  30-70N,  130-10W)):
ua_sub = ua.sel(plev=slice(700, 925), lat=slice(30, 70), lon=slice(230, 350))

# Run statistic:
bp15 = jsmetrics.jet_statistics.barnes_polvani_2015(ua_sub)

jsmetrics.metrics.jet_statistics.barnes_simpson_2017(data)

This method calculates two outputs: ‘jet_lat’ and ‘jet_speed’ which are defined as the latitude and speed of the 10-day-averaged maximum zonally-averaged wind speed.

This method was originally introduce in Barnes & Simpson 2017 https://doi.org/10.1175/JCLI-D-17-0299.1 and is described in Section 2b of that study.

Please see ‘Notes’ below for any additional information about the implementation of this method to this package.

Parameters:

dataxarray.Dataset

Data which should containing the variables: ‘ua’, and the coordinates: ‘lon’, ‘lat’, ‘plev’ and ‘time’.

Returns:

outputxarray.Dataset

Data containing the x outputs: ‘jet_lat’ and ‘jet_speed’

Notes

The original methodology was intended to work on one pressure level (700 hPa) and on daily data, for the implementation included in this package, we have included methods to automatically average any inputted pressure levels and to return the data without 10-day averaging if data above the 10-day resolution is inputted. Instead, warnings are returned to user to help them use the method the way it was originally intended.

Examples

import jsmetrics
import xarray as xr

# Load in dataset with u component wind:
ua_data = xr.open_dataset('path_to_u_data')

# Subset dataset to range used in original methodology (700 hPa & North Atlantic & North Pacific)):
ua_na = ua.sel(plev=700, lat=slice(0, 90), lon=slice(280, 350)) # North Atlantic
ua_np = ua.sel(plev=700, lat=slice(0, 90), lon=slice(120, 230)) # North Pacific

# Run statistic:
bp17_na = jsmetrics.jet_statistics.barnes_simpson_2017(ua_na)
bp17_np = jsmetrics.jet_statistics.barnes_simpson_2017(ua_np)

jsmetrics.metrics.jet_statistics.bracegirdle_et_al_2018(data)

This method calculates the seasonal and annual jet-stream position (‘JPOS’) and strength (‘JSTR’) by applying a 0.075 degrees cubic spline interpolation to a zonally-averaged wind climatology and selecting the maximum.

This method was originally introduce in Bracegirdle et al (2018) https://doi.org/10.1175/JCLI-D-17-0320.1 and is described in Section 2 of that study.

Please see ‘Notes’ below for any additional information about the implementation of this method to this package.

Parameters:

dataxarray.Dataset

Data which should containing the variables: ‘ua’, and the coordinates: ‘lon’, ‘lat’, ‘plev’ and ‘time’.

Returns:

outputxarray.Dataset

Data containing the four outputs: ‘annual_JPOS’, ‘seasonal_JPOS’, ‘annual_JSTR’ and ‘seasonal_JSTR’

Notes

This method was originally developed for the jet streams in the Southern Hemisphere.

Examples

import jsmetrics
import xarray as xr

# Load in dataset with u component wind:
ua_data = xr.open_dataset('path_to_u_data')

# Subset dataset to range used in original methodology (850 hPa &  -75--10N)):
ua_sub = ua.sel(plev=850, lat=slice(-75, -10))

# Run statistic:
b18 = jsmetrics.jet_statistics.bracegirdle_et_al_2018(ua_sub)

jsmetrics.metrics.jet_statistics.ceppi_et_al_2018(data, lon_resolution=None)

This method calculates the jet latitude (‘jet-lat’) as defined by selecting the centroid of a zonally-averaged wind profile.

The centroid is calculated by:

\[\phi_{jet} = \frac{\int_{30°}^{60°} \phi\bar{u}^2, d\phi}{\int_{30°}^{60°} \bar{u}^2, d\phi}\]

This method has been slightly adapted to include a ‘jet_speed’ extraction (provided for this method in Screen et al. (2022) https://doi.org/10.1029/2022GL100523).

Note: The implementation here does not explicit limit the centroid calculation to latitude between 20°-70°, instead this range is determined by the input data.

This method was originally introduce in Ceppi et al (2018) https://doi.org/10.1175/JCLI-D-17-0323.1 and is described in Section 2b of that study.

Please see ‘Notes’ below for any additional information about the implementation of this method to this package.

Parameters:

dataxarray.Dataset

Data which should containing the variables: ‘ua’, and the coordinates: ‘lon’, ‘lat’, ‘plev’ and ‘time’.

lon_resolutionnumeric

Resolution to use for longitude coord if size 1

Returns:

outputxarray.Dataset

Data containing the three outputs: ‘jet_lat’, ‘jet_speed’, ‘total_area_m2’

Notes

This method was improved by Zappa et al. (2018) https://doi.org/10.1029/2019GL083653, which includes exclusion of \(<0 m s^{-1}\) u-wind. This method is also available in jsmetrics.

Examples

import jsmetrics
import xarray as xr

# Load in dataset with u component wind:
ua_data = xr.open_dataset('path_to_u_data')

# Subset dataset to range used in original methodology (850 hPa &  30-60N/S)):
ua_na = ua.sel(plev=850, lat=slice(30, 60), lon=slice(300, 60)) # North Atlantic-European Sector
ua_na = ua.sel(plev=850, lat=slice(30, 60), lon=slice(140, 240)) # North Pacific
ua_sh = ua.sel(plev=850, lat=slice(-60, -30)) # Southern Hemisphere

# Run statistic:
c18_na = jsmetrics.jet_statistics.ceppi_et_al_2018(ua_na)
c18_np = jsmetrics.jet_statistics.ceppi_et_al_2018(ua_np)
c18_sh = jsmetrics.jet_statistics.ceppi_et_al_2018(ua_sh)

jsmetrics.metrics.jet_statistics.zappa_et_al_2018(data, lon_resolution=None)

This method calculates the jet latitude (‘jet-lat’) as defined by selecting the centroid of a zonally-averaged wind profile.

The centroid is calculated by:

\[\phi_{jet} = \frac{\int_{20°}^{70°} \phi\bar{u}^2_0, d\phi}{\int_{20°}^{70°} \bar{u}^2_0, d\phi}\]

\[u_0(\phi) = \max(0, u(\phi))\]

Note: The implementation here does not explicit limit the centroid calculation to latitude between 20°-70°, instead this range is determined by the input data.

This method has been slightly adapted to include a ‘jet_speed’ extraction (provided for this method in Screen et al. (2022) https://doi.org/10.1029/2022GL100523).

This method was originally introduced in Zappa et al. 2018 https://doi.org/10.1029/2019GL083653 and is described in Section 2.3 of that study.

Please see ‘Notes’ below for any additional information about the implementation of this method to this package.

Parameters:

dataxarray.Dataset

Data which should containing the variables: ‘ua’, and the coordinates: ‘lon’, ‘lat’, ‘plev’ and ‘time’.

lon_resolutionnumeric

Resolution to use for longitude coord if size 1

Returns:

outputxarray.Dataset

Data containing the three outputs: ‘jet_lat’, ‘jet_speed’, ‘total_area_m2’

Notes

This method was adapted from and very similar to Ceppi et al (2018) https://doi.org/10.1175/JCLI-D-17-0323.1.

This method is also used in Ayres & Screen, 2019 and Screen et al. 2022. Similar methods are used in: Chen et al. 2008; Ceppi et al. 2014, Ceppi et al. 2018.

Examples

import jsmetrics
import xarray as xr

# Load in dataset with u component wind:
ua_data = xr.open_dataset('path_to_u_data')

# Subset dataset to range used in original methodology (850 hPa &  20-70N, 140∘E-240∘E or 300-360 E))):
ua_na = ua.sel(plev=850, lat=slice(20, 70), lon=slice(140, 240))
ua_np = ua.sel(plev=850, lat=slice(20, 70), lon=slice(300, 360))

# Run statistic:
z18_na = jsmetrics.jet_statistics.zappa_et_al_2018(ua_na)
z18_np = jsmetrics.jet_statistics.zappa_et_al_2018(ua_np)

jsmetrics.metrics.jet_statistics.kerr_et_al_2020(data, width_of_pulse=10)

This method defines the latitude and speed of the jet-stream where the maximum zonal winds occur for each longitude and for each time unit (i.e. day). These values are then smoothed across the longitudes with a rectangular pulse (by default this has a width of 10 degrees).

This method was originally introduced in Kerr et al. (2020) https://onlinelibrary.wiley.com/doi/10.1029/2020JD032735 and is described in Section 2.4.2 of that study.

Please see ‘Notes’ below for any additional information about the implementation of this method to this package.

Parameters:

dataxarray.Dataset

Data which should containing the variables: ‘ua’, and the coordinates: ‘lon’, ‘lat’, ‘plev’ and ‘time’.

Returns:

outputxarray.Dataset

Data containing the two outputs: ‘jet_lat’ and ‘smoothed_jet_lat’

Notes

This method was based on the method from Barnes and Fiore (2013) https://doi.org/10.1002/grl.50411

The implementation here returns both smoothed and unsmoothed jet latitude outputs.

Examples

import jsmetrics
import xarray as xr

# Load in dataset with u component wind:
ua_data = xr.open_dataset('path_to_u_data')

# Subset dataset to range used in original methodology (700 hPa &  20-70N,  300-360W)):
ua_sub = ua.sel(plev=700, lat=slice(20, 70), lon=slice(300, 360))

# Run statistic:
k20 = jsmetrics.jet_statistics.kerr_et_al_2020(ua_sub)

Jet core algorithms

Methods that return a mask of coordinates related to the jet location, e.g., identifying the maximum wind speed throughout the horizontal and/or vertical plane within a given time window.

Overview table:

Algorithm	Status		Algorithm	Status
Koch et al. 2006	Complete		Archer & Caldiera 2008	Complete
Schiemann et al. 2009	Complete		Manney et al. 2011	To verify
Pena-Ortiz et al. 2013	To verify		Kuang et al. 2014	To verify

* == help needed

Documentation:

Methods that return a mask of coordinates related to the jet location, e.g., identifying the maximum wind speed throughout the horizontal and/or vertical plane within a given time window.

The following algorithms each return a xarray.Dataset and are ordered by paper publish year.

jsmetrics.metrics.jet_core_algorithms.sort_xarray_data_coords(coords, ascending=True)

Sort the coordinates of the input data (e.g. lat or lon). Sort in ascending order by default. Used in all jet statistics and jet core algorithms included in this package.

Parameters:

coordsarray-like

Coords to sort in xarray data

ascendingBoolean (default=True)

Whether the coords should be sorted in ascending order

jsmetrics.metrics.jet_core_algorithms.koch_et_al_2006(data, ws_threshold=30)

This method follows a two-step procedure used to detect jet-event occurences (here: ‘jet_events_ws’).

The weighted average windspeed (\(\alpha vel\)) for the jet events is calculated as follows:

\[\alpha vel = \frac{1}{p2-p1} \int_{p1}^{p2} (u^2+v^2)^{1/2} \,dp\]

where \(p1\), \(p2\) is min, max pressure level.

After calculating \(\alpha vel\), in a second step a windspeed threshold to isolate jet events (the default is \(30 m s^{-1}\)).

This method was first introduced in Koch et al (2006) (https://doi.org/10.1002/joc.1255) and is described in section 2.2.2 of that study. The original methodology provides a third step (to produce a climatology of jet events), but this has been ignored in this implementation. Instead, we have provided an example of how to calculate this after running this method in ‘Examples’ below.

Please see ‘Notes’ below for any additional information about the implementation of this method to this package.

Parameters:

dataxarray.Dataset

Data which should containing the variables: ‘ua’ and ‘va’, and the coordinates: ‘lon’, ‘lat’, ‘plev’ and ‘time’.

ws_thresholdint or float

Windspeed threshold used to extract from weighted average (default: 30 ms-1)

Returns:

outputxr.Dataset

Data containing the output variable ‘jet_events_ws’

Notes

The original methodology uses windspeed between 100-400 hPa to calculated the weighted average and 30 meters per second as the windspeed threshold.

This equation for this method is provided on pg 287 of the Koch et al. 2006 paper.

In the original paper, they accumulate the jet events into two-class jet typology (described in section 2.2.3 of Koch et al. 2006)

Examples

import jsmetrics
import xarray as xr

# Load in dataset with u and v components:
uv_data = xr.open_dataset('path_to_uv_data')

# Subset dataset to range used in original methodology (100-400 hPa)):
uv_sub = uv_data.sel(plev=slice(100, 400))

# Run algorithm:
koch_outputs = jsmetrics.jet_core_algorithms.koch_et_al_2006(uv_sub, ws_threshold=30)

# Produce climatology of jet occurence events for each season and each month:
koch_month_climatology = koch_outputs.groupby("time.month").mean("time")
koch_season_climatology = koch_outputs.groupby("time.season").mean("time")

jsmetrics.metrics.jet_core_algorithms.schiemann_et_al_2009(data, ws_threshold=30, u_threshold=0)

This method detects ‘jet occurrences’, whereby each jet occurence is detected based on three rules applied to inputted wind speed (\(V = [u, v]\)):

1. \(|V|\) is a local maxima in latitude and altitude plane

2. \(|V| \ge 30 m s^{-1}\)

3. \(u \ge 0 m s^{-1}\).

The implementation of this method here allows you to edit the \(|V|\) threshold (by changing ‘ws_threshold’), and \(u\) threshold (by changing ‘u_threshold’).

This method was originally introduce in Schiemann et al 2009 (https://doi.org/10.1175/2008JCLI2625.1) and is described in Section 2 of that study.

Please see ‘Notes’ below for any additional information about the implementation of this method to this package.

Parameters:

dataxarray.Dataset

Data which should containing the variables: ‘ua’ and ‘va’, and the coordinates: ‘lon’, ‘lat’, ‘plev’ and ‘time’.

ws_thresholdint or float

Windspeed threshold used to extract jet maxima from resultant windspeed (default: 30 m/s)

u_thresholdint or float

Windspeed threshold used to extract u-component wind speed (default: 0 m/s)

Returns:

outputxr.Dataset

Data containing the two output variables: ‘ws’ and ‘jet_occurence’

Notes

While the original method is built on a four dimension slice of wind speed (time, lat, lon, plev), This implementation will work where there is only one pressure level, so a 3-d slice (time, lat, lon).

Slow method: due to the nature of this method, it currently takes a moderately long time to run, i.e. 7.6 seconds per time unit on AMD Ryzen 5 3600 6-core processor.

Examples

import jsmetrics
import xarray as xr

# Load in dataset with u and v components:
uv_data = xr.open_dataset('path_to_uv_data')

# Subset dataset to range used in original methodology (100-500 hPa & 16.7-58.25 N, 42.5-220.5 E)):
uv_sub = uv_data.sel(plev=slice(100, 500), lat=slice(16.7, 58.25), lon=slice(42.5, 220.5))

# Run algorithm:
schiemann_outputs = jsmetrics.jet_core_algorithms.schiemann_et_al_2009(uv_sub, ws_threshold=30)

# Get jet occurence for first day in data
schiemann_jet_occurence_first_day = schiemann_outputs['jet_occurence'].isel(time=0).max('plev')

# Produce a jet occurence count across all pressure levels
schiemann_jet_counts_all_levels = schiemann_outputs['jet_occurence'].sum(('time', 'plev'))

# Use the jet occurence values as a mask to extract the jet windspeeds
schiemann_jet_ws = schiemann_outputs.where(schiemann_outputs['jet_occurence'] > 0)['ws']

jsmetrics.metrics.jet_core_algorithms.manney_et_al_2011(data, jet_core_plev_limit, jet_core_ws_threshold=40, jet_boundary_ws_threshold=30, ws_drop_threshold=25, jet_core_lat_distance=15, check_diagonals=False)

This method detects jet cores (within an altitude range see ‘jet_core_plev_limit’) and a boundary region surrounding those cores based on two windspeed thresholds. Two checks are applied after initial detection of cores to check whether boundaries with more then one core are part of the same feature (the default threshold for these boundaries is 30 m/s, see ‘jet_boundary_ws_threshold’). The two checks seperate cores based on whether the cores are more than a certain distance apart (default is 15 degrees, see ‘jet_core_lat_distance’) and whether the windspeed between two given cores does not drop below a windspeed threshold (default is 25 m/s, see ‘ws_drop_threshold’)

This method returns four outputs

1. jet_core_mask – Regions within each latitude-altitude slice that are local maxima and have windspeeds above the ‘jet_core_ws_threshold’

2. jet_region_mask – Regions above, below, left and right of any given jet core with windspeed above the ‘jet_boundary_ws_threshold’

3. jet_region_contour_mask – All contigious regions of windspeeds emcompassing a jet core above the ‘jet_boundary_ws_threshold’ (i.e. not just above, below, left and right)

4. ws – Resultant wind speed calculated from ‘ua’, ‘va’ inputs.

This method was originally introduce in Manney et al. (2011) (https://doi.org/10.5194/acp-11-6115-2011), and is described in Section 3.1 of that study. This method is also known as the JETPAC (Jet and Tropopause Products for Analysis and Characterization) software package, and available in its original form at request to NASA JPL.

Please see ‘Notes’ below for any additional information about the implementation of this method to this package.

Parameters:

dataxarray.Dataset

Data which should containing the variables: ‘ua’ and ‘va’, and the coordinates: ‘lon’, ‘lat’, ‘plev’ and ‘time’.

jet_core_plev_limit: tuple or array

Sequence of two values relating to the pressure level limit of the jet cores (original paper uses 100-400 hPa)

jet_core_ws_thresholdint or float

Threshold used for jet cores (default=40 m/s)

jet_boundary_ws_thresholdint or float

Threshold for jet boundaries (default=30 m/s)

ws_drop_thresholdint or float

Threshold for drop in windspeed along direct interpolated path between cores (default: 25 m/s)

jet_core_lat_distanceint or float

Threshold for maximum distance between cores to be counted the same (default: 15 degrees)

check_diagonalsbool

Whether to check the diagonal edges of each maxima in the latitude-altitude plane. The original method does not include this (default: False).

Returns:

outputxarray.Dataset

Data containing the four output variables: ‘ws’, ‘jet_region_mask’, ‘jet_region_contour_mask’, and ‘jet_core_mask’

Notes

The implementation of this method varies slightly from the original, because this method will return a mask rather than dynamical values. The intention of this, is to allow these masks to be used to subset other variables such as windspeed (see ‘Examples’ for demonstration of how to use these masks).

There is an update to this method introduced in Manney & Hegglin 2018 to include physically-based method to extract the subtropical jet is identified (and thus distinguished from polar jets). This update is not included in jsmetrics.

‘jet_region_above_ws_threshold_mask’ is provided here, but not explicitly in the original methodology. It is meant as an alternative to using a contour to check which regions encompass jet cores.

Examples

import jsmetrics
import xarray as xr

# Load in dataset with u and v components:
uv_data = xr.open_dataset('path_to_uv_data')

# Subset dataset to range appropriate for original methodology (100-1000 hPa)):
uv_sub = uv_data.sel(plev=slice(100, 1900))

# Run algorithm:
manney_outputs = jsmetrics.jet_core_algorithms.manney_et_al_2011(uv_sub, ws_core_threshold=40, ws_boundary_threshold=30, jet_core_plev_limit=(100, 400))

# Use the jet core mask to extract the jet windspeeds
manney_jet_ws = manney_outputs.where(manney_outputs['jet_core_mask'])['ws']

jsmetrics.metrics.jet_core_algorithms.penaortiz_et_al_2013(data, ws_threshold=30)

This method follows a two step procedure for calculating local wind maxima and then subcategorising the local maxima into two distinct jet masks: the Subtropical Jet (STJ) and Polar Front Jet (PFJ).

This method returns 4 outputs:

1. local_wind_maxima – Binary mask of local wind maxima

2. local_wind_maxima_by_monthyear – Same as above, but for monthyear frequency

3. polar_front_jet – Binary mask of the PFJ (by monthyear).

4. subtropical_jet – Binary ask of the STJ (by monthyear).

This method was first introduced in Pena-Ortiz et al. (2013) (https://doi.org/10.1002/jgrd.50305) and is described in section 2 of that study.

Please see ‘Notes’ below for any additional information about the implementation of this method to this package.

Parameters:

dataxarray.Dataset

Data which should containing the variables: ‘ua’ and ‘va’, and the coordinates: ‘lon’, ‘lat’, ‘plev’ and ‘time’.

ws_thresholdint or float

windspeed threshold to apply (default=30 ms-1)

Returns:

outputxarray.Dataset

Data containing the four output variables: ‘local_wind_maxima’, ‘local_wind_maxima_by_monthyear’, ‘polar_front_jet’, and ‘subtropical_jet’

Notes

See Table 1 in the respective paper for the categories used to seperate the STJ and PFJ. The STJ is only seperated in DJF for the Northern Hemisphere.

Slow method: currently takes a long time i.e. 1.3 seconds per time unit with 8 plevs (i.e. 1.3 seconds per day) on a AMD Ryzen 5 3600 6-core processor.

Examples

import jsmetrics
import xarray as xr

# Load in dataset with u and v components:
uv_data = xr.open_dataset('path_to_uv_data')

# Subset dataset to range used in original methodology (100-400 hPa)):
uv_sub = uv_data.sel(plev=slice(100, 400))

# Run algorithm:
pena_outputs = jsmetrics.penaortiz_et_al_2013(uv_sub)

# Produce a count of polar front jet values across all pressure levels
pena_pfj_counts_all_levels = pena_outputs['polar_front_jet'].sum(('monthyear', 'plev'))

# Use the polar front jet values as a mask to extract the jet windspeeds
uv_sub["ws_by_monthyear"] = (
    data["ws"]
    .resample(time="MS")
    .sum()
    .rename({"time": "monthyear"})
)
pena_pfj_ws = pena_outputs.where(pena_outputs['polar_front_jet'] > 0)['ws_by_monthyear']

jsmetrics.metrics.jet_core_algorithms.kuang_et_al_2014(data, occurence_ws_threshold=30)

This method produces an event-based jet occurrences and jet center occurrences of the jet stream in a given atmospheric column.

The outputs of this method will produce categorical values of three types:

0. is not determined to be part of the jet

1. is a jet occurence

2. is jet core (upgraded from a jet occurence)

This method was first introduced in Kuang et al (2014) (https://doi.org/10.1007/s00704-013-0994-x) and is described in section 2 of that study.

Please see ‘Notes’ below for any additional information about the implementation of this method to this package.

Parameters:

dataxarray.Dataset

Data which should containing the variables: ‘ua’ and ‘va’, and the coordinates: ‘lon’, ‘lat’, ‘plev’ and ‘time’.

occurence_ws_thresholdint or float

Threshold used to identify a jet-stream occurence point (default=30)

Returns:

outputxarray.Dataset

Data containing jet-occurence and jet-centres (1 is occurence, 2 is core, 0 is no jet)

Notes

Slow method: currently takes a long time i.e. 2 seconds per time unit with 1 plev (i.e. 2 seconds per day) on AMD Ryzen 5 3600 6-core processor

Examples

import jsmetrics
import xarray as xr

# Load in dataset with u and v components:
uv_data = xr.open_dataset('path_to_uv_data')

# Subset dataset to range used in original methodology (200 hPa)):
uv_sub = uv_data.sel(plev=slice(200, 200))

# Run algorithm:
kuang_outputs = jsmetrics.kuang_et_al_2014(uv_sub, occurence_ws_threshold=30)

# Extract only jet centers
kuang_jet_centers = kuang_outputs.where(kuang_outputs['jet_ocurrence1_jet_centre2']==2)

# Produce a count of jet centers across all pressure levels
kuang_jet_centers_counts_all_levels = kuang_jet_centers.sum(('time', 'plev'))

# Use the jet occurence values as a mask to extract the jet windspeeds
kuang_jet_ws = kuang_outputs.where(kuang_outputs['jet_ocurrence1_jet_centre2'] > 0)['ws']

jsmetrics.metrics.jet_core_algorithms.jet_core_identification_algorithm(data, ws_core_threshold=40, ws_boundary_threshold=30)

This method extract seperate jet cores based on boundary and core windspeed thresholds.

The output variable of this method includes two types:

• 0 – regions not determined to be part of the jet

• 1-n – Seperate jet core regions seperated by one condition: if two cores in the same region are more than 15 degrees of latitude away

This method is inspired by the method from Manney et al. (2011) (https://doi.org/10.5194/acp-11-6115-2011), which is described in Section 3.1 of that study.

Parameters:

dataxarray.Dataset

Data which should containing the variables: ‘ua’ and ‘va’, and the coordinates: ‘lon’, ‘lat’, ‘plev’ and ‘time’.

ws_core_thresholdint or float

Threshold used for jet-stream core point (default=40)

ws_boundary_thresholdint or float

Threshold for jet-stream boundary point (default=30)

Returns:

outputxarray.Dataset

Data containing the variable ‘jet-core_id’ (ID number relates to each unique core)

Notes

Slow method: currently takes a long time i.e. 7.6 seconds per time unit with 8 plevs (i.e. 7.6 seconds per day) on AMD Ryzen 5 3600 6-core processor

Examples

import jsmetrics
import xarray as xr

# Load in dataset with u and v components:
uv_data = xr.open_dataset('path_to_uv_data')

# Subset dataset to range appropriate for methodology (100-400 hPa)):
uv_sub = uv_data.sel(plev=slice(100, 400))

# Run algorithm:
jca_outputs = jsmetrics.jet_core_algorithms.jet_core_identification_algorithm(uv_sub, ws_core_threshold=40, ws_boundary_threshold=30)

Waviness metrics

Statistics and algorithms for determining the “waviness” of upper-level mean flow within a given time window. These metrics only have meaning at an integrated global scale.

Overview table:

Algorithm	Status		Algorithm	Status
Francis & Vavrus 2015	Complete		Cattiaux et al. 2009	To verify
Local Wave Activity	In progress*

* == help needed

Documentation:

Statistics and algorithms for determining the “waviness” of upper-level mean flow within a given time window. These metrics only have meaning at an integrated global scale.

The following metrics each return a xarray.Dataset and are ordered by paper publish year.

jsmetrics.metrics.waviness_metrics.sort_xarray_data_coords(coords, ascending=True)

Sort the coordinates of the input data (e.g. lat or lon). Sort in ascending order by default. Used in all jet statistics and jet core algorithms included in this package.

Parameters:

coordsarray-like

Coords to sort in xarray data

ascendingBoolean (default=True)

Whether the coords should be sorted in ascending order

jsmetrics.metrics.waviness_metrics.francis_vavrus_2015(data)

This method calculates a waviness metric: Meridional Circulation Index (MCI) from u- and v-components of wind.

MCI is calculated as follows:

\[MCI = \frac{v*|v|}{u^2+v^2}\]

When MCI = 0, the wind is purely zonal, and when MCI= 1 (-1), the flow is from the South (North).

This method was originally introduce in Francis & Vavrus (2015) https://doi.org/10.1088/1748-9326/10/1/014005 and is described in the Section entitled: ‘Meridional circulation Index (MCI)’.

Please see ‘Notes’ below for any additional information about the implementation of this method to this package.

Parameters:

dataxarray.Dataset

Data which should containing the variables: ‘ua’ and ‘va’, and the coordinates: ‘lon’, ‘lat’, ‘plev’ and ‘time’.

Returns:

outputxarray.Dataset

Data containing the one output: ‘mci’

Notes

In the original methodology, MCI is expressed in terms of seasonal anomaly, we show you how to do this in ‘Examples’

Examples

import jsmetrics
import xarray as xr

# Load in dataset with u and v components:
uv_data = xr.open_dataset('path_to_uv_data')

# Subset dataset to range used in original methodology (500 hPa & 20-80 N)):
uv_sub = uv_data.sel(plev=500, lat=slice(20, 80))

# Run statistic:
fv15 = jsmetrics.waviness_metrics.francis_vavrus_2015(uv_sub)

# Express MCI as a seasonal anomaly
fv15_seasonal_anomalies = (mci['mci'] - mci['mci'].groupby('time.season').mean())

jsmetrics.metrics.waviness_metrics.cattiaux_et_al_2016(data)

This method calculates a sinousity metric for upper-air flow using geopotential height. The value of sinuosity is selected using an isohypse which precisely corresponds to the Z500 average over 30-70∘N . Then this value is compared to the perimeter at 50∘N to calculate sinuosity.

This method was first introduce in Cattiaux et al (2016) https://doi.org/10.1002/2016GL070309 and is described in section 3.1 of that study.

Please see ‘Notes’ below for any additional information about the implementation of this method to this package.

Parameters:

dataxarray.Dataset

Data which should containing the variables: ‘zg’, and the coordinates: ‘lon’, ‘lat’, ‘plev’ and ‘time’.

Returns:

outputxarray.Dataset

Data containing the two outputs: ‘sinousity’ and ‘zonal_mean_zg_30Nto70N’

Notes

The original implementation used the R package ‘geosphere’ to calculate sinuosity. Here we use the Python package Shapely to calculate great circle distances between the two perimeters to calculate sinuosity.

Moderately slow method: currently takes a moderate amount of time i.e. 2 seconds per 100 time unitswith 1 plev on AMD Ryzen 5 3600 6-core processor

Examples

import jsmetrics
import xarray as xr

# Load in dataset with u and v components:
zg_data = xr.open_dataset('path_to_zg_data')

# Subset dataset to range used in original methodology (500 hPa & 00-90 N)):
zg_sub = zg_data.sel(plev=500, lat=slice(0, 90))

# Run statistic:
c16 = jsmetrics.waviness_metrics.cattiaux_et_al_2016(zg_sub)

See the Metric sub-components for more detail about the implementation of the metrics included in this package.