Tracking options

We explore the impact of some of the tracking options for two of the sample datasets.

import warnings
from string import ascii_lowercase

import cartopy
import cartopy.crs as ccrs
import matplotlib as mpl
import matplotlib.patheffects as path_effects
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import xarray as xr
from IPython.display import display, Markdown
from joblib import delayed, Parallel

import tams

warnings.filterwarnings("ignore", category=cartopy.io.DownloadWarning)

plt.rcParams.update(
    {
        "axes.formatter.use_mathtext": True,
    }
)

%matplotlib inline

xr.set_options(display_expand_data=False)

Show code cell output Hide code cell output

<xarray.core.options.set_options at 0x230061013f0>

For sample satellite data

tb = tams.load_example_tb()
tb

<xarray.DataArray 'ch9' (time: 6, y: 714, x: 2829)>
295.9 295.9 295.8 295.8 295.8 295.9 295.9 295.8 ... nan nan nan nan nan nan nan
Coordinates:
    lon      (time, y, x) float32 ...
    lat      (time, y, x) float32 ...
  * time     (time) datetime64[ns] 2006-09-01 ... 2006-09-01T10:00:00
Dimensions without coordinates: y, x
Attributes:
    units:      K
    long_name:  Brightness temperature

xarray.DataArray

'ch9'

• time: 6
• y: 714
• x: 2829

• 295.9 295.9 295.8 295.8 295.8 295.9 295.9 ... nan nan nan nan nan nan

  array([[[295.92957, 295.92957, 295.8033 , ..., 292.85956, 292.85956,
           292.7298 ],
          [295.8033 , 295.8033 , 295.92957, ..., 293.11865, 292.98917,
           292.46976],
          [295.8033 , 295.67694, 295.8033 , ..., 292.85956, 292.98917,
           292.59985],
          ...,
          [      nan,       nan,       nan, ...,       nan,       nan,
                 nan],
          [      nan,       nan,       nan, ...,       nan,       nan,
                 nan],
          [      nan,       nan,       nan, ...,       nan,       nan,
                 nan]],
  
         [[295.8033 , 295.8033 , 295.8033 , ..., 287.54327, 285.89505,
           284.92136],
          [296.05566, 295.8033 , 295.4237 , ..., 288.76266, 287.67944,
           286.0334 ],
          [295.92957, 295.8033 , 295.55038, ..., 290.1013 , 289.70148,
           288.08698],
  ...
          [      nan,       nan,       nan, ...,       nan,       nan,
                 nan],
          [      nan,       nan,       nan, ...,       nan,       nan,
                 nan],
          [      nan,       nan,       nan, ...,       nan,       nan,
                 nan]],
  
         [[293.63504, 293.76373, 293.24796, ..., 290.8965 , 289.3001 ,
           286.99673],
          [293.76373, 293.76373, 293.37717, ..., 289.43405, 286.17157,
           286.30954],
          [293.76373, 293.63504, 293.63504, ..., 288.62784, 285.20047,
           287.13364],
          ...,
          [      nan,       nan,       nan, ...,       nan,       nan,
                 nan],
          [      nan,       nan,       nan, ...,       nan,       nan,
                 nan],
          [      nan,       nan,       nan, ...,       nan,       nan,
                 nan]]], dtype=float32)

• Coordinates: (3)
  • lon
    (time, y, x)
    float32
    ...

    _CoordinateAxisType :

    Lon

    standard_name :

    longitude

    valid_max :

    180.0

    long_name :

    Longitude

    valid_min :

    -180.0

    units :

    degrees_east

    [12119436 values with dtype=float32]

  • lat
    (time, y, x)
    float32
    ...

    _CoordinateAxisType :

    Lat

    standard_name :

    latitude

    valid_max :

    90.0

    long_name :

    Latitude

    valid_min :

    -90.0

    units :

    degrees_north

    [12119436 values with dtype=float32]

  • time
    (time)
    datetime64[ns]
    2006-09-01 ... 2006-09-01T10:00:00

    array(['2006-09-01T00:00:00.000000000', '2006-09-01T02:00:00.000000000',
           '2006-09-01T04:00:00.000000000', '2006-09-01T06:00:00.000000000',
           '2006-09-01T08:00:00.000000000', '2006-09-01T10:00:00.000000000'],
          dtype='datetime64[ns]')

• Indexes: (1)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2006-09-01 00:00:00', '2006-09-01 02:00:00',
                   '2006-09-01 04:00:00', '2006-09-01 06:00:00',
                   '2006-09-01 08:00:00', '2006-09-01 10:00:00'],
                  dtype='datetime64[ns]', name='time', freq='2H'))

• Attributes: (2)

  units :

  K

  long_name :

  Brightness temperature

Identify cloud elements (CEs)

Our tracking options cases here are all in the tracking stage, so we only have to run tams.identify() once.

%%time

ces, _ = tams.identify(tb, parallel=True)

Show code cell output Hide code cell output

[Parallel(n_jobs=-2)]: Using backend LokyBackend with 11 concurrent workers.

CPU times: total: 172 ms
Wall time: 5.26 s

[Parallel(n_jobs=-2)]: Done   2 out of   6 | elapsed:    5.0s remaining:   10.1s
[Parallel(n_jobs=-2)]: Done   3 out of   6 | elapsed:    5.0s remaining:    5.0s
[Parallel(n_jobs=-2)]: Done   4 out of   6 | elapsed:    5.0s remaining:    2.5s
[Parallel(n_jobs=-2)]: Done   6 out of   6 | elapsed:    5.1s finished

Run cases

%%time

cases = {}
time = tb.time.values.tolist()

proj = -15

with warnings.catch_warnings():
    warnings.filterwarnings("ignore", message="forward `look` considered experimental")
    cases["default"] = tams.track(ces, time)
    cases["link=f, norm=min"] = tams.track(ces, time, look="f", overlap_norm="min")
    cases["link=f, norm=max"] = tams.track(ces, time, look="f", overlap_norm="max")
    cases[f"u_proj={proj}"] = tams.track(ces, time, u_projection=proj)
    cases[f"u_proj={proj}, link=f, norm=min"] = tams.track(ces, time, u_projection=proj, look="f", overlap_norm="min")
    cases[f"u_proj={proj}, link=f, norm=max"] = tams.track(ces, time, u_projection=proj, look="f", overlap_norm="max")

for key, ce in cases.items():
    cases[key] = tams.classify(ce)

CPU times: total: 5.83 s
Wall time: 5.83 s

Plot CEs

And check if the method produced a single track within the box, saving a table of these results.

m, n = 2, 3
assert len(cases) <= m * n

# Restrict range
# lon_min, lat_min, lon_max, lat_max = -10, 8, -3, 15
# lon_min, lat_min, lon_max, lat_max = 17, 1, 27, 8
# lon_min, lat_min, lon_max, lat_max = 20, 1, 27, 4
lon_min, lat_min, lon_max, lat_max = 30, 10, 38, 16

cases_full = []
for key in cases:
    d = {"u_proj": 0, "link": "b", "norm": "a"}
    updates = {}
    if key != "default":
        for kv in key.split(", "):
            k, v = kv.split("=")
            updates[k] = v
        if "norm" not in updates and updates.get("link", "b").startswith("f"):
            d["norm"] = "b"
    d.update(updates)
    cases_full.append(d)
tbl = pd.DataFrame(cases_full)
tbl = tbl.assign(continued=False)

for i, (_, ce) in enumerate(cases.items()):
    n = ce.cx[lon_min:lon_max, lat_min:lat_max].mcs_id.nunique()
    assert n in {1, 2}
    if n == 1:
        tbl.loc[i, "continued"] = True

tams.plot_tracked(ce, label="none", size=3.5, add_colorbar=True, cbar_kwargs=dict(fraction=0.05, shrink=0.3))

ax = plt.gca()
patch = mpl.patches.Polygon(
    [(lon_min, lat_min), (lon_min, lat_max), (lon_max, lat_max), (lon_max, lat_min)],
    ec="orangered",
    fill=False,
    transform=ccrs.PlateCarree(),
)
ax.add_patch(patch)

# proj = ccrs.PlateCarree()

# fig, axs = plt.subplots(m, n, figsize=(8, 5), sharex=True, sharey=True, subplot_kw=dict(projection=proj), constrained_layout=True)

# for i, ((key, ce), ax) in enumerate(zip(cases.items(), axs.flat)):
#     if i == 0:
#         assert key == "default"
#         tams.plot_tracked(ce)
    
#     tams.plot_tracked(ce.cx[lon_min:lon_max, lat_min:lat_max], ax=ax)
#     ax.text(0.03, 0.98, key, ha="left", va="top", transform=ax.transAxes, size=8)
#     gl = ax.gridlines(draw_labels=True, color="none")
#     if not i < n:
#         gl.top_labels = False
#     if not i % n == 0:
#         gl.left_labels = False

with open("sat-tracking-options-table.txt", "w") as f:
    f.write(tbl.assign(continued=tbl.continued.map({True: "y", False: "n"})).style.to_latex())

plt.savefig("sat-tracking-options-ces-box.pdf", bbox_inches="tight", pad_inches=0.05, transparent=False)

For sample MPAS lat/lon dataset

ds = tams.load_example_mpas().isel(time=slice(1, None)).rename_vars(tb="ctt", precip="pr")
ds

<xarray.Dataset>
Dimensions:  (time: 126, lon: 341, lat: 180)
Coordinates:
  * time     (time) datetime64[ns] 2006-09-08T13:00:00 ... 2006-09-13T18:00:00
  * lon      (lon) float64 85.0 85.25 85.5 85.75 ... 169.2 169.5 169.8 170.0
  * lat      (lat) float64 -4.875 -4.625 -4.375 -4.125 ... 39.38 39.62 39.88
Data variables:
    ctt      (time, lat, lon) float32 266.8 266.9 267.0 ... 266.0 265.9 265.7
    pr       (time, lat, lon) float32 ...

xarray.Dataset

• Dimensions:
  • time: 126
  • lon: 341
  • lat: 180
• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    2006-09-08T13:00:00 ... 2006-09-...

    array(['2006-09-08T13:00:00.000000000', '2006-09-08T14:00:00.000000000',
           '2006-09-08T15:00:00.000000000', '2006-09-08T16:00:00.000000000',
           '2006-09-08T17:00:00.000000000', '2006-09-08T18:00:00.000000000',
           '2006-09-08T19:00:00.000000000', '2006-09-08T20:00:00.000000000',
           '2006-09-08T21:00:00.000000000', '2006-09-08T22:00:00.000000000',
           '2006-09-08T23:00:00.000000000', '2006-09-09T00:00:00.000000000',
           '2006-09-09T01:00:00.000000000', '2006-09-09T02:00:00.000000000',
           '2006-09-09T03:00:00.000000000', '2006-09-09T04:00:00.000000000',
           '2006-09-09T05:00:00.000000000', '2006-09-09T06:00:00.000000000',
           '2006-09-09T07:00:00.000000000', '2006-09-09T08:00:00.000000000',
           '2006-09-09T09:00:00.000000000', '2006-09-09T10:00:00.000000000',
           '2006-09-09T11:00:00.000000000', '2006-09-09T12:00:00.000000000',
           '2006-09-09T13:00:00.000000000', '2006-09-09T14:00:00.000000000',
           '2006-09-09T15:00:00.000000000', '2006-09-09T16:00:00.000000000',
           '2006-09-09T17:00:00.000000000', '2006-09-09T18:00:00.000000000',
           '2006-09-09T19:00:00.000000000', '2006-09-09T20:00:00.000000000',
           '2006-09-09T21:00:00.000000000', '2006-09-09T22:00:00.000000000',
           '2006-09-09T23:00:00.000000000', '2006-09-10T00:00:00.000000000',
           '2006-09-10T01:00:00.000000000', '2006-09-10T02:00:00.000000000',
           '2006-09-10T03:00:00.000000000', '2006-09-10T04:00:00.000000000',
           '2006-09-10T05:00:00.000000000', '2006-09-10T06:00:00.000000000',
           '2006-09-10T07:00:00.000000000', '2006-09-10T08:00:00.000000000',
           '2006-09-10T09:00:00.000000000', '2006-09-10T10:00:00.000000000',
           '2006-09-10T11:00:00.000000000', '2006-09-10T12:00:00.000000000',
           '2006-09-10T13:00:00.000000000', '2006-09-10T14:00:00.000000000',
           '2006-09-10T15:00:00.000000000', '2006-09-10T16:00:00.000000000',
           '2006-09-10T17:00:00.000000000', '2006-09-10T18:00:00.000000000',
           '2006-09-10T19:00:00.000000000', '2006-09-10T20:00:00.000000000',
           '2006-09-10T21:00:00.000000000', '2006-09-10T22:00:00.000000000',
           '2006-09-10T23:00:00.000000000', '2006-09-11T00:00:00.000000000',
           '2006-09-11T01:00:00.000000000', '2006-09-11T02:00:00.000000000',
           '2006-09-11T03:00:00.000000000', '2006-09-11T04:00:00.000000000',
           '2006-09-11T05:00:00.000000000', '2006-09-11T06:00:00.000000000',
           '2006-09-11T07:00:00.000000000', '2006-09-11T08:00:00.000000000',
           '2006-09-11T09:00:00.000000000', '2006-09-11T10:00:00.000000000',
           '2006-09-11T11:00:00.000000000', '2006-09-11T12:00:00.000000000',
           '2006-09-11T13:00:00.000000000', '2006-09-11T14:00:00.000000000',
           '2006-09-11T15:00:00.000000000', '2006-09-11T16:00:00.000000000',
           '2006-09-11T17:00:00.000000000', '2006-09-11T18:00:00.000000000',
           '2006-09-11T19:00:00.000000000', '2006-09-11T20:00:00.000000000',
           '2006-09-11T21:00:00.000000000', '2006-09-11T22:00:00.000000000',
           '2006-09-11T23:00:00.000000000', '2006-09-12T00:00:00.000000000',
           '2006-09-12T01:00:00.000000000', '2006-09-12T02:00:00.000000000',
           '2006-09-12T03:00:00.000000000', '2006-09-12T04:00:00.000000000',
           '2006-09-12T05:00:00.000000000', '2006-09-12T06:00:00.000000000',
           '2006-09-12T07:00:00.000000000', '2006-09-12T08:00:00.000000000',
           '2006-09-12T09:00:00.000000000', '2006-09-12T10:00:00.000000000',
           '2006-09-12T11:00:00.000000000', '2006-09-12T12:00:00.000000000',
           '2006-09-12T13:00:00.000000000', '2006-09-12T14:00:00.000000000',
           '2006-09-12T15:00:00.000000000', '2006-09-12T16:00:00.000000000',
           '2006-09-12T17:00:00.000000000', '2006-09-12T18:00:00.000000000',
           '2006-09-12T19:00:00.000000000', '2006-09-12T20:00:00.000000000',
           '2006-09-12T21:00:00.000000000', '2006-09-12T22:00:00.000000000',
           '2006-09-12T23:00:00.000000000', '2006-09-13T00:00:00.000000000',
           '2006-09-13T01:00:00.000000000', '2006-09-13T02:00:00.000000000',
           '2006-09-13T03:00:00.000000000', '2006-09-13T04:00:00.000000000',
           '2006-09-13T05:00:00.000000000', '2006-09-13T06:00:00.000000000',
           '2006-09-13T07:00:00.000000000', '2006-09-13T08:00:00.000000000',
           '2006-09-13T09:00:00.000000000', '2006-09-13T10:00:00.000000000',
           '2006-09-13T11:00:00.000000000', '2006-09-13T12:00:00.000000000',
           '2006-09-13T13:00:00.000000000', '2006-09-13T14:00:00.000000000',
           '2006-09-13T15:00:00.000000000', '2006-09-13T16:00:00.000000000',
           '2006-09-13T17:00:00.000000000', '2006-09-13T18:00:00.000000000'],
          dtype='datetime64[ns]')

  • lon
    (lon)
    float64
    85.0 85.25 85.5 ... 169.8 170.0

    long_name :

    Longitude

    units :

    degrees_east

    array([ 85.  ,  85.25,  85.5 , ..., 169.5 , 169.75, 170.  ])

  • lat
    (lat)
    float64
    -4.875 -4.625 ... 39.62 39.88

    standard_name :

    latitude

    long_name :

    Latitude

    units :

    degrees_north

    axis :

    Y

    array([-4.875, -4.625, -4.375, -4.125, -3.875, -3.625, -3.375, -3.125, -2.875,
           -2.625, -2.375, -2.125, -1.875, -1.625, -1.375, -1.125, -0.875, -0.625,
           -0.375, -0.125,  0.125,  0.375,  0.625,  0.875,  1.125,  1.375,  1.625,
            1.875,  2.125,  2.375,  2.625,  2.875,  3.125,  3.375,  3.625,  3.875,
            4.125,  4.375,  4.625,  4.875,  5.125,  5.375,  5.625,  5.875,  6.125,
            6.375,  6.625,  6.875,  7.125,  7.375,  7.625,  7.875,  8.125,  8.375,
            8.625,  8.875,  9.125,  9.375,  9.625,  9.875, 10.125, 10.375, 10.625,
           10.875, 11.125, 11.375, 11.625, 11.875, 12.125, 12.375, 12.625, 12.875,
           13.125, 13.375, 13.625, 13.875, 14.125, 14.375, 14.625, 14.875, 15.125,
           15.375, 15.625, 15.875, 16.125, 16.375, 16.625, 16.875, 17.125, 17.375,
           17.625, 17.875, 18.125, 18.375, 18.625, 18.875, 19.125, 19.375, 19.625,
           19.875, 20.125, 20.375, 20.625, 20.875, 21.125, 21.375, 21.625, 21.875,
           22.125, 22.375, 22.625, 22.875, 23.125, 23.375, 23.625, 23.875, 24.125,
           24.375, 24.625, 24.875, 25.125, 25.375, 25.625, 25.875, 26.125, 26.375,
           26.625, 26.875, 27.125, 27.375, 27.625, 27.875, 28.125, 28.375, 28.625,
           28.875, 29.125, 29.375, 29.625, 29.875, 30.125, 30.375, 30.625, 30.875,
           31.125, 31.375, 31.625, 31.875, 32.125, 32.375, 32.625, 32.875, 33.125,
           33.375, 33.625, 33.875, 34.125, 34.375, 34.625, 34.875, 35.125, 35.375,
           35.625, 35.875, 36.125, 36.375, 36.625, 36.875, 37.125, 37.375, 37.625,
           37.875, 38.125, 38.375, 38.625, 38.875, 39.125, 39.375, 39.625, 39.875])

• Data variables: (2)
  • ctt
    (time, lat, lon)
    float32
    266.8 266.9 267.0 ... 265.9 265.7

    long_name :

    Brightness temperature

    units :

    K

    array([[[266.8233 , 266.8882 , 266.95795, ..., 268.88196, 269.06625,
             269.22308],
            [266.7498 , 266.91086, 266.94043, ..., 269.02643, 269.17902,
             269.28247],
            [266.86603, 267.0336 , 267.0265 , ..., 269.07828, 269.27002,
             269.40576],
            ...,
            [265.21896, 265.2382 , 265.0005 , ..., 265.66608, 265.7233 ,
             266.0509 ],
            [264.91156, 265.0018 , 264.59222, ..., 266.05298, 266.25204,
             266.3739 ],
            [264.79828, 264.50156, 263.52072, ..., 266.12006, 266.1328 ,
             266.4638 ]],
    
           [[266.18558, 266.23065, 265.9952 , ..., 268.7587 , 268.52835,
             269.0521 ],
            [265.47928, 266.09277, 266.58337, ..., 268.40442, 267.71396,
             268.22873],
            [266.18246, 266.36987, 266.6086 , ..., 268.21286, 269.1558 ,
             269.3137 ],
    ...
            [264.63684, 264.8857 , 265.37686, ..., 265.22318, 265.0068 ,
             264.94736],
            [264.40012, 264.61005, 265.11108, ..., 265.5867 , 265.48117,
             265.38724],
            [264.32895, 264.52274, 264.95453, ..., 266.00293, 265.82333,
             265.5836 ]],
    
           [[259.4579 , 255.84242, 254.06657, ..., 265.8721 , 265.70184,
             265.81763],
            [258.09866, 255.88245, 252.54802, ..., 265.5966 , 265.71997,
             265.9623 ],
            [253.25778, 257.53015, 257.4035 , ..., 265.59103, 265.8222 ,
             265.94907],
            ...,
            [264.2128 , 264.21365, 264.36285, ..., 265.51105, 265.2532 ,
             265.05634],
            [264.06696, 264.01376, 264.12125, ..., 265.6935 , 265.6174 ,
             265.53763],
            [263.98972, 263.95242, 264.03027, ..., 265.9807 , 265.9347 ,
             265.68768]]], dtype=float32)

  • pr
    (time, lat, lon)
    float32
    ...

    long_name :

    Precipitation rate

    units :

    mm h-1

    [7733880 values with dtype=float32]

• Indexes: (3)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2006-09-08 13:00:00', '2006-09-08 14:00:00',
                   '2006-09-08 15:00:00', '2006-09-08 16:00:00',
                   '2006-09-08 17:00:00', '2006-09-08 18:00:00',
                   '2006-09-08 19:00:00', '2006-09-08 20:00:00',
                   '2006-09-08 21:00:00', '2006-09-08 22:00:00',
                   ...
                   '2006-09-13 09:00:00', '2006-09-13 10:00:00',
                   '2006-09-13 11:00:00', '2006-09-13 12:00:00',
                   '2006-09-13 13:00:00', '2006-09-13 14:00:00',
                   '2006-09-13 15:00:00', '2006-09-13 16:00:00',
                   '2006-09-13 17:00:00', '2006-09-13 18:00:00'],
                  dtype='datetime64[ns]', name='time', length=126, freq=None))

  • lon
    PandasIndex

    PandasIndex(Float64Index([  85.0,  85.25,   85.5,  85.75,   86.0,  86.25,   86.5,  86.75,
                    87.0,  87.25,
                  ...
                  167.75,  168.0, 168.25,  168.5, 168.75,  169.0, 169.25,  169.5,
                  169.75,  170.0],
                 dtype='float64', name='lon', length=341))

  • lat
    PandasIndex

    PandasIndex(Float64Index([-4.875, -4.625, -4.375, -4.125, -3.875, -3.625, -3.375, -3.125,
                  -2.875, -2.625,
                  ...
                  37.625, 37.875, 38.125, 38.375, 38.625, 38.875, 39.125, 39.375,
                  39.625, 39.875],
                 dtype='float64', name='lat', length=180))

• Attributes: (0)

Identify CEs using different CTT thresholds

For CEs (not necessarily for MCSs), the CTT threshold only affects identification.

Here we use tams.data_in_contours() to add a few cloud-top temperature and precipitation rate statistics to our CE dataset.

%%time

ctt = ds.ctt
precip = ds.pr

thresh = 235
thresh_core = 219

cases_ces = {}
for delta in [-15, 0, 5]:
    if delta == 0:
        key = rf"default ($T = {thresh}\,\mathrm{{K}}$, $T_\mathrm{{core}} = {thresh_core}\,\mathrm{{K}}$)"
    else:
        s_delta = str(delta) if delta <= 0 else f"+{delta}"
        key = rf"$(T, T_\mathrm{{core}}) {s_delta}$"
    ces = tams.identify(ctt,
        parallel=True,
        ctt_threshold=thresh + delta,
        ctt_core_threshold=thresh_core + delta,
    )[0]

    assert len(ces) == ds.sizes["time"]
    
    # Add min Tb and avg precip -- parallel
    # NOTE: parallel attempt giving can't pickle threadlock error
    def fun(ds, ce):
        ce = tams.data_in_contours(ds.ctt, ce, agg=("mean", "min",), merge=True)
        ce = tams.data_in_contours(ds.pr, ce, agg=("mean",), merge=True)
        return ce
    
    ces_ = Parallel(n_jobs=-2, verbose=10)(
        delayed(fun)(ds.isel(time=i).copy(deep=False), ce.copy())
        for i, ce in enumerate(ces)
    )
    cases_ces[key] = ces_

    # # Add min Tb and avg precip -- serial
    # for i, ce in enumerate(ces):
    #     ce = tams.data_in_contours(ctt.isel(time=i), ce, agg=("min",), merge=True)
    #     ce = tams.data_in_contours(precip.isel(time=i), ce, agg=("mean",), merge=True)
    #     ces[i] = ce
    #     print(".", end="")
    # print()
    # cases_ces[key] = ces

Show code cell output Hide code cell output

[Parallel(n_jobs=-2)]: Using backend LokyBackend with 11 concurrent workers.
[Parallel(n_jobs=-2)]: Done   3 tasks      | elapsed:    0.8s
[Parallel(n_jobs=-2)]: Done  10 tasks      | elapsed:    1.8s
[Parallel(n_jobs=-2)]: Done  19 tasks      | elapsed:    3.6s
[Parallel(n_jobs=-2)]: Done  28 tasks      | elapsed:    4.5s
[Parallel(n_jobs=-2)]: Done  39 tasks      | elapsed:    5.2s
[Parallel(n_jobs=-2)]: Done  50 tasks      | elapsed:    6.2s
[Parallel(n_jobs=-2)]: Done  63 tasks      | elapsed:    7.2s
[Parallel(n_jobs=-2)]: Done  76 tasks      | elapsed:    8.5s
[Parallel(n_jobs=-2)]: Done  91 tasks      | elapsed:    9.6s
[Parallel(n_jobs=-2)]: Done 118 out of 126 | elapsed:   11.6s remaining:    0.7s
[Parallel(n_jobs=-2)]: Done 126 out of 126 | elapsed:   12.2s finished
[Parallel(n_jobs=-2)]: Using backend LokyBackend with 11 concurrent workers.
[Parallel(n_jobs=-2)]: Done   3 tasks      | elapsed:    1.3s
[Parallel(n_jobs=-2)]: Done  10 tasks      | elapsed:    1.4s
[Parallel(n_jobs=-2)]: Done  19 tasks      | elapsed:    2.6s
[Parallel(n_jobs=-2)]: Done  28 tasks      | elapsed:    3.7s
[Parallel(n_jobs=-2)]: Done  39 tasks      | elapsed:    4.9s
[Parallel(n_jobs=-2)]: Done  50 tasks      | elapsed:    6.0s
[Parallel(n_jobs=-2)]: Done  63 tasks      | elapsed:    7.3s
[Parallel(n_jobs=-2)]: Done  76 tasks      | elapsed:    8.6s
[Parallel(n_jobs=-2)]: Done  91 tasks      | elapsed:   10.6s
[Parallel(n_jobs=-2)]: Done 118 out of 126 | elapsed:   13.5s remaining:    0.8s
[Parallel(n_jobs=-2)]: Done 126 out of 126 | elapsed:   14.2s finished
[Parallel(n_jobs=-2)]: Using backend LokyBackend with 11 concurrent workers.
[Parallel(n_jobs=-2)]: Done   3 tasks      | elapsed:    1.1s
[Parallel(n_jobs=-2)]: Done  10 tasks      | elapsed:    1.3s
[Parallel(n_jobs=-2)]: Done  19 tasks      | elapsed:    2.8s
[Parallel(n_jobs=-2)]: Done  28 tasks      | elapsed:    4.3s
[Parallel(n_jobs=-2)]: Done  39 tasks      | elapsed:    5.9s
[Parallel(n_jobs=-2)]: Done  50 tasks      | elapsed:    7.6s
[Parallel(n_jobs=-2)]: Done  63 tasks      | elapsed:    9.5s
[Parallel(n_jobs=-2)]: Done  76 tasks      | elapsed:   11.5s
[Parallel(n_jobs=-2)]: Done  91 tasks      | elapsed:   13.8s
[Parallel(n_jobs=-2)]: Done 118 out of 126 | elapsed:   17.7s remaining:    1.1s
[Parallel(n_jobs=-2)]: Done 126 out of 126 | elapsed:   18.4s finished
[Parallel(n_jobs=-2)]: Using backend LokyBackend with 11 concurrent workers.
[Parallel(n_jobs=-2)]: Done   3 tasks      | elapsed:    1.6s
[Parallel(n_jobs=-2)]: Done  10 tasks      | elapsed:    1.9s
[Parallel(n_jobs=-2)]: Done  19 tasks      | elapsed:    3.6s
[Parallel(n_jobs=-2)]: Done  28 tasks      | elapsed:    4.9s
[Parallel(n_jobs=-2)]: Done  39 tasks      | elapsed:    6.3s
[Parallel(n_jobs=-2)]: Done  50 tasks      | elapsed:    8.1s
[Parallel(n_jobs=-2)]: Done  63 tasks      | elapsed:    9.8s
[Parallel(n_jobs=-2)]: Done  76 tasks      | elapsed:   11.7s
[Parallel(n_jobs=-2)]: Done  91 tasks      | elapsed:   14.6s
[Parallel(n_jobs=-2)]: Done 118 out of 126 | elapsed:   18.2s remaining:    1.1s
[Parallel(n_jobs=-2)]: Done 126 out of 126 | elapsed:   19.2s finished
[Parallel(n_jobs=-2)]: Using backend LokyBackend with 11 concurrent workers.
[Parallel(n_jobs=-2)]: Done   3 tasks      | elapsed:    1.2s
[Parallel(n_jobs=-2)]: Done  10 tasks      | elapsed:    1.3s
[Parallel(n_jobs=-2)]: Done  19 tasks      | elapsed:    2.9s
[Parallel(n_jobs=-2)]: Done  28 tasks      | elapsed:    4.7s
[Parallel(n_jobs=-2)]: Done  39 tasks      | elapsed:    6.4s
[Parallel(n_jobs=-2)]: Done  50 tasks      | elapsed:    8.3s
[Parallel(n_jobs=-2)]: Done  63 tasks      | elapsed:   10.3s
[Parallel(n_jobs=-2)]: Done  76 tasks      | elapsed:   12.3s
[Parallel(n_jobs=-2)]: Done  91 tasks      | elapsed:   14.5s
[Parallel(n_jobs=-2)]: Done 118 out of 126 | elapsed:   18.9s remaining:    1.2s
[Parallel(n_jobs=-2)]: Done 126 out of 126 | elapsed:   19.8s finished
[Parallel(n_jobs=-2)]: Using backend LokyBackend with 11 concurrent workers.
[Parallel(n_jobs=-2)]: Done   3 tasks      | elapsed:    1.9s
[Parallel(n_jobs=-2)]: Done  10 tasks      | elapsed:    2.1s
[Parallel(n_jobs=-2)]: Done  19 tasks      | elapsed:    3.8s
[Parallel(n_jobs=-2)]: Done  28 tasks      | elapsed:    5.3s
[Parallel(n_jobs=-2)]: Done  39 tasks      | elapsed:    7.1s
[Parallel(n_jobs=-2)]: Done  50 tasks      | elapsed:    8.9s
[Parallel(n_jobs=-2)]: Done  63 tasks      | elapsed:   10.5s
[Parallel(n_jobs=-2)]: Done  76 tasks      | elapsed:   12.7s
[Parallel(n_jobs=-2)]: Done  91 tasks      | elapsed:   15.2s
[Parallel(n_jobs=-2)]: Done 118 out of 126 | elapsed:   19.3s remaining:    1.2s

CPU times: total: 7.64 s
Wall time: 1min 44s

[Parallel(n_jobs=-2)]: Done 126 out of 126 | elapsed:   20.2s finished

Track for different \(u\) projection values

For each of our CTT threshold cases, we track using different \(u\) projection values, bringing the case count to 9.

%%time

cases = {}
time = ctt.time.values.tolist()

projs = [-5, -10, -12]
for proj in projs:
    for thresh_key, ces in cases_ces.items():
        key = rf"{thresh_key}, $u_\mathrm{{proj}} = {proj}\,\mathrm{{m}}\,\mathrm{{s}}^{{-1}}$"
        cases[key] = tams.track(ces, time, u_projection=proj)

CPU times: total: 1min 40s
Wall time: 1min 40s

Check the case labels

fig, ax = plt.subplots(figsize=(4, 3))

ax.axis("off")

lines = []
for i, k in enumerate(cases):
    lines.append(f"{i+1}. `{k}`")
    y = len(cases) - i - 0.5
    ax.text(0, y, i + 1, va="center", ha="left", color="orangered")
    ax.text(0.07, y, k, va="center", ha="left")

ax.set_ylim(ymax=len(cases))

display(Markdown("\n\n".join(lines)))

1. $(T, T_\mathrm{core}) -15$, $u_\mathrm{proj} = -5\,\mathrm{m}\,\mathrm{s}^{-1}$

2. default ($T = 235\,\mathrm{K}$, $T_\mathrm{core} = 219\,\mathrm{K}$), $u_\mathrm{proj} = -5\,\mathrm{m}\,\mathrm{s}^{-1}$

3. $(T, T_\mathrm{core}) +5$, $u_\mathrm{proj} = -5\,\mathrm{m}\,\mathrm{s}^{-1}$

4. $(T, T_\mathrm{core}) -15$, $u_\mathrm{proj} = -10\,\mathrm{m}\,\mathrm{s}^{-1}$

5. default ($T = 235\,\mathrm{K}$, $T_\mathrm{core} = 219\,\mathrm{K}$), $u_\mathrm{proj} = -10\,\mathrm{m}\,\mathrm{s}^{-1}$

6. $(T, T_\mathrm{core}) +5$, $u_\mathrm{proj} = -10\,\mathrm{m}\,\mathrm{s}^{-1}$

7. $(T, T_\mathrm{core}) -15$, $u_\mathrm{proj} = -12\,\mathrm{m}\,\mathrm{s}^{-1}$

8. default ($T = 235\,\mathrm{K}$, $T_\mathrm{core} = 219\,\mathrm{K}$), $u_\mathrm{proj} = -12\,\mathrm{m}\,\mathrm{s}^{-1}$

9. $(T, T_\mathrm{core}) +5$, $u_\mathrm{proj} = -12\,\mathrm{m}\,\mathrm{s}^{-1}$

Plot area dist

m = len(projs)
n = 3
assert m * n == len(cases)

fig, axs = plt.subplots(m, n, figsize=(8.2, 7), sharex=True, sharey=True, constrained_layout=True)

for a, (key, ce), ax in zip(ascii_lowercase, cases.items(), axs.flat):
    n = ce.mcs_id.nunique()
    # Dist of max area
    x = ce.groupby("mcs_id").area_km2.max()
    x.plot.kde(ax=ax, label="CE")
    
    mcs = ce.groupby(["mcs_id", "itime"]).area_km2.sum()
    x2 = mcs.groupby("mcs_id").max()
    x2.plot.kde(ax=ax, label="MCS")
    
    # ax.text(0.99, 0.98, f"$N={n}$", size=9, ha="right", va="top", transform=ax.transAxes)
    ax.text(0.99, 0.03, f"$N={n}$", size=9, ha="right", va="bottom", transform=ax.transAxes)
    
    # ax.set_title(key, size=9)
    i = key.index(", $u_")
    l1, l2 = key[:i], key[i+2:]
    if ax.get_subplotspec().is_first_col():
        xt, yt = 0.02, 0.22
    else:
        xt, yt = 0.1, 0.97
    ax.text(xt, yt, l1, size=10, ha="left", va="top", transform=ax.transAxes)
    ax.text(xt, yt - 0.09, l2, size=10, ha="left", va="top", transform=ax.transAxes)
    
    ax.set_ylabel("")
    ax.text(0.02, 0.97, f"{a}",
        size=12, weight="bold",
        ha="left", va="top",
        transform=ax.transAxes,
    )
    if a == "c":
        ax.legend(loc="center right")

fig.supxlabel("Area [km$^2$]")
fig.supylabel("Density", x=-0.03)

# ax.set_xlim(xmin=0)
ax.set_xlim(xmin=1000); ax.set_xscale("log")
ax.set_ylim(ymin=0)

fig.savefig("mpas-tracking-options-area-kde.pdf", bbox_inches="tight", pad_inches=0.05, transparent=False)

Plot 2-D duration and area dist

m = len(projs)
n = 3
assert m * n == len(cases)

fig, axs = plt.subplots(m, n, figsize=(8.2, 7), sharex=True, sharey=True, constrained_layout=True)

for a, (key, ce), ax in zip(ascii_lowercase, cases.items(), axs.flat):
    n = ce.mcs_id.nunique()
    gb = ce.groupby("mcs_id")
    area = gb.area_km2.max()
    dur = ((gb.time.max() - gb.time.min()).dt.total_seconds() / 3600).rename("duration_h")
    data = pd.concat([area, dur], axis="columns")
    
    # sns.jointplot(x="area_km2", y="duration_h", kind="kde", joint_kws=dict(fill=True), ax=ax, data=data)
    sns.kdeplot(x="area_km2", y="duration_h",
        fill=True,
        common_norm=True, clip=(0, None), hue_norm=(0, 0.5),
        ax=ax, data=data,
    )
    # TODO: ensure same levels?
    # sns.rugplot(x="area_km2", y="duration_h", ax=ax, data=data)
    
    #ax.text(0.99, 0.98, f"$N={n}$", size=9, ha="right", va="top", transform=ax.transAxes)
    ax.text(0.99, 0.03, f"$N={n}$",
        size=9,
        ha="right", va="bottom",
        path_effects=[path_effects.withStroke(linewidth=1.5, foreground="w")],
        transform=ax.transAxes,
    )
   
    # ax.set_title(key, size=10)
    i = key.index(", $u_")
    l1, l2 = key[:i], key[i+2:]
    xt, yt = 0.09, 0.975
    kws = dict(
        size=10,
        ha="left", va="top",
        path_effects=[path_effects.withStroke(linewidth=1.5, foreground="w")],
        transform=ax.transAxes,
    )
    ax.text(xt, yt, l1, **kws)
    ax.text(xt, yt - 0.09, l2, **kws)
    
    ax.set_xlabel("")
    ax.set_ylabel("")
    ax.text(0.02, 0.97, f"{a}",
        size=12, weight="bold",
        ha="left", va="top",
        path_effects=[path_effects.withStroke(linewidth=1.5, foreground="w")],
        transform=ax.transAxes,
    )

fig.supxlabel("Area [km$^2$]")
fig.supylabel("Duration [h]", x=-0.03)

ax.set_xlim(xmin=0, xmax=0.5e7)
ax.set_ylim(ymin=0, ymax=30)

fig.savefig("mpas-tracking-options-duration-area-2d-kde.pdf", bbox_inches="tight", pad_inches=0.05, transparent=False)

list(cases.values())[0].head(3)

	geometry	area_km2	inds219	area219_km2	cs219	mean_ctt	min_ctt	mean_pr	time	itime	dtime	mcs_id
0	POLYGON ((152.50000 8.93176, 154.09671 11.1250...	65189.730096	[1, 2, 3, 4]	7632.145729	MULTIPOLYGON (((154.00000 11.33169, 154.00332 ...	223.781916	199.467636	2.174953	2006-09-08 13:00:00	0	0 days 01:00:00	0
1	POLYGON ((105.00000 20.86816, 105.04808 20.875...	43468.589233	[8]	11043.783834	MULTIPOLYGON (((105.00000 21.04446, 105.25000 ...	217.004782	196.3992	3.237674	2006-09-08 13:00:00	0	0 days 01:00:00	1
2	POLYGON ((103.75000 0.07912, 103.83888 0.12500...	98642.412303	[9, 10, 11]	38699.492801	MULTIPOLYGON (((103.25000 1.66696, 103.28213 1...	216.819586	195.909149	3.098599	2006-09-08 13:00:00	0	0 days 01:00:00	2

Plot 2-D CTT and precip dist

m = len(projs)
n = 3
assert m * n == len(cases)

fig, axs = plt.subplots(m, n, figsize=(8.3, 6.8), sharex=True, sharey=True, constrained_layout=True)

for a, (key, ce), ax in zip(ascii_lowercase, cases.items(), axs.flat):    
    n = ce.mcs_id.nunique()
    gb = ce.groupby("mcs_id")
    tb = gb.min_ctt.min().astype(float)  # np.sqrt not supported for `Float64`
    pr = gb.mean_pr.mean().astype(float)  # TODO: area-weighted average of CE precip at each time?
    data = pd.concat([tb, pr], axis="columns")

    ax.plot(tb, pr, ".", ms=5, mec="none", mfc="0.35", alpha=0.9)
    
    sns.kdeplot(x="min_ctt", y="mean_pr",
        fill=False, color=sns.color_palette()[0], alpha=0.6,
        common_norm=True, clip=(0, None),
        ax=ax, data=data,
    )
    
    ax.text(0.99, 0.03, f"$N={n}$",
        size=9,
        ha="right", va="bottom",
        path_effects=[path_effects.withStroke(linewidth=1.5, foreground="w")],
        transform=ax.transAxes,
    )
   
    i = key.index(", $u_")
    l1, l2 = key[:i], key[i+2:]
    xt, yt = 0.09, 0.975
    kws = dict(
        size=10,
        ha="left", va="top",
        path_effects=[path_effects.withStroke(linewidth=1.5, foreground="w")],
        transform=ax.transAxes,
    )
    ax.text(xt, yt, l1, **kws)
    ax.text(xt, yt - 0.09, l2, **kws)
    
    ax.set_xlabel("")
    ax.set_ylabel("")
    ax.text(0.02, 0.97, f"{a}",
        size=12, weight="bold",
        ha="left", va="top",
        path_effects=[path_effects.withStroke(linewidth=1.5, foreground="w")],
        transform=ax.transAxes,
    )

fig.supxlabel("Minimum CTT [K]")
fig.supylabel("Average precipitation rate [mm h$^{-1}$]", x=-0.03)

ax.set_xlim(xmin=185, xmax=230)
ax.set_ylim(ymin=0, ymax=4.5)

fig.savefig("mpas-tracking-options-tb-pr-2d-kde.pdf", bbox_inches="tight", pad_inches=0.05, transparent=False)