Welcome to py-eddy-tracker’s documentation!¶
How do I get set up ?¶
Source are available on github https://github.com/AntSimi/py-eddy-tracker
Use python3. To avoid problems with installation, use of the virtualenv Python virtual environment is recommended or conda.
Then use pip to install all dependencies (numpy, scipy, matplotlib, netCDF4, …), e.g.:
pip install numpy scipy netCDF4 matplotlib opencv-python pyyaml pint polygon3
Then run the following to install the eddy tracker:
python setup.py install
Several executables are available in your PATH:
GridFiltering # Allow to apply a high frequency filter on a NetCDF grid
EddyId # Provide identification of eddies for one grid
EddySubSetter # Allow to apply sub setting on eddies dataset
EddyTracking # Allow to track Identification dataset
Py eddy tracker toolbox¶
All figures in this gallery, used an experimental dataset, compute with this dataset : cmems_product.
Eddy detection¶
Note
Click here to download the full example code or to run this example in your browser via Binder
Display contour & circle¶
from matplotlib import pyplot as plt
from py_eddy_tracker import data
from py_eddy_tracker.observations.observation import EddiesObservations
Load detection files
a = EddiesObservations.load_file(data.get_path("Anticyclonic_20190223.nc"))
Plot the speed and effective (dashed) contours
fig = plt.figure(figsize=(15, 8))
ax = fig.add_axes((0.05, 0.05, 0.9, 0.9))
ax.set_aspect("equal")
ax.set_xlim(10, 70)
ax.set_ylim(-50, -25)
a.display(ax, label="Anticyclonic contour", color="r", lw=1)
# Replace contours by circles using center and radius (effective is dashed)
a.circle_contour()
a.display(ax, label="Anticyclonic circle", color="g", lw=1)
ax.legend(loc="upper right")
Total running time of the script: ( 0 minutes 1.568 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Display identification¶
from matplotlib import pyplot as plt
from py_eddy_tracker import data
from py_eddy_tracker.observations.observation import EddiesObservations
Load detection files
a = EddiesObservations.load_file(data.get_path("Anticyclonic_20190223.nc"))
c = EddiesObservations.load_file(data.get_path("Cyclonic_20190223.nc"))
Fill effective contour with amplitude
fig = plt.figure(figsize=(15, 8))
ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])
ax.set_aspect("equal")
ax.set_xlim(0, 140)
ax.set_ylim(-80, 0)
kwargs = dict(extern_only=True, color="k", lw=1)
a.display(ax, **kwargs), c.display(ax, **kwargs)
a.filled(ax, "amplitude", cmap="magma_r", vmin=0, vmax=0.5)
m = c.filled(ax, "amplitude", cmap="magma_r", vmin=0, vmax=0.5)
colorbar = plt.colorbar(m, cax=ax.figure.add_axes([0.95, 0.03, 0.02, 0.94]))
colorbar.set_label("Amplitude (m)")
Draw speed contours
fig = plt.figure(figsize=(15, 8))
ax = fig.add_axes([0.03, 0.03, 0.94, 0.94])
ax.set_aspect("equal")
ax.set_xlim(0, 360)
ax.set_ylim(-80, 80)
a.display(ax, label="Anticyclonic ({nb_obs} eddies)", color="r", lw=1)
c.display(ax, label="Cyclonic ({nb_obs} eddies)", color="b", lw=1)
ax.legend(loc="upper right")
Get general informations
print(a)
Out:
| 3137 observations from 25255 to 25255 (1 days, ~3137 obs/day)
| Speed area : 32.98 Mkm²/day
| Effective area : 45.65 Mkm²/day
----Distribution in Amplitude:
| Amplitude bounds (cm) 0.00 1.00 2.00 3.00 4.00 5.00 10.00 500.00
| Percent of eddies : 19.35 22.73 15.40 10.30 6.18 15.91 10.14
----Distribution in Radius:
| Speed radius (km) 0.00 15.00 30.00 45.00 60.00 75.00 100.00 200.00 2000.00
| Percent of eddies : 0.00 9.47 34.56 24.55 13.29 11.67 6.34 0.13
| Effective radius (km) 0.00 15.00 30.00 45.00 60.00 75.00 100.00 200.00 2000.00
| Percent of eddies : 0.00 7.52 26.62 20.88 15.40 15.94 13.32 0.32
----Distribution in Latitude
Latitude bounds -90.00 -60.00 -15.00 15.00 60.00 90.00
Percent of eddies : 7.62 46.86 12.81 30.06 2.65
Percent of speed area : 4.69 41.94 26.90 25.30 1.17
Percent of effective area : 4.74 43.40 25.53 25.11 1.21
Mean speed radius (km) : 43.94 52.75 81.69 51.01 37.91
Mean effective radius (km): 52.14 62.43 94.14 59.44 44.81
Mean amplitude (cm) : 3.53 5.30 2.19 4.32 3.12
print(c)
Out:
| 3360 observations from 25255 to 25255 (1 days, ~3360 obs/day)
| Speed area : 32.89 Mkm²/day
| Effective area : 46.42 Mkm²/day
----Distribution in Amplitude:
| Amplitude bounds (cm) 0.00 1.00 2.00 3.00 4.00 5.00 10.00 500.00
| Percent of eddies : 18.81 24.02 14.11 10.89 5.98 16.19 10.00
----Distribution in Radius:
| Speed radius (km) 0.00 15.00 30.00 45.00 60.00 75.00 100.00 200.00 2000.00
| Percent of eddies : 0.03 10.15 35.03 25.15 14.40 10.09 5.12 0.03
| Effective radius (km) 0.00 15.00 30.00 45.00 60.00 75.00 100.00 200.00 2000.00
| Percent of eddies : 0.03 7.98 26.88 21.61 15.92 15.09 12.14 0.36
----Distribution in Latitude
Latitude bounds -90.00 -60.00 -15.00 15.00 60.00 90.00
Percent of eddies : 7.92 46.96 13.12 29.61 2.38
Percent of speed area : 4.80 41.08 27.30 25.87 0.93
Percent of effective area : 4.83 42.35 25.36 26.55 0.92
Mean speed radius (km) : 42.23 50.71 78.76 50.80 34.64
Mean effective radius (km): 49.25 60.50 89.91 59.96 40.20
Mean amplitude (cm) : 3.19 5.71 2.19 4.24 2.42
Total running time of the script: ( 0 minutes 2.586 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Radius vs area¶
from matplotlib import pyplot as plt
from numpy import array, pi
from py_eddy_tracker import data
from py_eddy_tracker.generic import coordinates_to_local
from py_eddy_tracker.observations.observation import EddiesObservations
from py_eddy_tracker.poly import poly_area
Load detection files
a = EddiesObservations.load_file(data.get_path("Anticyclonic_20190223.nc"))
areas = list()
# For each contour area will be compute in local reference
for i in a:
x, y = coordinates_to_local(
i["contour_lon_s"], i["contour_lat_s"], i["lon"], i["lat"]
)
areas.append(poly_area(x, y))
areas = array(areas)
Radius provided by eddy detection is computed with fit_circle() method.
This radius will be compared with an equivalent radius deduced from polygon area.
ax = plt.subplot(111)
ax.set_aspect("equal")
ax.grid()
ax.set_xlabel("Speed radius computed with fit_circle")
ax.set_ylabel("Radius deduced from area\nof contour_lon_s/contour_lat_s")
ax.set_title("Area vs radius")
ax.plot(a["radius_s"] / 1000.0, (areas / pi) ** 0.5 / 1000.0, ".")
ax.plot((0, 250), (0, 250), "r")
Fit circle give a radius bigger than polygon area
When error is tiny, radius are very close.
ax = plt.subplot(111)
ax.grid()
ax.set_xlabel("Radius ratio")
ax.set_ylabel("Shape error")
ax.set_title("err = f(radius_ratio)")
ax.plot(a["radius_s"] / (areas / pi) ** 0.5, a["shape_error_s"], ".")
Total running time of the script: ( 0 minutes 3.475 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Shape error gallery¶
Gallery of contours with shape error
from matplotlib import pyplot as plt
from numpy import arange, cos, linspace, radians, sin
from py_eddy_tracker import data
from py_eddy_tracker.dataset.grid import RegularGridDataset
from py_eddy_tracker.eddy_feature import Contours
from py_eddy_tracker.generic import local_to_coordinates
Method to built circle from center coordinates
def build_circle(x0, y0, r):
angle = radians(linspace(0, 360, 50))
x_norm, y_norm = cos(angle), sin(angle)
return local_to_coordinates(x_norm * r, y_norm * r, x0, y0)
We iterate over closed contours and sort with regards of shape error
g = RegularGridDataset(
data.get_path("dt_med_allsat_phy_l4_20160515_20190101.nc"), "longitude", "latitude"
)
c = Contours(g.x_c, g.y_c, g.grid("adt") * 100, arange(-50, 50, 0.2))
contours = dict()
for coll in c.iter():
for current_contour in coll.get_paths():
_, _, _, aerr = current_contour.fit_circle()
i = int(aerr // 4) + 1
if i not in contours:
contours[i] = list()
contours[i].append(current_contour)
Out:
We assume pixel position of grid is center for /home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/pyEddyTracker-3.3.0-py3.7.egg/py_eddy_tracker/data/dt_med_allsat_phy_l4_20160515_20190101.nc
Shape error gallery¶
For each contour display, we display circle fitted, we work at different latitude circle could have distorsion
fig = plt.figure(figsize=(12, 12))
for i in range(1, 26):
e_min, e_max = (i - 1) * 4, i * 4
ax = plt.subplot(5, 5, i, title=f" {e_min} < err < {e_max}")
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
ax.set_aspect("equal")
ax.grid()
if i in contours:
for contour in contours[i]:
x, y = contour.lon, contour.lat
x0, y0, radius, _ = contour.fit_circle()
if x.shape[0] > 30 and 30000 < radius < 70000:
# Plot only first contour found
m = ax.plot(x, y, "r")[0]
ax.plot(*build_circle(x0, y0, radius), "g--")
ax.plot(x0, y0, "k.")
break
plt.tight_layout()
Total running time of the script: ( 0 minutes 9.358 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Get mean of grid in each eddies¶
from matplotlib import pyplot as plt
from py_eddy_tracker import data
from py_eddy_tracker.dataset.grid import RegularGridDataset
from py_eddy_tracker.observations.observation import EddiesObservations
def start_axes(title):
fig = plt.figure(figsize=(13, 5))
ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])
ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
ax.set_aspect("equal")
ax.set_title(title)
return ax
def update_axes(ax, mappable=None):
ax.grid()
ax.legend()
if mappable:
plt.colorbar(mappable, cax=ax.figure.add_axes([0.95, 0.05, 0.01, 0.9]))
Load detection files and data to interp
a = EddiesObservations.load_file(data.get_path("Anticyclonic_20160515.nc"))
c = EddiesObservations.load_file(data.get_path("Cyclonic_20160515.nc"))
aviso_map = RegularGridDataset(
data.get_path("dt_med_allsat_phy_l4_20160515_20190101.nc"), "longitude", "latitude"
)
aviso_map.add_uv("adt")
Out:
We assume pixel position of grid is center for /home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/pyEddyTracker-3.3.0-py3.7.egg/py_eddy_tracker/data/dt_med_allsat_phy_l4_20160515_20190101.nc
Compute and store eke in cm²/s²
aviso_map.add_grid(
"eke", (aviso_map.grid("u") ** 2 + aviso_map.grid("v") ** 2) * 0.5 * (100 ** 2)
)
eke_kwargs = dict(vmin=1, vmax=1000, cmap="magma_r")
ax = start_axes("EKE (cm²/s²)")
m = aviso_map.display(ax, "eke", **eke_kwargs)
a.display(ax, color="r", linewidth=0.5, label="Anticyclonic", ref=-10)
c.display(ax, color="b", linewidth=0.5, label="Cyclonic", ref=-10)
update_axes(ax, m)
Get mean of eke in each effective contour
ax = start_axes("EKE mean (cm²/s²)")
a.display(ax, color="r", linewidth=0.5, label="Anticyclonic ({nb_obs} eddies)", ref=-10)
c.display(ax, color="b", linewidth=0.5, label="Cyclonic ({nb_obs} eddies)", ref=-10)
eke = a.interp_grid(aviso_map, "eke", method="mean", intern=False)
a.filled(ax, eke, ref=-10, **eke_kwargs)
eke = c.interp_grid(aviso_map, "eke", method="mean", intern=False)
m = c.filled(ax, eke, ref=-10, **eke_kwargs)
update_axes(ax, m)
Total running time of the script: ( 0 minutes 5.625 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Eddy detection : Med¶
Script will detect eddies on adt field, and compute u,v with method add_uv(which could use, only if equator is avoid)
Figures will show different step to detect eddies.
from datetime import datetime
from matplotlib import pyplot as plt
from numpy import arange
from py_eddy_tracker import data
from py_eddy_tracker.dataset.grid import RegularGridDataset
def start_axes(title):
fig = plt.figure(figsize=(13, 5))
ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])
ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
ax.set_aspect("equal")
ax.set_title(title, weight="bold")
return ax
def update_axes(ax, mappable=None):
ax.grid()
if mappable:
plt.colorbar(mappable, cax=ax.figure.add_axes([0.94, 0.05, 0.01, 0.9]))
Load Input grid, ADT is used to detect eddies
g = RegularGridDataset(
data.get_path("dt_med_allsat_phy_l4_20160515_20190101.nc"), "longitude", "latitude"
)
ax = start_axes("ADT (m)")
m = g.display(ax, "adt", vmin=-0.15, vmax=0.15, cmap="RdBu_r")
update_axes(ax, m)
Out:
We assume pixel position of grid is center for /home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/pyEddyTracker-3.3.0-py3.7.egg/py_eddy_tracker/data/dt_med_allsat_phy_l4_20160515_20190101.nc
Get geostrophic speed u,v¶
U/V are deduced from ADT, this algortihm is not ok near the equator (~+- 2°)
g.add_uv("adt")
ax = start_axes("U/V deduce from ADT (m)")
ax.set_xlim(2.5, 9), ax.set_ylim(37.5, 40)
m = g.display(ax, "adt", vmin=-0.15, vmax=0.15, cmap="RdBu_r")
u, v = g.grid("u").T, g.grid("v").T
ax.quiver(g.x_c, g.y_c, u, v, scale=10)
update_axes(ax, m)
Pre-processings¶
Apply a high-pass filter to remove the large scale and highlight the mesoscale
g.bessel_high_filter("adt", 500)
ax = start_axes("ADT (m) filtered (500km)")
m = g.display(ax, "adt", vmin=-0.15, vmax=0.15, cmap="RdBu_r")
update_axes(ax, m)
Identification¶
Run the identification step with slices of 2 mm
date = datetime(2016, 5, 15)
a, c = g.eddy_identification("adt", "u", "v", date, 0.002, shape_error=55)
Display of all closed contours found in the grid (only 1 contour every 4)
ax = start_axes("ADT closed contours (only 1 / 4 levels)")
g.contours.display(ax, step=4)
update_axes(ax)
Contours included in eddies
ax = start_axes("ADT contours used as eddies")
g.contours.display(ax, only_used=True)
update_axes(ax)
Post analysis¶
Contours can be rejected for several reasons (shape error to high, several extremum in contour, …)
ax = start_axes("ADT rejected contours")
g.contours.display(ax, only_unused=True)
update_axes(ax)
Creteria for rejecting a contour 0. - Accepted (green) 1. - Rejection for shape error (red) 2. - Masked value within contour (blue) 3. - Under or over the pixel limit bounds (black) 4. - Amplitude criterion (yellow)
ax = start_axes("Contours' rejection criteria")
g.contours.display(ax, only_unused=True, lw=0.5, display_criterion=True)
update_axes(ax)
Display the shape error of each tested contour, the limit of shape error is set to 55 %
ax = start_axes("Contour shape error")
m = g.contours.display(
ax, lw=0.5, field="shape_error", bins=arange(20, 90.1, 5), cmap="PRGn_r"
)
update_axes(ax, m)
Some closed contours contains several eddies (aka, more than one extremum)
ax = start_axes("ADT rejected contours containing eddies")
g.contours.label_contour_unused_which_contain_eddies(a)
g.contours.label_contour_unused_which_contain_eddies(c)
g.contours.display(
ax,
only_contain_eddies=True,
color="k",
lw=1,
label="Could be a contour of interaction",
)
a.display(ax, color="r", linewidth=0.75, label="Anticyclonic", ref=-10)
c.display(ax, color="b", linewidth=0.75, label="Cyclonic", ref=-10)
ax.legend()
update_axes(ax)
Output¶
When displaying the detected eddies, dashed lines are for effective contour, solide lines for the contour of the maximum mean speed. See figure 1 of https://doi.org/10.1175/JTECH-D-14-00019.1
ax = start_axes("Detected Eddies")
a.display(
ax, color="r", linewidth=0.75, label="Anticyclonic ({nb_obs} eddies)", ref=-10
)
c.display(ax, color="b", linewidth=0.75, label="Cyclonic ({nb_obs} eddies)", ref=-10)
ax.legend()
update_axes(ax)
Display the speed radius of the detected eddies
ax = start_axes("Speed Radius (km)")
a.scatter(ax, "radius_s", vmin=10, vmax=50, s=80, ref=-10, cmap="magma_r", factor=0.001)
m = c.scatter(
ax, "radius_s", vmin=10, vmax=50, s=80, ref=-10, cmap="magma_r", factor=0.001
)
update_axes(ax, m)
Filling the effective radius contours with the effective radius values
ax = start_axes("Effective Radius (km)")
kwargs = dict(vmin=10, vmax=80, cmap="magma_r", factor=0.001, lut=14, ref=-10)
a.filled(ax, "effective_radius", **kwargs)
m = c.filled(
ax, "radius_e", vmin=10, vmax=80, cmap="magma_r", factor=0.001, lut=14, ref=-10
)
update_axes(ax, m)
Total running time of the script: ( 0 minutes 19.600 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Eddy detection : Gulf stream¶
Script will detect eddies on adt field, and compute u,v with method add_uv(which could use, only if equator is avoid)
Figures will show different step to detect eddies.
from datetime import datetime
from matplotlib import pyplot as plt
from numpy import arange
from py_eddy_tracker import data
from py_eddy_tracker.dataset.grid import RegularGridDataset
from py_eddy_tracker.eddy_feature import Contours
def start_axes(title):
fig = plt.figure(figsize=(13, 8))
ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])
ax.set_xlim(279, 304), ax.set_ylim(29, 44)
ax.set_aspect("equal")
ax.set_title(title, weight="bold")
return ax
def update_axes(ax, mappable=None):
ax.grid()
if mappable:
plt.colorbar(mappable, cax=ax.figure.add_axes([0.94, 0.05, 0.01, 0.9]))
Load Input grid, ADT is used to detect eddies
margin = 30
g = RegularGridDataset(
data.get_path("nrt_global_allsat_phy_l4_20190223_20190226.nc"),
"longitude",
"latitude",
# Manual area subset
indexs=dict(
longitude=slice(1116 - margin, 1216 + margin),
latitude=slice(476 - margin, 536 + margin),
),
)
ax = start_axes("ADT (m)")
m = g.display(ax, "adt", vmin=-1, vmax=1, cmap="RdBu_r")
# Draw line on the gulf stream front
great_current = Contours(g.x_c, g.y_c, g.grid("adt"), levels=(0.35,), keep_unclose=True)
great_current.display(ax, color="k")
update_axes(ax, m)
Out:
We assume pixel position of grid is center for /home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/pyEddyTracker-3.3.0-py3.7.egg/py_eddy_tracker/data/nrt_global_allsat_phy_l4_20190223_20190226.nc
Get geostrophic speed u,v¶
U/V are deduced from ADT, this algortihm is not ok near the equator (~+- 2°)
g.add_uv("adt")
Pre-processings¶
Apply a high-pass filter to remove the large scale and highlight the mesoscale
g.bessel_high_filter("adt", 700)
ax = start_axes("ADT (m) filtered (700km)")
m = g.display(ax, "adt", vmin=-0.4, vmax=0.4, cmap="RdBu_r")
great_current.display(ax, color="k")
update_axes(ax, m)
Identification¶
Run the identification step with slices of 2 mm
date = datetime(2016, 5, 15)
a, c = g.eddy_identification("adt", "u", "v", date, 0.002, shape_error=55)
Display of all closed contours found in the grid (only 1 contour every 5)
ax = start_axes("ADT closed contours (only 1 / 5 levels)")
g.contours.display(ax, step=5, lw=1)
great_current.display(ax, color="k")
update_axes(ax)
Contours included in eddies
ax = start_axes("ADT contours used as eddies")
g.contours.display(ax, only_used=True, lw=0.25)
great_current.display(ax, color="k")
update_axes(ax)
Post analysis¶
Contours can be rejected for several reasons (shape error to high, several extremum in contour, …)
ax = start_axes("ADT rejected contours")
g.contours.display(ax, only_unused=True, lw=0.25)
great_current.display(ax, color="k")
update_axes(ax)
- Criteria for rejecting a contour :
Accepted (green)
Rejection for shape error (red)
Masked value within contour (blue)
Under or over the pixel limit bounds (black)
Amplitude criterion (yellow)
ax = start_axes("Contours' rejection criteria")
g.contours.display(ax, only_unused=True, lw=0.5, display_criterion=True)
update_axes(ax)
Display the shape error of each tested contour, the limit of shape error is set to 55 %
ax = start_axes("Contour shape error")
m = g.contours.display(
ax, lw=0.5, field="shape_error", bins=arange(20, 90.1, 5), cmap="PRGn_r"
)
update_axes(ax, m)
Some closed contours contains several eddies (aka, more than one extremum)
ax = start_axes("ADT rejected contours containing eddies")
g.contours.label_contour_unused_which_contain_eddies(a)
g.contours.label_contour_unused_which_contain_eddies(c)
g.contours.display(
ax,
only_contain_eddies=True,
color="k",
lw=1,
label="Could be a contour of interaction",
)
a.display(ax, color="r", linewidth=0.75, label="Anticyclonic", ref=-10)
c.display(ax, color="b", linewidth=0.75, label="Cyclonic", ref=-10)
ax.legend()
update_axes(ax)
Output¶
When displaying the detected eddies, dashed lines are for effective contour, solide lines for the contour of the maximum mean speed. See figure 1 of https://doi.org/10.1175/JTECH-D-14-00019.1
ax = start_axes("Eddies detected")
a.display(
ax, color="r", linewidth=0.75, label="Anticyclonic ({nb_obs} eddies)", ref=-10
)
c.display(ax, color="b", linewidth=0.75, label="Cyclonic ({nb_obs} eddies)", ref=-10)
ax.legend()
great_current.display(ax, color="k")
update_axes(ax)
Display the effective radius of the detected eddies
ax = start_axes("Effective radius (km)")
a.filled(ax, "radius_e", vmin=10, vmax=150, cmap="magma_r", factor=0.001, lut=14)
m = c.filled(ax, "radius_e", vmin=10, vmax=150, cmap="magma_r", factor=0.001, lut=14)
great_current.display(ax, color="k")
update_axes(ax, m)
Total running time of the script: ( 0 minutes 10.348 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Eddy detection and filter¶
from datetime import datetime
from matplotlib import pyplot as plt
from numpy import arange
from py_eddy_tracker import data
from py_eddy_tracker.dataset.grid import RegularGridDataset
def start_axes(title):
fig = plt.figure(figsize=(13, 5))
ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])
ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
ax.set_aspect("equal")
ax.set_title(title, weight="bold")
return ax
def update_axes(ax, mappable=None):
ax.grid()
if mappable:
plt.colorbar(mappable, cax=ax.figure.add_axes([0.94, 0.05, 0.01, 0.9]))
Load Input grid, ADT is used to detect eddies. Add a new filed to store the high-pass filtered ADT
g = RegularGridDataset(
data.get_path("dt_med_allsat_phy_l4_20160515_20190101.nc"), "longitude", "latitude"
)
g.add_uv("adt")
g.copy("adt", "adt_high")
wavelength = 800
g.bessel_high_filter("adt_high", wavelength)
date = datetime(2016, 5, 15)
Out:
We assume pixel position of grid is center for /home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/pyEddyTracker-3.3.0-py3.7.egg/py_eddy_tracker/data/dt_med_allsat_phy_l4_20160515_20190101.nc
Run the detection for the total grid and the filtered grid
a_filtered, c_filtered = g.eddy_identification("adt_high", "u", "v", date, 0.002)
merge_f = a_filtered.merge(c_filtered)
a_tot, c_tot = g.eddy_identification("adt", "u", "v", date, 0.002)
merge_t = a_tot.merge(c_tot)
Display the two detections
ax = start_axes("Eddies detected over ADT")
m = g.display(ax, "adt", vmin=-0.15, vmax=0.15)
merge_f.display(
ax,
lw=0.75,
label="Eddies in the filtered grid ({nb_obs} eddies)",
ref=-10,
color="k",
)
merge_t.display(
ax, lw=0.75, label="Eddies without filter ({nb_obs} eddies)", ref=-10, color="r"
)
ax.legend()
update_axes(ax, m)
Amplitude and Speed Radius distributions¶
fig = plt.figure(figsize=(12, 5))
ax_a = fig.add_subplot(121, xlabel="Amplitude (cm)")
ax_r = fig.add_subplot(122, xlabel="Speed Radius (km)")
ax_a.hist(
merge_f.amplitude * 100,
bins=arange(0.0005, 100, 1),
label="Eddies in the filtered grid",
histtype="step",
)
ax_a.hist(
merge_t.amplitude * 100,
bins=arange(0.0005, 100, 1),
label="Eddies without filter",
histtype="step",
)
ax_a.set_xlim(0, 10)
ax_r.hist(merge_f.radius_s / 1000.0, bins=arange(0, 300, 5), histtype="step")
ax_r.hist(merge_t.radius_s / 1000.0, bins=arange(0, 300, 5), histtype="step")
ax_r.set_xlim(0, 100)
ax_a.legend()
Match detection and compare¶
i_, j_, c = merge_f.match(merge_t, cmin=0.1)
Where are the lonely eddies?
kwargs_f = dict(lw=1.5, label="Lonely eddies in the filtered grid", ref=-10, color="k")
kwargs_t = dict(lw=1.5, label="Lonely eddies without filter", ref=-10, color="r")
ax = start_axes("Eddies with no match, over filtered ADT")
mappable = g.display(ax, "adt_high", vmin=-0.15, vmax=0.15)
merge_f.index(i_, reverse=True).display(ax, **kwargs_f)
merge_t.index(j_, reverse=True).display(ax, **kwargs_t)
ax.legend()
update_axes(ax, mappable)
ax = start_axes("Eddies with no match, over filtered ADT (zoom)")
ax.set_xlim(25, 36), ax.set_ylim(31, 35.25)
mappable = g.display(ax, "adt_high", vmin=-0.15, vmax=0.15)
u, v = g.grid("u").T, g.grid("v").T
ax.quiver(g.x_c, g.y_c, u, v, scale=10, pivot="mid", color="gray")
merge_f.index(i_, reverse=True).display(ax, **kwargs_f)
merge_t.index(j_, reverse=True).display(ax, **kwargs_t)
ax.legend()
update_axes(ax, mappable)
fig = plt.figure(figsize=(12, 12))
fig.suptitle(f"Scatter plot ({i_.shape[0]} matches)", weight="bold")
for i, (label, field, factor, stop) in enumerate(
(
("Speed radius (km)", "radius_s", 0.001, 80),
("Effective radius (km)", "radius_e", 0.001, 120),
("Amplitude (cm)", "amplitude", 100, 25),
("Maximum Speed (cm/s)", "speed_average", 100, 25),
)
):
ax = fig.add_subplot(
2, 2, i + 1, xlabel="Filtered grid", ylabel="Without filter", title=label
)
ax.plot(merge_f[field][i_] * factor, merge_t[field][j_] * factor, ".")
ax.set_aspect("equal"), ax.grid()
ax.plot((0, 1000), (0, 1000), "r")
ax.set_xlim(0, stop), ax.set_ylim(0, stop)
Total running time of the script: ( 0 minutes 8.243 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Eddy detection on SLA and ADT¶
from datetime import datetime
from matplotlib import pyplot as plt
from py_eddy_tracker import data
from py_eddy_tracker.dataset.grid import RegularGridDataset
def start_axes(title):
fig = plt.figure(figsize=(13, 5))
ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])
ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
ax.set_aspect("equal")
ax.set_title(title)
return ax
def update_axes(ax, mappable=None):
ax.grid()
if mappable:
plt.colorbar(mappable, cax=ax.figure.add_axes([0.95, 0.05, 0.01, 0.9]))
Load Input grid, ADT will be used to detect eddies
g = RegularGridDataset(
data.get_path("dt_med_allsat_phy_l4_20160515_20190101.nc"), "longitude", "latitude"
)
g.add_uv("adt", "ugos", "vgos")
g.add_uv("sla", "ugosa", "vgosa")
wavelength = 400
g.copy("adt", "adt_raw")
g.copy("sla", "sla_raw")
g.bessel_high_filter("adt", wavelength)
g.bessel_high_filter("sla", wavelength)
date = datetime(2016, 5, 15)
Out:
We assume pixel position of grid is center for /home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/pyEddyTracker-3.3.0-py3.7.egg/py_eddy_tracker/data/dt_med_allsat_phy_l4_20160515_20190101.nc
kwargs_a_adt = dict(
lw=0.5, label="Anticyclonic ADT ({nb_obs} eddies)", ref=-10, color="k"
)
kwargs_c_adt = dict(lw=0.5, label="Cyclonic ADT ({nb_obs} eddies)", ref=-10, color="r")
kwargs_a_sla = dict(
lw=0.5, label="Anticyclonic SLA ({nb_obs} eddies)", ref=-10, color="g"
)
kwargs_c_sla = dict(lw=0.5, label="Cyclonic SLA ({nb_obs} eddies)", ref=-10, color="b")
Run algorithm of detection
a_adt, c_adt = g.eddy_identification("adt", "ugos", "vgos", date, 0.002)
a_sla, c_sla = g.eddy_identification("sla", "ugosa", "vgosa", date, 0.002)
over filtered
ax = start_axes(f"ADT (m) filtered ({wavelength}km)")
m = g.display(ax, "adt", vmin=-0.15, vmax=0.15)
a_adt.display(ax, **kwargs_a_adt), c_adt.display(ax, **kwargs_c_adt)
ax.legend(), update_axes(ax, m)
ax = start_axes(f"SLA (m) filtered ({wavelength}km)")
m = g.display(ax, "sla", vmin=-0.15, vmax=0.15)
a_sla.display(ax, **kwargs_a_sla), c_sla.display(ax, **kwargs_c_sla)
ax.legend(), update_axes(ax, m)
over raw
ax = start_axes("ADT (m)")
m = g.display(ax, "adt_raw", vmin=-0.15, vmax=0.15)
a_adt.display(ax, **kwargs_a_adt), c_adt.display(ax, **kwargs_c_adt)
ax.legend(), update_axes(ax, m)
ax = start_axes("SLA (m)")
m = g.display(ax, "sla_raw", vmin=-0.15, vmax=0.15)
a_sla.display(ax, **kwargs_a_sla), c_sla.display(ax, **kwargs_c_sla)
ax.legend(), update_axes(ax, m)
Display detection
ax = start_axes("Eddies detected")
a_adt.display(ax, **kwargs_a_adt)
a_sla.display(ax, **kwargs_a_sla)
c_adt.display(ax, **kwargs_c_adt)
c_sla.display(ax, **kwargs_c_sla)
ax.legend()
update_axes(ax)
Match¶
Where cyclone meet anticyclone
i_c_adt, i_a_sla, c = c_adt.match(a_sla, cmin=0.01)
i_a_adt, i_c_sla, c = a_adt.match(c_sla, cmin=0.01)
ax = start_axes("Cyclone share area with anticyclone")
a_adt.index(i_a_adt).display(ax, **kwargs_a_adt)
c_adt.index(i_c_adt).display(ax, **kwargs_c_adt)
a_sla.index(i_a_sla).display(ax, **kwargs_a_sla)
c_sla.index(i_c_sla).display(ax, **kwargs_c_sla)
ax.legend()
update_axes(ax)
Scatter plot¶
i_a_adt, i_a_sla, c = a_adt.match(a_sla, cmin=0.1)
i_c_adt, i_c_sla, c = c_adt.match(c_sla, cmin=0.1)
where is lonely eddies
ax = start_axes("Eddies with no match")
a_adt.index(i_a_adt, reverse=True).display(ax, **kwargs_a_adt)
c_adt.index(i_c_adt, reverse=True).display(ax, **kwargs_c_adt)
a_sla.index(i_a_sla, reverse=True).display(ax, **kwargs_a_sla)
c_sla.index(i_c_sla, reverse=True).display(ax, **kwargs_c_sla)
ax.legend()
update_axes(ax)
fig = plt.figure(figsize=(12, 12))
fig.suptitle(f"Scatter plot (A : {i_a_adt.shape[0]}, C : {i_c_adt.shape[0]} matches)")
for i, (label, field, factor, stop) in enumerate(
(
("speed radius (km)", "radius_s", 0.001, 80),
("outter radius (km)", "radius_e", 0.001, 120),
("amplitude (cm)", "amplitude", 100, 25),
("speed max (cm/s)", "speed_average", 100, 25),
)
):
ax = fig.add_subplot(2, 2, i + 1, title=label)
ax.set_xlabel("Absolute Dynamic Topography")
ax.set_ylabel("Sea Level Anomaly")
ax.plot(
a_adt[field][i_a_adt] * factor,
a_sla[field][i_a_sla] * factor,
"r.",
label="Anticyclonic",
)
ax.plot(
c_adt[field][i_c_adt] * factor,
c_sla[field][i_c_sla] * factor,
"b.",
label="Cyclonic",
)
ax.set_aspect("equal"), ax.grid()
ax.plot((0, 1000), (0, 1000), "g")
ax.set_xlim(0, stop), ax.set_ylim(0, stop)
ax.legend()
Total running time of the script: ( 0 minutes 6.418 seconds)
Grid Manipulation¶
Note
Click here to download the full example code or to run this example in your browser via Binder
Select pixel in eddies¶
from matplotlib import pyplot as plt
from matplotlib.path import Path
from numpy import ones
from py_eddy_tracker import data
from py_eddy_tracker.dataset.grid import RegularGridDataset
from py_eddy_tracker.observations.observation import EddiesObservations
from py_eddy_tracker.poly import create_vertice
Load an eddy file which contains contours
a = EddiesObservations.load_file(data.get_path("Anticyclonic_20190223.nc"))
Load a grid where we want found pixels in eddies or out
g = RegularGridDataset(
data.get_path("nrt_global_allsat_phy_l4_20190223_20190226.nc"),
"longitude",
"latitude",
)
Out:
We assume pixel position of grid is center for /home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/pyEddyTracker-3.3.0-py3.7.egg/py_eddy_tracker/data/nrt_global_allsat_phy_l4_20190223_20190226.nc
For each contours, we will get pixels indice in contour.
fig = plt.figure(figsize=(12, 6))
ax = fig.add_axes((0.05, 0.05, 0.9, 0.9))
ax.set_aspect("equal")
ax.set_xlim(10, 70)
ax.set_ylim(-50, -25)
# We will used the outter contour
x_name, y_name = a.intern(False)
adt = g.grid("adt")
mask = ones(adt.shape, dtype="bool")
for eddy in a:
i, j = Path(create_vertice(eddy[x_name], eddy[y_name])).pixels_in(g)
mask[i, j] = False
adt.mask[:] += ~mask
g.display(ax, "adt")
a.display(ax, label="Anticyclonic", color="g", lw=1, extern_only=True)
fig = plt.figure(figsize=(12, 6))
ax = fig.add_axes((0.05, 0.05, 0.9, 0.9))
ax.set_aspect("equal")
ax.set_xlim(10, 70)
ax.set_ylim(-50, -25)
adt.mask[:] = mask
g.display(ax, "adt")
a.display(ax, label="Anticyclonic", color="g", lw=1, extern_only=True)
Total running time of the script: ( 0 minutes 1.355 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Grid filtering in PET¶
How filter work in py eddy tracker. This implementation maybe doesn’t respect state art, but …
We code a specific filter in order to filter grid with same wavelength at each pixel.
from matplotlib import pyplot as plt
from numpy import arange
from py_eddy_tracker import data
from py_eddy_tracker.dataset.grid import RegularGridDataset
def start_axes(title):
fig = plt.figure(figsize=(13, 5))
ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])
ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
ax.set_aspect("equal")
ax.set_title(title)
return ax
def update_axes(ax, mappable=None):
ax.grid()
if mappable:
plt.colorbar(mappable, cax=ax.figure.add_axes([0.95, 0.05, 0.01, 0.9]))
All information will be for regular grid
g = RegularGridDataset(
data.get_path("dt_med_allsat_phy_l4_20160515_20190101.nc"), "longitude", "latitude"
)
Out:
We assume pixel position of grid is center for /home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/pyEddyTracker-3.3.0-py3.7.egg/py_eddy_tracker/data/dt_med_allsat_phy_l4_20160515_20190101.nc
Kernel¶
Shape of kernel will increase in x, when latitude increase
fig = plt.figure(figsize=(12, 8))
for i, latitude in enumerate((15, 35, 55, 75)):
k = g.kernel_bessel(latitude, 500, order=3).T
ax0 = fig.add_subplot(
2,
2,
i + 1,
title=f"Kernel at {latitude}° of latitude\nfor 1/8° grid, shape : {k.shape}",
aspect="equal",
)
m = ax0.pcolormesh(k, vmin=-0.5, vmax=2, cmap="viridis_r")
plt.colorbar(m, cax=fig.add_axes((0.92, 0.05, 0.01, 0.9)))
Kernel along latitude
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(
111,
ylabel="Kernel weight",
xlabel="Latitude in °",
title="Kernel in latitude, centered at 0° of latitude ",
)
k = g.kernel_bessel(0, 500, order=3)
k_lat = k[k.shape[0] // 2 + 1]
nb = k_lat.shape[0] // 2
ax.plot(
arange(-nb * g.xstep, (nb + 0.5) * g.xstep, g.xstep), k_lat, label="Bessel kernel"
)
ax.legend()
ax.grid()
Kernel applying¶
Original grid
ax = start_axes("ADT")
m = g.display(ax, "adt", vmin=-0.15, vmax=0.15)
update_axes(ax, m)
We will select wavelength of 300 km
Low frequency
ax = start_axes("ADT low frequency")
g.copy("adt", "adt_low_300")
g.bessel_low_filter("adt_low_300", 300, order=3)
m = g.display(ax, "adt_low_300", vmin=-0.15, vmax=0.15)
update_axes(ax, m)
High frequency
ax = start_axes("ADT high frequency")
g.copy("adt", "adt_high_300")
g.bessel_high_filter("adt_high_300", 300, order=3)
m = g.display(ax, "adt_high_300", vmin=-0.15, vmax=0.15)
update_axes(ax, m)
Clues¶
wavelength : 80km
g.copy("adt", "adt_high_bessel")
g.bessel_high_filter("adt_high_bessel", 80, order=3)
g.copy("adt", "adt_low_bessel")
g.bessel_low_filter("adt_low_bessel", 80, order=3)
area = dict(llcrnrlon=11.75, urcrnrlon=21, llcrnrlat=33, urcrnrlat=36.75)
Spectrum
fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(111)
ax.set_title("Spectrum")
ax.set_xlabel("km")
for label in ("adt_high_bessel", "adt_low_bessel", "adt"):
lon_spec, lat_spec = g.spectrum_lonlat(label, area=area)
mappable = ax.loglog(*lat_spec, label=f"lat {label}")[0]
ax.loglog(
*lon_spec, label=f"lon {label}", color=mappable.get_color(), linestyle="--"
)
ax.set_xlim(10, 1000)
ax.set_ylim(1e-6, 1)
ax.set_xscale("log")
ax.legend()
ax.grid()
Out:
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/scipy/signal/spectral.py:1963: UserWarning: nperseg = 256 is greater than input length = 30, using nperseg = 30
.format(nperseg, input_length))
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/scipy/signal/spectral.py:1963: UserWarning: nperseg = 256 is greater than input length = 74, using nperseg = 74
.format(nperseg, input_length))
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/scipy/signal/spectral.py:1963: UserWarning: nperseg = 256 is greater than input length = 30, using nperseg = 30
.format(nperseg, input_length))
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/scipy/signal/spectral.py:1963: UserWarning: nperseg = 256 is greater than input length = 74, using nperseg = 74
.format(nperseg, input_length))
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/scipy/signal/spectral.py:1963: UserWarning: nperseg = 256 is greater than input length = 30, using nperseg = 30
.format(nperseg, input_length))
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/scipy/signal/spectral.py:1963: UserWarning: nperseg = 256 is greater than input length = 74, using nperseg = 74
.format(nperseg, input_length))
Spectrum ratio
fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(111)
ax.set_title("Spectrum ratio")
ax.set_xlabel("km")
for label in ("adt_high_bessel", "adt_low_bessel"):
lon_spec, lat_spec = g.spectrum_lonlat(label, area=area, ref=g, ref_grid_name="adt")
mappable = ax.plot(*lat_spec, label=f"lat {label}")[0]
ax.plot(*lon_spec, label=f"lon {label}", color=mappable.get_color(), linestyle="--")
ax.set_xlim(10, 1000)
ax.set_ylim(0, 1)
ax.set_xscale("log")
ax.legend()
ax.grid()
Out:
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/scipy/signal/spectral.py:1963: UserWarning: nperseg = 256 is greater than input length = 30, using nperseg = 30
.format(nperseg, input_length))
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/scipy/signal/spectral.py:1963: UserWarning: nperseg = 256 is greater than input length = 74, using nperseg = 74
.format(nperseg, input_length))
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/scipy/signal/spectral.py:1963: UserWarning: nperseg = 256 is greater than input length = 30, using nperseg = 30
.format(nperseg, input_length))
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/scipy/signal/spectral.py:1963: UserWarning: nperseg = 256 is greater than input length = 74, using nperseg = 74
.format(nperseg, input_length))
Note
Click here to download the full example code or to run this example in your browser via Binder
Get Okubo Weis¶
\[OW = S_n^2 + S_s^2 + \omega^2\]
with normal strain (\(S_n\)), shear strain (\(S_s\)) and vorticity (\(\omega\))
\[S_n = \frac{\partial u}{\partial x} - \frac{\partial v}{\partial y},
S_s = \frac{\partial v}{\partial x} + \frac{\partial u}{\partial y},
\omega = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}\]
from matplotlib import pyplot as plt
from numpy import arange, ma, where
from py_eddy_tracker import data
from py_eddy_tracker.dataset.grid import RegularGridDataset
from py_eddy_tracker.observations.observation import EddiesObservations
def start_axes(title, zoom=False):
fig = plt.figure(figsize=(12, 6))
axes = fig.add_axes([0.03, 0.03, 0.90, 0.94])
axes.set_xlim(0, 360), axes.set_ylim(-80, 80)
if zoom:
axes.set_xlim(270, 340), axes.set_ylim(20, 50)
axes.set_aspect("equal")
axes.set_title(title)
return axes
def update_axes(axes, mappable=None):
axes.grid()
if mappable:
plt.colorbar(mappable, cax=axes.figure.add_axes([0.94, 0.05, 0.01, 0.9]))
Load detection files
a = EddiesObservations.load_file(data.get_path("Anticyclonic_20190223.nc"))
c = EddiesObservations.load_file(data.get_path("Cyclonic_20190223.nc"))
Load Input grid, ADT will be used to detect eddies
g = RegularGridDataset(
data.get_path("nrt_global_allsat_phy_l4_20190223_20190226.nc"),
"longitude",
"latitude",
)
ax = start_axes("ADT (cm)")
m = g.display(ax, "adt", vmin=-120, vmax=120, factor=100)
update_axes(ax, m)
Out:
We assume pixel position of grid is center for /home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/pyEddyTracker-3.3.0-py3.7.egg/py_eddy_tracker/data/nrt_global_allsat_phy_l4_20190223_20190226.nc
Get parameter for ow
u_x = g.compute_stencil(g.grid("ugos"))
u_y = g.compute_stencil(g.grid("ugos"), vertical=True)
v_x = g.compute_stencil(g.grid("vgos"))
v_y = g.compute_stencil(g.grid("vgos"), vertical=True)
ow = g.vars["ow"] = (u_x - v_y) ** 2 + (v_x + u_y) ** 2 - (v_x - u_y) ** 2
ax = start_axes("Okubo weis")
m = g.display(ax, "ow", vmin=-1e-10, vmax=1e-10, cmap="bwr")
update_axes(ax, m)
Gulf stream zoom
ax = start_axes("Okubo weis, Gulf stream", zoom=True)
m = g.display(ax, "ow", vmin=-1e-10, vmax=1e-10, cmap="bwr")
kw_ed = dict(intern_only=True, color="k", lw=1)
a.display(ax, **kw_ed), c.display(ax, **kw_ed)
update_axes(ax, m)
only negative OW
ax = start_axes("Okubo weis, Gulf stream", zoom=True)
threshold = ow.std() * -0.2
ow = ma.array(ow, mask=ow > threshold)
m = g.display(ax, ow, vmin=-1e-10, vmax=1e-10, cmap="bwr")
a.display(ax, **kw_ed), c.display(ax, **kw_ed)
update_axes(ax, m)
Get okubo-weiss mean/min/center in eddies
plt.figure(figsize=(8, 6))
ax = plt.subplot(111)
ax.set_xlabel("Okubo-Weiss parameter")
kw_hist = dict(bins=arange(-20e-10, 20e-10, 50e-12), histtype="step")
for method in ("mean", "center", "min"):
kw_interp = dict(grid_object=g, varname="ow", method=method, intern=True)
_, _, m = ax.hist(
a.interp_grid(**kw_interp), label=f"Anticyclonic - OW {method}", **kw_hist
)
ax.hist(
c.interp_grid(**kw_interp),
label=f"Cyclonic - OW {method}",
color=m[0].get_edgecolor(),
ls="--",
**kw_hist,
)
ax.axvline(threshold, color="r")
ax.set_yscale("log")
ax.grid()
ax.set_ylim(1, 1e4)
ax.set_xlim(-15e-10, 15e-10)
ax.legend()
Catch eddies with bad OW
ax = start_axes("Eddies with a min OW in speed contour over threshold")
ow_min = a.interp_grid(**kw_interp)
a_bad_ow = a.index(where(ow_min > threshold)[0])
a_bad_ow.display(ax, color="r", label="Anticyclonic")
ow_min = c.interp_grid(**kw_interp)
c_bad_ow = c.index(where(ow_min > threshold)[0])
c_bad_ow.display(ax, color="b", label="Cyclonic")
ax.legend()
Display Radius and amplitude of eddies
fig = plt.figure(figsize=(12, 5))
fig.suptitle(
"Parameter distribution (solid line) and cumulative distribution (dashed line)"
)
ax_amp, ax_rad = fig.add_subplot(121), fig.add_subplot(122)
ax_amp_c, ax_rad_c = ax_amp.twinx(), ax_rad.twinx()
ax_amp_c.set_ylim(0, 1), ax_rad_c.set_ylim(0, 1)
kw_a = dict(xname="amplitude", bins=arange(0, 2, 0.002).astype("f4"))
kw_r = dict(xname="radius_s", bins=arange(0, 500e6, 2e3).astype("f4"))
for d, label, color in (
(a, "Anticyclonic all", "r"),
(a_bad_ow, "Anticyclonic bad OW", "orange"),
(c, "Cyclonic all", "blue"),
(c_bad_ow, "Cyclonic bad OW", "lightblue"),
):
x, y = d.bins_stat(**kw_a)
ax_amp.plot(x * 100, y, label=label, color=color)
ax_amp_c.plot(
x * 100, y.cumsum() / y.sum(), label=label, color=color, ls="-.", lw=0.5
)
x, y = d.bins_stat(**kw_r)
ax_rad.plot(x * 1e-3, y, label=label, color=color)
ax_rad_c.plot(
x * 1e-3, y.cumsum() / y.sum(), label=label, color=color, ls="-.", lw=0.5
)
ax_amp.set_xlim(0, 12.5), ax_amp.grid(), ax_amp.set_ylim(0), ax_amp.legend()
ax_rad.set_xlim(0, 120), ax_rad.grid(), ax_rad.set_ylim(0)
ax_amp.set_xlabel("Amplitude (cm)"), ax_amp.set_ylabel("Nb eddies")
ax_rad.set_xlabel("Speed radius (km)")
Total running time of the script: ( 0 minutes 7.616 seconds)
Tracking Manipulation¶
Note
Click here to download the full example code or to run this example in your browser via Binder
Track animation¶
Run in a terminal this script, which allow to watch eddy evolution
import py_eddy_tracker_sample
from py_eddy_tracker.appli.gui import Anim
from py_eddy_tracker.observations.tracking import TrackEddiesObservations
Load experimental atlas, and we select one eddy
a = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Anticyclonic.zarr")
)
# We get only 300 first step to save time of documentation builder
eddy = a.extract_ids([9672]).index(slice(0, 300))
Run animation Key shortcut :
Escape => exit
SpaceBar => pause
left arrow => t - 1
right arrow => t + 1
+ => speed increase of 10 %
- => speed decrease of 10 %
a = Anim(eddy, sleep_event=1e-10, intern=True, figsize=(8, 3.5), cmap="viridis")
a.txt.set_position((17, 34.6))
a.ax.set_xlim(16.5, 23)
a.ax.set_ylim(34.5, 37)
a.show(infinity_loop=False)
Total running time of the script: ( 0 minutes 6.810 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Display fields¶
import py_eddy_tracker_sample
from matplotlib import pyplot as plt
from py_eddy_tracker.observations.tracking import TrackEddiesObservations
Load an experimental cyclonic atlas, we keep only eddies which are follow more than 180 days
c = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Cyclonic.zarr")
)
c = c.extract_with_length((180, -1))
Plot amplitude field
fig = plt.figure(figsize=(12, 6))
ax = fig.add_axes((0.05, 0.1, 0.9, 0.9))
ax.set_aspect("equal")
ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
m = c.scatter(ax, "amplitude", ref=-10, vmin=0, vmax=0.1)
ax.grid()
cb = plt.colorbar(
m, cax=fig.add_axes([0.05, 0.07, 0.9, 0.01]), orientation="horizontal"
)
cb.set_label("Amplitude (m)")
Total running time of the script: ( 0 minutes 2.271 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Track animation with standard matplotlib¶
Run in a terminal this script, which allow to watch eddy evolution
import py_eddy_tracker_sample
from matplotlib.animation import FuncAnimation
from numpy import arange
from py_eddy_tracker.appli.gui import Anim
from py_eddy_tracker.observations.tracking import TrackEddiesObservations
Load experimental atlas, and we select one eddy
a = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Anticyclonic.zarr")
)
eddy = a.extract_ids([9672])
Run animation
a = Anim(eddy, intern=True, figsize=(8, 3.5), cmap="magma_r", nb_step=6)
a.txt.set_position((17, 34.6))
a.ax.set_xlim(16.5, 23)
a.ax.set_ylim(34.5, 37)
# arguments to get full animation
# kwargs = dict(frames=arange(*a.period), interval=50)
# arguments to reduce compute cost for doucmentation, we display only every 10 days
kwargs = dict(frames=arange(*a.period)[200:800:10], save_count=60, interval=200)
ani = FuncAnimation(a.fig, a.func_animation, **kwargs)
Total running time of the script: ( 0 minutes 16.154 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Display Tracks¶
import py_eddy_tracker_sample
from matplotlib import pyplot as plt
from py_eddy_tracker.observations.tracking import TrackEddiesObservations
Load experimental atlas
a = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Anticyclonic.zarr")
)
c = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Cyclonic.zarr")
)
print(a)
Out:
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/pyEddyTracker-3.3.0-py3.7.egg/py_eddy_tracker/observations/observation.py:220: RuntimeWarning: invalid value encountered in true_divide
v = v.astype("f4") / v.sum() * 100
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/pyEddyTracker-3.3.0-py3.7.egg/py_eddy_tracker/observations/observation.py:222: RuntimeWarning: invalid value encountered in true_divide
v /= hist_numba(self[x], bins=bins)[0]
| 473554 observations from 15706 to 25334 (9629 days, ~49 obs/day)
| Speed area : 0.00 Mkm²/day
| Effective area : 0.00 Mkm²/day
----Distribution in Amplitude:
| Amplitude bounds (cm) 0.00 1.00 2.00 3.00 4.00 5.00 10.00 500.00
| Percent of eddies : 21.69 25.72 15.32 9.95 6.20 14.80 6.32
----Distribution in Radius:
| Speed radius (km) 0.00 15.00 30.00 45.00 60.00 75.00 100.00 200.00 2000.00
| Percent of eddies : 0.85 48.58 36.33 10.97 2.59 0.64 0.03 0.00
| Effective radius (km) 0.00 15.00 30.00 45.00 60.00 75.00 100.00 200.00 2000.00
| Percent of eddies : 100.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
----Distribution in Latitude
Latitude bounds -90.00 -60.00 -15.00 15.00 60.00 90.00
Percent of eddies : 0.00 0.00 0.00 100.00 0.00
Percent of speed area : nan nan nan nan nan
Percent of effective area : nan nan nan nan nan
Mean speed radius (km) : nan nan nan 32.28 nan
Mean effective radius (km): nan nan nan 0.00 nan
Mean amplitude (cm) : nan nan nan 3.51 nan
| 12224 tracks (38.74 obs/tracks, shorter 10 obs, longer 1183 obs)
| 27246 filled observations (2.23 obs/tracks, 5.75 % of total)
| Intepolated speed area : 0.00 Mkm²/day
| Intepolated effective area : 0.00 Mkm²/day
| Distance by day : Mean 3.68 , Median 2.53 km/day
| Distance by track : Mean 139.01 , Median 77.28 km/track
----Distribution in lifetime:
| Lifetime (days ) 1.00 30.00 90.00 180.00 270.00 365.00 1000.00 10000.00
| Percent of tracks : 65.03 26.87 5.72 1.20 0.51 0.65 0.03
| Percent of eddies : 27.88 33.88 17.90 6.74 4.02 8.65 0.93
keep only eddies longer than 20 weeks, use -1 to have no upper limit
a = a.extract_with_length((7 * 20, -1))
c = c.extract_with_length((7 * 20, -1))
print(a)
Out:
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/pyEddyTracker-3.3.0-py3.7.egg/py_eddy_tracker/observations/observation.py:220: RuntimeWarning: invalid value encountered in true_divide
v = v.astype("f4") / v.sum() * 100
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/pyEddyTracker-3.3.0-py3.7.egg/py_eddy_tracker/observations/observation.py:222: RuntimeWarning: invalid value encountered in true_divide
v /= hist_numba(self[x], bins=bins)[0]
| 124198 observations from 15706 to 25334 (9629 days, ~13 obs/day)
| Speed area : 0.00 Mkm²/day
| Effective area : 0.00 Mkm²/day
----Distribution in Amplitude:
| Amplitude bounds (cm) 0.00 1.00 2.00 3.00 4.00 5.00 10.00 500.00
| Percent of eddies : 5.11 9.99 12.06 11.79 9.32 31.51 20.21
----Distribution in Radius:
| Speed radius (km) 0.00 15.00 30.00 45.00 60.00 75.00 100.00 200.00 2000.00
| Percent of eddies : 0.19 22.20 46.88 22.96 6.31 1.39 0.07 0.00
| Effective radius (km) 0.00 15.00 30.00 45.00 60.00 75.00 100.00 200.00 2000.00
| Percent of eddies : 100.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
----Distribution in Latitude
Latitude bounds -90.00 -60.00 -15.00 15.00 60.00 90.00
Percent of eddies : 0.00 0.00 0.00 100.00 0.00
Percent of speed area : nan nan nan nan nan
Percent of effective area : nan nan nan nan nan
Mean speed radius (km) : nan nan nan 40.03 nan
Mean effective radius (km): nan nan nan 0.00 nan
Mean amplitude (cm) : nan nan nan 6.62 nan
| 469 tracks (264.81 obs/tracks, shorter 140 obs, longer 1183 obs)
| 3348 filled observations (7.14 obs/tracks, 2.70 % of total)
| Intepolated speed area : 0.00 Mkm²/day
| Intepolated effective area : 0.00 Mkm²/day
| Distance by day : Mean 3.13 , Median 2.00 km/day
| Distance by track : Mean 825.11 , Median 695.99 km/track
----Distribution in lifetime:
| Lifetime (days ) 1.00 30.00 90.00 180.00 270.00 365.00 1000.00 10000.00
| Percent of tracks : 0.00 0.00 37.74 31.34 13.22 16.84 0.85
| Percent of eddies : 0.00 0.00 22.44 25.72 15.34 32.97 3.53
Position filtering for nice display
a.position_filter(median_half_window=1, loess_half_window=5)
c.position_filter(median_half_window=1, loess_half_window=5)
Plot
fig = plt.figure(figsize=(12, 5))
ax = fig.add_axes((0.05, 0.1, 0.9, 0.9))
ax.set_aspect("equal")
ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
a.plot(ax, ref=-10, label="Anticyclonic ({nb_tracks} tracks)", color="r", lw=0.1)
c.plot(ax, ref=-10, label="Cyclonic ({nb_tracks} tracks)", color="b", lw=0.1)
ax.legend()
ax.grid()
Total running time of the script: ( 0 minutes 6.547 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
One Track¶
import py_eddy_tracker_sample
from matplotlib import pyplot as plt
from py_eddy_tracker.observations.tracking import TrackEddiesObservations
Load experimental atlas, and we select one eddy
a = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Anticyclonic.zarr")
)
eddy = a.extract_ids([9672])
eddy_f = a.extract_ids([9672])
eddy_f.position_filter(median_half_window=1, loess_half_window=5)
plot
fig = plt.figure(figsize=(12, 5))
ax = fig.add_axes((0.05, 0.05, 0.9, 0.9))
ax.set_xlim(17.5, 22)
ax.set_ylim(35, 36.5)
ax.set_aspect("equal")
ax.grid()
eddy.plot(ax, color="r", lw=0.5)
eddy_f.scatter(ax, "n", cmap="jet", s=80)
fig = plt.figure(figsize=(12, 5))
ax = fig.add_axes((0.05, 0.05, 0.9, 0.9))
ax.set_xlim(17, 23)
ax.set_ylim(34.5, 37)
ax.set_aspect("equal")
ax.grid()
eddy.plot(ax, color="r", lw=0.5, label="track")
eddy.index(range(0, len(eddy), 40)).display(
ax, intern_only=True, label="observations every 40 days"
)
ax.legend()
Total running time of the script: ( 0 minutes 1.384 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Tracks which go through area¶
import py_eddy_tracker_sample
from matplotlib import pyplot as plt
from py_eddy_tracker.observations.tracking import TrackEddiesObservations
Load experimental atlas, we filter position to have nice display
c = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Cyclonic.zarr")
)
c.position_filter(median_half_window=1, loess_half_window=5)
We extract eddies in the area set below, but we ask to keep full_path
x0, x1, y0, y1 = 3, 4, 37, 38
area = dict(llcrnrlon=x0, llcrnrlat=y0, urcrnrlon=x1, urcrnrlat=y1)
c_subset = c.extract_with_area(area, full_path=True)
Plot
fig = plt.figure(figsize=(12, 5))
ax = fig.add_axes((0.05, 0.05, 0.9, 0.9))
ax.set_xlim(-1, 9)
ax.set_ylim(36, 40)
ax.set_aspect("equal")
ax.grid()
c.plot(ax, color="gray", lw=0.1, ref=-10, label="All tracks ({nb_tracks} tracks)")
c_subset.plot(
ax, color="red", lw=0.2, ref=-10, label="selected tracks ({nb_tracks} tracks)"
)
ax.plot(
(x0, x0, x1, x1, x0),
(y0, y1, y1, y0, y0),
color="green",
lw=1.5,
label="Box of selection",
)
ax.legend()
Total running time of the script: ( 0 minutes 2.784 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Track in python¶
This example didn’t replace EddyTracking, we remove check that application do and also postprocessing step.
from py_eddy_tracker.data import get_remote_sample
from py_eddy_tracker.featured_tracking.area_tracker import AreaTracker
from py_eddy_tracker.gui import GUI
from py_eddy_tracker.tracking import Correspondances
Get remote data, we will keep only 180 first days, get_remote_sample function is only to get demo dataset, in your own case give a list of identification filename and don’t mix cyclonic and anticyclonic files.
file_objects = get_remote_sample(
"eddies_med_adt_allsat_dt2018/Anticyclonic_2010_2011_2012"
)[:180]
We run a traking with a tracker which use contour overlap
c = Correspondances(datasets=file_objects, class_method=AreaTracker, virtual=3)
c.track()
c.prepare_merging()
# We have now an eddy object
eddies_area_tracker = c.merge(raw_data=False)
eddies_area_tracker.virtual[:] = eddies_area_tracker.time == 0
eddies_area_tracker.filled_by_interpolation(eddies_area_tracker.virtual == 1)
We run a traking with default tracker
c = Correspondances(datasets=file_objects, virtual=3)
c.track()
c.prepare_merging()
eddies_default_tracker = c.merge(raw_data=False)
eddies_default_tracker.virtual[:] = eddies_default_tracker.time == 0
eddies_default_tracker.filled_by_interpolation(eddies_default_tracker.virtual == 1)
Out:
High number of conflict : 56 (nb_conflict)
High number of conflict : 46 (nb_conflict)
High number of conflict : 49 (nb_conflict)
High number of conflict : 50 (nb_conflict)
High number of conflict : 58 (nb_conflict)
High number of conflict : 62 (nb_conflict)
High number of conflict : 67 (nb_conflict)
High number of conflict : 67 (nb_conflict)
High number of conflict : 51 (nb_conflict)
High number of conflict : 50 (nb_conflict)
High number of conflict : 54 (nb_conflict)
High number of conflict : 60 (nb_conflict)
High number of conflict : 59 (nb_conflict)
High number of conflict : 61 (nb_conflict)
High number of conflict : 68 (nb_conflict)
High number of conflict : 74 (nb_conflict)
High number of conflict : 64 (nb_conflict)
High number of conflict : 71 (nb_conflict)
High number of conflict : 67 (nb_conflict)
High number of conflict : 67 (nb_conflict)
High number of conflict : 71 (nb_conflict)
High number of conflict : 78 (nb_conflict)
High number of conflict : 76 (nb_conflict)
High number of conflict : 78 (nb_conflict)
High number of conflict : 71 (nb_conflict)
High number of conflict : 66 (nb_conflict)
High number of conflict : 55 (nb_conflict)
High number of conflict : 60 (nb_conflict)
High number of conflict : 59 (nb_conflict)
High number of conflict : 67 (nb_conflict)
High number of conflict : 57 (nb_conflict)
High number of conflict : 57 (nb_conflict)
High number of conflict : 37 (nb_conflict)
High number of conflict : 72 (nb_conflict)
High number of conflict : 75 (nb_conflict)
High number of conflict : 73 (nb_conflict)
High number of conflict : 77 (nb_conflict)
High number of conflict : 90 (nb_conflict)
High number of conflict : 89 (nb_conflict)
High number of conflict : 88 (nb_conflict)
High number of conflict : 95 (nb_conflict)
High number of conflict : 87 (nb_conflict)
High number of conflict : 78 (nb_conflict)
High number of conflict : 73 (nb_conflict)
High number of conflict : 83 (nb_conflict)
High number of conflict : 86 (nb_conflict)
High number of conflict : 89 (nb_conflict)
High number of conflict : 72 (nb_conflict)
High number of conflict : 65 (nb_conflict)
High number of conflict : 48 (nb_conflict)
High number of conflict : 83 (nb_conflict)
High number of conflict : 81 (nb_conflict)
High number of conflict : 77 (nb_conflict)
High number of conflict : 89 (nb_conflict)
High number of conflict : 84 (nb_conflict)
High number of conflict : 89 (nb_conflict)
High number of conflict : 99 (nb_conflict)
High number of conflict : 79 (nb_conflict)
High number of conflict : 79 (nb_conflict)
High number of conflict : 75 (nb_conflict)
High number of conflict : 75 (nb_conflict)
High number of conflict : 65 (nb_conflict)
High number of conflict : 80 (nb_conflict)
High number of conflict : 50 (nb_conflict)
High number of conflict : 73 (nb_conflict)
High number of conflict : 73 (nb_conflict)
High number of conflict : 72 (nb_conflict)
High number of conflict : 75 (nb_conflict)
High number of conflict : 62 (nb_conflict)
High number of conflict : 53 (nb_conflict)
High number of conflict : 45 (nb_conflict)
High number of conflict : 49 (nb_conflict)
High number of conflict : 41 (nb_conflict)
High number of conflict : 50 (nb_conflict)
High number of conflict : 46 (nb_conflict)
High number of conflict : 37 (nb_conflict)
High number of conflict : 52 (nb_conflict)
Start GUI to compare tracking
g = GUI(
Acyc_area_tracker=eddies_area_tracker, Acyc_default_tracker=eddies_default_tracker
)
g.now = 22000
g.bbox = 0, 9, 36, 40
g.adjust()
g.show()
Start GUI with area tracker
g = GUI(Acyc_area_tracker=eddies_area_tracker)
g.now = 22000
g.bbox = 0, 9, 36, 40
g.adjust()
g.show()
Start GUI with default one
g = GUI(Acyc_default_tracker=eddies_default_tracker)
g.now = 22000
g.bbox = 0, 9, 36, 40
g.adjust()
g.show()
Total running time of the script: ( 0 minutes 39.290 seconds)
Tracking diagnostics¶
Note
Click here to download the full example code or to run this example in your browser via Binder
Geographical statistics¶
import py_eddy_tracker_sample
from matplotlib import pyplot as plt
from py_eddy_tracker.observations.tracking import TrackEddiesObservations
def start_axes(title):
fig = plt.figure(figsize=(13.5, 5))
ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])
ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
ax.set_aspect("equal")
ax.set_title(title)
return ax
Load an experimental med atlas over a period of 26 years (1993-2019), we merge the 2 datasets
a = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Anticyclonic.zarr")
)
c = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Cyclonic.zarr")
)
a = a.merge(c)
step = 0.1
Mean of amplitude in each box
ax = start_axes("Amplitude mean by box of %s°" % step)
g = a.grid_stat(((-7, 37, step), (30, 46, step)), "amplitude")
m = g.display(ax, name="amplitude", vmin=0, vmax=10, factor=100)
ax.grid()
cb = plt.colorbar(m, cax=ax.figure.add_axes([0.94, 0.05, 0.01, 0.9]))
cb.set_label("Amplitude (cm)")
Out:
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/pyEddyTracker-3.3.0-py3.7.egg/py_eddy_tracker/dataset/grid.py:1896: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3. Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading']. This will become an error two minor releases later.
return ax.pcolormesh(x, self.y_bounds, data.T * factor, **kwargs)
Mean of speed radius in each box
ax = start_axes("Speed radius mean by box of %s°" % step)
g = a.grid_stat(((-7, 37, step), (30, 46, step)), "radius_s")
m = g.display(ax, name="radius_s", vmin=10, vmax=50, factor=0.001)
ax.grid()
cb = plt.colorbar(m, cax=ax.figure.add_axes([0.94, 0.05, 0.01, 0.9]))
cb.set_label("Speed radius (km)")
Out:
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/pyEddyTracker-3.3.0-py3.7.egg/py_eddy_tracker/dataset/grid.py:1896: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3. Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading']. This will become an error two minor releases later.
return ax.pcolormesh(x, self.y_bounds, data.T * factor, **kwargs)
Percent of virtual on the whole obs in each box
ax = start_axes("Percent of virtual by box of %s°" % step)
g = a.grid_stat(((-7, 37, step), (30, 46, step)), "virtual")
g.vars["virtual"] *= 100
m = g.display(ax, name="virtual", vmin=0, vmax=15)
ax.grid()
cb = plt.colorbar(m, cax=ax.figure.add_axes([0.94, 0.05, 0.01, 0.9]))
cb.set_label("Percent of virtual (%)")
Out:
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/pyEddyTracker-3.3.0-py3.7.egg/py_eddy_tracker/dataset/grid.py:1896: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3. Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading']. This will become an error two minor releases later.
return ax.pcolormesh(x, self.y_bounds, data.T * factor, **kwargs)
Total running time of the script: ( 0 minutes 4.987 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Birth and death¶
Following figures are based on https://doi.org/10.1016/j.pocean.2011.01.002
import py_eddy_tracker_sample
from matplotlib import pyplot as plt
from py_eddy_tracker.observations.tracking import TrackEddiesObservations
def start_axes(title):
fig = plt.figure(figsize=(13, 5))
ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])
ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
ax.set_aspect("equal")
ax.set_title(title)
return ax
def update_axes(ax, mappable=None):
ax.grid()
if mappable:
plt.colorbar(mappable, cax=ax.figure.add_axes([0.95, 0.05, 0.01, 0.9]))
Load an experimental med atlas over a period of 26 years (1993-2019)
kwargs_load = dict(
include_vars=(
"longitude",
"latitude",
"observation_number",
"track",
"time",
"speed_contour_longitude",
"speed_contour_latitude",
)
)
a = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Anticyclonic.zarr")
)
c = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Cyclonic.zarr")
)
t0, t1 = a.period
step = 0.125
bins = ((-10, 37, step), (30, 46, step))
kwargs = dict(cmap="terrain_r", factor=100 / (t1 - t0), name="count", vmin=0, vmax=1)
Cyclonic¶
ax = start_axes("Birth cyclonic frenquency (%)")
g_c_first = c.first_obs().grid_count(bins, intern=True)
m = g_c_first.display(ax, **kwargs)
update_axes(ax, m)
ax = start_axes("Death cyclonic frenquency (%)")
g_c_last = c.last_obs().grid_count(bins, intern=True)
m = g_c_last.display(ax, **kwargs)
update_axes(ax, m)
Anticyclonic¶
ax = start_axes("Birth anticyclonic frequency (%)")
g_a_first = a.first_obs().grid_count(bins, intern=True)
m = g_a_first.display(ax, **kwargs)
update_axes(ax, m)
ax = start_axes("Death anticyclonic frequency (%)")
g_a_last = a.last_obs().grid_count(bins, intern=True)
m = g_a_last.display(ax, **kwargs)
update_axes(ax, m)
Total running time of the script: ( 0 minutes 5.135 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Lifetime Histogram¶
import py_eddy_tracker_sample
from matplotlib import pyplot as plt
from numpy import arange, ones
from py_eddy_tracker.observations.tracking import TrackEddiesObservations
Load an experimental med atlas over a period of 26 years (1993-2019)
a = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Anticyclonic.zarr")
)
c = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Cyclonic.zarr")
)
nb_year = (a.period[1] - a.period[0] + 1) / 365.25
Setup axes
figure = plt.figure(figsize=(12, 8))
ax_ratio_cum = figure.add_axes([0.55, 0.06, 0.42, 0.34])
ax_ratio = figure.add_axes([0.07, 0.06, 0.46, 0.34])
ax_cum = figure.add_axes([0.55, 0.43, 0.42, 0.54])
ax = figure.add_axes([0.07, 0.43, 0.46, 0.54])
ax.set_ylabel("Eddies by year")
ax_ratio.set_ylabel("Ratio Cyclonic/Anticyclonic")
for ax_ in (ax, ax_cum, ax_ratio_cum, ax_ratio):
ax_.set_xlim(0, 400)
if ax_ in (ax, ax_cum):
ax_.set_ylim(1e-1, 1e4), ax_.set_yscale("log")
else:
ax_.set_xlabel("Lifetime in days (by week bins)")
ax_.set_ylim(0, 2)
ax_.axhline(1, color="g", lw=2)
ax_.grid()
ax_cum.xaxis.set_ticklabels([]), ax_cum.yaxis.set_ticklabels([])
ax.xaxis.set_ticklabels([]), ax_ratio_cum.yaxis.set_ticklabels([])
# plot data
bin_hist = arange(7, 2000, 7)
x = (bin_hist[1:] + bin_hist[:-1]) / 2.0
a_nb, c_nb = a.nb_obs_by_track, c.nb_obs_by_track
a_nb, c_nb = a_nb[a_nb != 0], c_nb[c_nb != 0]
w_a, w_c = ones(a_nb.shape) / nb_year, ones(c_nb.shape) / nb_year
kwargs_a = dict(histtype="step", bins=bin_hist, x=a_nb, color="r", weights=w_a)
kwargs_c = dict(histtype="step", bins=bin_hist, x=c_nb, color="b", weights=w_c)
cum_a, _, _ = ax_cum.hist(cumulative=-1, **kwargs_a)
cum_c, _, _ = ax_cum.hist(cumulative=-1, **kwargs_c)
nb_a, _, _ = ax.hist(label="Anticyclonic", **kwargs_a)
nb_c, _, _ = ax.hist(label="Cyclonic", **kwargs_c)
ax_ratio_cum.plot(x, cum_c / cum_a)
ax_ratio.plot(x, nb_c / nb_a)
ax.legend()
Out:
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/checkouts/v3.3.0/examples/10_tracking_diagnostics/pet_lifetime.py:55: RuntimeWarning: invalid value encountered in true_divide
ax_ratio_cum.plot(x, cum_c / cum_a)
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/checkouts/v3.3.0/examples/10_tracking_diagnostics/pet_lifetime.py:56: RuntimeWarning: divide by zero encountered in true_divide
ax_ratio.plot(x, nb_c / nb_a)
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/checkouts/v3.3.0/examples/10_tracking_diagnostics/pet_lifetime.py:56: RuntimeWarning: invalid value encountered in true_divide
ax_ratio.plot(x, nb_c / nb_a)
Total running time of the script: ( 0 minutes 3.156 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Parameter Histogram¶
import py_eddy_tracker_sample
from matplotlib import pyplot as plt
from numpy import arange
from py_eddy_tracker.observations.tracking import TrackEddiesObservations
Load an experimental med atlas over a period of 26 years (1993-2019)
a = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Anticyclonic.zarr")
)
c = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Cyclonic.zarr")
)
kwargs_a = dict(label="Anticyclonic", color="r", histtype="step", density=True)
kwargs_c = dict(label="Cyclonic", color="b", histtype="step", density=True)
Plot
fig = plt.figure(figsize=(12, 7))
for x0, name, title, xmax, factor, bins in zip(
(0.4, 0.72, 0.08),
("speed_radius", "speed_average", "amplitude"),
("Speed radius (km)", "Speed average (cm/s)", "Amplitude (cm)"),
(100, 50, 20),
(0.001, 100, 100),
(arange(0, 2000, 1), arange(0, 1000, 0.5), arange(0.0005, 1000, 0.2)),
):
ax_hist = fig.add_axes((x0, 0.24, 0.27, 0.35))
nb_a, _, _ = ax_hist.hist(a[name] * factor, bins=bins, **kwargs_a)
nb_c, _, _ = ax_hist.hist(c[name] * factor, bins=bins, **kwargs_c)
ax_hist.set_xticklabels([])
ax_hist.set_xlim(0, xmax)
ax_hist.grid()
ax_cum = fig.add_axes((x0, 0.62, 0.27, 0.35))
ax_cum.hist(a[name] * factor, bins=bins, cumulative=-1, **kwargs_a)
ax_cum.hist(c[name] * factor, bins=bins, cumulative=-1, **kwargs_c)
ax_cum.set_xticklabels([])
ax_cum.set_title(title)
ax_cum.set_xlim(0, xmax)
ax_cum.set_ylim(0, 1)
ax_cum.grid()
ax_ratio = fig.add_axes((x0, 0.06, 0.27, 0.15))
ax_ratio.set_xlim(0, xmax)
ax_ratio.set_ylim(0, 2)
ax_ratio.plot((bins[1:] + bins[:-1]) / 2, nb_c / nb_a)
ax_ratio.axhline(1, color="k")
ax_ratio.grid()
ax_ratio.set_xlabel(title)
ax_cum.set_ylabel("Cumulative\npercent distribution")
ax_hist.set_ylabel("Percent of observations")
ax_ratio.set_ylabel("Ratio percent\nCyc/Acyc")
ax_cum.legend()
Out:
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/checkouts/v3.3.0/examples/10_tracking_diagnostics/pet_histo.py:54: RuntimeWarning: divide by zero encountered in true_divide
ax_ratio.plot((bins[1:] + bins[:-1]) / 2, nb_c / nb_a)
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/checkouts/v3.3.0/examples/10_tracking_diagnostics/pet_histo.py:54: RuntimeWarning: invalid value encountered in true_divide
ax_ratio.plot((bins[1:] + bins[:-1]) / 2, nb_c / nb_a)
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/checkouts/v3.3.0/examples/10_tracking_diagnostics/pet_histo.py:54: RuntimeWarning: divide by zero encountered in true_divide
ax_ratio.plot((bins[1:] + bins[:-1]) / 2, nb_c / nb_a)
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/checkouts/v3.3.0/examples/10_tracking_diagnostics/pet_histo.py:54: RuntimeWarning: invalid value encountered in true_divide
ax_ratio.plot((bins[1:] + bins[:-1]) / 2, nb_c / nb_a)
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/checkouts/v3.3.0/examples/10_tracking_diagnostics/pet_histo.py:54: RuntimeWarning: invalid value encountered in true_divide
ax_ratio.plot((bins[1:] + bins[:-1]) / 2, nb_c / nb_a)
Total running time of the script: ( 0 minutes 3.634 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Groups distribution¶
import py_eddy_tracker_sample
from matplotlib import pyplot as plt
from numpy import arange, ones, percentile
from py_eddy_tracker.observations.tracking import TrackEddiesObservations
Load an experimental med atlas over a period of 26 years (1993-2019)
a = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Anticyclonic.zarr")
)
Group distribution
groups = dict()
bins_time = [10, 20, 30, 60, 90, 180, 360, 100000]
for t0, t1 in zip(bins_time[:-1], bins_time[1:]):
groups[f"lifetime_{t0}_{t1}"] = lambda dataset, t0=t0, t1=t1: (
dataset.lifetime >= t0
) * (dataset.lifetime < t1)
bins_percentile = arange(0, 100.0001, 5)
Function to build stats
def stats_compilation(dataset, groups, field, bins, filter=None):
datas = dict(ref=dataset.bins_stat(field, bins=bins, mask=filter)[1], y=dict())
for k, index in groups.items():
i = dataset.merge_filters(filter, index)
x, datas["y"][k] = dataset.bins_stat(field, bins=bins, mask=i)
datas["x"], datas["bins"] = x, bins
return datas
def plot_stats(ax, bins, x, y, ref, box=False, cmap=None, percentiles=None, **kw):
base, ref = ones(x.shape) * 100.0, ref / 100.0
x = arange(bins.shape[0]).repeat(2)[1:-1] if box else x
y0 = base
if cmap is not None:
cmap, nb_groups = plt.get_cmap(cmap), len(y)
keys = tuple(y.keys())
for i, k in enumerate(keys[::-1]):
y1 = y0 - y[k] / ref
args = (y0.repeat(2), y1.repeat(2)) if box else (y0, y1)
if cmap is not None:
kw["color"] = cmap(1 - i / (nb_groups - 1))
ax.fill_between(x, *args, label=k, **kw)
y0 = y1
if percentiles:
for b in bins:
ax.axvline(b, **percentiles)
Speed radius by track period
stats = stats_compilation(
a, groups, "radius_s", percentile(a.radius_s, bins_percentile)
)
fig = plt.figure()
ax = fig.add_subplot(111)
plot_stats(ax, **stats, cmap="magma", percentiles=dict(color="gray", ls="-.", lw=0.4))
ax.set_xlabel("Speed radius (m)"), ax.set_ylabel("% of class"), ax.set_ylim(0, 100)
ax.grid(), ax.legend()
Amplitude by track period
stats = stats_compilation(
a, groups, "amplitude", percentile(a.amplitude, bins_percentile)
)
fig = plt.figure()
ax = fig.add_subplot(111)
plot_stats(ax, **stats, cmap="magma")
ax.set_xlabel("Amplitude (m)"), ax.set_ylabel("% of class"), ax.set_ylim(0, 100)
ax.grid(), ax.legend()
Total running time of the script: ( 0 minutes 2.031 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Propagation Histogram¶
import py_eddy_tracker_sample
from matplotlib import pyplot as plt
from numpy import arange, ones
from py_eddy_tracker.generic import cumsum_by_track
from py_eddy_tracker.observations.tracking import TrackEddiesObservations
Load an experimental med atlas over a period of 26 years (1993-2019)
a = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Anticyclonic.zarr")
)
c = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Cyclonic.zarr")
)
nb_year = (a.period[1] - a.period[0] + 1) / 365.25
Filtering position to remove noisy position
a.position_filter(median_half_window=1, loess_half_window=5)
c.position_filter(median_half_window=1, loess_half_window=5)
Compute curvilign distance
i0, nb = a.index_from_track, a.nb_obs_by_track
d_a = cumsum_by_track(a.distance_to_next(), a.tracks)[(i0 - 1 + nb)[nb != 0]] / 1000.0
i0, nb = c.index_from_track, c.nb_obs_by_track
d_c = cumsum_by_track(c.distance_to_next(), c.tracks)[(i0 - 1 + nb)[nb != 0]] / 1000.0
Setup axes
figure = plt.figure(figsize=(12, 8))
ax_ratio_cum = figure.add_axes([0.55, 0.06, 0.42, 0.34])
ax_ratio = figure.add_axes([0.07, 0.06, 0.46, 0.34])
ax_cum = figure.add_axes([0.55, 0.43, 0.42, 0.54])
ax = figure.add_axes([0.07, 0.43, 0.46, 0.54])
ax.set_ylabel("Eddies by year")
ax_ratio.set_ylabel("Ratio Cyclonic/Anticyclonic")
for ax_ in (ax, ax_cum, ax_ratio_cum, ax_ratio):
ax_.set_xlim(0, 1000)
if ax_ in (ax, ax_cum):
ax_.set_ylim(1e-1, 1e4), ax_.set_yscale("log")
else:
ax_.set_xlabel("Propagation in km (with bins of 20 km)")
ax_.set_ylim(0, 2)
ax_.axhline(1, color="g", lw=2)
ax_.grid()
ax_cum.xaxis.set_ticklabels([]), ax_cum.yaxis.set_ticklabels([])
ax.xaxis.set_ticklabels([]), ax_ratio_cum.yaxis.set_ticklabels([])
# plot data
bin_hist = arange(0, 2000, 20)
x = (bin_hist[1:] + bin_hist[:-1]) / 2.0
w_a, w_c = ones(d_a.shape) / nb_year, ones(d_c.shape) / nb_year
kwargs_a = dict(histtype="step", bins=bin_hist, x=d_a, color="r", weights=w_a)
kwargs_c = dict(histtype="step", bins=bin_hist, x=d_c, color="b", weights=w_c)
cum_a, _, _ = ax_cum.hist(cumulative=-1, **kwargs_a)
cum_c, _, _ = ax_cum.hist(cumulative=-1, **kwargs_c)
nb_a, _, _ = ax.hist(label="Anticyclonic", **kwargs_a)
nb_c, _, _ = ax.hist(label="Cyclonic", **kwargs_c)
ax_ratio_cum.plot(x, cum_c / cum_a)
ax_ratio.plot(x, nb_c / nb_a)
ax.legend()
Out:
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/checkouts/v3.3.0/examples/10_tracking_diagnostics/pet_propagation.py:66: RuntimeWarning: invalid value encountered in true_divide
ax_ratio_cum.plot(x, cum_c / cum_a)
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/checkouts/v3.3.0/examples/10_tracking_diagnostics/pet_propagation.py:67: RuntimeWarning: divide by zero encountered in true_divide
ax_ratio.plot(x, nb_c / nb_a)
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/checkouts/v3.3.0/examples/10_tracking_diagnostics/pet_propagation.py:67: RuntimeWarning: invalid value encountered in true_divide
ax_ratio.plot(x, nb_c / nb_a)
Total running time of the script: ( 0 minutes 4.234 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Count pixel used¶
Do Geo stat with frequency and compare with center count method: Count center
import py_eddy_tracker_sample
from matplotlib import pyplot as plt
from matplotlib.colors import LogNorm
from py_eddy_tracker.observations.tracking import TrackEddiesObservations
Load an experimental med atlas over a period of 26 years (1993-2019)
a = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Anticyclonic.zarr")
)
c = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Cyclonic.zarr")
)
Parameters
step = 0.125
bins = ((-10, 37, step), (30, 46, step))
kwargs_pcolormesh = dict(
cmap="terrain_r", vmin=0, vmax=0.75, factor=1 / a.nb_days, name="count"
)
Plot
fig = plt.figure(figsize=(12, 18.5))
ax_a = fig.add_axes([0.03, 0.75, 0.90, 0.25])
ax_a.set_title("Anticyclonic frequency")
ax_c = fig.add_axes([0.03, 0.5, 0.90, 0.25])
ax_c.set_title("Cyclonic frequency")
ax_all = fig.add_axes([0.03, 0.25, 0.90, 0.25])
ax_all.set_title("All eddies frequency")
ax_ratio = fig.add_axes([0.03, 0.0, 0.90, 0.25])
ax_ratio.set_title("Ratio cyclonic / Anticyclonic")
# Count pixel used for each contour
g_a = a.grid_count(bins, intern=True)
g_a.display(ax_a, **kwargs_pcolormesh)
g_c = c.grid_count(bins, intern=True)
g_c.display(ax_c, **kwargs_pcolormesh)
# Compute a ratio Cyclonic / Anticyclonic
ratio = g_c.vars["count"] / g_a.vars["count"]
# Mask manipulation to be able to sum the 2 grids
m_c = g_c.vars["count"].mask
m = m_c & g_a.vars["count"].mask
g_c.vars["count"][m_c] = 0
g_c.vars["count"] += g_a.vars["count"]
g_c.vars["count"].mask = m
m = g_c.display(ax_all, **kwargs_pcolormesh)
plt.colorbar(m, cax=fig.add_axes([0.95, 0.27, 0.01, 0.7]))
g_c.vars["count"] = ratio
m = g_c.display(
ax_ratio, name="count", vmin=0.1, vmax=10, norm=LogNorm(), cmap="coolwarm_r"
)
plt.colorbar(m, cax=fig.add_axes([0.95, 0.02, 0.01, 0.2]))
for ax in (ax_a, ax_c, ax_all, ax_ratio):
ax.set_aspect("equal")
ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
ax.grid()
Out:
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/pyEddyTracker-3.3.0-py3.7.egg/py_eddy_tracker/dataset/grid.py:1896: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it.
return ax.pcolormesh(x, self.y_bounds, data.T * factor, **kwargs)
Count Anticyclones as a function of lifetime¶
fig = plt.figure(figsize=(12, 10))
mask = a.lifetime >= 60
ax_long = fig.add_axes([0.03, 0.53, 0.90, 0.45])
g_a = a.grid_count(bins, intern=True, filter=mask)
g_a.display(ax_long, **kwargs_pcolormesh)
ax_long.set_title(f"Anticyclones with lifetime >= 60 days ({mask.sum()} Obs)")
ax_short = fig.add_axes([0.03, 0.03, 0.90, 0.45])
g_a = a.grid_count(bins, intern=True, filter=~mask)
m = g_a.display(ax_short, **kwargs_pcolormesh)
ax_short.set_title(f"Anticyclones with lifetime < 60 days ({(~mask).sum()} Obs)")
for ax in (ax_short, ax_long):
ax.set_aspect("equal"), ax.grid()
ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
cb = plt.colorbar(m, cax=fig.add_axes([0.94, 0.05, 0.015, 0.9]))
Total running time of the script: ( 0 minutes 27.104 seconds)
Note
Click here to download the full example code or to run this example in your browser via Binder
Count center¶
Do Geo stat with center and compare with frequency method show: Count pixel used
import py_eddy_tracker_sample
from matplotlib import pyplot as plt
from matplotlib.colors import LogNorm
from py_eddy_tracker.observations.tracking import TrackEddiesObservations
Load an experimental med atlas over a period of 26 years (1993-2019)
a = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Anticyclonic.zarr")
)
c = TrackEddiesObservations.load_file(
py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Cyclonic.zarr")
)
Parameters
step = 0.125
bins = ((-10, 37, step), (30, 46, step))
kwargs_pcolormesh = dict(
cmap="terrain_r", vmin=0, vmax=2, factor=1 / (a.nb_days * step ** 2), name="count"
)
Plot
fig = plt.figure(figsize=(12, 18.5))
ax_a = fig.add_axes([0.03, 0.75, 0.90, 0.25])
ax_a.set_title("Anticyclonic center frequency")
ax_c = fig.add_axes([0.03, 0.5, 0.90, 0.25])
ax_c.set_title("Cyclonic center frequency")
ax_all = fig.add_axes([0.03, 0.25, 0.90, 0.25])
ax_all.set_title("All eddies center frequency")
ax_ratio = fig.add_axes([0.03, 0.0, 0.90, 0.25])
ax_ratio.set_title("Ratio cyclonic / Anticyclonic")
# Count pixel used for each center
g_a = a.grid_count(bins, intern=True, center=True)
g_a.display(ax_a, **kwargs_pcolormesh)
g_c = c.grid_count(bins, intern=True, center=True)
g_c.display(ax_c, **kwargs_pcolormesh)
# Compute a ratio Cyclonic / Anticyclonic
ratio = g_c.vars["count"] / g_a.vars["count"]
# Mask manipulation to be able to sum the 2 grids
m_c = g_c.vars["count"].mask
m = m_c & g_a.vars["count"].mask
g_c.vars["count"][m_c] = 0
g_c.vars["count"] += g_a.vars["count"]
g_c.vars["count"].mask = m
m = g_c.display(ax_all, **kwargs_pcolormesh)
cb = plt.colorbar(m, cax=fig.add_axes([0.94, 0.27, 0.01, 0.7]))
cb.set_label("Eddies by 1°^2 by day")
g_c.vars["count"] = ratio
m = g_c.display(
ax_ratio, name="count", vmin=0.1, vmax=10, norm=LogNorm(), cmap="coolwarm_r"
)
plt.colorbar(m, cax=fig.add_axes([0.94, 0.02, 0.01, 0.2]))
for ax in (ax_a, ax_c, ax_all, ax_ratio):
ax.set_aspect("equal")
ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
ax.grid()
Out:
/home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/pyEddyTracker-3.3.0-py3.7.egg/py_eddy_tracker/dataset/grid.py:1896: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it.
return ax.pcolormesh(x, self.y_bounds, data.T * factor, **kwargs)
Count Anticyclones as a function of lifetime¶
Count at the center’s position
fig = plt.figure(figsize=(12, 10))
mask = a.lifetime >= 60
ax_long = fig.add_axes([0.03, 0.53, 0.90, 0.45])
g_a = a.grid_count(bins, center=True, filter=mask)
g_a.display(ax_long, **kwargs_pcolormesh)
ax_long.set_title(f"Anticyclones with lifetime >= 60 days ({mask.sum()} Obs)")
ax_short = fig.add_axes([0.03, 0.03, 0.90, 0.45])
g_a = a.grid_count(bins, center=True, filter=~mask)
m = g_a.display(ax_short, **kwargs_pcolormesh)
ax_short.set_title(f"Anticyclones with lifetime < 60 days ({(~mask).sum()} Obs)")
for ax in (ax_short, ax_long):
ax.set_aspect("equal"), ax.grid()
ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
cb = plt.colorbar(m, cax=fig.add_axes([0.94, 0.05, 0.015, 0.9]))
Total running time of the script: ( 0 minutes 3.789 seconds)
External data¶
Note
Click here to download the full example code or to run this example in your browser via Binder
Collocating external data¶
Script will use py-eddy-tracker methods to upload external data (sea surface temperature, SST) in a common structure with altimetry.
Figures higlights the different steps.
from datetime import datetime
from matplotlib import pyplot as plt
from py_eddy_tracker import data
from py_eddy_tracker.dataset.grid import RegularGridDataset
date = datetime(2016, 7, 7)
filename_alt = data.get_path(f"dt_blacksea_allsat_phy_l4_{date:%Y%m%d}_20200801.nc")
filename_sst = data.get_path(
f"{date:%Y%m%d}000000-GOS-L4_GHRSST-SSTfnd-OISST_HR_REP-BLK-v02.0-fv01.0.nc"
)
var_name_sst = "analysed_sst"
extent = [27, 42, 40.5, 47]
Loading data¶
sst = RegularGridDataset(filename=filename_sst, x_name="lon", y_name="lat")
alti = RegularGridDataset(
data.get_path(filename_alt), x_name="longitude", y_name="latitude"
)
# We can use `Grid` tools to interpolate ADT on the sst grid
sst.regrid(alti, "sla")
sst.add_uv("sla")
Out:
We assume pixel position of grid is center for /home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/pyEddyTracker-3.3.0-py3.7.egg/py_eddy_tracker/data/20160707000000-GOS-L4_GHRSST-SSTfnd-OISST_HR_REP-BLK-v02.0-fv01.0.nc
We assume pixel position of grid is center for /home/docs/checkouts/readthedocs.org/user_builds/py-eddy-tracker/envs/v3.3.0/lib/python3.7/site-packages/pyEddyTracker-3.3.0-py3.7.egg/py_eddy_tracker/data/dt_blacksea_allsat_phy_l4_20160707_20200801.nc
Functions to initiate figure axes
def start_axes(title, extent=extent):
fig = plt.figure(figsize=(13, 6), dpi=120)
ax = fig.add_axes([0.03, 0.05, 0.89, 0.91])
ax.set_xlim(extent[0], extent[1])
ax.set_ylim(extent[2], extent[3])
ax.set_title(title)
ax.set_aspect("equal")
return ax
def update_axes(ax, mappable=None, unit=""):
ax.grid()
if mappable:
cax = ax.figure.add_axes([0.93, 0.05, 0.01, 0.9], title=unit)
plt.colorbar(mappable, cax=cax)
ADT first display¶
ax = start_axes("SLA", extent=extent)
m = sst.display(ax, "sla", vmin=0.05, vmax=0.35)
update_axes(ax, m, unit="[m]")
![SLA, [m]](_images/sphx_glr_pet_SST_collocation_001.png)
SST first display¶
We can now plot SST from sst
ax = start_axes("SST")
m = sst.display(ax, "analysed_sst", vmin=295, vmax=300)
update_axes(ax, m, unit="[°K]")
![SST, [°K]](_images/sphx_glr_pet_SST_collocation_002.png)
ax = start_axes("SST")
m = sst.display(ax, "analysed_sst", vmin=295, vmax=300)
u, v = sst.grid("u").T, sst.grid("v").T
ax.quiver(sst.x_c[::3], sst.y_c[::3], u[::3, ::3], v[::3, ::3], scale=10)
update_axes(ax, m, unit="[°K]")
![SST, [°K]](_images/sphx_glr_pet_SST_collocation_003.png)
Now, with eddy contours, and displaying SST anomaly
sst.bessel_high_filter("analysed_sst", 400)
Eddy detection
sst.bessel_high_filter("sla", 400)
# ADT filtered
ax = start_axes("SLA", extent=extent)
m = sst.display(ax, "sla", vmin=-0.1, vmax=0.1)
update_axes(ax, m, unit="[m]")
a, c = sst.eddy_identification("sla", "u", "v", date, 0.002)
![SLA, [m]](_images/sphx_glr_pet_SST_collocation_004.png)
kwargs_a = dict(lw=2, label="Anticyclonic", ref=-10, color="b")
kwargs_c = dict(lw=2, label="Cyclonic", ref=-10, color="r")
ax = start_axes("SST anomaly")
m = sst.display(ax, "analysed_sst", vmin=-1, vmax=1)
a.display(ax, **kwargs_a), c.display(ax, **kwargs_c)
ax.legend()
update_axes(ax, m, unit="[°K]")
![SST anomaly, [°K]](_images/sphx_glr_pet_SST_collocation_005.png)
Example of post-processing¶
Get mean of sst anomaly_high in each internal contour
anom_a = a.interp_grid(sst, "analysed_sst", method="mean", intern=True)
anom_c = c.interp_grid(sst, "analysed_sst", method="mean", intern=True)
Are cyclonic (resp. anticyclonic) eddies generally associated with positive (resp. negative) SST anomaly ?
fig = plt.figure(figsize=(7, 5))
ax = fig.add_axes([0.05, 0.05, 0.90, 0.90])
ax.set_xlabel("SST anomaly")
ax.set_xlim([-1, 1])
ax.set_title("Histograms of SST anomalies")
ax.hist(
anom_a, 5, alpha=0.5, color="b", label="Anticyclonic (mean:%s)" % (anom_a.mean())
)
ax.hist(anom_c, 5, alpha=0.5, color="r", label="Cyclonic (mean:%s)" % (anom_c.mean()))
ax.legend()
Not clearly so in that case ..
Total running time of the script: ( 0 minutes 10.201 seconds)
Eddy identification¶
Run the identification process for a single day
Shell/bash command¶
Bash command will allow to process one grid, it will apply a filter and an identification.
EddyId share/nrt_global_allsat_phy_l4_20190223_20190226.nc 20190223 \
adt ugos vgos longitude latitude \
out_directory -v DEBUG
Filter could be modify with options –cut_wavelength and –filter_order. You could also defined height between two isolines with –isoline_step, which could improve speed profile quality and detect accurately tiny eddies. You could also use –fit_errmax to manage acceptable shape of eddies.
An eddy identification will produce two files in the output directory, one for anticyclonic eddies and the other one for cyclonic.
In regional area which are away from the equator, current could be deduce from height, juste write None None inplace of ugos vgos
In case of datacube, you need to specify index for each layer (time, depth, …) wiht –indexs option like:
EddyId share/nrt_global_allsat_phy_l4_20190223_20190226.nc 20190223 \
adt ugos vgos longitude latitude \
out_directory -v DEBUG --indexs time=0
Warning
If no index are specified, you will apply identification only on dataset first layer, which could be a problem for datacube. Date set in command is used only for output storage.
Python code¶
If we want customize eddies identification, python module is here.
Activate verbose
from py_eddy_tracker import start_logger
start_logger().setLevel('DEBUG') # Available options: ERROR, WARNING, INFO, DEBUG
Run identification
from datetime import datetime
h = RegularGridDataset(grid_name, lon_name, lat_name)
h.bessel_high_filter('adt', 500, order=3)
date = datetime(2019, 2, 23)
a, c = h.eddy_identification(
'adt', 'ugos', 'vgos', # Variables used for identification
date, # Date of identification
0.002, # step between two isolines of detection (m)
pixel_limit=(5, 2000), # Min and max pixel count for valid contour
shape_error=55, # Error max (%) between ratio of circle fit and contour
)
Plot the resulting identification
fig = plt.figure(figsize=(15,7))
ax = fig.add_axes([.03,.03,.94,.94])
ax.set_title('Eddies detected -- Cyclonic(red) and Anticyclonic(blue)')
ax.set_ylim(-75,75)
ax.set_xlim(0,360)
ax.set_aspect('equal')
a.display(ax, color='b', linewidth=.5)
c.display(ax, color='r', linewidth=.5)
ax.grid()
fig.savefig('share/png/eddies.png')
Save identification data
from netCDF import Dataset
with Dataset(date.strftime('share/Anticyclonic_%Y%m%d.nc'), 'w') as h:
a.to_netcdf(h)
with Dataset(date.strftime('share/Cyclonic_%Y%m%d.nc'), 'w') as h:
c.to_netcdf(h)
Load, Display and Filtering¶
Loading grid
from py_eddy_tracker.dataset.grid import RegularGridDataset
grid_name, lon_name, lat_name = 'share/nrt_global_allsat_phy_l4_20190223_20190226.nc', 'longitude', 'latitude'
h = RegularGridDataset(grid_name, lon_name, lat_name)
Plotting grid
from matplotlib import pyplot as plt
fig = plt.figure(figsize=(14, 12))
ax = fig.add_axes([.02, .51, .9, .45])
ax.set_title('ADT (m)')
ax.set_ylim(-75, 75)
ax.set_aspect('equal')
m = h.display(ax, name='adt', vmin=-1, vmax=1)
ax.grid(True)
plt.colorbar(m, cax=fig.add_axes([.94, .51, .01, .45]))
Filtering
h = RegularGridDataset(grid_name, lon_name, lat_name)
h.bessel_high_filter('adt', 500, order=3)
Save grid
h.write('/tmp/grid.nc')
Add second plot
ax = fig.add_axes([.02, .02, .9, .45])
ax.set_title('ADT Filtered (m)')
ax.set_aspect('equal')
ax.set_ylim(-75, 75)
m = h.display(ax, name='adt', vmin=-.1, vmax=.1)
ax.grid(True)
plt.colorbar(m, cax=fig.add_axes([.94, .02, .01, .45]))
fig.savefig('share/png/filter.png')
Spectrum¶
Compute spectrum and spectrum ratio on some area¶
Load data
raw = RegularGridDataset(grid_name, lon_name, lat_name)
filtered = RegularGridDataset(grid_name, lon_name, lat_name)
filtered.bessel_low_filter('adt', 150, order=3)
areas = dict(
sud_pacific=dict(llcrnrlon=188, urcrnrlon=280, llcrnrlat=-64, urcrnrlat=-7),
atlantic_nord=dict(llcrnrlon=290, urcrnrlon=340, llcrnrlat=19.5, urcrnrlat=43),
indien_sud=dict(llcrnrlon=35, urcrnrlon=110, llcrnrlat=-49, urcrnrlat=-26),
)
Compute and display spectrum
fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111)
ax.set_title('Spectrum')
ax.set_xlabel('km')
for name_area, area in areas.items():
lon_spec, lat_spec = raw.spectrum_lonlat('adt', area=area)
mappable = ax.loglog(*lat_spec, label='lat %s raw' % name_area)[0]
ax.loglog(*lon_spec, label='lon %s raw' % name_area, color=mappable.get_color(), linestyle='--')
lon_spec, lat_spec = filtered.spectrum_lonlat('adt', area=area)
mappable = ax.loglog(*lat_spec, label='lat %s high' % name_area)[0]
ax.loglog(*lon_spec, label='lon %s high' % name_area, color=mappable.get_color(), linestyle='--')
ax.set_xscale('log')
ax.legend()
ax.grid()
fig.savefig('share/png/spectrum.png')
Compute and display spectrum ratio
fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111)
ax.set_title('Spectrum ratio')
ax.set_xlabel('km')
for name_area, area in areas.items():
lon_spec, lat_spec = filtered.spectrum_lonlat('adt', area=area, ref=raw)
mappable = ax.plot(*lat_spec, label='lat %s high' % name_area)[0]
ax.plot(*lon_spec, label='lon %s high' % name_area, color=mappable.get_color(), linestyle='--')
ax.set_xscale('log')
ax.legend()
ax.grid()
fig.savefig('share/png/spectrum_ratio.png')
Tracking¶
Requirements¶
Before to run tracking, you will need to run identification on every time step of the period(Period of your study).
Advice : Before to run tracking, display some identification file allow to learn a lot
Default method¶
To run a tracking just create an yaml file with minimal specification (FILES_PATTERN and SAVE_DIR). You will run tracking separately between Cyclonic eddies and Anticyclonic eddies.
Example of conf.yaml
PATHS:
# Files produces with EddyIdentification
FILES_PATTERN: MY_IDENTIFICATION_PATH/Anticyclonic*.nc
SAVE_DIR: MY_OUTPUT_PATH
# Number of timestep for missing detection
VIRTUAL_LENGTH_MAX: 3
# Minimal time to consider as a full track
TRACK_DURATION_MIN: 10
To run:
EddyTracking conf.yaml -v DEBUG
It will use default tracker:
No travel longer than 125 km between two observation
Amplitude and speed radius must be close to previous observation
In case of several candidate only closest is kept
It will produce 4 files by run:
A file of correspondances which will contains all the information to merge all identifications file
A file which will contains all the observations which are alone
A file which will contains all the short track which are shorter than TRACK_DURATION_MIN
A file which will contains all the long track which are longer than TRACK_DURATION_MIN
Use python module¶
An example of tracking with python module is available in the gallery: Track in python
Choose a tracker¶
With yaml you could also select another tracker:
PATHS:
# Files produces with EddyIdentification
FILES_PATTERN: MY/IDENTIFICATION_PATH/Anticyclonic*.nc
SAVE_DIR: MY_OUTPUT_PATH
# Number of timestep for missing detection
VIRTUAL_LENGTH_MAX: 3
# Minimal time to consider as a full track
TRACK_DURATION_MIN: 10
CLASS:
# Give the module to import,
# must be available when you do "import module" in python
MODULE: py_eddy_tracker.featured_tracking.old_tracker_reference
# Give class name which must be inherit from
# py_eddy_tracker.observations.observation.EddiesObservations
CLASS: CheltonTracker
This tracker is like described in CHELTON11[https://doi.org/10.1016/j.pocean.2011.01.002].
Code is here py_eddy_tracker.featured_tracking.old_tracker_reference()
Customize tracking¶
Code my own tracking¶
To use your own tracking method, you just need to create a class which inherit
from py_eddy_tracker.observations.observation.EddiesObservations() and set this class in yaml file like
we see in the previous topic.
API reference¶
Entry point |
|
Class to load and manipulate RegularGrid and UnRegularGrid |
|
Class to create network of observations |
|
Base class to manage eddy observation |
|
Class to manage observations gathered in track |
|
Class to compute Amplitude and average speed profile |
|
Tool method which use mostly numba |
|
GUI class |
|
Method for polygon |
|
Class to store link between observations |
Changelog¶
All notable changes to this project will be documented in this file.
The format is based on Keep a Changelog and this project adheres to Semantic Versioning.
[Unreleased]¶
[3.3.0] - 2020-12-03¶
Added¶
Add an implementation of visvalingam algorithm to simplify polygons with low modification
Add method to found close tracks in an other atlas
Allow to give a x reference when we display grid to be able to change xlim
Add option to EddyId to select data index like –indexs time=5 depth=2
Add a method to merge several indexs type for eddy obs
Get dataset variable like attribute, and lifetime/age are available for all observations
Add EddyInfos application to get general information about eddies dataset
Add method to inspect contour rejection (which are not in eddies)
Grid interp could be “nearest” or “bilinear”
Changed¶
Now to have object informations in plot label used python
`format`style, several key are available :“t0”
“t1”
“nb_obs”
“nb_tracks” (only for tracked eddies)