Basic use for assessing storm climatology in a dataset

Here, we examplify usage of huracanpy with the dataset of TC in ERA-20C detected by the TRACK algorithm. This is meant to show an example of workflow. Please refer to specific parts of the documentation to learn about each part (e.g. loading, plotting, etc.) in more detail.

import huracanpy

import matplotlib.pyplot as plt
import seaborn as sns

Read the file

huracanpy’s load function can handle different track file types. Here, the data is available as a netcdf file.

data = huracanpy.load(huracanpy.example_ERA20C_file)
data.psl.attrs["units"] = "hPa"  # Fixing misspelled pressure unit
data

/home/docs/checkouts/readthedocs.org/user_builds/huracanpy/envs/v1.4.0/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

<xarray.Dataset> Size: 2MB
Dimensions:         (record: 24865)
Coordinates:
    time            (record) datetime64[ns] 199kB ...
Dimensions without coordinates: record
Data variables:
    longitude_psl   (record) float32 99kB ...
    latitude_psl    (record) float32 99kB ...
    psl             (record) float64 199kB ...
    wind_speed_10m  (record) float32 99kB ...
    lon             (record) float32 99kB ...
    lat             (record) float32 99kB ...
    track_id        (record) <U7 696kB '1900-12' '1900-12' ... '2010-9' '2010-9'
Attributes:
    realm:        atmos
    history:      testing
    Conventions:  CF-1.7
    featureType:  trajectory

xarray.Dataset

• Dimensions:
  • record: 24865
• Coordinates: (1)
  • time
    (record)
    datetime64[ns]
    ...

    standard_name :

    time

    long_name :

    Time

    time_calendar :

    gregorian

    start :

    1979010100

    step :

    6

    [24865 values with dtype=datetime64[ns]]

• Data variables: (7)
  • longitude_psl
    (record)
    float32
    ...

    standard_name :

    longitude

    long_name :

    Longitude

    units :

    degrees_east

    [24865 values with dtype=float32]

  • latitude_psl
    (record)
    float32
    ...

    standard_name :

    latitude

    long_name :

    Latitude

    units :

    degrees_north

    [24865 values with dtype=float32]

  • psl
    (record)
    float64
    ...

    standard_name :

    psl

    long_name :

    air pressure at mean sea level

    units :

    hPa

    [24865 values with dtype=float64]

  • wind_speed_10m
    (record)
    float32
    ...

    standard_name :

    wind_speed_10m

    long_name :

    Wind Speed at 10m

    units :

    m s**-1

    [24865 values with dtype=float32]

  • lon
    (record)
    float32
    ...

    standard_name :

    longitude

    long_name :

    Longitude

    units :

    degrees_east

    [24865 values with dtype=float32]

  • lat
    (record)
    float32
    ...

    standard_name :

    latitude

    long_name :

    Latitude

    units :

    degrees_north

    [24865 values with dtype=float32]

  • track_id
    (record)
    <U7
    '1900-12' '1900-12' ... '2010-9'

    cf_role :

    trajectory_id

    array(['1900-12', '1900-12', '1900-12', ..., '2010-9', '2010-9', '2010-9'],
          shape=(24865,), dtype='<U7')

• Attributes: (4)

  realm :

  atmos

  history :

  testing

  Conventions :

  CF-1.7

  featureType :

  trajectory

Add useful information

After loading, you can add various useful information for the analysis (basin, season, category…)

# Use the accessor's add_ function to add the info you want
## SSHS category
data = data.hrcn.add_saffir_simpson_category(wind_name="wind_speed_10m")
## Presure category
data = data.hrcn.add_pressure_category(slp_name="psl")
## Season
data = data.hrcn.add_season()
# More info are available, in particular geographical ones, but we do not need them for
# this example.

data

<xarray.Dataset> Size: 2MB
Dimensions:                  (record: 24865)
Coordinates:
    time                     (record) datetime64[ns] 199kB 1900-07-23T18:00:0...
Dimensions without coordinates: record
Data variables:
    longitude_psl            (record) float32 99kB ...
    latitude_psl             (record) float32 99kB ...
    psl                      (record) float64 199kB 1.012e+03 ... 996.0
    wind_speed_10m           (record) float32 99kB 7.798 7.301 ... 12.17 19.06
    lon                      (record) float32 99kB ...
    lat                      (record) float32 99kB 9.978 10.07 ... 61.93 62.89
    track_id                 (record) <U7 696kB '1900-12' '1900-12' ... '2010-9'
    saffir_simpson_category  (record) int64 199kB -1 -1 -1 -1 -1 ... -1 -1 -1 0
    pressure_category        (record) int64 199kB -1 -1 -1 -1 -1 ... 0 0 0 0 0
    season                   (record) float64 199kB 1.9e+03 1.9e+03 ... 2.01e+03
Attributes:
    realm:        atmos
    history:      testing
    Conventions:  CF-1.7
    featureType:  trajectory

xarray.Dataset

• Dimensions:
  • record: 24865
• Coordinates: (1)
  • time
    (record)
    datetime64[ns]
    1900-07-23T18:00:00 ... 2010-09-...

    standard_name :

    time

    long_name :

    Time

    time_calendar :

    gregorian

    start :

    1979010100

    step :

    6

    array(['1900-07-23T18:00:00.000000000', '1900-07-24T00:00:00.000000000',
           '1900-07-24T06:00:00.000000000', ..., '2010-09-04T00:00:00.000000000',
           '2010-09-04T06:00:00.000000000', '2010-09-04T12:00:00.000000000'],
          shape=(24865,), dtype='datetime64[ns]')

• Data variables: (10)
  • longitude_psl
    (record)
    float32
    ...

    standard_name :

    longitude

    long_name :

    Longitude

    units :

    degrees_east

    [24865 values with dtype=float32]

  • latitude_psl
    (record)
    float32
    ...

    standard_name :

    latitude

    long_name :

    Latitude

    units :

    degrees_north

    [24865 values with dtype=float32]

  • psl
    (record)
    float64
    1.012e+03 1.015e+03 ... 996.7 996.0

    standard_name :

    psl

    long_name :

    air pressure at mean sea level

    units :

    hPa

    array([1012.330017, 1015.466003, 1011.830994, ...,  993.535217,  996.716919,
            995.981079], shape=(24865,))

  • wind_speed_10m
    (record)
    float32
    7.798 7.301 8.169 ... 12.17 19.06

    standard_name :

    wind_speed_10m

    long_name :

    Wind Speed at 10m

    units :

    m s**-1

    array([ 7.798183,  7.301081,  8.16876 , ..., 13.69531 , 12.16668 , 19.05791 ],
          shape=(24865,), dtype=float32)

  • lon
    (record)
    float32
    ...

    standard_name :

    longitude

    long_name :

    Longitude

    units :

    degrees_east

    [24865 values with dtype=float32]

  • lat
    (record)
    float32
    9.978 10.07 9.752 ... 61.93 62.89

    standard_name :

    latitude

    long_name :

    Latitude

    units :

    degrees_north

    array([ 9.977631, 10.066322,  9.751551, ..., 62.86075 , 61.92685 , 62.894787],
          shape=(24865,), dtype=float32)

  • track_id
    (record)
    <U7
    '1900-12' '1900-12' ... '2010-9'

    cf_role :

    trajectory_id

    array(['1900-12', '1900-12', '1900-12', ..., '2010-9', '2010-9', '2010-9'],
          shape=(24865,), dtype='<U7')

  • saffir_simpson_category
    (record)
    int64
    -1 -1 -1 -1 -1 -1 ... -1 -1 -1 -1 0

    array([-1, -1, -1, ..., -1, -1,  0], shape=(24865,))

  • pressure_category
    (record)
    int64
    -1 -1 -1 -1 -1 -1 ... 0 0 0 0 0 0

    units :

    dimensionless

    array([-1, -1, -1, ...,  0,  0,  0], shape=(24865,))

  • season
    (record)
    float64
    1.9e+03 1.9e+03 ... 2.01e+03

    array([1900., 1900., 1900., ..., 2010., 2010., 2010.], shape=(24865,))

• Attributes: (4)

  realm :

  atmos

  history :

  testing

  Conventions :

  CF-1.7

  featureType :

  trajectory

Check the content of the file

huracanpy provide a coarse plotting function that you can use for checking what is in your data.

# Basic plot of the data points
data.hrcn.plot_tracks(intensity_var_name="wind_speed_10m")

(<Figure size 1000x1000 with 1 Axes>, <GeoAxes: xlabel='lon', ylabel='lat'>)

/home/docs/checkouts/readthedocs.org/user_builds/huracanpy/envs/v1.4.0/lib/python3.12/site-packages/cartopy/io/__init__.py:242: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/110m_physical/ne_110m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

Climatological metrics

You can compute basic statistics: frequency, TC days, ACE. Here shown as yearly averages.

# Frequency (Number of track per year)
data.track_id.hrcn.nunique() / data.season.hrcn.nunique()
# number of unique tracks / number of unique season

5.611111111111111

# TCD (Accumulated duration of storms per year)
data.hrcn.get_track_duration().sum() / 24 / data.season.hrcn.nunique()
# Compute duration per track, convert to days / number of unique season

<xarray.DataArray 'duration' ()> Size: 8B
array(56.71990741)
Attributes:
    standard_name:  time
    long_name:      Time
    time_calendar:   gregorian
    start:          1979010100
    step:           6
    units:          hours

xarray.DataArray

'duration'

• 56.72

  array(56.71990741)

• Attributes: (6)

  standard_name :

  time

  long_name :

  Time

  time_calendar :

  gregorian

  start :

  1979010100

  step :

  6

  units :

  hours

# ACE per year
data = data.hrcn.add_ace(wind_name="wind_speed_10m")
data.ace.groupby(data.season).sum().mean()
# Get ace for each point, sum by season and average over the seasons
# NB: By default, huracanpy computes ACE only for points with wind above 34 knots

<xarray.DataArray 'ace' ()> Size: 4B
array(5.252404, dtype=float32)
Attributes:
    units:    knot ** 2

xarray.DataArray

'ace'

• 5.252

  array(5.252404, dtype=float32)

• Attributes: (1)

  units :

  knot ** 2

Variability

With xarray’s grouping functionnalities, you can show variations of these statistics.

## Interannual

fig, axs = plt.subplots(3, sharex=True)
# Frequency
# In this case, it is easier to go through a dataframe
data.to_dataframe().groupby("season").track_id.nunique().plot(ax=axs[0])
axs[0].set_ylabel("Number of tracks")
# TCD
data.groupby("season").apply(lambda s: s.hrcn.get_track_duration().sum()).plot(
    ax=axs[1]
)
axs[1].set_ylabel("TC days")
# ACE
data.groupby("season").sum().ace.plot(ax=axs[2])
axs[2].set_ylabel("ACE")

for ax in axs:
    ax.set_ylim(0)
    ax.set_xlabel("")

/home/docs/checkouts/readthedocs.org/user_builds/huracanpy/envs/v1.4.0/lib/python3.12/site-packages/xarray/structure/concat.py:674: UserWarning: No index created for dimension season because variable season is not a coordinate. To create an index for season, please first call `.set_coords('season')` on this object.
  ds.expand_dims(dim_name, create_index_for_new_dim=create_index_for_new_dim)

## Seasonal

gen = data.hrcn.get_gen_vals()  # Extract the point of genesis for each track
(
    gen.groupby("time.month").count().lon  # compute number of genesis points per month
    / gen.season.hrcn.nunique()  # Normalize by number of season
).plot(marker="o")  # plot

[<matplotlib.lines.Line2D at 0x7a947c7a2a50>]

Track statistics

You can also compute track-level statistics such as duration and lifetime maximum intensity.

fig, axs = plt.subplots(1, 3, sharey=True, figsize=[15, 5])

# Duration
data.hrcn.get_track_duration().plot.hist(ax=axs[0])

# Maximum wind speed
data.wind_speed_10m.groupby(data.track_id).max().plot.hist(ax=axs[1])

# Minimum SLP
data.psl.groupby(data.track_id).min().plot.hist(ax=axs[2])

(array([  6.,  17.,  34.,  51.,  58.,  77.,  89., 235.,  37.,   2.]),
 array([ 941.6126709 ,  950.1151001 ,  958.6175293 ,  967.1199585 ,
         975.6223877 ,  984.12481689,  992.62724609, 1001.12967529,
        1009.63210449, 1018.13453369, 1026.63696289]),
 <BarContainer object of 10 artists>)

Lifecycles

# Compute times from extremum
data["time_from_slp_min"] = data.hrcn.get_time_from_apex(
    intensity_var_name="psl", stat="min"
)
data["time_from_wind_max"] = data.hrcn.get_time_from_apex(
    intensity_var_name="wind_speed_10m",
)

# Plot with seaborn
fig, axs = plt.subplots(2, sharex=True)
# SLP lifecycle
sns.lineplot(data=data, x="time_from_slp_min", y="psl", ax=axs[0])
# Wind lifecycle
sns.lineplot(data=data, x="time_from_wind_max", y="wind_speed_10m", ax=axs[1])

<Axes: xlabel='time_from_wind_max', ylabel='wind_speed_10m'>

Track density

# Density of all points
d = data.hrcn.get_density(bin_size=5) / data.season.hrcn.nunique()
huracanpy.plot.density(
    d,
    cbar_kwargs=dict(label="Number of points per 5x5° box per year"),
)

/home/docs/checkouts/readthedocs.org/user_builds/huracanpy/envs/v1.4.0/lib/python3.12/site-packages/huracanpy/calc/_density.py:88: UserWarning: By default density does not take into account the spherical geometry ofthe Earth. Set spherical=True to account for this
  warnings.warn(

(<Figure size 640x480 with 2 Axes>, <GeoAxes: xlabel='lon', ylabel='lat'>)

# Density of all points
d = gen.hrcn.get_density(bin_size=10) / data.season.hrcn.nunique()
huracanpy.plot.density(
    d,
    cbar_kwargs=dict(label="Number of genesis per 10x10° box per year"),
)

/home/docs/checkouts/readthedocs.org/user_builds/huracanpy/envs/v1.4.0/lib/python3.12/site-packages/huracanpy/calc/_density.py:88: UserWarning: By default density does not take into account the spherical geometry ofthe Earth. Set spherical=True to account for this
  warnings.warn(

(<Figure size 640x480 with 2 Axes>, <GeoAxes: xlabel='lon', ylabel='lat'>)