Particles drift

README.md

@@ -17,15 +17,17 @@ Method was described in :
 ### Use case ###
 
 Method is used in :
- 
+
 [Mason, E., A. Pascual, P. Gaube, S.Ruiz, J. Pelegrí, A. Delepoulle, 2017: Subregional characterization of mesoscale eddies across the Brazil-Malvinas Confluence](https://doi.org/10.1002/2016JC012611)
 
 ### How do I get set up? ###
 
 #### Short story ####
+
 ```bash
 pip install pyeddytracker
 ```
+
 #### Long story ####
 
 To avoid problems with installation, use of the virtualenv Python virtual environment is recommended.
@@ -36,12 +38,20 @@ Then use pip to install all dependencies (numpy, scipy, matplotlib, netCDF4, ...
 pip install numpy scipy netCDF4 matplotlib opencv-python pyyaml pint polygon3
 ```
 
-Then run the following to install the eddy tracker:
+Clone :
+
+```bash
+git clone https://github.com/AntSimi/py-eddy-tracker
+```
+
+Then run the following to install the eddy tracker :
 
 ```bash
 python setup.py install
 ```
+
 ### Tools gallery ###
+
 Several examples based on py eddy tracker module are [here](https://py-eddy-tracker.readthedocs.io/en/latest/python_module/index.html).
 
 [![](https://py-eddy-tracker.readthedocs.io/en/latest/_static/logo.png)](https://py-eddy-tracker.readthedocs.io/en/latest/python_module/index.html)

examples/02_eddy_identification/pet_statistics_on_identification.py

@@ -0,0 +1,105 @@
+"""
+Stastics on identification files
+================================
+
+Some statistics on raw identification without any tracking
+"""
+import numpy as np
+from matplotlib import pyplot as plt
+from matplotlib.dates import date2num
+
+from py_eddy_tracker import start_logger
+from py_eddy_tracker.data import get_remote_demo_sample
+from py_eddy_tracker.observations.observation import EddiesObservations
+
+start_logger().setLevel("ERROR")
+
+
+# %%
+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]))
+
+
+# %%
+# We load demo sample and take only first year.
+#
+# Replace by a list of filename to apply on your own dataset.
+file_objects = get_remote_demo_sample(
+    "eddies_med_adt_allsat_dt2018/Anticyclonic_2010_2011_2012"
+)[:365]
+
+# %%
+# Merge all identification datasets in one object
+all_a = EddiesObservations.concatenate(
+    [EddiesObservations.load_file(i) for i in file_objects]
+)
+
+# %%
+# We define polygon bound
+x0, x1, y0, y1 = 15, 20, 33, 38
+xs = np.array([[x0, x1, x1, x0, x0]], dtype="f8")
+ys = np.array([[y0, y0, y1, y1, y0]], dtype="f8")
+# Polygon object created for the match function use.
+polygon = dict(contour_lon_e=xs, contour_lat_e=ys, contour_lon_s=xs, contour_lat_s=ys)
+
+# %%
+# Geographic frequency of eddies
+step = 0.125
+ax = start_axes("")
+# Count pixel encompassed in each contour
+g_a = all_a.grid_count(bins=((-10, 37, step), (30, 46, step)), intern=True)
+m = g_a.display(
+    ax, cmap="terrain_r", vmin=0, vmax=0.75, factor=1 / all_a.nb_days, name="count"
+)
+ax.plot(polygon["contour_lon_e"][0], polygon["contour_lat_e"][0], "r")
+update_axes(ax, m)
+
+# %%
+# We use the match function to count the number of eddies that intersect the polygon defined previously
+# `p1_area` option allow to get in c_e/c_s output, precentage of area occupy by eddies in the polygon.
+i_e, j_e, c_e = all_a.match(polygon, p1_area=True, intern=False)
+i_s, j_s, c_s = all_a.match(polygon, p1_area=True, intern=True)
+
+# %%
+dt = np.datetime64("1970-01-01") - np.datetime64("1950-01-01")
+kw_hist = dict(
+    bins=date2num(np.arange(21900, 22300).astype("datetime64") - dt), histtype="step"
+)
+# translate julian day in datetime64
+t = all_a.time.astype("datetime64") - dt
+# %%
+# Number of eddies within a polygon
+ax = plt.figure(figsize=(12, 6)).add_subplot(111)
+ax.set_title("Different ways to count eddies within a polygon")
+ax.set_ylabel("Count")
+m = all_a.mask_from_polygons(((xs, ys),))
+ax.hist(t[m], label="Eddy Center in polygon", **kw_hist)
+ax.hist(t[i_s[c_s > 0]], label="Intersection Speed contour and polygon", **kw_hist)
+ax.hist(t[i_e[c_e > 0]], label="Intersection Effective contour and polygon", **kw_hist)
+ax.legend()
+ax.set_xlim(np.datetime64("2010"), np.datetime64("2011"))
+ax.grid()
+
+# %%
+# Percent of the area of interest occupied by eddies.
+ax = plt.figure(figsize=(12, 6)).add_subplot(111)
+ax.set_title("Percent of polygon occupied by an anticyclonic eddy")
+ax.set_ylabel("Percent of polygon")
+ax.hist(t[i_s[c_s > 0]], weights=c_s[c_s > 0] * 100.0, label="speed contour", **kw_hist)
+ax.hist(
+    t[i_e[c_e > 0]], weights=c_e[c_e > 0] * 100.0, label="effective contour", **kw_hist
+)
+ax.legend(), ax.set_ylim(0, 25)
+ax.set_xlim(np.datetime64("2010"), np.datetime64("2011"))
+ax.grid()

examples/07_cube_manipulation/README.rst

@@ -1,6 +1,2 @@
 Time grid computation
 =====================
-
-.. warning::
-
-    Time grid is under development, API could move quickly!

examples/07_cube_manipulation/pet_particles_drift.py

@@ -0,0 +1,46 @@
+"""
+Build path of particle drifting
+===============================
+
+"""
+
+from matplotlib import pyplot as plt
+from numpy import arange, meshgrid
+
+from py_eddy_tracker import start_logger
+from py_eddy_tracker.data import get_demo_path
+from py_eddy_tracker.dataset.grid import GridCollection
+
+start_logger().setLevel("ERROR")
+
+# %%
+# Load data cube
+c = GridCollection.from_netcdf_cube(
+    get_demo_path("dt_med_allsat_phy_l4_2005T2.nc"),
+    "longitude",
+    "latitude",
+    "time",
+    heigth="adt",
+)
+
+# %%
+# Advection properties
+nb_days, step_by_day = 10, 6
+nb_time = step_by_day * nb_days
+kw_p = dict(nb_step=1, time_step=86400 / step_by_day)
+t0 = 20210
+
+# %%
+# Get paths
+x0, y0 = meshgrid(arange(32, 35, 0.5), arange(32.5, 34.5, 0.5))
+x0, y0 = x0.reshape(-1), y0.reshape(-1)
+t, x, y = c.path(x0, y0, "u", "v", t_init=t0, **kw_p, nb_time=nb_time)
+
+# %%
+# Plot paths
+ax = plt.figure(figsize=(9, 6)).add_subplot(111, aspect="equal")
+ax.plot(x0, y0, "k.", ms=20)
+ax.plot(x, y, lw=3)
+ax.set_title("10 days particle paths")
+ax.set_xlim(31, 35), ax.set_ylim(32, 34.5)
+ax.grid()