Calculation

Calculation#

The Skyborn calculation module provides statistical, atmospheric, and mathematical functions for climate data analysis.

Atmospheric Physics Functions#

Tropopause Calculations#

trop_wmo(temperature, pressure, xdim, ydim, levdim, timedim=None, pressure_unit='hPa', lapse_criterion=2.0, missing_value=-999.0, check_pressure_order=True)[source]#

Calculate WMO tropopause properties for multi-dimensional atmospheric data.

This function processes gridded atmospheric data to find the tropopause following the WMO definition (lapse rate < 2 K/km). Designed for isobaric (constant pressure level) data.

Parameters:

  • temperature (array_like) – Atmospheric temperature data [K] on isobaric levels. Must be ordered to correspond with the pressure levels.

  • pressure (array_like) –

    Atmospheric pressure data [hPa or Pa]. Can be either:

    • 1D array: Pressure levels (recommended for isobaric data). Length must match the level dimension of temperature.

    • Multi-dimensional array: Same shape as temperature (legacy support).

    CRITICAL: Pressure levels must be sorted in ascending order (from low pressure/high altitude to high pressure/low altitude). This is required by the underlying WMO tropopause algorithm.

  • xdim (int) – Longitude dimension index

  • ydim (int) – Latitude dimension index

  • levdim (int) – Vertical level dimension index

  • timedim (int, optional) – Time dimension index for 4D data

  • pressure_unit (str, default 'hPa') – Unit of pressure data (‘hPa’ or ‘Pa’)

  • lapse_criterion (float, default 2.0) – WMO lapse rate criterion [K/km] for tropopause definition

  • missing_value (float, default -999.0) – Value to use for missing data

  • check_pressure_order (bool, default True) – Whether to check and enforce ascending pressure order along the level dimension

Returns:

Dictionary containing: - ‘pressure’: Tropopause pressure [hPa] - ‘height’: Tropopause height [m] - ‘level_index’: Tropopause level index (0-based) - ‘lapse_rate’: Tropopause lapse rate [K/km] - ‘success’: Success flag for each grid point

Return type:

dict

Examples

Calculate tropopause for 4D isobaric data (time, level, lat, lon):

>>> # Define isobaric pressure levels (ascending order: high altitude to surface)
>>> pressure_levels = np.array([10, 20, 50, 100, 200, 300, 500, 700, 850, 1000])  # hPa
>>> temperature = 300 - np.random.rand(12, 10, 180, 360) * 80  # K
>>> result = trop_wmo(temperature, pressure_levels,
...                   xdim=3, ydim=2, levdim=1, timedim=0)
>>> print(f"Pressure shape: {result['pressure'].shape}")  # (12, 180, 360)
>>> print(f"Height shape: {result['height'].shape}")      # (12, 180, 360)

Calculate tropopause for 3D isobaric data (level, lat, lon):

>>> pressure_levels = np.array([50, 100, 200, 300, 500, 700, 850, 1000])  # hPa
>>> temperature_3d = 300 - np.random.rand(8, 180, 360) * 80  # K
>>> result = trop_wmo(temperature_3d, pressure_levels,
...                   xdim=2, ydim=1, levdim=0)
>>> print(f"Pressure shape: {result['pressure'].shape}")  # (180, 360)

Notes

  • This function is optimized for isobaric data (constant pressure levels).

  • For model level data, first interpolate to pressure levels.

  • Requires compiled Fortran extensions. Install with: pip install skyborn[fortran]

  • The underlying algorithm follows the WMO (1957) tropopause definition.

trop_wmo_profile(temperature, pressure, pressure_unit='hPa', lapse_criterion=2.0, missing_value=-999.0)[source]#

Calculate WMO tropopause properties for a single vertical profile.

This function processes a single atmospheric vertical profile to find the tropopause following the WMO definition (lapse rate < 2 K/km). Optimized for profile data analysis with isobaric level data.

Parameters:

  • temperature (array_like) – Atmospheric temperature profile [K] on isobaric levels. 1D array.

  • pressure (array_like) – Atmospheric pressure profile [hPa or Pa]. Must have same length as temperature. CRITICAL: Pressure levels must be sorted in ascending order (from low pressure/high altitude to high pressure/low altitude). This is required by the underlying WMO tropopause algorithm.

  • pressure_unit (str, default 'hPa') – Unit of pressure data (‘hPa’ or ‘Pa’)

  • lapse_criterion (float, default 2.0) – WMO lapse rate criterion [K/km] for tropopause definition

  • missing_value (float, default -999.0) – Value to use for missing data

Returns:

Dictionary containing: - ‘pressure’: Tropopause pressure [hPa] (float) - ‘height’: Tropopause height [m] (float) - ‘level_index’: Tropopause level index (int, 0-based) - ‘lapse_rate’: Tropopause lapse rate [K/km] (float) - ‘success’: Success flag (bool)

Return type:

dict

Examples

Basic usage with isobaric profile data:

>>> import numpy as np
>>> from skyborn.calc.troposphere.core import trop_wmo_profile
>>>
>>> # Create atmospheric profile on isobaric levels (ascending order)
>>> pressure = np.array([10, 20, 50, 100, 200, 300, 500, 700, 850, 1000])  # hPa
>>> temperature = np.array([190, 200, 210, 220, 230, 250, 270, 280, 285, 288])  # K
>>>
>>> result = trop_wmo_profile(temperature, pressure)
>>> print(f"Tropopause pressure: {result['pressure']:.1f} hPa")
>>> print(f"Tropopause height: {result['height']:.0f} m")

Notes

  • This function is optimized for isobaric data (constant pressure levels).

  • For model level data, first interpolate to pressure levels.

  • Uses dedicated Fortran routine for single profile analysis.

  • For gridded data, use trop_wmo() instead.

trop_wmo(temperature, pressure=None, xdim=None, ydim=None, levdim=None, timedim=None, pressure_unit='hPa', lapse_criterion=2.0, missing_value=-999.0, keep_attrs=True, auto_sort_levels=True)[source]#

Calculate WMO tropopause properties for xarray DataArrays.

This function processes gridded atmospheric data to find the tropopause following the WMO definition (lapse rate < 2 K/km). It automatically detects coordinate dimensions and preserves all metadata.

* DESIGNED FOR ISOBARIC DATA *

Parameters:

  • temperature (xarray.DataArray) – Atmospheric temperature data [K] on isobaric (constant pressure) levels. Can be a 1D profile, a 2D (level, lat) / (level, lon) cross-section, a 3D (level, lat, lon) grid or (time, level, lat) / (time, level, lon) time-varying cross-section, or a 4D (time, level, lat, lon) array. Must have level coordinate (for pressure generation) if pressure is not provided.

  • pressure (xarray.DataArray, optional) – Atmospheric pressure data [hPa or Pa] on isobaric levels. If None, will be automatically generated from temperature’s level coordinate. CRITICAL: Pressure levels MUST be sorted in ASCENDING order (from low pressure/high altitude to high pressure/low altitude). This is required by the underlying WMO tropopause algorithm.

  • xdim (int or str, optional) – Longitude dimension index/name. Auto-detected if None. Not required for 1D profiles. For (time, level, lat) cross-sections, the missing longitude axis may be left as None or passed explicitly as -1.

  • ydim (int or str, optional) – Latitude dimension index/name. Auto-detected if None. Not required for 1D profiles. For (time, level, lon) cross-sections, the missing latitude axis may be left as None or passed explicitly as -1.

  • levdim (int or str, optional) – Vertical level dimension index/name. Auto-detected if None.

  • timedim (int or str, optional) – Time dimension index/name. Auto-detected if None.

  • pressure_unit (str, default 'hPa') – Unit of pressure data (‘hPa’ or ‘Pa’)

  • lapse_criterion (float, default 2.0) – WMO lapse rate criterion [K/km] for tropopause definition

  • missing_value (float, default -999.0) – Value to use for missing data

  • keep_attrs (bool, default True) – Preserve input DataArray attributes in output

  • auto_sort_levels (bool, default True) – Automatically sort pressure levels in ascending order along the level dimension

Returns:

Dataset containing tropopause properties:

For multi-dimensional data (2D, 3D, 4D): - ‘pressure’: Tropopause pressure [hPa] with spatial/temporal coordinates - ‘height’: Tropopause height [m] with spatial/temporal coordinates - ‘level_index’: Tropopause level index (0-based) - ‘lapse_rate’: Tropopause lapse rate [K/km] - ‘success’: Success flag for each grid point

For 1D profile data: - Same variables but as scalar values (0D DataArrays)

Return type:

xarray.Dataset

Examples

1D Profile Analysis (Isobaric Data):

>>> import xarray as xr
>>> import numpy as np
>>> from skyborn.calc.troposphere.xarray import trop_wmo
>>>
>>> # Create 1D isobaric profile (ascending pressure order)
>>> temp_profile = xr.DataArray(
...     [210, 230, 250, 270, 280, 288],  # Temperature decreasing with altitude
...     dims=['level'],
...     coords={'level': [100, 300, 500, 700, 850, 1000]}  # hPa - ASCENDING order
... )
>>> result = trop_wmo(temp_profile)
>>> print(f"Tropopause: {float(result.pressure)} hPa, {float(result.height)} m")

Simplified Interface (Auto-pressure generation from isobaric levels):

>>> # Load ERA5 pressure level data (already isobaric)
>>> ds = xr.open_dataset('era5_pressure_levels.nc')  # Has 'level' coordinate in hPa
>>> result = trop_wmo(ds.temperature)  # Pressure auto-generated from level coordinate
>>> print(f"Tropopause pressure shape: {result.pressure.shape}")

2D Spatial Analysis (Isobaric Cross-sections):

>>> # Analyze latitude or longitude cross-sections on isobaric levels
>>> temp_2d = ds.temperature.isel(time=0, lon=0)  # (level, lat) - isobaric levels
>>> result = trop_wmo(temp_2d)
>>> # Result shape: (lat,)

3D Time-Varying Cross-section (Isobaric Data):

>>> temp_time_lat = ds.temperature.isel(lon=0)  # (time, level, lat)
>>> result = trop_wmo(temp_time_lat)
>>> # Result shape: (time, lat)
>>>
>>> # Explicitly mark the missing horizontal axis if desired
>>> result = trop_wmo(
...     temp_time_lat,
...     xdim=None,  # xdim=-1 is also accepted in the xarray wrapper
...     ydim='lat',
...     levdim='level',
...     timedim='time',
... )

Advanced usage with explicit isobaric pressure:

>>> # Ensure pressure levels are in ascending order
>>> result = trop_wmo(
...     temperature_data,  # On isobaric levels
...     pressure=pressure_data,  # Corresponding isobaric pressure levels
...     xdim='longitude', ydim='latitude', levdim='level',
...     lapse_criterion=2.5  # Custom WMO criterion
... )

4D Time Series (Isobaric Data):

>>> # Multi-year isobaric climate data
>>> result = trop_wmo(temperature_4d)  # (time, level, lat, lon) - isobaric levels
>>> # Result preserves time dimension: (time, lat, lon)
>>> seasonal_mean = result.height.groupby('time.season').mean()

Notes

  • This function is optimized for ISOBARIC data (constant pressure levels).

  • For model level data, first interpolate to pressure levels before using this function.

  • Requires compiled Fortran extensions. Install with: pip install skyborn[fortran]

  • The underlying algorithm follows the WMO (1957) tropopause definition.

The function automatically: - Detects spatial and temporal coordinates using metadata - Handles missing values (NaN or masked arrays) - Preserves all coordinate information and attributes - Works with multi-dimensional isobaric data

See also

skyborn.calc.troposphere.core.trop_wmo

Lower-level function for numpy arrays

Growth-rate Diagnostics#

baroc_growth_rate(u, temperature, pressure, *, lat=None, lat_bounds=None, lon=None, wavenumber_mode=DEFAULT_WAVENUMBER_MODE, method='log', tropopause_pressure=None, solver_levels=DEFAULT_SOLVER_LEVELS, smooth_window=DEFAULT_SMOOTH_WINDOW)[source]#

Compute the maximum baroclinic normal-mode growth rate.

Parameters:

  • u (xarray.DataArray, numpy.ndarray) – One-dimensional zonal-wind profile in m s^-1 or a two-dimensional stack of profiles that shares one vertical axis with pressure.

  • temperature (xarray.DataArray, numpy.ndarray) – One-dimensional temperature profile in Kelvin or a two-dimensional stack of profiles that shares one vertical axis with pressure.

  • pressure (xarray.DataArray, numpy.ndarray) – One-dimensional pressure coordinate in Pa or hPa. Values must be strictly monotonic and strictly positive.

  • lat (scalar, xarray.DataArray, numpy.ndarray, optional) –

    Latitude in degrees north. This may be a scalar or, for batched profile input, a one-dimensional vector matching the profile dimension. The compiled solver uses this latitude to build the local Coriolis parameter and meridional gradient of planetary vorticity, unless lat_bounds is provided.

    When lat_bounds is given and raw NumPy inputs still contain an explicit latitude axis, lat may instead be the one-dimensional latitude coordinate used to select and cosine-weight that latitude band before the remaining profile solve.

  • lat_bounds (tuple of float, optional) –

    Two-element latitude range in degrees north. When provided, the baroclinic solver uses the Chemke-style band-mean definitions

    \[\sin_{\mathrm{avg}} = \frac{\sin(\phi_1) + \sin(\phi_2)}{2}, \qquad \cos_{\mathrm{avg}} = \frac{\cos(\phi_1) + \cos(\phi_2)}{2},\]

    together with

    \[f = 2\Omega \sin_{\mathrm{avg}}, \qquad \beta = \frac{2\Omega\cos_{\mathrm{avg}}}{a},\]

    to match the latitude-band formulation used in the reference Chemke scripts.

    If u and temperature still carry a latitude axis, this option also reduces that axis with cosine-latitude weighting before the vertical-profile solve. For xarray inputs, the latitude coordinate is taken from the lat or latitude dimension. For NumPy inputs with three dimensions, pass the one-dimensional latitude coordinate through lat=... so the requested latitude band can be selected explicitly. The recommended NumPy layout is (time, level, lat). More generally, the wrapper requires exactly one axis matching the pressure-coordinate length and one non-pressure axis matching the latitude-coordinate length.

  • lon (xarray.DataArray, numpy.ndarray, optional) – One-dimensional longitude coordinate in degrees used only when wavenumber_mode="low". The coordinate must be strictly monotonic. If omitted in low-resolution mode, the wrapper uses the original Chemke Python-share default np.arange(0.0, 361.0, 1.5).

  • wavenumber_mode ({"high", "low"}, optional) –

    Zonal-wavenumber grid used by the compiled instability solver. "high" uses the fixed high-resolution Chemke/MATLAB-style grid k = 0, 1\times10^{-7}, \dots, 1.99\times10^{-5}\,\mathrm{m}^{-1}. "low" uses the lower-resolution Chemke Python-share formula

    \[k = \frac{2\pi w}{\left|\lambda_0 - \lambda_1\right|\pi a \cos\phi / 180},\]

    where \(w = 0, 1, \dots, \mathrm{round}(N_{\lambda}/3)-1\). When lon is omitted, the low-resolution mode uses the original Chemke Python-share longitude grid np.arange(0.0, 361.0, 1.5). Defaults to "high".

  • method ({"log", "linear"}, optional) – Vertical interpolation method used to remap the input profiles onto the solver pressure grid. "log" means linear interpolation in log-pressure, not logarithmic interpolation of the field itself. Defaults to "log". The legacy alias "logp" is still accepted.

  • tropopause_pressure (scalar, xarray.DataArray, numpy.ndarray, optional) – Explicit tropopause pressure in Pa or hPa. If provided, the automatic WMO tropopause diagnosis is skipped and the solver grid is built from this pressure to the lower-tropospheric bound. For batched profile input, this may be a scalar climatology or a one-dimensional vector matching the profile dimension.

  • solver_levels (int, optional) – Number of pressure levels used when the solver grid is built automatically from the diagnosed or explicit tropopause pressure. Defaults to DEFAULT_SOLVER_LEVELS.

  • smooth_window (int, optional) – Optional centered running-mean window applied to the growth-rate spectrum over zonal wavenumber inside the compiled Fortran backend before the final maximum-growth diagnostic is taken. 1 disables smoothing. If given, this must be a positive odd integer.

Returns:

Maximum baroclinic growth rate in s^-1. One-dimensional input returns a scalar float. Batched profile input returns a one-dimensional NumPy array or DataArray over the non-pressure dimension. Missing or unusable profiles return NaN.

Return type:

float, numpy.ndarray, xarray.DataArray

Notes

This function diagnoses the WMO tropopause pressure by default, or uses an explicit tropopause_pressure when provided, and then interpolates u and temperature to a fixed pressure grid between that tropopause and the lower troposphere. The thermodynamic profiles passed to the compiled kernel are

\[\theta = T\left(\frac{p_0}{p}\right)^{\kappa}, \qquad \rho = \frac{p}{R_d T}, \qquad N_z = -\frac{1}{\rho\theta}\frac{\partial \theta}{\partial p}.\]

The compiled solver then forms the finite-difference generalized eigenproblem

\[A(k)\psi = \lambda B(k)\psi,\]

where the vertical quasi-geostrophic operator is represented in pressure coordinates by

\[A(k) \sim k\beta - k^3 U + kUf^2 \partial_p \left(N_z^{-1}\partial_p\right) - kf^2 \left(\partial_p U\right) N_z^{-1}\partial_p,\]

and

\[B(k) \sim f^2 \partial_p \left(N_z^{-1}\partial_p\right) - k^2,\]

and the returned value is \(\max_k \operatorname{Im}(\lambda_k)\). If smooth_window is greater than 1, the compiled backend first applies a centered running mean along the discrete zonal-wavenumber spectrum and then evaluates the final Chemke-style turning-point maximum on that smoothed spectrum. When u or temperature contains multiple profiles, one-dimensional climatologies are broadcast across the profile dimension in Python, while the per-profile pressure interpolation and normal-mode solve loop run in the compiled Fortran backend.

The default DEFAULT_SOLVER_LEVELS = 45 is a practical resolution choice, not a formal optimum. Increasing solver_levels can improve vertical-grid convergence up to a point, but it also increases the per-profile eigenvalue cost and does not guarantee a better diagnostic once the solution is already converged on the chosen pressure interval.

References

Chemke, R., and Ming, Y. (2020):

Large Atmospheric Waves Will Get Stronger, While Small Waves Will Get Weaker by the End of the 21st Century. Geophysical Research Letters, 47, e2020GL090441. https://doi.org/10.1029/2020GL090441

Chemke, R., Zanna, L., Orbe, C., Sentman, L. T., and Polvani, L. M. (2022):

The Future Intensification of the North Atlantic Winter Storm Track: The Key Role of Dynamic Ocean Coupling. Journal of Climate, 35(8), 2407-2421. https://doi.org/10.1175/JCLI-D-21-0407.1

Chemke, R. (2022):

The future poleward shift of Southern Hemisphere summer mid-latitude storm tracks stems from ocean coupling. Nature Communications, 13, 2531. https://doi.org/10.1038/s41467-022-29392-4

Chemke, R., Ming, Y., and Yuval, J. (2022):

The intensification of winter mid-latitude storm tracks in the Southern Hemisphere. Nature Climate Change, 12, 553-557. https://doi.org/10.1038/s41558-022-01368-8

barot_growth_rate(u_barotropic, lat, *, radius=6_371_000.0, omega=7.292e-5)[source]#

Compute the maximum barotropic normal-mode growth rate from U(lat).

Parameters:

  • u_barotropic (xarray.DataArray, numpy.ndarray) – One-dimensional vertically averaged zonal-wind profile in m s^-1.

  • lat (xarray.DataArray, numpy.ndarray) – One-dimensional latitude coordinate in degrees north. Values must be strictly monotonic.

  • radius (float, optional) – Earth radius in meters. Only the default value is currently supported by the compiled kernel.

  • omega (float, optional) – Planetary rotation rate in s^-1. Only the default value is currently supported by the compiled kernel.

Returns:

Maximum barotropic growth rate in s^-1. Missing or unusable profiles return NaN.

Return type:

float

Notes

This diagnostic solves the finite-difference generalized eigenproblem

\[A(k)\psi = \lambda B(k)\psi\]

for each zonal wavenumber \(k\), where the continuous barotropic operator is

\[\left[k\beta - k^3 U + kU\partial_{yy} - kU_{yy}\right]\psi = \lambda \left(\partial_{yy} - k^2\right)\psi,\]

with

\[\beta = \frac{2\Omega\cos(\phi)}{a}.\]

The returned value is \(\max_k \operatorname{Im}(\lambda_k)\).

References

Chemke, R., Ming, Y., and Yuval, J. (2022):

The intensification of winter mid-latitude storm tracks in the Southern Hemisphere. Nature Climate Change, 12, 553-557. https://doi.org/10.1038/s41558-022-01368-8

Note

skyborn.calc.baroc_growth_rate uses method="log" by default. Here "log" means interpolation that is linear in log-pressure, not a logarithmic transform of the field values.

The public growth-rate wrappers return s^-1 directly. Missing or unusable profiles return NaN instead of a user-supplied missing-value marker.

By default, the compiled solver uses the fixed high-resolution Chemke/MATLAB-style zonal-wavenumber grid

\[k = 0, 1\times10^{-7}, \dots, 1.99\times10^{-5}\ \mathrm{m}^{-1}.\]

If you want the lower-resolution Chemke Python-share behavior instead, pass wavenumber_mode="low". When lon is omitted, Skyborn uses the original Chemke Python-share longitude grid np.arange(0.0, 361.0, 1.5) by default; you may still pass a custom one-dimensional longitude coordinate lon if needed. In that mode the solver uses

\[k = \frac{2\pi w}{\left|\lambda_0 - \lambda_1\right|\pi a \cos\phi / 180},\]

with \(w = 0, 1, \dots, \mathrm{round}(N_{\lambda}/3)-1\).

When a latitude band is more appropriate than a single representative latitude, pass the band endpoints and the compiled solver uses the same Chemke-style averages as the original reference scripts:

\[\sin_{\mathrm{avg}} = \frac{\sin(\phi_1) + \sin(\phi_2)}{2}, \qquad \cos_{\mathrm{avg}} = \frac{\cos(\phi_1) + \cos(\phi_2)}{2}\]
\[f = 2\Omega \sin_{\mathrm{avg}}, \qquad \beta = \frac{2\Omega\cos_{\mathrm{avg}}}{a}\]

If the input fields still carry a latitude axis, lat_bounds can also drive the latitude-band reduction step directly. Xarray inputs such as (time, level, lat) are selected and cosine-weighted over that band automatically before the compiled profile solver runs. For raw NumPy arrays with an explicit latitude axis, provide the one-dimensional latitude coordinate through lat=... together with lat_bounds so the requested band can be selected unambiguously. The recommended NumPy layout is (time, level, lat). More generally, the wrapper expects exactly one axis matching the pressure-coordinate length and one non-pressure axis matching the supplied latitude-coordinate length.

The function diagnoses the WMO tropopause pressure by default, or uses an explicit tropopause_pressure when given, and then builds a fixed solver grid with solver_levels=45 by default between that tropopause and the lower troposphere. This is a practical resolution choice for the compiled eigenvalue problem, not a universal optimum: increasing solver_levels can improve vertical-grid convergence up to a point, but it also raises the per-profile cost and is not automatically better once the solution is converged.

In the original Chemke-style analysis scripts, a parameter named window_size is often used for a centered running mean applied to the growth-rate spectrum over zonal wavenumber. In Skyborn this role is represented by smooth_window in baroc_growth_rate. That smoothing width is not the same thing as solver_levels. Even for a single atmospheric column, the compiled solver still sweeps many zonal wavenumbers and solves a generalized eigenvalue problem for each one, so the dominant cost is the repeated linear-algebra solve rather than tropopause diagnosis or pressure interpolation. The public default is smooth_window=1; values greater than 1 apply the centered running mean inside the Fortran backend before the final Chemke-style maximum-growth diagnostic is taken.

References

The current growth-rate implementation is documented against the following Chemke publications that are directly related to baroclinic wave growth and storm-track changes:

  • Chemke, R., and Ming, Y. (2020): Large Atmospheric Waves Will Get Stronger, While Small Waves Will Get Weaker by the End of the 21st Century. Geophysical Research Letters, 47, e2020GL090441. https://doi.org/10.1029/2020GL090441

  • Chemke, R., Zanna, L., Orbe, C., Sentman, L. T., and Polvani, L. M. (2022): The Future Intensification of the North Atlantic Winter Storm Track: The Key Role of Dynamic Ocean Coupling. Journal of Climate, 35(8), 2407-2421. https://doi.org/10.1175/JCLI-D-21-0407.1

  • Chemke, R. (2022): The future poleward shift of Southern Hemisphere summer mid-latitude storm tracks stems from ocean coupling. Nature Communications, 13, 2531. https://doi.org/10.1038/s41467-022-29392-4

  • Chemke, R., Ming, Y., and Yuval, J. (2022): The intensification of winter mid-latitude storm tracks in the Southern Hemisphere. Nature Climate Change, 12, 553-557. https://doi.org/10.1038/s41558-022-01368-8

Geostrophic Wind Calculations#

geostrophic_wind(z, glon, glat, dim_order, missing_value=-999.0)[source]#

Calculate geostrophic wind components from geopotential height.

This function can handle various input shapes by using windspharm’s prep_data and recover_data utilities to reshape data for batch processing.

Parameters:

  • z (ndarray) – Geopotential height data [gpm]. Can be 2D, 3D, or 4D. Must contain latitude (‘y’) and longitude (‘x’) dimensions.

  • glon (ndarray, shape (nlon,)) – Longitude coordinates in degrees

  • glat (ndarray, shape (nlat,)) – Latitude coordinates in degrees (south to north)

  • dim_order (str) – String specifying dimension order using: - ‘x’ for longitude - ‘y’ for latitude - ‘t’ for time - ‘z’ for level Example: ‘tzyx’ for (time, level, lat, lon)

  • missing_value (float, optional) – Missing value identifier (default: -999.0)

Returns:

  • ug (ndarray) – Zonal geostrophic wind component [m/s] (same shape as input z)

  • vg (ndarray) – Meridional geostrophic wind component [m/s] (same shape as input z)

Return type:

Tuple[ndarray, ndarray]

Examples

# 2D case: single time/level >>> z2d = np.random.randn(73, 144) # (lat, lon) >>> ug, vg = geostrophic_wind(z2d, glon, glat, ‘yx’)

# 3D case: multiple times >>> z3d = np.random.randn(73, 144, 12) # (lat, lon, time) >>> ug, vg = geostrophic_wind(z3d, glon, glat, ‘yxt’)

# 4D case: multiple levels and times >>> z4d = np.random.randn(12, 17, 73, 144) # (time, level, lat, lon) >>> ug, vg = geostrophic_wind(z4d, glon, glat, ‘tzyx’)

# Alternative 4D ordering >>> z4d_alt = np.random.randn(73, 144, 17, 12) # (lat, lon, level, time) >>> ug, vg = geostrophic_wind(z4d_alt, glon, glat, ‘yxzt’)

class GeostrophicWind(z, glon, glat, dim_order, missing_value=-999.0)[source]#

Class-based interface for geostrophic wind calculations.

This class provides a high-level interface similar to windspharm’s VectorWind, allowing for easy calculation of various geostrophic wind quantities.

Parameters:

  • z (ndarray) – Geopotential height data [gpm]

  • glon (ndarray) – Longitude coordinates [degrees]

  • glat (ndarray) – Latitude coordinates [degrees] (south to north)

  • dim_order (str) – Dimension ordering specification

  • missing_value (float, optional) – Missing value identifier (default: -999.0)

Examples

>>> # Create GeostrophicWind instance (longitude cyclicity auto-detected)
>>> gw = GeostrophicWind(z, glon, glat, 'tzyx')
>>>
>>> # Get wind components
>>> ug, vg = gw.uv_components()
>>>
>>> # Calculate derived quantities
>>> speed = gw.speed()
>>>
>>> # Access original data
>>> z_orig = gw.geopotential_height
__init__(z, glon, glat, dim_order, missing_value=-999.0)[source]#
Parameters:
property geopotential_height: ndarray#

Original geopotential height data.

property longitude: ndarray#

Longitude coordinates.

property latitude: ndarray#

Latitude coordinates.

uv_components()[source]#

Return zonal and meridional wind components.

Returns:

  • ug (ndarray) – Zonal (eastward) geostrophic wind component [m/s]

  • vg (ndarray) – Meridional (northward) geostrophic wind component [m/s]

Return type:

Tuple[ndarray, ndarray]

speed()[source]#

Calculate geostrophic wind speed.

Returns:

speed – Geostrophic wind speed [m/s]

Return type:

ndarray

geostrophic_wind(z, missing_value=-999.0, keep_attrs=True)[source]#

Calculate geostrophic wind components for xarray DataArrays.

This function processes geopotential height data to calculate geostrophic wind components (ug, vg). It automatically detects coordinate dimensions and preserves all metadata.

Parameters:

  • z (xarray.DataArray) – Geopotential height data [gpm]. Can be 2D, 3D, or 4D. Must contain longitude and latitude dimensions with coordinate information.

  • missing_value (float, optional) – Missing value identifier (default: -999.0)

  • keep_attrs (bool, optional) – Preserve input DataArray attributes in output (default: True)

Returns:

Dataset containing geostrophic wind components: - ‘ug’: Zonal geostrophic wind component [m/s] with spatial/temporal coordinates - ‘vg’: Meridional geostrophic wind component [m/s] with spatial/temporal coordinates

Return type:

xarray.Dataset

Examples

2D Geopotential Height Analysis:

>>> import xarray as xr
>>> import numpy as np
>>> from skyborn.calc.geostrophic.xarray import geostrophic_wind
>>>
>>> # Load 500 hPa geopotential height
>>> z = xr.open_dataarray('z500_era5.nc')  # Shape: (lat, lon)
>>> result = geostrophic_wind(z)
>>> print(f"Wind components: ug{result.ug.shape}, vg{result.vg.shape}")

3D Time Series Analysis:

>>> # Multi-time geopotential height data
>>> z_3d = xr.open_dataarray('z500_timeseries.nc')  # Shape: (time, lat, lon)
>>> result = geostrophic_wind(z_3d)
>>> # Result preserves time dimension: (time, lat, lon)
>>> monthly_mean = result.ug.groupby('time.month').mean()

4D Multi-level Analysis:

>>> # Multi-level, multi-time data
>>> z_4d = xr.open_dataarray('z_multilevel.nc')  # Shape: (time, level, lat, lon)
>>> result = geostrophic_wind(z_4d)
>>> # Result shape: (time, level, lat, lon)
>>> surface_winds = result.sel(level=1000)  # 1000 hPa level

Simplified Interface (No coordinate specification needed):

>>> # Automatic coordinate detection
>>> result = geostrophic_wind(z_data)  # Longitude cyclicity auto-detected
>>> print(f"Longitude cyclicity auto-detected: {result.attrs['longitude_cyclic']}")
>>> print(f"Latitude ordering: {result.attrs['latitude_ordering']}")

Notes

  • Longitude cyclicity is automatically detected from coordinate spacing

  • Latitude ordering is automatically ensured to be south-to-north as required

  • Requires compiled Fortran extensions for optimal performance

  • All coordinate information and attributes are preserved

The function automatically: - Detects longitude and latitude coordinates using metadata - Handles missing values (NaN or masked arrays) - Preserves all coordinate information and attributes - Works with multi-dimensional data of any supported shape

See also

skyborn.calc.geostrophic.core.geostrophic_wind

Lower-level function for numpy arrays

GeostrophicWind

Class-based interface for xarray DataArrays

class GeostrophicWind(z, missing_value=-999.0, keep_attrs=True)[source]#

Class-based geostrophic wind analysis using xarray DataArrays.

This class provides a high-level interface for geostrophic wind calculations that preserves xarray coordinate information and metadata. It wraps the standard interface implementation while maintaining CF-compliant attributes.

Parameters:

  • z (xarray.DataArray) – Geopotential height data [gpm]. Must contain longitude and latitude dimensions with appropriate coordinate information.

  • missing_value (float, optional) – Missing value identifier (default: -999.0)

  • keep_attrs (bool, optional) – Preserve input DataArray attributes in output (default: True)

_z_original#

Original geopotential height data

Type:

xarray.DataArray

_result#

Computed geostrophic wind components

Type:

xarray.Dataset

_glon#

Longitude coordinates

Type:

np.ndarray

_glat#

Latitude coordinates

Type:

np.ndarray

Examples

>>> import xarray as xr
>>> from skyborn.calc.geostrophic.xarray import GeostrophicWind
>>>
>>> # Load geopotential height
>>> z = xr.open_dataarray('z500.nc')
>>>
>>> # Create GeostrophicWind instance
>>> gw = GeostrophicWind(z)
>>>
>>> # Get wind components with preserved metadata
>>> ug, vg = gw.uv_components()
>>> print(ug.attrs)  # CF-compliant attributes
>>>
>>> # Calculate derived quantities
>>> speed = gw.speed()
>>> print(f"Max wind speed: {float(speed.max()):.1f} m/s")
>>>
>>> # Access original data
>>> z_orig = gw.geopotential_height
__init__(z, missing_value=-999.0, keep_attrs=True)[source]#

Initialize GeostrophicWind with xarray DataArray.

Parameters:
property geopotential_height: DataArray#

Original geopotential height data.

property longitude: ndarray#

Longitude coordinates.

property latitude: ndarray#

Latitude coordinates.

uv_components()[source]#

Return zonal and meridional wind components.

Returns:

  • ug (xarray.DataArray) – Zonal (eastward) geostrophic wind component [m/s]

  • vg (xarray.DataArray) – Meridional (northward) geostrophic wind component [m/s]

Return type:

Tuple[DataArray, DataArray]

speed()[source]#

Calculate geostrophic wind speed.

Returns:

speed – Geostrophic wind speed [m/s]

Return type:

xarray.DataArray

Genesis Potential Index (GPI) / Tropical Cyclone Potential Intensity#

potential_intensity(sst, psl, pressure_levels, temperature, mixing_ratio, actual_levels=None)[source]#

Calculate tropical cyclone potential intensity with automatic dimension detection.

This is the primary NumPy entry point. It automatically detects whether the thermodynamic input describes a single profile, a 3D gridded field, or a 4D time series and dispatches to the compiled backend directly.

Parameters:

  • sst (float or ndarray) – Sea surface temperature [K]

  • psl (float or ndarray) – Sea level pressure [Pa]

  • pressure_levels (ndarray) – Atmospheric pressure levels [mb]

  • temperature (ndarray) – Temperature data [K]. Supported shapes are (level,), (level, y, x), and (time, level, y, x).

  • mixing_ratio (ndarray) – Water vapor mixing ratio [kg/kg]

  • actual_levels (int, optional) – Number of valid vertical levels to use for a 1D profile input. Ignored for gridded 3D and 4D inputs.

Returns:

  • min_pressure (float or ndarray) – Minimum central pressure [mb]

  • max_wind (float or ndarray) – Maximum sustained wind speed [m/s]

  • error_flag (int) – Error status (1 = success, other values indicate non-convergence or invalid input)

Return type:

Tuple[Any, Any, int]

log_decompose_pi(pi, sst, t0, CKCD=0.9, *, sst_units='K')[source]#

Log-decompose potential intensity into efficiency, disequilibrium, and Ck/Cd.

Parameters:

  • pi (Any) – Potential intensity wind speed [m/s].

  • sst (Any) – Sea surface temperature in units given by sst_units.

  • t0 (Any) – Outflow temperature [K].

  • CKCD (float, default: 0.9) – Ratio of exchange coefficients.

  • sst_units ({"K", "C"}, default: "K") – Units of sst.

Returns:

(lnpi, lneff, lndiseq, lnCKCD) where lnpi = ln(V^2).

Return type:

tuple

pi_log_decomposition(sst, psl, pressure_levels, temperature, mixing_ratio, CKCD=0.9, *, outflow_source='cape_star')[source]#

Calculate PI plus the Wing et al. (2015) logarithmic decomposition.

Parameters:

  • sst (Any) – Sea surface temperature [K].

  • psl (Any) – Sea level pressure [Pa].

  • pressure_levels (ndarray) – Atmospheric pressure levels [mb].

  • temperature (ndarray) – Thermodynamic input field. Supported shapes are (level,), (level, y, x), and (time, level, y, x).

  • mixing_ratio (ndarray) – Water vapor mixing ratio [kg/kg] with the same shape as temperature.

  • CKCD (float, default: 0.9) – Ratio of exchange coefficients.

  • outflow_source ({"cape_star", "cape_env"}, default: "cape_star") – Outflow branch used for the backend diagnostics.

Returns:

max_wind, min_pressure, error_flag, t0, otl, lnpi, lneff, lndiseq, and lnCKCD.

Return type:

dict

potential_intensity(sst, psl, pressure_levels, temperature, mixing_ratio)[source]#

Calculate tropical cyclone potential intensity with automatic dimension detection.

This function automatically detects the input data dimensions and calls the appropriate specific function (profile, 3D, or 4D).

Parameters:

  • sst (xr.DataArray) – Sea surface temperature [K]

  • psl (xr.DataArray) – Sea level pressure [Pa]

  • pressure_levels (xr.DataArray) – Atmospheric pressure levels [mb]

  • temperature (xr.DataArray) – Temperature data [K]

  • mixing_ratio (xr.DataArray) – Water vapor mixing ratio [kg/kg]

Returns:

result – Dataset containing potential intensity results with appropriate dimensions.

Return type:

xr.Dataset

Notes

The function determines the calculation type based on temperature array dimensions: - 1D (vertical only): Single profile calculation - 3D (vertical + 2 spatial): Gridded data calculation - 4D (time + vertical + 2 spatial): Time series calculation

pi_log_decomposition(sst, psl, pressure_levels, temperature, mixing_ratio, CKCD=0.9, *, outflow_source='cape_star')[source]#

Return PI plus outflow and logarithmic decomposition diagnostics.

Parameters:
Return type:

Dataset

Note

error_flag == 1 indicates a successful PI solve.

potential_intensity(...) remains the pure PI entry point.

Use pi_log_decomposition(...) for supported 1D, 3D, and 4D inputs when you need tcpyPI-style outflow diagnostics (outflow_source={"cape_star", "cape_env"}, t0, otl) and the Wing et al. (2015) logarithmic decomposition where lnpi = ln(V^2).

Statistical Functions#

linear_regression(data, predictor)[source]#

Perform linear regression between a 3D data array and a predictor sequence. Handles both numpy arrays and xarray DataArrays with NaN handling.

Parameters:

  • data (np.ndarray or xr.DataArray) – A 3D array of shape (n_samples, dim1, dim2) containing dependent variables. Missing values should be represented as NaN.

  • predictor (np.ndarray or xr.DataArray) – A 1D array of shape (n_samples,) containing the independent variable. Missing values should be represented as NaN.

Returns:

A tuple containing: - regression_coefficients: The slope of the regression line with shape (dim1, dim2) - p_values: The p-values of the regression with shape (dim1, dim2)

Return type:

Tuple[np.ndarray, np.ndarray]

Raises:

ValueError – If the number of samples in data doesn’t match the length of the predictor.

spatial_correlation(data, predictor)[source]#

Perform fast spatial correlation between a 3D data array and a predictor time series. Optimized for vectorized operations to avoid slow loops over lat/lon grid points.

Parameters:

  • data (np.ndarray or xr.DataArray) – A 3D array of shape (time, lat, lon) containing spatial data. Missing values should be represented as NaN.

  • predictor (np.ndarray or xr.DataArray) – A 1D array of shape (time,) containing the predictor time series. Missing values should be represented as NaN.

Returns:

A tuple containing: - correlation_coefficients: Pearson correlation coefficients with shape (lat, lon) - p_values: The p-values of the correlations with shape (lat, lon)

Return type:

Tuple[np.ndarray, np.ndarray]

Raises:

ValueError – If the time dimension of data doesn’t match the length of the predictor.

pearson_correlation(x, y)[source]#

Calculate Pearson correlation coefficient.

Parameters:

  • x (array-like) – Input data arrays.

  • y (array-like) – Input data arrays.

Returns:

Pearson correlation coefficient.

Return type:

float

spearman_correlation(x, y)[source]#

Calculate Spearman rank correlation coefficient.

Parameters:

  • x (array-like) – Input data arrays.

  • y (array-like) – Input data arrays.

Returns:

Spearman correlation coefficient.

Return type:

float

kendall_correlation(x, y)[source]#

Calculate Kendall’s tau correlation coefficient.

Parameters:

  • x (array-like) – Input data arrays.

  • y (array-like) – Input data arrays.

Returns:

Kendall’s tau correlation coefficient.

Return type:

float

Emergent Constraint Functions#

gaussian_pdf(mu, sigma, x)[source]#

Calculate Gaussian probability density function.

Parameters:

  • mu (float) – Mean of the Gaussian distribution.

  • sigma (float) – Standard deviation of the Gaussian distribution.

  • x (Union[np.ndarray, float]) – Input values at which to evaluate the PDF.

Returns:

Probability density function values.

Return type:

Union[np.ndarray, float]

References

Adapted from: blackcata/Emergent_Constraints

emergent_constraint_posterior(constraint_data, target_data, constraint_grid, target_grid, obs_pdf)[source]#

Calculate posterior PDF using emergent constraint method.

This function applies the emergent constraint method to reduce uncertainty in climate projections by utilizing observational constraints.

Parameters:

  • constraint_data (xr.DataArray) – Inter-model spread data for the constraint variable (e.g., model sensitivity).

  • target_data (xr.DataArray) – Inter-model spread data for the target variable (e.g., future projection).

  • constraint_grid (np.ndarray) – Grid values for the constraint variable.

  • target_grid (np.ndarray) – Grid values for the target variable.

  • obs_pdf (np.ndarray) – Observational PDF of the constraint variable.

Returns:

A tuple containing: - posterior_pdf : np.ndarray - Posterior PDF of the target variable - posterior_std : float - Standard deviation of the target variable - posterior_mean : float - Mean of the target variable

Return type:

Tuple[np.ndarray, float, float]

References

Adapted from: blackcata/Emergent_Constraints Cox, P. M., et al. (2013). Nature, 494(7437), 341-344.

emergent_constraint_prior(constraint_data, target_data, constraint_grid, target_grid)[source]#

Calculate prior probability distribution for emergent constraints.

Parameters:

  • constraint_data (xr.DataArray) – Inter-model spread data for the constraint variable.

  • target_data (xr.DataArray) – Inter-model spread data for the target variable.

  • constraint_grid (np.ndarray) – Grid values for the constraint variable.

  • target_grid (np.ndarray) – Grid values for the target variable.

Returns:

A tuple containing: - prior_pdf : np.ndarray - Prior PDF - prediction_error : np.ndarray - Prediction error array - regression_line : np.ndarray - Linear regression values

Return type:

Tuple[np.ndarray, np.ndarray, np.ndarray]

References

Adapted from: blackcata/Emergent_Constraints

Legacy Functions (for backward compatibility)#

calc_GAUSSIAN_PDF(mu, sigma, x)[source]#

Legacy function name. Use gaussian_pdf() instead.

Parameters:
Return type:

ndarray | float

calc_PDF_EC(tmp_x, tmp_y, x, y, PDF_x)[source]#

Legacy function name. Use emergent_constraint_posterior() instead.

Parameters:
Return type:

Tuple[ndarray, float, float]

find_std_from_PDF(thres, y, PDF)[source]#

Legacy function name. Use _calculate_std_from_pdf() instead.

Parameters:
Return type:

float

calc_PDF_EC_PRIOR(tmp_x, tmp_y, x, y)[source]#

Legacy function name. Use emergent_constraint_prior() instead.

Parameters:
Return type:

Tuple[ndarray, ndarray, ndarray]

Utility Functions#

calculate_potential_temperature(temperature, pressure, reference_pressure=1000.0)[source]#

Calculate potential temperature using fast numpy operations.

This implementation uses lazy imports and avoids heavy metpy dependencies for simple potential temperature calculations.

Parameters:

  • temperature (array-like) – Temperature values in Kelvin.

  • pressure (array-like) – Pressure values in hPa.

  • reference_pressure (float, optional) – Reference pressure in hPa, default is 1000.0.

Returns:

Potential temperature values in Kelvin.

Return type:

array-like

Notes

Uses the standard formula: theta = T * (P0/P)^(R/cp) where R/cp = 0.286 for dry air

convert_longitude_range(data, lon='lon', center_on_180=True)[source]#

Wrap longitude coordinates of DataArray or Dataset to either -180..179 or 0..359.

Parameters:

  • data (xr.DataArray or xr.Dataset) – An xarray DataArray or Dataset object containing longitude coordinates.

  • lon (str, optional) – The name of the longitude coordinate, default is ‘lon’.

  • center_on_180 (bool, optional) – If True, wrap longitude from 0..359 to -180..179; If False, wrap longitude from -180..179 to 0..359.

Returns:

The DataArray or Dataset with wrapped longitude coordinates.

Return type:

xr.DataArray or xr.Dataset

Example Usage#

Tropopause Calculation Examples

import skyborn as skb
import numpy as np
import xarray as xr

# === 1D Profile Analysis ===
# Single atmospheric profile
pressure = np.array([50, 100, 200, 300, 500, 700, 850, 1000])  # hPa
temperature = np.array([200, 210, 220, 230, 250, 270, 280, 288])  # K

result = skb.calc.trop_wmo_profile(temperature, pressure)
print(f"Tropopause: {result['pressure']:.1f} hPa at {result['height']:.0f} m")

# === xarray Interface - Simplified ===
# Load atmospheric data with level coordinate
ds = xr.open_dataset('era5_data.nc')  # Temperature with level coordinate in hPa

# Auto-generate pressure from level coordinate - no pressure array needed!
result = skb.calc.troposphere.xarray.trop_wmo(ds.temperature)
print(f"Global tropopause calculated for {result.pressure.shape}")

# === 2D Cross-section Analysis ===
# Meridional cross-section (level, lat)
temp_meridional = ds.temperature.isel(time=0, lon=0)  # (level, lat)
result_2d = skb.calc.troposphere.xarray.trop_wmo(temp_meridional)
# Result shape: (lat,) - tropopause height at each latitude

# === 4D Climate Analysis ===
# Multi-year dataset (time, level, lat, lon)
result_4d = skb.calc.troposphere.xarray.trop_wmo(ds.temperature)
# Result shape: (time, lat, lon) - preserves time and spatial dimensions

# Seasonal analysis
seasonal_mean = result_4d.height.groupby('time.season').mean()

# === Advanced Usage ===
# Custom pressure field and WMO criterion
result = skb.calc.troposphere.xarray.trop_wmo(
    temperature=ds.temperature,
    pressure=ds.pressure,  # Custom pressure field
    lapse_criterion=2.5,   # Custom WMO threshold (K/km)
    pressure_unit='Pa'     # If pressure is in Pascals
)

Geostrophic Wind Examples

import skyborn as skb
import numpy as np
import xarray as xr

# === NumPy Interface ===
# Traditional interface with manual parameter specification
nlat, nlon = 73, 144
z_data = np.random.randn(nlat, nlon) * 100 + 5500  # Geopotential height [gpm]
lat = np.linspace(-90, 90, nlat)  # Degrees north
lon = np.linspace(0, 360, nlon)[:-1]  # Degrees east (0-357.5)

# Calculate geostrophic wind components
ug, vg = skb.calc.geostrophic_wind(z_data, lon, lat, 'yx')
print(f"Wind components: ug{ug.shape}, vg{vg.shape}")

# Class interface for derived quantities
gw = skb.calc.GeostrophicWind(z_data, lon, lat, 'yx')
speed = gw.speed()
print(f"Max wind speed: {speed.max():.1f} m/s")

# === xarray Interface - Simplified ===
# Load geopotential height data
z = xr.DataArray(
    z_data,
    dims=['lat', 'lon'],
    coords={
        'lat': (['lat'], lat, {'units': 'degrees_north'}),
        'lon': (['lon'], lon, {'units': 'degrees_east'})
    },
    attrs={'long_name': '500 hPa geopotential height', 'units': 'gpm'}
)

# Automatic coordinate detection and parameter inference
from skyborn.calc.geostrophic.xarray import geostrophic_wind
result = geostrophic_wind(z)  # That's it!

print(f"Automatic features:")
print(f"  Longitude cyclic: {result.attrs['longitude_cyclic']}")
print(f"  Latitude ordering: {result.attrs['latitude_ordering']}")
print(f"  Output: ug{result.ug.shape}, vg{result.vg.shape}")

# === Multi-dimensional Examples ===
# 3D time series (time, lat, lon)
nt = 12
z_3d_data = np.random.randn(nt, nlat, nlon) * 100 + 5500
z_3d = xr.DataArray(
    z_3d_data,
    dims=['time', 'lat', 'lon'],
    coords={
        'time': pd.date_range('2023-01-01', periods=nt, freq='MS'),
        'lat': lat, 'lon': lon
    }
)

result_3d = geostrophic_wind(z_3d)
print(f"3D result: ug{result_3d.ug.shape}, vg{result_3d.vg.shape}")

# Seasonal analysis
seasonal_winds = result_3d.ug.groupby('time.season').mean()

# 4D multi-level data (time, level, lat, lon)
nz = 17
levels = [50, 100, 150, 200, 250, 300, 400, 500, 600, 700, 775, 850, 925, 1000]
z_4d_data = np.random.randn(nt, len(levels), nlat, nlon) * 100 + 5500
z_4d = xr.DataArray(
    z_4d_data,
    dims=['time', 'level', 'lat', 'lon'],
    coords={
        'time': pd.date_range('2023-01-01', periods=nt, freq='MS'),
        'level': (['level'], levels, {'units': 'hPa'}),
        'lat': lat, 'lon': lon
    }
)

result_4d = geostrophic_wind(z_4d)
print(f"4D result: ug{result_4d.ug.shape}, vg{result_4d.vg.shape}")

# Class interface with xarray
from skyborn.calc.geostrophic.xarray import GeostrophicWind
gw_xr = GeostrophicWind(z)
speed_xr = gw_xr.speed()
print(f"xarray class speed: {speed_xr.attrs['standard_name']}")

Tropical Cyclone Potential Intensity (GPI) Examples

import skyborn as skb
import numpy as np
import xarray as xr

# === Single Profile Calculation ===
# Atmospheric profile data
pressure_levels = np.array([1000, 850, 700, 500, 300, 200, 100])  # mb
temperature = np.array([300, 290, 280, 270, 250, 230, 210])  # K
mixing_ratio = np.array([0.015, 0.010, 0.005, 0.002, 0.0005, 0.0001, 0.00005])  # kg/kg

# Surface conditions
sst = 302.0  # K (29°C - warm tropical ocean)
psl = 101325.0  # Pa (sea level pressure)

# NumPy interface
from skyborn.calc.GPI.core import potential_intensity
min_p, pi, err = potential_intensity(
    sst, psl, pressure_levels, temperature, mixing_ratio
)
print(f"Potential Intensity: {pi:.1f} m/s, Min Pressure: {min_p:.1f} mb")
print(f"Error flag: {err} (1 = success)")

# === xarray Interface with Automatic Unit Conversion ===
from skyborn.calc.GPI.xarray import potential_intensity

# Create xarray data with various units (will be auto-converted)
temp_xr = xr.DataArray(
    temperature - 273.15,  # Convert to Celsius for demo
    dims=['level'],
    attrs={'units': '°C'}  # Will be auto-converted to K
)

mixr_xr = xr.DataArray(
    mixing_ratio * 1000,  # Convert to g/kg for demo
    dims=['level'],
    attrs={'units': 'g/kg'}  # Will be auto-converted to kg/kg
)

pres_xr = xr.DataArray(
    pressure_levels,
    dims=['level'],
    attrs={'units': 'mb'}
)

# Automatic unit conversion and dimension detection
result = potential_intensity(sst, psl, pres_xr, temp_xr, mixr_xr)
print(f"xarray PI: {result.pi.values:.1f} m/s")
print(f"Minimum pressure: {result.min_pressure.values:.1f} mb")

# === 3D Gridded Data (Spatial Analysis) ===
nlat, nlon = 20, 30
nlevs = len(pressure_levels)

# Create 3D atmospheric data
temp_3d = xr.DataArray(
    np.random.randn(nlevs, nlat, nlon) * 5 + 280,
    dims=['level', 'lat', 'lon'],
    coords={
        'level': pressure_levels,
        'lat': np.linspace(-30, 30, nlat),
        'lon': np.linspace(0, 360, nlon, endpoint=False)
    },
    attrs={'units': 'K', 'long_name': 'Temperature'}
)

mixr_3d = xr.DataArray(
    np.random.rand(nlevs, nlat, nlon) * 0.01 + 0.005,
    dims=['level', 'lat', 'lon'],
    coords=temp_3d.coords,
    attrs={'units': 'kg/kg', 'long_name': 'Water vapor mixing ratio'}
)

sst_3d = xr.DataArray(
    np.random.randn(nlat, nlon) * 2 + 300,
    dims=['lat', 'lon'],
    coords={'lat': temp_3d.lat, 'lon': temp_3d.lon},
    attrs={'units': 'K', 'long_name': 'Sea surface temperature'}
)

psl_3d = xr.DataArray(
    np.random.randn(nlat, nlon) * 500 + 101325,
    dims=['lat', 'lon'],
    coords={'lat': temp_3d.lat, 'lon': temp_3d.lon},
    attrs={'units': 'Pa', 'long_name': 'Sea level pressure'}
)

# Calculate potential intensity for entire grid
result_3d = potential_intensity(sst_3d, psl_3d, pres_xr, temp_3d, mixr_3d)
print(f"3D PI shape: {result_3d.pi.shape}")
print(f"Max PI: {result_3d.pi.max().values:.1f} m/s")
print(f"Min central pressure: {result_3d.min_pressure.min().values:.1f} mb")

# Analyze tropical regions only
tropical_mask = np.abs(result_3d.lat) <= 25
tropical_pi = result_3d.pi.where(tropical_mask)
print(f"Mean tropical PI: {tropical_pi.mean().values:.1f} m/s")

# === 4D Time Series Analysis ===
ntimes = 12  # Monthly data

# Create 4D data (time, level, lat, lon)
temp_4d = xr.DataArray(
    np.random.randn(ntimes, nlevs, nlat, nlon) * 5 + 280,
    dims=['time', 'level', 'lat', 'lon'],
    coords={
        'time': pd.date_range('2023-01', periods=ntimes, freq='MS'),
        'level': pressure_levels,
        'lat': np.linspace(-30, 30, nlat),
        'lon': np.linspace(0, 360, nlon, endpoint=False)
    },
    attrs={'units': 'K'}
)

mixr_4d = xr.DataArray(
    np.random.rand(ntimes, nlevs, nlat, nlon) * 0.01 + 0.005,
    dims=['time', 'level', 'lat', 'lon'],
    coords=temp_4d.coords,
    attrs={'units': 'kg/kg'}
)

sst_4d = xr.DataArray(
    np.random.randn(ntimes, nlat, nlon) * 2 + 300,
    dims=['time', 'lat', 'lon'],
    coords={'time': temp_4d.time, 'lat': temp_4d.lat, 'lon': temp_4d.lon},
    attrs={'units': 'K'}
)

psl_4d = xr.DataArray(
    np.random.randn(ntimes, nlat, nlon) * 500 + 101325,
    dims=['time', 'lat', 'lon'],
    coords={'time': temp_4d.time, 'lat': temp_4d.lat, 'lon': temp_4d.lon},
    attrs={'units': 'Pa'}
)

# Calculate monthly potential intensity
result_4d = potential_intensity(sst_4d, psl_4d, pres_xr, temp_4d, mixr_4d)
print(f"4D result shape: {result_4d.pi.shape}")  # (time, lat, lon)

# Seasonal analysis
seasonal_pi = result_4d.pi.groupby('time.season').mean()
print(f"Summer mean PI: {seasonal_pi.sel(season='JJA').mean().values:.1f} m/s")
print(f"Winter mean PI: {seasonal_pi.sel(season='DJF').mean().values:.1f} m/s")

# === PI log decomposition diagnostics ===
sst_profile = xr.DataArray(sst, attrs={'units': 'K'})
psl_profile = xr.DataArray(psl, attrs={'units': 'Pa'})
temp_profile_xr = xr.DataArray(temperature, dims=['level'], coords={'level': pressure_levels}, attrs={'units': 'K'})
mixr_profile_xr = xr.DataArray(mixing_ratio, dims=['level'], coords={'level': pressure_levels}, attrs={'units': 'kg/kg'})
pres_profile_xr = xr.DataArray(pressure_levels, dims=['level'], attrs={'units': 'mb'})

from skyborn.calc.GPI.xarray import pi_log_decomposition

profile_full = pi_log_decomposition(
    sst_profile,
    psl_profile,
    pres_profile_xr,
    temp_profile_xr,
    mixr_profile_xr,
    outflow_source='cape_env'
)
print(profile_full[['pi', 't0', 'otl', 'lnpi', 'lneff', 'lndiseq', 'lnCKCD']])

Statistical Analysis Examples

# Statistical analysis
x = np.random.randn(100)
y = 2 * x + np.random.randn(100) * 0.5

slope, p_value = skb.linear_regression(x, y)
print(f"Linear regression: slope={slope:.4f}, p_value={p_value:.6f}")

correlation = skb.pearson_correlation(x, y)
print(f"Pearson correlation: {correlation:.4f}")

# Spatial correlation analysis
# Create sample spatial data (time, lat, lon)
n_time, n_lat, n_lon = 120, 36, 72
spatial_data = np.random.randn(n_time, n_lat, n_lon)
time_series = np.random.randn(n_time)

# Calculate spatial correlations efficiently
corr_map, p_values = skb.spatial_correlation(spatial_data, time_series)

# Works with xarray too
data_xr = xr.DataArray(spatial_data, dims=['time', 'lat', 'lon'])
predictor_xr = xr.DataArray(time_series, dims=['time'])
corr_xr, p_xr = skb.spatial_correlation(data_xr, predictor_xr)

# Emergent constraint analysis
x_values = np.linspace(-3, 3, 100)
pdf = skb.gaussian_pdf(mu=0, sigma=1, x=x_values)

# Apply emergent constraint
posterior_mean, posterior_std = skb.emergent_constraint_posterior(
    prior_mean=3.0, prior_std=1.5,
    obs_mean=0.5, obs_std=0.2,
    relationship_slope=2.0, relationship_intercept=0.1
)

Emergent Constraints#

The emergent constraint module implements methods for reducing uncertainty in climate projections by leveraging relationships between observable present-day quantities and uncertain future projections.

For a complete interactive example, see Emergent Constraints Analysis.