Seed grids and flow maps¶
This notebook is about the structures and the API, not about advection.
A SeedGrid holds a grid of release positions. SeedGrid.to_parcels_pset()
emits them as flat lon, lat lists. SeedGrid.pset_to_flowmap(...) ingests the
advected positions and returns a FlowMap, which computes the deformation
gradient, the Cauchy–Green tensor and the FTLE. No Parcels is imported below;
the advection is replaced by particles that do not move.
There are two stencils, and each is a pair of classes. NeighborSeedGrid releases
one particle per grid point and differences against the neighbouring grid
points. AuxiliarySeedGrid releases four arms around each grid point and
differences across the arms.
import numpy as np
from lcs_parcels import AuxiliarySeedGrid, NeighborSeedGrid
Grid and window¶
The axes set the grid. The times are not attached to anything until ingest.
t1 > t0 makes the window \(T = t_1 - t_0\) positive, so the FTLE below is the
forward one.
lon_axis = np.linspace(-25.0, -20.0, 6)
lat_axis = np.linspace(15.0, 20.0, 5)
t0 = np.datetime64("2020-01-01")
t1 = np.datetime64("2020-01-11")
The neighbour stencil¶
One particle per grid point: the release positions lon_0/lat_0 are the
diagnostic grid points lon_grid/lat_grid, stored twice rather than left to
be reconstructed. The outermost ring of grid points has no neighbour to
difference against on one side, so the FTLE is undefined there.
seed = NeighborSeedGrid.from_axes(lon=lon_axis, lat=lat_axis)
print(seed)
seed.ds
<NeighborSeedGrid 6x5 grid, lon -25.00..-20.00, lat 15.00..20.00>
<xarray.Dataset> Size: 1kB
Dimensions: (i: 6, j: 5)
Coordinates:
* i (i) int64 48B 0 1 2 3 4 5
* j (j) int64 40B 0 1 2 3 4
lon_grid (i, j) float64 240B -25.0 -25.0 -25.0 -25.0 ... -20.0 -20.0 -20.0
lat_grid (i, j) float64 240B 15.0 16.25 17.5 18.75 ... 17.5 18.75 20.0
lon_0 (i, j) float64 240B -25.0 -25.0 -25.0 -25.0 ... -20.0 -20.0 -20.0
lat_0 (i, j) float64 240B 15.0 16.25 17.5 18.75 ... 17.5 18.75 20.0
Data variables:
*empty*The auxiliary stencil¶
Four arms (east, north, west, south) at aux_separation_m around each
grid point, so lon_0/lat_0 gain a displacement dim and there are four
particles per grid point. The finite-difference step is the arm separation
rather than the grid spacing, and the FTLE is defined at every grid point,
including the boundary.
aux_seed = AuxiliarySeedGrid.from_axes(
lon=lon_axis, lat=lat_axis, aux_separation_m=1_000.0
)
print(aux_seed)
aux_seed.ds
<AuxiliarySeedGrid 6x5 grid, lon -25.00..-20.00, lat 15.00..20.00>
<xarray.Dataset> Size: 3kB
Dimensions: (displacement: 4, i: 6, j: 5)
Coordinates:
* displacement (displacement) <U5 80B 'east' 'north' 'west' 'south'
lon_0 (i, j, displacement) float64 960B -24.99 -25.0 ... -20.0
lat_0 (i, j, displacement) float64 960B 15.0 15.01 ... 20.0 19.99
* i (i) int64 48B 0 1 2 3 4 5
* j (j) int64 40B 0 1 2 3 4
lon_grid (i, j) float64 240B -25.0 -25.0 -25.0 ... -20.0 -20.0 -20.0
lat_grid (i, j) float64 240B 15.0 16.25 17.5 18.75 ... 17.5 18.75 20.0
Data variables:
*empty*Emit the particle sets¶
to_parcels_pset() flattens the release positions to two plain lists, one
entry per particle: 30 grid points for the neighbour stencil, 4 arms each for
the auxiliary one. The two lists are everything Parcels needs, and ingest
reattaches the advected positions to the grid by their order.
lon_0, lat_0 = seed.to_parcels_pset()
aux_lon_0, aux_lat_0 = aux_seed.to_parcels_pset()
print("neighbour particles:", len(lon_0))
print("auxiliary particles:", len(aux_lon_0))
neighbour particles: 30
auxiliary particles: 120
Advect, without Parcels¶
The particles sit still: \(u = v = 0\) for ten days. The flow map is the identity, \(\nabla F = I\), so the FTLE is zero wherever it is defined. A non-zero window keeps \(1/T\) finite.
lon_1, lat_1 = lon_0, lat_0
aux_lon_1, aux_lat_1 = aux_lon_0, aux_lat_0
Ingest the advected positions¶
pset_to_flowmap puts the flat positions back on the grid as lon/lat
beside the release positions lon_0/lat_0, and records t0 and the signed
window T it derives from t0 and t1. Both stencils ingest through the same
call; the datasets differ only where their release positions did.
flowmap = seed.pset_to_flowmap(lon=lon_1, lat=lat_1, t0=t0, t1=t1)
print(flowmap)
flowmap.ds
<NeighborFlowMap 6x5 grid, lon -25.00..-20.00, lat 15.00..20.00,
t0 2020-01-01T00:00:00, T +10.0 days>
<xarray.Dataset> Size: 2kB
Dimensions: (i: 6, j: 5)
Coordinates:
* i (i) int64 48B 0 1 2 3 4 5
* j (j) int64 40B 0 1 2 3 4
lon_grid (i, j) float64 240B -25.0 -25.0 -25.0 -25.0 ... -20.0 -20.0 -20.0
lat_grid (i, j) float64 240B 15.0 16.25 17.5 18.75 ... 17.5 18.75 20.0
lon_0 (i, j) float64 240B -25.0 -25.0 -25.0 -25.0 ... -20.0 -20.0 -20.0
lat_0 (i, j) float64 240B 15.0 16.25 17.5 18.75 ... 17.5 18.75 20.0
t0 datetime64[s] 8B 2020-01-01
T timedelta64[s] 8B 10 days
Data variables:
lon (i, j) float64 240B -25.0 -25.0 -25.0 -25.0 ... -20.0 -20.0 -20.0
lat (i, j) float64 240B 15.0 16.25 17.5 18.75 ... 17.5 18.75 20.0aux_flowmap = aux_seed.pset_to_flowmap(lon=aux_lon_1, lat=aux_lat_1, t0=t0, t1=t1)
print(aux_flowmap)
aux_flowmap.ds
<AuxiliaryFlowMap 6x5 grid, lon -25.00..-20.00, lat 15.00..20.00,
t0 2020-01-01T00:00:00, T +10.0 days>
<xarray.Dataset> Size: 5kB
Dimensions: (displacement: 4, i: 6, j: 5)
Coordinates:
* displacement (displacement) <U5 80B 'east' 'north' 'west' 'south'
lon_0 (i, j, displacement) float64 960B -24.99 -25.0 ... -20.0
lat_0 (i, j, displacement) float64 960B 15.0 15.01 ... 20.0 19.99
* i (i) int64 48B 0 1 2 3 4 5
* j (j) int64 40B 0 1 2 3 4
lon_grid (i, j) float64 240B -25.0 -25.0 -25.0 ... -20.0 -20.0 -20.0
lat_grid (i, j) float64 240B 15.0 16.25 17.5 18.75 ... 17.5 18.75 20.0
t0 datetime64[s] 8B 2020-01-01
T timedelta64[s] 8B 10 days
Data variables:
lon (i, j, displacement) float64 960B -24.99 -25.0 ... -20.0
lat (i, j, displacement) float64 960B 15.0 15.01 ... 20.0 19.99grid_image is the one place the two datasets are reconciled: the advected
position of the grid point, on (i, j). The neighbour flow map passes its
single advected position through; the auxiliary one takes the centroid of its
four advected arms.
aux_flowmap.grid_image
<xarray.Dataset> Size: 1kB
Dimensions: (i: 6, j: 5)
Coordinates:
* i (i) int64 48B 0 1 2 3 4 5
* j (j) int64 40B 0 1 2 3 4
lon_grid (i, j) float64 240B -25.0 -25.0 -25.0 -25.0 ... -20.0 -20.0 -20.0
lat_grid (i, j) float64 240B 15.0 16.25 17.5 18.75 ... 17.5 18.75 20.0
t0 datetime64[s] 8B 2020-01-01
T timedelta64[s] 8B 10 days
Data variables:
lon (i, j) float64 240B -25.0 -25.0 -25.0 -25.0 ... -20.0 -20.0 -20.0
lat (i, j) float64 240B 15.0 16.25 17.5 18.75 ... 17.5 18.75 20.0FTLE from both stencils¶
ftle() returns 1/s on the diagnostic grid. Both stencils give zero; they
differ in where the FTLE is defined.
ftle = flowmap.ftle()
ftle
<xarray.DataArray 'ftle' (i: 6, j: 5)> Size: 240B
array([[nan, nan, nan, nan, nan],
[nan, 0., 0., 0., nan],
[nan, 0., 0., 0., nan],
[nan, 0., 0., 0., nan],
[nan, 0., 0., 0., nan],
[nan, nan, nan, nan, nan]])
Coordinates:
* i (i) int64 48B 0 1 2 3 4 5
* j (j) int64 40B 0 1 2 3 4
lon_grid (i, j) float64 240B -25.0 -25.0 -25.0 -25.0 ... -20.0 -20.0 -20.0
lat_grid (i, j) float64 240B 15.0 16.25 17.5 18.75 ... 17.5 18.75 20.0
t0 datetime64[s] 8B 2020-01-01
T timedelta64[s] 8B 10 days
Attributes:
long_name: finite-time Lyapunov exponent
units: 1/saux_ftle = aux_flowmap.ftle()
aux_ftle
<xarray.DataArray 'ftle' (i: 6, j: 5)> Size: 240B
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.]])
Coordinates:
* i (i) int64 48B 0 1 2 3 4 5
* j (j) int64 40B 0 1 2 3 4
lon_grid (i, j) float64 240B -25.0 -25.0 -25.0 -25.0 ... -20.0 -20.0 -20.0
lat_grid (i, j) float64 240B 15.0 16.25 17.5 18.75 ... 17.5 18.75 20.0
t0 datetime64[s] 8B 2020-01-01
T timedelta64[s] 8B 10 days
Attributes:
long_name: finite-time Lyapunov exponent
units: 1/sprint("neighbour: max |ftle| =", float(abs(ftle).max()))
print("auxiliary: max |ftle| =", float(abs(aux_ftle).max()))
print("neighbour NaN points:", int(ftle.isnull().sum()))
print("auxiliary NaN points:", int(aux_ftle.isnull().sum()))
neighbour: max |ftle| = 0.0
auxiliary: max |ftle| = 0.0
neighbour NaN points: 18
auxiliary NaN points: 0