data_file = "../Simulations/toy_data_01.nc"

import hextraj

import xarray as xr
from matplotlib import pyplot as plt
import cartopy

Load data and check raw data

ds_traj = xr.open_dataset(data_file)
ds_traj

<xarray.Dataset>
Dimensions:     (traj: 144, obs: 121)
Dimensions without coordinates: traj, obs
Data variables:
    trajectory  (traj, obs) float64 ...
    time        (traj, obs) datetime64[ns] ...
    lat         (traj, obs) float32 ...
    lon         (traj, obs) float32 ...
    z           (traj, obs) float32 ...
    U           (traj, obs) float32 ...
    V           (traj, obs) float32 ...
Attributes:
    feature_type:           trajectory
    Conventions:            CF-1.6/CF-1.7
    ncei_template_version:  NCEI_NetCDF_Trajectory_Template_v2.0
    parcels_version:        2.3.1.dev20+g92f2fb90
    parcels_mesh:           spherical

xarray.Dataset

• Dimensions:
  • traj: 144
  • obs: 121
• Coordinates: (0)
• Data variables: (7)
  • trajectory
    (traj, obs)
    float64
    ...

    long_name :

    Unique identifier for each particle

    cf_role :

    trajectory_id

    [17424 values with dtype=float64]

  • time
    (traj, obs)
    datetime64[ns]
    ...

    long_name :

    standard_name :

    time

    axis :

    T

    [17424 values with dtype=datetime64[ns]]

  • lat
    (traj, obs)
    float32
    ...

    long_name :

    standard_name :

    latitude

    units :

    degrees_north

    axis :

    Y

    [17424 values with dtype=float32]

  • lon
    (traj, obs)
    float32
    ...

    long_name :

    standard_name :

    longitude

    units :

    degrees_east

    axis :

    X

    [17424 values with dtype=float32]

  • z
    (traj, obs)
    float32
    ...

    long_name :

    standard_name :

    depth

    units :

    m

    positive :

    down

    [17424 values with dtype=float32]

  • U
    (traj, obs)
    float32
    ...

    long_name :

    standard_name :

    U

    units :

    unknown

    [17424 values with dtype=float32]

  • V
    (traj, obs)
    float32
    ...

    long_name :

    standard_name :

    V

    units :

    unknown

    [17424 values with dtype=float32]

• Indexes: (0)
• Attributes: (5)

  feature_type :

  trajectory

  Conventions :

  CF-1.6/CF-1.7

  ncei_template_version :

  NCEI_NetCDF_Trajectory_Template_v2.0

  parcels_version :

  2.3.1.dev20+g92f2fb90

  parcels_mesh :

  spherical

fig, ax = plt.subplots(
    1, 1,
    subplot_kw={"projection": cartopy.crs.Stereographic(
        central_latitude=ds_traj.lat.mean().data[()],
        central_longitude=ds_traj.lon.mean().data[()]
    )},
    figsize=(4, 4),
)
ds_traj.to_dataframe().plot.scatter(x="lon", y="lat", s=2, alpha=0.1, ax=ax, transform=cartopy.crs.PlateCarree())
ax.coastlines()
ax.gridlines(draw_labels=True);

Create hex labeller

This uses a Lamber azimuthal equal-area projection which assures equal area of all hexagons and preserves angles between all vectors pointing away from the center of the projection.

We center the projection about the middle point of all trajectories.

lon_origin = ((ds_traj.lon.max() + ds_traj.lon.min()) / 2.0).data[()]
lat_origin = ((ds_traj.lat.max() + ds_traj.lat.min()) / 2.0).data[()]
hex_size_meters = 200_000

hex_proj = hextraj.HexProj(
    lon_origin=lon_origin,
    lat_origin=lat_origin,
    hex_size_meters=hex_size_meters,
)

hex_proj

/srv/conda/envs/notebook/lib/python3.10/site-packages/pyproj/crs/crs.py:141: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
  in_crs_string = _prepare_from_proj_string(in_crs_string)

HexProj(projection_name='laea', lon_origin=29.791340827941895, lat_origin=-38.089969635009766, hex_size_meters=200000, hex_orientation='flat', )

Compute labels

Labels are three element tuples.

ds_traj["hex_labels"] = xr.apply_ufunc(
    hex_proj.lon_lat_to_hex_AoS,
    ds_traj.lon,
    ds_traj.lat,
)

display(ds_traj)

<xarray.Dataset>
Dimensions:     (traj: 144, obs: 121)
Dimensions without coordinates: traj, obs
Data variables:
    trajectory  (traj, obs) float64 0.0 0.0 0.0 0.0 ... 143.0 143.0 143.0 143.0
    time        (traj, obs) datetime64[ns] 2021-01-01T00:30:00 ... 2021-01-31...
    lat         (traj, obs) float32 -39.17 -39.02 -38.9 ... -39.19 -39.21 -39.25
    lon         (traj, obs) float32 22.5 22.55 22.61 22.66 ... 17.16 17.18 17.17
    z           (traj, obs) float32 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    U           (traj, obs) float32 nan 2.598e-06 ... 1.845e-07 -7.776e-07
    V           (traj, obs) float32 nan 5.796e-06 ... -1.211e-06 -2.669e-06
    hex_labels  (traj, obs) object (-2, 1, 1) (-2, 1, 1) ... (-4, 1, 3)
Attributes:
    feature_type:           trajectory
    Conventions:            CF-1.6/CF-1.7
    ncei_template_version:  NCEI_NetCDF_Trajectory_Template_v2.0
    parcels_version:        2.3.1.dev20+g92f2fb90
    parcels_mesh:           spherical

xarray.Dataset

• Dimensions:
  • traj: 144
  • obs: 121
• Coordinates: (0)
• Data variables: (8)
  • trajectory
    (traj, obs)
    float64
    0.0 0.0 0.0 ... 143.0 143.0 143.0

    long_name :

    Unique identifier for each particle

    cf_role :

    trajectory_id

    array([[  0.,   0.,   0., ...,   0.,   0.,   0.],
           [  1.,   1.,   1., ...,   1.,   1.,   1.],
           [  2.,   2.,   2., ...,   2.,   2.,   2.],
           ...,
           [141., 141., 141., ..., 141., 141., 141.],
           [142., 142., 142., ..., 142., 142., 142.],
           [143., 143., 143., ..., 143., 143., 143.]])

  • time
    (traj, obs)
    datetime64[ns]
    2021-01-01T00:30:00 ... 2021-01-...

    long_name :

    standard_name :

    time

    axis :

    T

    array([['2021-01-01T00:30:00.000000000', '2021-01-01T06:30:00.000000000',
            '2021-01-01T12:30:00.000000000', ..., '2021-01-30T12:30:00.000000000',
            '2021-01-30T18:30:00.000000000', '2021-01-31T00:30:00.000000000'],
           ['2021-01-01T00:30:00.000000000', '2021-01-01T06:30:00.000000000',
            '2021-01-01T12:30:00.000000000', ..., '2021-01-30T12:30:00.000000000',
            '2021-01-30T18:30:00.000000000', '2021-01-31T00:30:00.000000000'],
           ['2021-01-01T00:30:00.000000000', '2021-01-01T06:30:00.000000000',
            '2021-01-01T12:30:00.000000000', ..., '2021-01-30T12:30:00.000000000',
            '2021-01-30T18:30:00.000000000', '2021-01-31T00:30:00.000000000'],
           ...,
           ['2021-01-01T00:30:00.000000000', '2021-01-01T06:30:00.000000000',
            '2021-01-01T12:30:00.000000000', ..., '2021-01-30T12:30:00.000000000',
            '2021-01-30T18:30:00.000000000', '2021-01-31T00:30:00.000000000'],
           ['2021-01-01T00:30:00.000000000', '2021-01-01T06:30:00.000000000',
            '2021-01-01T12:30:00.000000000', ..., '2021-01-30T12:30:00.000000000',
            '2021-01-30T18:30:00.000000000', '2021-01-31T00:30:00.000000000'],
           ['2021-01-01T00:30:00.000000000', '2021-01-01T06:30:00.000000000',
            '2021-01-01T12:30:00.000000000', ..., '2021-01-30T12:30:00.000000000',
            '2021-01-30T18:30:00.000000000', '2021-01-31T00:30:00.000000000']],
          dtype='datetime64[ns]')

  • lat
    (traj, obs)
    float32
    -39.17 -39.02 ... -39.21 -39.25

    long_name :

    standard_name :

    latitude

    units :

    degrees_north

    axis :

    Y

    array([[-39.166668, -39.023876, -38.904987, ..., -37.456646, -37.43389 ,
            -37.488907],
           [-39.166668, -39.055485, -39.003204, ..., -38.98379 , -38.98158 ,
            -38.997955],
           [-39.166668, -39.209393, -39.376793, ..., -40.00204 , -40.010117,
            -40.02809 ],
           ...,
           [-34.583332, -34.773434, -35.00076 , ..., -34.65815 , -34.588226,
            -34.52226 ],
           [-34.583332, -34.66703 , -34.7913  , ..., -39.15673 , -39.334305,
            -39.513103],
           [-34.583332, -34.587185, -34.60994 , ..., -39.193733, -39.211082,
            -39.25398 ]], dtype=float32)

  • lon
    (traj, obs)
    float32
    22.5 22.55 22.61 ... 17.18 17.17

    long_name :

    standard_name :

    longitude

    units :

    degrees_east

    axis :

    X

    array([[22.5     , 22.547514, 22.607336, ..., 26.906063, 26.797806, 26.698011],
           [22.916666, 23.089357, 23.296543, ..., 27.901531, 28.012968, 28.164171],
           [23.333334, 23.61503 , 23.908384, ..., 27.667067, 27.645962, 27.643606],
           ...,
           [26.25    , 25.866007, 25.449963, ..., 17.211655, 17.284988, 17.336267],
           [26.666666, 26.344868, 25.925909, ..., 13.803558, 13.742613, 13.677705],
           [27.083334, 26.926388, 26.71817 , ..., 17.15514 , 17.17541 , 17.166965]],
          dtype=float32)

  • z
    (traj, obs)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    long_name :

    standard_name :

    depth

    units :

    m

    positive :

    down

    array([[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.]], dtype=float32)

  • U
    (traj, obs)
    float32
    nan 2.598e-06 ... -7.776e-07

    long_name :

    standard_name :

    U

    units :

    unknown

    array([[          nan,  2.597646e-06,  2.767847e-06, ..., -3.317469e-06,
            -5.994304e-06, -2.201192e-06],
           [          nan,  8.876400e-06,  1.002607e-05, ...,  5.008096e-06,
             5.821934e-06,  7.830603e-06],
           [          nan,  1.393304e-05,  1.235174e-05, ..., -8.949883e-07,
            -6.879567e-07,  1.717140e-07],
           ...,
           [          nan, -1.905466e-05, -1.861986e-05, ...,  3.830198e-06,
             2.962539e-06,  1.925747e-06],
           [          nan, -1.713974e-05, -2.094000e-05, ..., -2.493254e-06,
            -3.214069e-06, -2.799593e-06],
           [          nan, -8.068844e-06, -1.109233e-05, ...,  1.658384e-06,
             1.844593e-07, -7.776144e-07]], dtype=float32)

  • V
    (traj, obs)
    float32
    nan 5.796e-06 ... -2.669e-06

    long_name :

    standard_name :

    V

    units :

    unknown

    array([[          nan,  5.795511e-06,  5.306074e-06, ...,  1.302789e-06,
            -1.119043e-07, -4.836137e-06],
           [          nan,  3.527368e-06,  1.323293e-06, ...,  9.175749e-07,
            -6.233992e-07, -4.962131e-07],
           [          nan, -4.892519e-06, -1.079000e-05, ...,  1.360431e-07,
            -8.286600e-07, -6.589117e-07],
           ...,
           [          nan, -8.928485e-06, -1.232287e-05, ...,  3.249278e-06,
             3.223366e-06,  2.770371e-06],
           [          nan, -4.034645e-06, -8.001747e-06, ..., -7.918349e-06,
            -8.593543e-06, -7.880011e-06],
           [          nan,  5.917896e-08, -2.591159e-06, ..., -6.002495e-07,
            -1.211249e-06, -2.668733e-06]], dtype=float32)

  • hex_labels
    (traj, obs)
    object
    (-2, 1, 1) ... (-4, 1, 3)

    array([[(-2, 1, 1), (-2, 1, 1), (-2, 1, 1), ..., (-1, 1, 0), (-1, 1, 0),
            (-1, 1, 0)],
           [(-2, 1, 1), (-2, 1, 1), (-2, 1, 1), ..., (-1, 0, 1), (-1, 0, 1),
            (0, 0, 0)],
           [(-2, 1, 1), (-2, 1, 1), (-2, 1, 1), ..., (-1, 0, 1), (-1, 0, 1),
            (-1, 0, 1)],
           ...,
           [(-1, 2, -1), (-1, 2, -1), (-1, 1, 0), ..., (-4, 3, 1),
            (-4, 3, 1), (-4, 3, 1)],
           [(-1, 2, -1), (-1, 2, -1), (-1, 2, -1), ..., (-5, 2, 3),
            (-5, 2, 3), (-5, 2, 3)],
           [(-1, 2, -1), (-1, 2, -1), (-1, 2, -1), ..., (-4, 1, 3),
            (-3, 1, 2), (-4, 1, 3)]], dtype=object)

• Indexes: (0)
• Attributes: (5)

  feature_type :

  trajectory

  Conventions :

  CF-1.6/CF-1.7

  ncei_template_version :

  NCEI_NetCDF_Trajectory_Template_v2.0

  parcels_version :

  2.3.1.dev20+g92f2fb90

  parcels_mesh :

  spherical

Aggregate data into hexagons

Create empty dataset:

data_agg_hex = xr.Dataset()
data_agg_hex

<xarray.Dataset>
Dimensions:  ()
Data variables:
    *empty*

xarray.Dataset

• Dimensions:
• Coordinates: (0)
• Data variables: (0)
• Indexes: (0)
• Attributes: (0)

Aggregate particle position count:

data_agg_hex["particle_obs_count"] = (
    ds_traj.lon
    .groupby(ds_traj["hex_labels"]).count()
)
data_agg_hex["particle_obs_count"].attrs["long_name"] = "number of positions recorded"
data_agg_hex["particle_obs_count"].attrs["units"] = None

Aggregate U, V and squared velocities:

data_agg_hex["U"] = (
    ds_traj.U
    .groupby(ds_traj["hex_labels"]).mean()
)

data_agg_hex["U2"] = (
    (ds_traj.U ** 2)
    .groupby(ds_traj["hex_labels"]).mean()
)

data_agg_hex["V"] = (
    ds_traj.V
    .groupby(ds_traj["hex_labels"]).mean()
)

data_agg_hex["V2"] = (
    (ds_traj.V ** 2)
    .groupby(ds_traj["hex_labels"]).mean()
)

Add positions of the hex labels:

_lon, _lat = hex_proj.hex_to_lon_lat_SoA(hextraj.aux.hex_AoS_to_SoA(data_agg_hex["hex_labels"].data))

data_agg_hex = data_agg_hex.assign_coords(lon=("hex_labels", _lon))
data_agg_hex = data_agg_hex.assign_coords(lat=("hex_labels", _lat))

data_agg_hex

<xarray.Dataset>
Dimensions:             (hex_labels: 28)
Coordinates:
  * hex_labels          (hex_labels) object (-5, 2, 3) ... (5, -3, -2)
    lon                 (hex_labels) float64 12.48 13.18 15.59 ... 43.39 47.1
    lat                 (hex_labels) float64 -38.38 -35.34 ... -37.3 -38.38
Data variables:
    particle_obs_count  (hex_labels) int64 82 27 206 635 311 ... 67 69 26 61 36
    U                   (hex_labels) float32 -1.309e-06 -1.814e-06 ... 2.058e-06
    U2                  (hex_labels) float32 6.788e-12 8.051e-12 ... 7.148e-12
    V                   (hex_labels) float32 -6.445e-06 -5.26e-07 ... -3.678e-06
    V2                  (hex_labels) float32 5.543e-11 1.045e-11 ... 1.404e-11

xarray.Dataset

• Dimensions:
  • hex_labels: 28
• Coordinates: (3)
  • hex_labels
    (hex_labels)
    object
    (-5, 2, 3) ... (5, -3, -2)

    array([(-5, 2, 3), (-5, 3, 2), (-4, 1, 3), (-4, 2, 2), (-4, 3, 1), (-3, 0, 3),
           (-3, 1, 2), (-3, 2, 1), (-2, 0, 2), (-2, 1, 1), (-2, 2, 0), (-1, -1, 2),
           (-1, 0, 1), (-1, 1, 0), (-1, 2, -1), (0, -1, 1), (0, 0, 0), (0, 1, -1),
           (1, -1, 0), (1, 0, -1), (1, 1, -2), (2, -2, 0), (2, -1, -1),
           (3, -2, -1), (3, -1, -2), (4, -3, -1), (4, -2, -2), (5, -3, -2)],
          dtype=object)

  • lon
    (hex_labels)
    float64
    12.48 13.18 15.59 ... 43.39 47.1

    array([12.48211915, 13.17564421, 15.58522063, 16.18882019, 16.71664805,
           18.85085713, 19.34302822, 19.77168759, 22.65244452, 22.96150975,
           23.23112543, 26.13090548, 26.29776831, 26.44283289, 26.56961643,
           29.79134083, 29.79134083, 29.79134083, 33.28491335, 33.13984877,
           33.01306523, 36.93023714, 36.62117191, 40.23965344, 39.81099406,
           43.99746102, 43.39386147, 47.1005625 ])

  • lat
    (hex_labels)
    float64
    -38.38 -35.34 ... -37.3 -38.38

    array([-38.38466406, -35.33513829, -40.37239933, -37.3013092 ,
           -34.22263832, -42.28263703, -39.19207996, -36.09734867,
           -40.99998013, -37.89202807, -34.78102815, -42.71644941,
           -39.59920826, -36.48119369, -33.36020696, -41.21047358,
           -38.08996964, -34.96798938, -39.59920826, -36.48119369,
           -33.36020696, -40.99998013, -37.89202807, -39.19207996,
           -36.09734867, -40.37239933, -37.3013092 , -38.38466406])

• Data variables: (5)
  • particle_obs_count
    (hex_labels)
    int64
    82 27 206 635 311 ... 69 26 61 36

    long_name :

    number of positions recorded

    standard_name :

    longitude

    units :

    None

    axis :

    X

    array([  82,   27,  206,  635,  311,   23, 1511, 1367,  625, 3716,  857,
            213, 1087, 2776,   25,   82, 1098,  299,  623,  821,   31,  137,
            613,   67,   69,   26,   61,   36])

  • U
    (hex_labels)
    float32
    -1.309e-06 -1.814e-06 ... 2.058e-06

    long_name :

    standard_name :

    U

    units :

    unknown

    array([-1.3086674e-06, -1.8142053e-06,  3.1030866e-06, -3.4642781e-06,
           -1.3312723e-06,  2.3717862e-06, -1.6678820e-06, -1.6510704e-06,
            7.8102657e-06, -1.4152438e-06, -4.8500997e-06,  2.8390896e-06,
            7.9041536e-07,  1.2421006e-06, -7.8585263e-06,  4.8709571e-06,
            6.2151644e-06,  1.6521816e-06,  1.6391774e-06,  5.1245765e-06,
           -3.0891438e-06,  4.8839220e-06,  2.7710628e-06,  1.1500605e-05,
            8.6671453e-06,  3.7185387e-06,  1.0577953e-05,  2.0579164e-06],
          dtype=float32)

  • U2
    (hex_labels)
    float32
    6.788e-12 8.051e-12 ... 7.148e-12

    array([6.7883620e-12, 8.0507337e-12, 3.2911462e-11, 3.8786616e-11,
           1.3031889e-11, 5.0984768e-11, 9.8177154e-11, 1.7780407e-11,
           9.7843601e-11, 2.0067465e-11, 4.3955558e-11, 1.8952912e-11,
           1.7511721e-11, 1.7034789e-11, 9.6907149e-11, 4.7802044e-11,
           6.5141670e-11, 5.9736946e-12, 2.1065517e-11, 4.5758477e-11,
           1.0288130e-11, 2.7704405e-11, 1.3757575e-11, 1.4102113e-10,
           8.0379085e-11, 3.4171214e-11, 1.2145748e-10, 7.1475542e-12],
          dtype=float32)

  • V
    (hex_labels)
    float32
    -6.445e-06 -5.26e-07 ... -3.678e-06

    long_name :

    standard_name :

    V

    units :

    unknown

    array([-6.4445285e-06, -5.2600450e-07, -1.9238626e-06,  1.3216653e-06,
            8.4195398e-07, -1.2551959e-06, -3.3487247e-06, -2.2633100e-07,
            1.3501857e-06, -7.0820477e-08, -1.9294000e-06, -1.6830599e-07,
           -3.8451827e-07,  1.2026364e-07, -1.4148388e-06, -1.2328758e-06,
           -6.6794433e-07, -5.6924381e-07, -4.0447514e-07,  1.0407711e-06,
            1.5914960e-06,  2.3512037e-07, -7.2813577e-08, -5.8157211e-06,
           -2.9786800e-06, -1.4632735e-06,  2.0550453e-06, -3.6782803e-06],
          dtype=float32)

  • V2
    (hex_labels)
    float32
    5.543e-11 1.045e-11 ... 1.404e-11

    array([5.5430365e-11, 1.0450220e-11, 2.3424446e-11, 1.6535875e-11,
           5.4810436e-12, 2.7552101e-11, 2.3563654e-11, 6.9954433e-12,
           2.8270305e-11, 1.7667829e-11, 1.0027489e-11, 5.0818893e-12,
           1.7434745e-11, 8.6656680e-12, 1.3146766e-11, 9.0282365e-12,
           2.5095318e-11, 3.2847839e-12, 1.4440381e-11, 1.5384017e-11,
           3.2383501e-12, 3.2777459e-12, 3.7399920e-11, 5.4875892e-11,
           4.8804148e-11, 1.0631654e-11, 1.6860277e-11, 1.4040002e-11],
          dtype=float32)

• Indexes: (1)
  • hex_labels
    PandasIndex

    PandasIndex(Index([ (-5, 2, 3),  (-5, 3, 2),  (-4, 1, 3),  (-4, 2, 2),  (-4, 3, 1),
            (-3, 0, 3),  (-3, 1, 2),  (-3, 2, 1),  (-2, 0, 2),  (-2, 1, 1),
            (-2, 2, 0), (-1, -1, 2),  (-1, 0, 1),  (-1, 1, 0), (-1, 2, -1),
            (0, -1, 1),   (0, 0, 0),  (0, 1, -1),  (1, -1, 0),  (1, 0, -1),
            (1, 1, -2),  (2, -2, 0), (2, -1, -1), (3, -2, -1), (3, -1, -2),
           (4, -3, -1), (4, -2, -2), (5, -3, -2)],
          dtype='object', name='hex_labels'))

• Attributes: (0)

Visualize aggregated data

fig, ax = plt.subplots(
    1, 3,
    subplot_kw={"projection": cartopy.crs.Stereographic(
        central_latitude=lat_origin,
        central_longitude=lon_origin,
    )},
    figsize=(13, 3),
)

data_agg_hex.plot.scatter(
    x="lon", y="lat", hue="particle_obs_count", s=190,
    ax=ax[0], transform=cartopy.crs.PlateCarree(),
)

data_agg_hex.plot.scatter(
    x="lon", y="lat", hue="U", s=190,
    ax=ax[1], transform=cartopy.crs.PlateCarree(),
)
data_agg_hex.plot.quiver(
    u="U", v="V", x="lon", y="lat",
    ax=ax[1], transform=cartopy.crs.PlateCarree(),
)

data_agg_hex.plot.scatter(
    x="lon", y="lat", hue="V", s=190,
    ax=ax[2], transform=cartopy.crs.PlateCarree(),
)
data_agg_hex.plot.quiver(
    u="U", v="V", x="lon", y="lat",
    ax=ax[2], transform=cartopy.crs.PlateCarree(),
)


for _ax in ax.flat:
    _ax.coastlines()
    _ax.gridlines(draw_labels=True)
    _ax.set_extent([10, 55, -45, -22])

fig.tight_layout()

Connectivity

connection_count = (
    xr.merge(
        (
            ds_traj.isel(obs=0).hex_labels.rename("hex_start"),
            ds_traj.isel(obs=-1).hex_labels.rename("hex_end"),
        )
    )
    .to_dataframe().reset_index()
    .groupby(["hex_start", "hex_end"])["traj"].count()
    .rename("connection_count")
    .to_frame()
)

display(connection_count)

		connection_count
hex_start	hex_end
(-2, 1, 1)	(-5, 3, 2)	2
	(-4, 1, 3)	4
	(-2, 0, 2)	2
	(-2, 1, 1)	14
	(-1, -1, 2)	3
	(-1, 0, 1)	1
	(-1, 1, 0)	5
	(0, 0, 0)	1
	(1, -1, 0)	1
	(1, 0, -1)	2
	(1, 1, -2)	1
	(2, -2, 0)	1
	(2, -1, -1)	3
	(3, -1, -2)	1
(-2, 2, 0)	(-4, 1, 3)	4
	(-4, 2, 2)	1
	(-4, 3, 1)	9
	(-3, 1, 2)	2
	(-3, 2, 1)	5
	(-2, 0, 2)	2
	(-2, 1, 1)	4
	(-1, -1, 2)	1
	(-1, 0, 1)	2
	(0, -1, 1)	1
(-1, 0, 1)	(-1, 1, 0)	1
	(0, -1, 1)	2
	(0, 0, 0)	1
	(1, -1, 0)	4
	(1, 0, -1)	3
	(2, -2, 0)	2
	(2, -1, -1)	2
	(3, -1, -2)	3
	(4, -2, -2)	1
(-1, 1, 0)	(-4, 1, 3)	1
	(-3, 1, 2)	2
	(-3, 2, 1)	2
	(-2, 0, 2)	2
	(-2, 1, 1)	13
	(-2, 2, 0)	1
	(-1, 0, 1)	1
	(-1, 1, 0)	1
	(0, 0, 0)	8
	(0, 1, -1)	2
	(1, -1, 0)	1
	(1, 0, -1)	1
	(2, -2, 0)	2
	(2, -1, -1)	5
	(3, -2, -1)	1
	(4, -3, -1)	3
	(5, -3, -2)	2
(-1, 2, -1)	(-5, 2, 3)	1
	(-4, 1, 3)	1
	(-4, 3, 1)	2
	(-3, 2, 1)	1

import seaborn as sns

sns.heatmap(connection_count.unstack(level=-1).fillna(0))

<AxesSubplot: xlabel='None-hex_end', ylabel='hex_start'>