w22 package

Submodules

w22.constants module

CLE15 python package constants.

w22.plot module

Plot the new potential size calculation results for PIPS chapter/paper

w22.plot.dataframe_to_latex_table(df)

Converts a pandas DataFrame to a publication-quality LaTeX table string.

This function automatically generates formal LaTeX headers and formats numerical values into a human-readable scientific notation suitable for academic papers.

Key features: - Dynamic Header Generation: Parses column names like ‘rho_vmax’ or

‘fit_r0_sst’ into LaTeX expressions, e.g., ‘(rho(V_p))’.

  • Advanced Number Formatting:
    • Single numbers are formatted as (m times 10^{e}).

    • Value ± error pairs are formatted as ((textit{m}_n pm textit{m}_e) times 10^{E}), where the exponent is factored out and values are rounded systematically based on the error’s magnitude.

  • Custom Table Style: Uses ‘topline’, ‘midline’, ‘bottomline’ for table rules.

Parameters:

df (pd.DataFrame) – The input DataFrame. Column names are expected to follow conventions like ‘rho_{…}’ and ‘fit_{…}’.

Returns:

A string containing the fully formatted LaTeX table.

Return type:

str

Doctest:
>>> # Example DataFrame mimicking a typical analysis output
>>> data = {
...     'place': ['Atlantic'],
...     'member': ['member_01'],
...     'rho_vmax': [0.891],
...     'fit_vmax': [43.2],
...     'fit_vmax_err': [5.73],
...     'fit_r0_sst': [-1.54],
...     'fit_r0_sst_err': [0.11],
... }
>>> df_test = pd.DataFrame(data)
w22.plot.figure_two(place='new_orleans')

Plot the solution for the GOM bbox for potential size and intensity.

Then plot New Orleans timeseries.

Return type:

None

w22.plot.get_cmip6_timeseries(place='new_orleans', pressure_assumption='isothermal', model='CESM2', member=4, pi_version=4)

Get the CMIP6 timeseries dataset for the given place and potential intensity version.

Parameters:
  • place (str) – The place to get the timeseries for (default is “new_orleans”).

  • pressure_assumption (str) – The pressure assumption to use (default is “isothermal”).

  • member (int) – The ensemble member to use (default is 4).

  • pi_version (int) – The potential intensity version to use (default is 4).

Returns:

The timeseries dataset.

Return type:

xr.Dataset

w22.plot.get_era5_timeseries(place='new_orleans', pi_version=4)

Get the ERA5 timeseries dataset for the given place and potential intensity version.

Parameters:
  • place (str) – The place to get the timeseries for (default is “new_orleans”).

  • pi_version (int) – The potential intensity version to use (default is 4).

Returns:

The timeseries dataset.

Return type:

xr.Dataset

w22.plot.get_timeseries(model='CESM2', place='new_orleans', pi_version=4, pressure_assumption='isothermal', member=4)
Return type:

Dataset

w22.plot.multipanel(place='new_orleans', vars=('vmax_3', 'r0_3', 'rmax_3', 'rmax_1'), models={'CESM2', 'ERA5', 'HADGEM3-GC31-MM', 'MIROC6'})

Plot the multipanel figure for the CMIP6 timeseries.

Parameters:

vars (tuple) – The variables to plot (default is (“vmax_3”, “r0_3”, “rmax_3”, “rmax_1”)).

w22.plot.place_to_position(place)

Convert a place name to a latitude and longitude position.

Parameters:

place (str) – The name of the place.

Returns:

A tuple containing the latitude and longitude.

Return type:

tuple

w22.plot.plot_era5_timeseries(axs, place='new_orleans', color='black', pi_version=4)

Plot the potential intensity and size timeseries for the point near New Orleans for ERA5 data.

Parameters:
  • axs (np.ndarray) – The axes to plot on.

  • color (str) – The color of the line (default is “black”).

  • place (str) – The place to plot (default is “new_orleans”).

  • pi_version (int) – The potential intensity version to use (default is 4).

Return type:

None

w22.plot.plot_panels()
Return type:

None

w22.plot.plot_seasonal_profiles()

Plot the seasonal profiles for New Orleans.

w22.plot.plot_spatials(axs, place='new_orleans', vars=('vmax_3', 'r0_3', 'rmax_3', 'rmax_1'), labels=('Potential intensity, $V_{p}$ [m s$^{-1}$]', 'PI potential outer size, $r_{a3}$ [km]', 'PI potential inner size, $r_{3}$ [km]', 'Cat1 potential inner size, $r_{1}$ [km]'), pi_version=4, trial=1, pressure_assumption='isothermal', model='HADGEM3-GC31-MM', member='r1i1p1f3')
Return type:

None

w22.plot.plot_timeserii(axs, vars=('vmax_3', 'r0_3', 'rmax_3', 'rmax_1'), labels=('Potential intensity, $V_{p}$ [m s$^{-1}$]', 'PI potential outer size, $r_{a3}$ [km]', 'PI potential inner size, $r_{3}$ [km]', 'Cat1 potential inner size, $r_{1}$ [km]'), linewidth=0.5, color='black', member=4, model='CESM2', pressure_assumption='isothermal', place='new_orleans', pi_version=4, year_min=2014, year_max=2100)
Return type:

None

w22.plot.plot_two_spatial(axs, place='new_orleans', pi_version=4, trial=1, pressure_assumption='isothermal')

Plot the potential intensity and size for the Gulf of Mexico area.

Parameters:

axs (np.ndarray) – axes to plot on. Should be two.

Return type:

None

w22.plot.plot_two_timeseries(axs, text=True, color='black', member=4, model='CESM2', pressure_assumption='isothermal', place='new_orleans', pi_version=4, year_min=2014, year_max=2100)

Plot the potential intensity and size timeseries for the point near New Orleans for different ensemble members.

Parameters:
  • axs (np.ndarray) – The axes to plot on.

  • text (bool) – Whether to add text to the plot (default is True).

  • color (str) – The color of the line (default is “black”).

  • member (int) – The ensemble member to plot (default is 4).

Return type:

None

w22.plot.safe_corr(xt, yt)

Calculate the correlation of the data with error handling.

Parameters:
  • xt (np.ndarray) – The x data.

  • yt (np.ndarray) – The y data.

Returns:

The correlation.

Return type:

float

w22.plot.safe_grad(xt, yt)

Calculate the gradient of the data with error handling.

Parameters:
  • xt (np.ndarray) – The x data.

  • yt (np.ndarray) – The y data.

Returns:

The gradient.

Return type:

ufloat

w22.plot.temporal_relationship_data(place='new_orleans', pi_version=4)

Get the temporal relationships data for the given place and potential intensity version.

Parameters:
  • place (str) – The place to get the data for (default is “new_orleans”).

  • pi_version (int) – The potential intensity version to use (default is 4).

Return type:

None

w22.plot.timeseries_plot(name='new_orleans', plot_name='New Orleans', years=[2025, 2097], members=[4, 10, 11])

Plot the timeseries for the given location and members.

Parameters:
  • name (str) – The name of the location (default is “new_orleans”).

  • plot_name (str) – The name of the location to use in the plot title (default is “New Orleans”).

  • years (List[int]) – The years to plot (default is [2025, 2097]).

  • members (List[int]) – The members to plot (default is [4, 10, 11]).

Return type:

None

w22.plot.timeseries_relationships(timeseries_ds, place, member, year_min=2014, year_max=2100)
Return type:

DataFrame

w22.ps module

Functions to calculate potential size from xarray inputs in parallel.

TODO: this seems like a unit nightmare, perhaps we should have a unit checking function, and also make sure that everything is in SI units rather than some in hPa and some in Pa, some in kelvin, some in celsius.

w22.ps.calculate_ps13_ufunc(vmax_1, vmax_3, msl, sst, t0, lat, rh, ck_cd, cd, w_cool, supergradient_factor, pressure_assumption='isothermal')

Core computation function designed for xr.apply_ufunc.

This calculates both the potential size assuming the intensity is the potential intensity (vmax) and the potential size assuming the intensity is the minimum required for category 1 hurricane (normally 33 m s-1 at 10m height) (vmin).

TODO: This function is not yet implemented.

Parameters:
  • vmax_1 (float) – Minimum intensity for category 1 hurricane in m/s. Assume it is as the gradient wind level.

  • vmax_3 (float) – Potential intensity in m/s. Assume it is as the gradient wind level.

  • msl (float) – Ambient surface pressure in mbar.

  • sst (float) – Sea surface temperature in degC.

  • t0 (float) – Outflow temperature in degK.

  • lat (float) – Latitude in degrees North.

  • rh (float) – Environmental relative humidity (0-1).

  • ck_cd (float) – Ratio of exchange coefficients [dimensionless].

  • cd (float) – Drag coefficient [dimensionless].

  • w_cool (float) – Cooling coefficient in m/s.

  • supergradient_factor (float) – Supergradient factor [dimensionless].

  • pressure_assumption (str, optional) – Assumption for pressure calculation. Defaults to “isothermal”. Alternative is “isopycnal”.

Returns:

(r0_1, pm_1, pc_1, rmax_1, r0_3, pm_3, pc_3, rmax_3) where

r0_1 (float): Potential size for category 1 hurricane in m. pm_1 (float): Pressure at maximum winds for category 1 hurricane in Pa. pc_1 (float): Central pressure for CLE15 profile for category 1 hurricane in Pa. rmax_1 (float): Radius of maximum winds for category 1 hurricane in m. r0_3 (float): Potential size for potential intensity in m. pm_3 (float): Pressure at maximum winds for potential intensity in Pa. pc_3 (float): Central pressure for CLE15 profile for potential intensity in Pa. rmax_3 (float): Radius of maximum winds for potential intensity in m.

Return type:

tuple

Examples

>>> _ = calculate_ps13_ufunc(
...    vmax_1=33.0,
...    vmax_3=50.0,
...    msl=1015.0,
...    sst=29.0,
...    t0=200.0,
...    lat=20.0,
...    rh=0.9,
...    ck_cd=0.9,
...    cd=0.0015,
...    w_cool=0.002,
...    supergradient_factor=1.2,
...    )  
w22.ps.calculate_ps_ufunc(vmax, msl, sst, t0, lat, rh, ck_cd, cd, w_cool, supergradient_factor, pressure_assumption='isothermal')

Core computation function designed for xr.apply_ufunc. Accepts scalar NumPy values and returns a tuple of scalar results.

Parameters:
  • vmax (float) – Potential intensity in m/s. Assume it is as the gradient wind level.

  • msl (float) – Ambient surface pressure in mbar.

  • sst (float) – Sea surface temperature in degC.

  • t0 (float) – Outflow temperature in degK.

  • lat (float) – Latitude in degrees North.

  • rh (float) – Environmental relative humidity (0-1).

  • ck_cd (float) – Ratio of exchange coefficients [dimensionless].

  • cd (float) – Drag coefficient [dimensionless].

  • w_cool (float) – Cooling coefficient in m/s.

  • supergradient_factor (float) – Supergradient factor [dimensionless].

  • pressure_assumption (str, optional) – Assumption for pressure calculation. Defaults to “isothermal”. Alternative is “isopycnal”.

Returns:

(r0, pm, pc, rmax) where

r0 (float): Potential size in m. pm (float): Pressure at maximum winds in Pa. pc (float): Central pressure for CLE15 profile in Pa. rmax (float): Radius of maximum winds in m.

Return type:

tuple

w22.ps.create_test_dataset(shape, dim_names)

Helper function to generate a generic test dataset.

w22.ps.multi_point_example_1d()

Example of a multi-point solution for potential size.

Return type:

None

w22.ps.multi_point_example_2d(autofail=True)

Example of a multi-point solution for potential size. This function creates a dataset with multiple points and applies the parallelized_ps function to it.

Parameters:

autofail (bool, optional) – If True, will fail on any error. Defaults to True.

Return type:

None

w22.ps.parallelized_ps(ds, jobs=-1, pressure_assumption='isothermal', dryrun=False, autofail=False)

Apply point solution to all of the points in the dataset, using joblib to paralelize.

Parameters:
  • ds (xr.Dataset) – contains msl, vmax, sst, t0, rh, and lat.

  • jobs (int, optional) – Number of threads, defaults to -1.

  • pressure_assumption (str, optional) – Assumption for pressure calculation. Defaults to “isothermal”. Alternative is “isopycnal”.

  • dryrun (bool, optional) – If True, perform a dry run without actual calculations. Defaults to False.

  • autofail (bool, optional) – If True, fail on any error. Defaults to False.

Returns:

additionally contains r0, pm, pc, and rmax.

Return type:

xr.Dataset

Example

>>> in_ds = xr.Dataset(data_vars={
...     "msl": (("y", "x"), [[1016.7, 1016.7], [1016.7, 1016.7]]),  # mbar or hPa
...     "vmax": (("y", "x"), [[50, 51], [49, 49.5]]),  # m/s, potential intensity
...     "sst": (("y", "x"), [[29, 30], [28, 28]]),  # degC
...     "t0": (("y", "x"), [[200, 200], [200, 200]]),  # degK
...     "rh": (("y", "x"), [[0.9, 0.9], [0.9, 0.9]]),  # [dimensionless], relative humidity
...     "ck_cd": (("y", "x"), [[0.95, 0.95], [0.95, 0.95]]),  # [dimensionless], ck_cd
...     "cd": (("y", "x"), [[0.0015, 0.0015], [0.0015, 0.0015]]),  # [dimensionless], cd
...     "w_cool": (("y", "x"), [[0.002, 0.002], [0.002, 0.002]]),  # mbar or hPa
...     "supergradient_factor": (("y", "x"), [[1.2, 1.2], [1.2, 1.2]]),  # mbar or hPa
...     "pressure_assumption": (("y", "x"), [["isothermal", "isothermal"], ["isothermal", "isothermal"]]),  # mbar or hPa
...     },
...     coords={"lat": (("y", "x"), [[30, 30], [25, 25]])},  # degNorth
... )
>>> out_ds = parallelized_ps(in_ds, jobs=1, autofail=True)  
About to conduct 4 jobs in parallel
Automatically failed in autofail mode.
Automatically failed in autofail mode.
Automatically failed in autofail mode.
Automatically failed in autofail mode.
'parallelized_ps' ... s
w22.ps.parallelized_ps13_dask(ds)

Calculate potential size at category 1 intensity, and at the potential intensity.

Parameters:

ds (xr.Dataset) – Must constain: - “vmax_1”: Minimum intensity for category 1 hurricane in m/s. Assume it is as the gradient wind level. - “vmax_3”: Potential intensity in m/s. Assume it is as the gradient wind level. - “msl”: Ambient surface pressure in mbar. - “sst”: Sea surface temperature in degrees Celsius. - “t0”: Outflow temperature in degK. - “lat”: Latitude in degrees North.

Returns:

Additionally contains:
  • ”r0_1”: Potential outer size for category 1 hurricane in m.

  • ”pm_1”: Pressure at maximum winds for minimum category 1 hurricane intensity in Pa.

  • ”pc_1”: Central pressure for CLE15 profile for category 1 hurricane in Pa.

  • ”rmax_1”: Potential inner size for category 1 hurricane in m.

  • ”r0_3”: Potential outer size for potential intensity in m.

  • ”pm_3”: Pressure at maximum winds for the potential intensity in Pa.

  • ”pc_3”: Central pressure for CLE15 profile for potential intensity in Pa.

  • ”rmax_3”: Potential inner size for the potential intensity in m.

Return type:

xr.Dataset

w22.ps.parallelized_ps_dask(ds)

Apply point solution to all points using xr.apply_ufunc and Dask.

Parameters:

ds (xr.Dataset) – contains msl, vmax, sst, t0, rh, and lat.

Returns:

additionally contains r0, pm, pc, and rmax.

Return type:

xr.Dataset

Doctest::
>>> in_ds = xr.Dataset(
...    data_vars={
...        "msl": (("y", "x"), [[1015.0, 1016.0], [1012.0, 1014.0]]),
...        "vmax": (("y", "x"), [[50.0, 0.0], [45.0, 60.0]]), # Includes a zero vmax
...        "sst": (("y", "x"), [[29.0, 30.0], [28.0, 28.5]]),
...        "t0": (("y", "x"), [[200.0, 201.0], [199.0, 200.0]]),
...        "rh": (("y", "x"), [[0.9, 0.85], [0.92, 0.88]]),
...        "ck_cd": 0.9,
...        "cd": 0.0015,
...        "w_cool": 0.002,
...        "supergradient_factor": 1.2,
...    },
...    coords={"lat": (("y", "x"), [[30.0, 25.0], [20.0, 15.0]])},
... )
>>> result_joblib = parallelized_ps(in_ds.copy(), jobs=2, autofail=False) 
About to conduct 4 jobs in parallel
'parallelized_ps' ... s
>>> result_dask = parallelized_ps_dask(in_ds.copy()).compute() 
'parallelized_ps_dask' ... s
>>> output_vars = ["r0", "pm", "pc", "rmax"]
>>> print(result_joblib, result_dask) 
<xarray.Dataset> ... <xarray.Dataset> ...
>>> if "potential_size_calculated_at_time" in result_joblib: del result_joblib["potential_size_calculated_at_time"]
>>> if "potential_size_calculated_at_time" in result_dask: del result_dask["potential_size_calculated_at_time"]
>>> for var in output_vars:
...    xr.testing.assert_allclose(result_joblib[var], result_dask[var])
w22.ps.point_solution_ps(ds, supergradient_factor=1.2, include_profile=False, pressure_assumption='isothermal')

Find the solution for a given point in the grid.

Parameters:
  • ds (xr.Dataset) – Dataset with the input values.

  • supergradient_factor (float, optional) – Supergradient. Defaults to 1.2.

  • include_profile (bool, optional) – If True, include the full profile in the output dataset. Defaults to False.

  • pressure_assumption (str, optional) – Assumption for pressure calculation. Defaults to “isothermal”. Alternative is “isopycnal”.

Returns:

Find the solution dataset.

Return type:

xr.Dataset

Example

>>> in_ds = xr.Dataset(data_vars={
...     "msl": 1016.7, # mbar or hPa
...     "vmax": 49.5, # m/s, potential intensity
...     "sst": 28, # degC
...     "t0": 200, # degK
...     "rh": 0.9, # [dimensionless], relative humidity
...     },
...     coords={"lat":28})
>>> out_ds = point_solution_ps(in_ds) 
    'point_solution_ps' ... s
w22.ps.single_point_example()

Example of a single point solution for potential size. This function creates a dataset with a single point and applies the point_solution_ps function to it.

Example: >>> single_point_example() # doctest: +ELLIPSIS ‘point_solution_ps’ … s

Return type:

None

w22.ps.test_equivalence_across_dimensionalities(shape, dim_names)

Tests that Joblib and Dask outputs are structurally identical for 1D, 2D, and 3D inputs with arbitrary dimension names.

w22.ps_runs module

Run potential size calculations on CMIP6 data to make specific subsets (very expensive computation).

w22.ps_runs.ex_data_path(pi_version=2, member=4)

Get the path to the example data.

Parameters:
  • pi_version (int, optional) – pi_version of the pi calculation. Defaults to 2.

  • member (int, optional) – member number. Defaults to 4.

Returns:

path to the example data.

Return type:

str

w22.ps_runs.get_processed_data_for_point(place='new_orleans', member='r10i1p1f1', model='CESM2', exp='ssp585')

Get processed data for a specific point.

Parameters:
  • place (str, optional) – place to get data for. Defaults to “new_orleans”.

  • member (str, optional) – member number. Defaults to “r10i1p1f1”.

  • model (str, optional) – model name. Defaults to “CESM2”.

  • exp (str, optional) – experiment name. Defaults to “ssp585”.

w22.ps_runs.get_regional_processed_data(place='new_orleans')

Get processed data for a specific region.

Parameters:

place (str, optional) – place to get data for. Defaults to “new_orleans

Return type:

Dataset

w22.ps_runs.global_cmip6(part='nw', year=2015, pi_version=2)

Run potential size calculations on CMIP6 data to make specific subsets. :type part: :param part: segment of the world to calculate. Defaults to “nw”. :type part: str, optional :type year: int :param year: year to calculate. Defaults to 2015. :type year: int, optional :type pi_version: :param pi_version: pi_version of the data. Defaults to 2. :type pi_version: int, optional

Return type:

None

w22.ps_runs.load_global(year=2015)

Load all parts of the global data and stitch together.

Parameters:

year (int, optional) – year to load. Defaults to 2015.

Returns:

stitched together dataset.

Return type:

xr.Dataset

w22.ps_runs.new_orleans_10year(pi_version=2, pressure_assumption='isothermal')

Run potential size calculations on CMIP6 data to get New Orleans data.

Parameters:
  • pi_version (int, optional) – pi_version of the pi calculation. Defaults to 2.

  • pressure_assumption (str, optional) – pressure assumption. Defaults to “isothermal”.

Return type:

None

w22.ps_runs.point_era5_timeseries(place='new_orleans', pressure_assumption='isothermal')

Point timeseries for ERA5 data.

Parameters:
  • member (int, optional) – member number. Defaults to 10.

  • place (str, optional) – location. Defaults to “new_orleans”.

  • pressure_assumption (str, optional) – pressure assumption. Defaults to “isothermal”.

Return type:

None

w22.ps_runs.point_timeseries(member=10, model='CESM2', recalculate_pi=True, exp='ssp585', place='new_orleans', pressure_assumption='isothermal', pi_version=3)

Point timeseries.

Parameters:
  • member (int, optional) – member number. Defaults to 10.

  • model (str, optional) – model name. Defaults to “CESM2”.

  • recalculate_pi (bool, optional) – whether to recalculate pi. Defaults to True.

  • exp (str, optional) – experiment name. Defaults to “ssp585”.

  • pi_version (int, optional) – pi version. Defaults to 3.xxeexexe

  • place (str, optional) – location. Defaults to “new_orleans”.

  • pressure_assumption (str, optional) – pressure assumption. Defaults to “isothermal”.

Return type:

None

w22.ps_runs.spatial_example(pressure_assumption='isothermal', model='CESM2', trial=1, pi_version=2, member='r4i1p1f1', recalculate_pi=True, place='new_orleans')

Run potential size calculations on CMIP6 data to get a region.

Parameters:
  • pressure_assumption (str, optional) – pressure assumption. Defaults to “isothermal”.

  • model (str, optional) – model name. Defaults to “CESM2”.

  • member (str, optional) – member number. Defaults to “r4i1p1f1”.

  • trial (int, optional) – trial number. Defaults to 1.

  • pi_version (int, optional) – pi_version of the pi calculation. Defaults to 2.

  • recalculate_pi (bool, optional) – whether to recalculate pi. Defaults to True.

  • place (str, optional) – place to get data for. Defaults to “new_orleans

Return type:

None

w22.solve module

Solve module for numerical methods.

w22.solve.bisection(f, left, right, tol)

Bisection numerical method.

https://en.wikipedia.org/wiki/Root-finding_algorithms#Bisection_method

Parameters:
  • f (Callable) – Function to find root of.

  • left (float) – Left boundary.

  • right (float) – Right boundary.

  • tol (float) – Tolerance for convergence.

Returns:

x such that |f(x)| < tol.

Return type:

float

Example::
>>> f = lambda x: x - 2
>>> np.isclose(bisection(f, 1, 3, 1e-6), 2, rtol=1e-2, atol=1e-2)
True
>>> f = lambda x: x**2 - 4
>>> np.isclose(bisection(f, 0, 3, 1e-3), 2, rtol=1e-3, atol=1e-3)
True
>>> f = lambda x: - x**2 + 4
>>> np.isclose(bisection(f, 0, 3, 1e-3), 2, rtol=1e-3, atol=1e-3)
True

w22.stats module

Generate tables of results.

TODO: Add in ability to process HadGEM3-GC31-MM r1i1p1f3, r2i1p1f3, r3i1p1f3, and MIROC6 r1i1p1f1, r2i1p1f1, r3i1p1f1.

TODO: Add in tables to assess biases, variability etc. for the different models compared to ERA5.

TODO: Add in ability to process new variables and for it to be more extensible.
New variables include:
  • rmax_3 (should be labeled (r_{3}))

  • rmax_1 (should be labeled (r_{1}))

  • r0_1 (should be labeled (r_{a1}))

  • r0_3 (should be labeled (r_{a3}))

  • vmax_3 (should be labeled (V_{p}))

w22.stats.analyze_bias(ground_truth, model)

Analyzes the bias between a model and ground truth data.

Calculates the mean of the bias and the trend in the bias, and determines if they are statistically significant from zero.

Parameters:
  • ground_truth (ndarray) – A 1D numpy array of the true values.

  • model (ndarray) – A 1D numpy array of the modeled values.

Return type:

dict

Returns:

A dictionary containing the mean bias, the p-value for the mean bias, the trend of the bias (slope), and the p-value for the trend.

Doctests: >>> truth = np.linspace(0, 10, 10) >>> # Model has a constant bias of +2 and a slight upward trend in bias >>> model_data = truth + 2 + np.linspace(0, 0.5, 10) >>> results = analyze_bias(truth, model_data) >>> round(results[‘mean_bias’], 4) 2.25 >>> # p-value should be very small, indicating significance >>> results[‘mean_bias_p_value’] < 0.001 True >>> round(results[‘bias_trend’], 4) 0.0556 >>> # p-value should be very small, indicating significance >>> results[‘bias_trend_p_value’] < 0.001 True

w22.stats.calculate_detrended_cv(timeseries_data)

Calculates the detrended standard deviation and coefficient of variation.

This function performs a linear regression on the input time series to remove the trend. It then calculates the standard deviation of the residuals (the detrended data) and the coefficient of variation (CV) relative to the mean of the original, non-detrended data.

Parameters:

timeseries_data (ndarray) – A 1-dimensional numpy array of numerical data representing the time series.

Returns:

  • residuals (np.ndarray): The detrended time series.

  • std_residuals (float): The standard deviation of the residuals.

  • detrended_cv (float): The coefficient of variation of the detrended

    data, normalized by the original mean.

Return type:

A tuple containing

Doctests: >>> # Test with a perfect linear trend and no noise >>> perfect_series = np.arange(10) * 2 + 50 >>> res, std_res, cv = calculate_detrended_cv(perfect_series) >>> np.allclose(std_res, 0.0) True >>> np.allclose(cv, 0.0) True >>> # Test with some noise - CORRECTED VALUE >>> noisy_series = np.array([51, 51, 54, 55, 56, 59, 60, 62, 63, 65]) >>> res, std_res, cv = calculate_detrended_cv(noisy_series) >>> round(std_res, 4) 0.5865 >>> round(cv, 4) 0.0102

w22.stats.dataframe_to_latex_table(df)

Converts a pandas DataFrame to a publication-quality LaTeX table string.

This function automatically generates formal LaTeX headers and formats numerical values into a human-readable scientific notation suitable for academic papers.

Key features: - Dynamic Header Generation: Parses column names like ‘rho_vmax’ or

‘fit_r0_sst’ into LaTeX expressions, e.g., ‘(rho(V_p))’.

  • Advanced Number Formatting:
    • Single numbers are formatted as (m times 10^{e}).

    • Value ± error pairs are formatted as ((textit{m}_n pm textit{m}_e) times 10^{E}), where the exponent is factored out and values are rounded systematically based on the error’s magnitude.

  • Custom Table Style: Uses ‘topline’, ‘midline’, ‘bottomline’ for table rules.

Parameters:

df (pd.DataFrame) – The input DataFrame. Column names are expected to follow conventions like ‘rho_{…}’ and ‘fit_{…}’.

Returns:

A string containing the fully formatted LaTeX table.

Return type:

str

Doctest:
>>> # Example DataFrame mimicking a typical analysis output
>>> data = {
...     'place': ['Atlantic'],
...     'member': ['member_01'],
...     'rho_vmax': [0.891],
...     'fit_vmax': [43.2],
...     'fit_vmax_err': [5.73],
...     'fit_r0_sst': [-1.54],
...     'fit_r0_sst_err': [0.11],
... }
>>> df_test = pd.DataFrame(data)
w22.stats.format_error_latex_sci(nominal, error)

Formats a nominal ± error pair with a common exponent.

Parameters:
  • nominal (float) – The nominal value.

  • error (float) – The error value.

Returns:

The formatted string in LaTeX scientific notation.

Return type:

str

Doctest:
>>> print(format_error_latex_sci(12345, 67))
\(\left(1.234 \pm 0.007\right)\times 10^{4}\)
>>> print(format_error_latex_sci(0.0012345, float('inf'))) # previous behaviour\(1.2 \times 10^{-3}\)\( \pm \infty \)
--
>>> print(format_error_latex_sci(0, 0))
\(0.0\)
>>> print(format_error_latex_sci(float('nan'), 1))
--
w22.stats.format_single_latex_sci(value, sig_figs=2)

Formats a single number as m x 10^e.

Parameters:
  • value (float) – The value to format.

  • sig_figs (int) – The number of significant figures to use (default is 2 which gives 1 decimal place for the mantissa).

Returns:

The formatted string in LaTeX scientific notation.

Return type:

str

Doctest:
>>> print(format_single_latex_sci(12345))
\(1.2 \times 10^{4}\)
>>> print(format_single_latex_sci(0.0012345))
\(1.2 \times 10^{-3}\)
>>> print(format_single_latex_sci(0))
\(0.0\)
>>> print(format_single_latex_sci(float('inf')))
\(\infty\)
w22.stats.safe_corr(xt, yt)

Calculate the correlation of the data with error handling.

Parameters:
  • xt (np.ndarray) – The x data.

  • yt (np.ndarray) – The y data.

Returns:

The correlation.

Return type:

float

w22.stats.safe_grad(xt, yt)

Calculate the gradient of the data with error handling.

Parameters:
  • xt (np.ndarray) – The x data.

  • yt (np.ndarray) – The y data.

Returns:

The gradient.

Return type:

ufloat

w22.stats.temporal_relationship_data(place='new_orleans', pi_version=4)

Get the temporal relationships data for the given place and potential intensity version.

Parameters:
  • place (str) – The place to get the data for (default is “new_orleans”).

  • pi_version (int) – The potential intensity version to use (default is 4).

Return type:

None

w22.stats.timeseries_relationships(timeseries_ds, place, member, year_min=2014, year_max=2100)

Timeseries relationships function.

Parameters:
  • timeseries_ds (xr.Dataset) – _description_

  • place (str) – place for description.

  • member (int) – _description_

  • year_min (int, optional) – _description_. Defaults to 2014.

  • year_max (int, optional) – _description_. Defaults to 2100.

Returns:

_description_

Return type:

pd.DataFrame

w22.stats2 module

Generate tables of results.

w22.stats2.analyze_bias(ground_truth, model)

Analyzes the bias between a model and ground truth data.

Calculates the mean of the bias and the trend in the bias, and determines if they are statistically significant from zero.

Parameters:
  • ground_truth (ndarray) – A 1D numpy array of the true values.

  • model (ndarray) – A 1D numpy array of the modeled values.

Return type:

dict

Returns:

A dictionary containing the mean bias, the p-value for the mean bias, the trend of the bias (slope), and the p-value for the trend.

Doctests: >>> truth = np.linspace(0, 10, 10) >>> # Model has a constant bias of +2 and a slight upward trend in bias >>> model_data = truth + 2 + np.linspace(0, 0.5, 10) >>> results = analyze_bias(truth, model_data) >>> round(results[‘mean_bias’], 4) 2.25 >>> results[‘mean_bias_p_value’] < 0.001 True >>> round(results[‘bias_trend’], 4) 0.0556 >>> results[‘bias_trend_p_value’] < 0.001 True

w22.stats2.analyze_bias_mixed_model(model_df)

Analyzes bias using a linear mixed-effects model.

This treats ‘member’ as a random effect to account for non-independence, providing a robust estimate of the overall mean bias (fixed intercept) and the trend in bias (fixed slope for time).

Parameters:

model_df (DataFrame) – DataFrame with columns [‘time’, ‘bias’, ‘member’].

Return type:

dict

Returns:

A dictionary with the mean bias and trend, their standard errors, p-values, and the variance of the inter-member variability.

w22.stats2.analyze_bias_newey_west_ensemble(model_aligned_ds, era5_aligned_ds, var)

Analyzes bias using OLS with Newey-West standard errors on the ensemble mean.

This function first calculates the bias for a given variable, then computes the mean across the ensemble members, and finally calculates the trend of this ensemble-mean bias time series.

Parameters:
  • model_aligned_ds (xr.Dataset) – Dataset with multiple members, aligned to ERA5.

  • era5_aligned_ds (xr.Dataset) – ERA5 dataset, aligned to the model.

  • var (str) – The variable to analyze.

Return type:

dict

Returns:

A dictionary containing the mean bias, bias trend, and the trend’s p-value.

w22.stats2.calculate_detrended_cv(timeseries_data)

Calculates the detrended standard deviation and coefficient of variation.

This function performs a linear regression on the input time series to remove the trend. It then calculates the standard deviation of the residuals (the detrended data) and the coefficient of variation (CV) relative to the mean of the original, non-detrended data.

Parameters:

timeseries_data (ndarray) – A 1-dimensional numpy array of numerical data representing the time series.

Returns:

  • residuals (np.ndarray): The detrended time series.

  • std_residuals (float): The standard deviation of the residuals.

  • detrended_cv (float): The coefficient of variation of the detrended

    data, normalized by the original mean.

Return type:

A tuple containing

Doctests: >>> # Test with a perfect linear trend and no noise >>> perfect_series = np.arange(10) * 2 + 50 >>> res, std_res, cv = calculate_detrended_cv(perfect_series) >>> np.allclose(std_res, 0.0) True >>> np.allclose(cv, 0.0) True >>> # Test with some noise - CORRECTED VALUE >>> noisy_series = np.array([51, 51, 54, 55, 56, 59, 60, 62, 63, 65]) >>> res, std_res, cv = calculate_detrended_cv(noisy_series) >>> round(std_res, 4) 0.5865 >>> round(cv, 4) 0.0102

w22.stats2.dataframe_to_latex_table(df)

Converts a pandas DataFrame to a publication-quality LaTeX table string.

Return type:

str

w22.stats2.format_error_latex_sci(nominal, error)
Return type:

str

w22.stats2.format_single_latex_sci(value, sig_figs=2)
Return type:

str

w22.stats2.generate_assessment_tables(place='new_orleans', pi_version=4, comparison_period=(1980, 2014))

Generates tables to assess biases and variability for models compared to ERA5.

Parameters:
  • place (str) – The place for which to generate the assessment.

  • pi_version (int) – The potential intensity version to use.

  • comparison_period (Tuple[int, int]) – The (start, end) year for the comparison.

Return type:

None

w22.stats2.generate_future_trend_table(place='new_orleans', pi_version=4, future_period=(2014, 2100))

Generates a table of future trends for CMIP6 models.

This function calculates the ensemble-mean trend and its significance using Newey-West standard errors. It also calculates the range of trends across individual ensemble members to quantify internal variability.

Parameters:
  • place (str) – The place for which to generate the trend table.

  • pi_version (int) – The potential intensity version to use.

  • future_period (Tuple[int, int]) – The (start, end) year for the trend.

Return type:

None

w22.stats2.generate_mixed_model_assessment_tables(place='new_orleans', pi_version=4, comparison_period=(1980, 2024))

Generates assessment tables using a robust mixed-effects model approach.

Return type:

None

w22.stats2.generate_newey_west_assessment_tables(place='new_orleans', pi_version=4, comparison_period=(1980, 2024))

Generates assessment tables using an ensemble-mean Newey-West approach.

Return type:

None

w22.stats2.generate_wide_assessment_tables(place='new_orleans', pi_version=4, comparison_period=(1980, 2024))

Generates wide-format assessment tables with aggregated model statistics.

This function calculates bias, RMSE, and variability metrics for CMIP models against ERA5. It then aggregates results across ensemble members for each model and pivots the data to create publication-ready ‘wide’ format tables for mean state assessment and variability assessment.

Parameters:
  • place (str) – The place for which to generate the assessment.

  • pi_version (int) – The potential intensity version to use.

  • comparison_period (Tuple[int, int]) – The (start, end) year for comparison.

Return type:

None

w22.stats2.p_to_str(p_val)

Formats a p-value into a LaTeX string for inclusion in tables.

Parameters:

p_val (float) – The p-value to format.

Returns:

A LaTeX-formatted string, e.g., “($p=1.3 times 10^{-3}$)”.

Return type:

str

Doctests: >>> p_to_str(0.000126) ‘ ($p=1.3 \times 10^{-4}$)’ >>> p_to_str(0.045) ‘ ($p=0.045$)’ >>> p_to_str(0.2314) ‘ ($p=0.231$)’ >>> p_to_str(np.nan) ‘’

w22.stats2.safe_corr(xt, yt)

Calculate the correlation of the data with error handling.

Parameters:
  • xt (np.ndarray) – The x data.

  • yt (np.ndarray) – The y data.

Returns:

The correlation.

Return type:

float

w22.stats2.safe_grad(xt, yt)

Calculate the gradient of the data with error handling.

Parameters:
  • xt (np.ndarray) – The x data.

  • yt (np.ndarray) – The y data.

Returns:

The gradient.

Return type:

ufloat

w22.stats2.temporal_relationship_data(place='new_orleans', pi_version=4)

Generates and saves temporal relationship statistics for multiple models.

Parameters:
  • place (str) – The place to get the data for.

  • pi_version (int) – The potential intensity version to use.

Return type:

None

w22.stats2.timeseries_relationships(timeseries_ds, place, member, year_min, year_max, variables)

Calculates trend and correlation statistics for a set of timeseries variables.

Parameters:
  • timeseries_ds (xr.Dataset) – Dataset containing the time series.

  • place (str) – Name of the location.

  • member (str) – Identifier for the model member (e.g., ‘r1i1p1f1’).

  • year_min (int) – Start year of the analysis period.

  • year_max (int) – End year of the analysis period.

  • variables (List[str]) – A list of variable names to analyze.

Returns:

A single-row DataFrame with calculated statistics.

Return type:

pd.DataFrame

w22.stats2.trend_with_neweywest_full(y)

Calculates the linear trend (slope, intercept) and its p-value using OLS with Newey-West standard errors.

Return type:

tuple[float, float, float]

w22.stats2.wide_dataframe_to_latex(df, caption, label, filename)

Converts a wide-format pandas DataFrame to a publication-quality LaTeX table.

This function handles DataFrames with multi-level columns (metric, variable) and formats them into a clean LaTeX table with grouped headers.

Parameters:
  • df (pd.DataFrame) – The wide-format DataFrame to convert.

  • caption (str) – The LaTeX table caption.

  • label (str) – The LaTeX label for cross-referencing.

  • filename (str) – The full path to save the .tex file.

Return type:

None

w22.test module

Test solutions against Wang 2022 figure and paper data (make sure that everything works as expected)

w22.test.make_matlab_v_python_plot(add_name='')
w22.test.octave_vs_python(add_name='')
w22.test.q_gibbs(near_surface_air_temperature=299, dry_pressure_at_inflow=98500, dry_pressure_at_maximum_winds=94400, gas_constant=287)
w22.test.q_s(efficiency_relative_to_carnot=0.5, near_surface_air_temperature=299, outflow_temperature=200, dry_pressure_at_inflow=98500, dry_pressure_at_maximum_winds=94400, gas_constant=287, latent_heat_of_vaporization=2500000, gas_constant_for_water_vapor=461.5)
w22.test.test_figure_4()
w22.test.test_figure_5()
w22.test.v_carnot(efficiency_relative_to_carnot=0.5, near_surface_air_temperature=299, outflow_temperature=200, latent_heat_of_vaporization=2500000, gas_constant=287, gas_constant_for_water_vapor=461.5, pressure_dry_at_inflow=98500)

Calculate the carnot engine velocity approximated in Wang 2022.

Parameters:
  • efficiency_relative_to_carnot (float) – The efficiency relative to Carnot.

  • near_surface_air_temperature (float) – The near surface air temperature in Kelvin.

  • outflow_temperature (float) – The outflow temperature in Kelvin.

  • latent_heat_of_vaporization (float) – The latent heat of vaporization in J/kg

  • gas_constant (float) – The gas constant for dry air in J/(kg*K).

  • gas_constant_for_water_vapor (float) – The gas constant for water vapor in J/(kg*K).

  • pressure_dry_at_inflow (float) – The dry pressure at inflow in Pascals (Pa).

Returns:

The velocity at which the Carnot efficiency is achieved in m/s.

Return type:

float

w22.test.w_out(coriolis_parameter=5e-05, outer_radius=2193000, maximum_velocity=83, radius_of_maximum_winds=64000, radius_of_outflow=inf)
w22.test.w_p(coriolis_parameter=5e-05, gas_constant=287, beta_lift_parameterization=1.25, near_surface_air_temperature=299, outer_radius=2193000, maximum_velocity=83, dry_pressure_at_inflow=98500, radius_of_maximum_winds=64000, dry_pressure_at_maximum_winds=94400)
w22.test.w_pbl(gas_constant=287, near_surface_air_temperature=299, dry_pressure_at_maximum_winds=94400, dry_pressure_at_inflow=98500, maximum_velocity=83)

w22.utils module

Utilities for idealized tropical cyclone calculations.

w22.utils.absolute_angular_momentum(v, r, f)

Calculate absolute angular momentum.

Parameters:
  • v (float) – Azimuthal wind speed [m/s].

  • r (float) – Radius from storm centre [m].

  • f (float) – Coriolis parameter [s-1].

Returns:

Absolute angular momentum [m2/s].

Return type:

float

Example::
>>> absolute_angular_momentum(0, 0, 0)
0.0
>>> absolute_angular_momentum(1, 1, 1)
1.5
w22.utils.buck_sat_vap_pressure(temp)

Arden buck_sat_vap_pressure equation.

https://en.wikipedia.org/wiki/Arden_buck_sat_vap_pressure_equation

Parameters:

temp (float) – temperature in Kelvin.

Returns:

saturation vapour pressure in Pa.

Return type:

float

w22.utils.carnot_efficiency(temp_hot, temp_cold)

Calculate Carnot efficiency.

Parameters:
  • temp_hot (float) – Temperature of hot reservoir [K].

  • temp_cold (float) – Temperature of cold reservoir [K].

Returns:

Carnot efficiency [dimensionless].

Return type:

float

w22.utils.coriolis_parameter_from_lat(lat)

Calculate the Coriolis parameter from latitude.

Parameters:

lat (np.ndarray) – Latitude [degrees].

Returns:

Coriolis parameter [s-1].

Return type:

np.ndarray

Example::
>>> cp_out = coriolis_parameter_from_lat(30)
>>> cp_manual = 2 * 2 * np.pi / 24 / 60 / 60 * 1/2  ## 2 * omega * sin(30deg)
>>> np.isclose(cp_out, cp_manual, rtol=1e-3, atol=1e-6)
True
>>> np.isclose(coriolis_parameter_from_lat(0), 0, rtol=1e-3, atol=1e-6)
True
w22.utils.curveintersect(x1, y1, x2, y2)

Find intersection points of two curves (x1, y1) and (x2, y2).

Parameters: x1, y1 : array-like, coordinates of the first curve x2, y2 : array-like, coordinates of the second curve

Returns: Tuple[List[float], List[float]]: intersections x and y coordinates.

Return type:

Tuple[List[float], List[float]]

w22.utils.pressure_from_wind(rr, vv, p0=101600, rho0=1.225, fcor=5e-05, assumption='isopycnal')

Pressure from wind, assuming cyclogeostrophic balance and either constant temperature or pressure.

We assume cyclogeostrophic balance, with pressure gradient balancing the circular acceleration and coriolis force.

frac{dp}{dr}=rho(r)left(frac{v(r)^2}{r} +fv(r)right)

where p is pressure, r is radius, rho is density, f is coriolis parameter and v is velocity. If ideal gas following isothermal expansion then rho1/p1 = rho2/p2, rho(r) = rho0/p0*p(r)

Where rho0 and p0 are the background density and pressure respectively. Substituting in we find that,

frac{dp}{dr}=p(r) frac{rho0}{p0}left(frac{v(r)^2}{r} +fv(r)right)

and then bringing p(r) to the left hand side yields

frac{1}{p}frac{dp}{dr}=frac{dlog(p)}{dr} = frac{rho0}{p0}left(frac{v(r)^2}{r} +fv(r)right),

allowing us to express p(r) by integrating from infinity to r,

p(r) = p0*expleft(-frac{rho0}{p0}int^{inf}_{r}left(frac{v(r)^2}{r} +fv(r)right)drright).

For an exact solution, let fcor=0, v=a0*exp(-lambda r)*r**0.5

Then we see,

p(r) = p0*expleft(-frac{rho0}{p0}int^{inf}_{r}left(a0**2*exp(-2lambda r)right)right)

= p0*expleft(frac{a0**2 * rho0}{2* p0 lambda}left[*exp(-2lambda r)right]^{inf}_{r}right) = p0expleft(-frac{a0**2 * rho0}{2 * p0 * lambda}*a0**2*exp(- 2 lambda r)right).

Parameters:
  • rr (np.ndarray) – radii array [m].

  • vv (np.ndarray) – velocity array [m/s]

  • p0 (float) – ambient pressure [Pa].

  • rho0 (float, optional) – Air density at ambient pressure [kg/m3].

  • fcor (float, optional) – Coriolis force [s-1].

  • assumption (Literal["isopycnal", "isothermal"], optional)

Returns:

Pressure array [Pa].

Return type:

np.ndarray

Example::
>>> rr = np.array([0, 1, 2, 3, 4, 5])
>>> vv = np.array([0] * 6)
>>> p = pressure_from_wind(rr, vv, assumption="isopycnal")
>>> np.allclose(p, np.array([101500] * 6), rtol=1e-3, atol=1e-6) # zero velocity -> no change.
True
>>> p = pressure_from_wind(rr, vv, assumption="isothermal")
>>> np.allclose(p, np.array([101500] * 6), rtol=1e-3, atol=1e-6) # zero velocity -> no change.
True
>>> r0 = 1_000_000
>>> ds = 1/1_000 # decay scale [m-1]
>>> rho0 = 1.0 # [kg m-3]
>>> p0 = 1000_00.0 # [Pa]
>>> rr = np.linspace(0, r0, 1_000_000)
>>> a0 = 1.0 # [m0.5 s-1]
>>> vv = a0 * np.exp(-ds * rr) * (rr)**0.5
>>> pest = pressure_from_wind(rr, vv, p0=p0, rho0=rho0, fcor=0.0, assumption="isopycnal")
>>> pcalc = p0 - rho0 * a0**2 / (2*ds) * (np.exp(-2*ds*rr) - np.exp(-2*ds*r0))
>>> if not np.allclose(pest, pcalc, rtol=1e-4, atol=1e-6):
...    print("pest", pest[:10], pest[::-1][:10][::-1])
...    print("pcalc", pcalc[:10], pcalc[::-1][:10][::-1])
>>> pest = pressure_from_wind(rr, vv, p0=p0, rho0=rho0, fcor=0.0, assumption="isothermal")
>>> integral = -(a0**2*rho0)/(2*p0*ds) * (np.exp(-2*ds*rr) - np.exp(-2*ds*r0))
>>> assert np.all(integral <= 0)
>>> pcalc = p0 * np.exp(integral)
>>> if not np.allclose(pest, pcalc, rtol=1e-3, atol=1e-6):
...    print("pest", pest[:10], pest[::-1][:10][::-1])
...    print("pcalc", pcalc[:10], pcalc[::-1][:10][::-1])
w22.utils.qair2rh(qair, temp, press=101600)

Convert specific humidity to relative humidity.

Inspired by: https://earthscience.stackexchange.com/a/2385

Parameters:
  • qair (Union[xr.DataArray, float]) – Specific humidity [dimensionless].

  • temp (Union[xr.DataArray, float]) – Temperature [K].

  • press (Union[xr.DataArray, float]) – Pressure [hPa].

Returns:

Relative humidity [dimensionless].

Return type:

Union[xr.DataArray, float]

Example::
>>> print(f'{qair2rh(0.01, 300, 1013.25):.3f}')
0.458
>>> qair = xr.DataArray([0.001, 0.002, 0.003], dims=["x"], attrs={"units": "dimensionless"})
>>> temp = xr.DataArray([300, 300, 300], dims=["x"], attrs={"units": "K"})
>>> press = xr.DataArray([1013.25, 1013.25, 1013.25], dims=["x"], attrs={"units": "hPa"})
>>> rh = qair2rh(qair, temp, press) #>>> np.allclose(rh, np.array([0.0458*i for i in range(1, 4)]), rtol=1e-3, atol=1e-6)
>>> assert rh.attrs["units"] == "dimensionless"
>>> assert rh.attrs["long_name"] == "Relative Humdity"
w22.utils.qtp2rh(qa, ta, msl)

Generate surface relative humidity from the specific humidity and temperature volumes, and pressure. The pressure level coordinate is “p” in Pa.

Parameters:
  • qa (xr.DataArray) – Specific humidity [dimensionless] at many pressure levels.

  • ta (xr.DataArray) – Air temperature [K] at the same pressure levels.

  • msl (xr.DataArray) – The mean sea level pressure [Pa].

Returns:

Relative humidity [dimensionless] at the surface.

Return type:

xr.DataArray

Example::
>>> pressure_levels = np.array([1000, 900])
>>> qa = xr.DataArray(np.array([0.01, 0.02]), dims=["p"], coords={"p": ("p", pressure_levels, {"units": "hPa"})}, attrs={"units": "dimensionless"})
>>> ta = xr.DataArray(np.array([300, 300]), dims=["p"], coords={"p": ("p", pressure_levels, {"units": "hPa"})}, attrs={"units": "K"})
>>> msl = xr.DataArray(np.array([1000]), attrs={"units": "hPa"})
>>> rh = qtp2rh(qa, ta, msl)
>>> assert rh.attrs["units"] == "dimensionless"
>>> assert rh.attrs["long_name"] == "Relative Humdity"
w22.utils.rho_air_f(total_pressure_hpa, air_temperature_k, water_vapour_pressure_pa)

Calculate air density given total pressure, temperature, and water vapour pressure.

Parameters:
  • total_pressure_hpa (float) – Total pressure in hPa.

  • air_temperature (float) – Air temperature in K.

  • water_vapour_pressure (float) – Water vapour pressure in Pa.

Returns:

Air density in kg/m^3.

Return type:

float

Doctests: >>> # Test with standard sea level conditions (approx) >>> # P = 1013.25 hPa, T = 15 C = 288.15 K, RH = 50% >>> # Saturated vapor pressure at 15 C is ~1705 Pa. 50% RH -> p_v = 852.5 Pa >>> round(rho_air_f(1013.25, 288.15, 852.5), 3) 1.221 >>> # Test for dry air (vapor pressure = 0) >>> round(rho_air_f(1013.25, 288.15, 0.0), 3) 1.225

w22.w22_carnot module

Some functions for the Wang 2022 “Tropical Cyclone Potential Size” sub-Carnot engine model.

w22.w22_carnot.wang_carnot_velocity(efficiency_relative_to_carnot=0.5, near_surface_air_temperature=299, outflow_temperature=200, latent_heat_of_vaporization=2500000, gas_constant_for_water_vapor=461.5, gas_constant=287, pressure_dry_at_inflow=98500)

Calculate the sub-Carnot velocity from the Wang 2022 sub-Carnot engine model.

In equation 27 they make a series of assumptions, including that intensity does not signicantly change with potential size to get the formula.

V_carnot**2 = (eta * epsilon_c * L_v - R_v * T_s) * q_va^{*}

q_va^{star} = R/R_v (e*s)

They interpret this as the largest possible wind speed that can be reached in the outflow at radius r_a.

Parameters:
  • efficiency_relative_to_carnot (float, optional) – Efficiency relative to Carnot. Defaults to EFFICIENCY_RELATIVE_TO_CARNOT_DEFAULT.

  • near_surface_air_temperature (float, optional) – Near surface air temperature in Kelvin. Defaults to NEAR_SURFACE_AIR_TEMPERATURE_DEFAULT.

  • outflow_temperature (float, optional) – Outflow temperature in Kelvin. Defaults to OUTFLOW_TEMPERATURE_DEFAULT.

  • latent_heat_of_vaporization (float, optional) – Latent heat of vaporization in J/kg. Defaults to LATENT_HEAT_OF_VAPORIZATION.

  • gas_constant_for_water_vapor (float, optional) – Gas constant for water vapor in J/kg/K. Defaults to GAS_CONSTANT_FOR_WATER_VAPOR.

  • gas_constant (float, optional) – Gas constant in J/kg/K. Defaults to GAS_CONSTANT.

  • pressure_dry_at_inflow (float, optional) – Dry air pressure at inflow in Pa. Defaults to PRESSURE_DRY_AT_INFLOW_DEFAULT.

Return type:

float

w22.w22_carnot.wang_consts(near_surface_air_temperature=299, outflow_temperature=200, latent_heat_of_vaporization=2500000, gas_constant_for_water_vapor=461.5, gas_constant=287, beta_lift_parameterization=1.25, efficiency_relative_to_carnot=0.5, pressure_dry_at_inflow=98500, coriolis_parameter=5e-05, maximum_wind_speed=83, radius_of_inflow=2193000, radius_of_max_wind=64000)

Wang 2022 sub-Carnot engine model parameters. This function calculates the constants that are then used to calculate the difference function whose root we want to find.

Parameters:
  • near_surface_air_temperature (float, optional) – Defaults to 299 [K].

  • outflow_temperature (float, optional) – Defaults to 200 [K].

  • latent_heat_of_vaporization (float, optional) – Defaults to 2.27e6 [J/kg].

  • gas_constant_for_water_vapor (float, optional) – Defaults to 461 [J/kg/K].

  • gas_constant (float, optional) – Defaults to 287 [J/kg/K].

  • beta_lift_parameterization (float, optional) – Defaults to 1.25 [dimesionless].

  • efficiency_relative_to_carnot (float, optional) – Defaults to 0.5 [dimensionless].

  • pressure_dry_at_inflow (float, optional) – Defaults to 985 * 100 [Pa].

  • coriolis_parameter (float, optional) – Defaults to 5e-5 [s-1].

  • maximum_wind_speed (float, optional) – Defaults to 83 [m/s].

  • radius_of_inflow (float, optional) – Defaults to 2193 * 1000 [m].

  • radius_of_max_wind (float, optional) – Defaults to 64 * 1000 [m].

Returns:

a, b, c (equations 17–19 in Wang et al. 2022).

Return type:

Tuple[float, float, float]

# only works if you use dodgy value for constant of latent heat of vaporization of water (2_500_000 J kg-1) # and not the more standard value (2_268_000 J kg-1) # this is because at a point in the derivation they ignore temperature variations, and so pick to use values at 0C # 2_500_000 J kg-1 is the value at 0C (https://met.nps.edu/~bcreasey/mr3222/files/helpful/UnitsandConstantsUsefulInMeteorology-PSU.html) >>> a, b, c = wang_consts(near_surface_air_temperature=299, outflow_temperature=200, latent_heat_of_vaporization=2_500_000+ 0*2_268_000, gas_constant_for_water_vapor=461, gas_constant=287, beta_lift_parameterization=5/4, efficiency_relative_to_carnot=0.5, pressure_dry_at_inflow=985_00, coriolis_parameter=5e-5, maximum_wind_speed=83, radius_of_inflow=2193_000, radius_of_max_wind=64_000) >>> f”{c:.3f}” ‘0.008’ >>> f”{b:.3f}” ‘0.031’ >>> f”{a:.3f}” ‘0.062’

w22.w22_carnot.wang_diff(a=0.062, b=0.031, c=0.008)

Wang et al. 2022 difference function to find roots of.

Parameters:
  • a (float, optional) –

    1. Defaults to 0.062.

  • b (float, optional) –

    1. Defaults to 0.031.

  • c (float, optional) –

    1. Defaults to 0.008.

Returns:

Function to find root of.

Return type:

Callable[[float], float]

Example::
>>> f = wang_diff(a=0.062, b=0.031, c=0.008)
>>> f"{f(1.081):.3f}"
'0.000'
>>> f"{f(19.1829):.3f}" # root said to be around 18 in paper. Is 19.18 close enough?
'0.000'
w22.w22_carnot.wang_outer_radius_approx(efficiency_relative_to_carnot=0.5, near_surface_air_temperature=299, outflow_temperature=200, latent_heat_of_vaporization=2500000, gas_constant_for_water_vapor=461.5, gas_constant=287, pressure_dry_at_inflow=98500, coriolis_parameter=5e-05)

Calculate the length scale from the Wang 2022 sub-Carnot engine model. The factor of 2 is from Appendix B.

Most of the time the length scale is referred to as V_carnot/f_cor, but there is a factor of 2 in derivation in Appendix B, so we include that here for completeness.

Parameters:
  • efficiency_relative_to_carnot (float, optional) – Efficiency relative to Carnot. Defaults to EFFICIENCY_RELATIVE_TO_CARNOT_DEFAULT.

  • near_surface_air_temperature (float, optional) – Near surface air temperature in Kelvin. Defaults to NEAR_SURFACE_AIR_TEMPERATURE_DEFAULT.

  • outflow_temperature (float, optional) – Outflow temperature in Kelvin. Defaults to OUTFLOW_TEMPERATURE_DEFAULT.

  • latent_heat_of_vaporization (float, optional) – Latent heat of vaporization in J/kg. Defaults to LATENT_HEAT_OF_VAPORIZATION.

  • gas_constant_for_water_vapor (float, optional) – Gas constant for water vapor in J/kg/K. Defaults to GAS_CONSTANT_FOR_WATER_VAPOR.

  • gas_constant (float, optional) – Gas constant in J/kg/K. Defaults to GAS_CONSTANT.

  • pressure_dry_at_inflow (float, optional) – Dry air pressure at inflow in Pa. Defaults to PRESSURE_DRY_AT_INFLOW_DEFAULT.

  • coriolis_parameter (float, optional) – Coriolis parameter in s-1. Defaults to F_COR_DEFAULT.

Module contents