add examples

AntSimi · AntSimi · commit a237a29f0474 · 2020-06-15T15:59:53.000+02:00
diff --git a/examples/00_tracking_manipulation/README.rst b/examples/00_tracking_manipulation/README.rst
@@ -0,0 +1,2 @@
+Tracking Manipulation
+=======================
diff --git a/examples/00_tracking_manipulation/pet_display_field.py b/examples/00_tracking_manipulation/pet_display_field.py
@@ -0,0 +1,29 @@
+"""
+Display fields
+==============
+
+"""
+
+from matplotlib import pyplot as plt
+from py_eddy_tracker.observations.tracking import TrackEddiesObservations
+import py_eddy_tracker_sample
+
+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
+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(-5, 37)
+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)")
diff --git a/examples/00_tracking_manipulation/pet_display_track.py b/examples/00_tracking_manipulation/pet_display_track.py
@@ -24,11 +24,12 @@
 
 
 # Plot
-fig = plt.figure(figsize=(15, 8))
-ax = fig.add_subplot(111)
+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(-5, 37)
 ax.set_ylim(30, 46)
 a.plot(ax, ref=-10, label="Anticyclonic", color="r", lw=0.1)
 c.plot(ax, ref=-10, label="Cyclonic", color="b", lw=0.1)
 ax.legend()
+ax.grid()
diff --git a/examples/00_tracking_manipulation/pet_one_track.py b/examples/00_tracking_manipulation/pet_one_track.py
@@ -15,10 +15,11 @@
 eddy_f = a.extract_ids([9672])
 eddy_f.median_filter(1, "time", "lon").loess_filter(5, "time", "lon")
 eddy_f.median_filter(1, "time", "lat").loess_filter(5, "time", "lat")
-fig = plt.figure(figsize=(20, 8))
-ax = fig.add_subplot(111)
-ax.set_xlim(17, 22.5)
+fig = plt.figure(figsize=(12, 5))
+ax = fig.add_axes((0.05, 0.1, 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.plot(ax, color="g", lw=1)
diff --git a/examples/00_tracking_manipulation/pet_select_track_across_area.py b/examples/00_tracking_manipulation/pet_select_track_across_area.py
@@ -0,0 +1,33 @@
+"""
+Tracks which go through area
+============================
+
+"""
+from matplotlib import pyplot as plt
+from py_eddy_tracker.observations.tracking import TrackEddiesObservations
+import py_eddy_tracker_sample
+
+c = TrackEddiesObservations.load_file(
+    py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Cyclonic.zarr")
+)
+c.median_filter(1, "time", "lon").loess_filter(5, "time", "lon")
+c.median_filter(1, "time", "lat").loess_filter(5, "time", "lat")
+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)
+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")
+c_subset.plot(ax, color="red", lw=0.2, ref=-10, label="selected tracks")
+ax.plot(
+    (x0, x0, x1, x1, x0,),
+    (y0, y1, y1, y0, y0,),
+    color="green",
+    lw=1.5,
+    label="Box of selection",
+)
+ax.legend()
diff --git a/examples/10_tracking_diagnostics/README.rst b/examples/10_tracking_diagnostics/README.rst
diff --git a/examples/10_tracking_diagnostics/pet_center_count.py b/examples/10_tracking_diagnostics/pet_center_count.py
@@ -21,11 +21,11 @@
 ax = fig.add_subplot(111)
 ax.set_xlim(-5, 37)
 ax.set_ylim(30, 46)
+ax.set_aspect("equal")
 step = 0.1
 t0, t1 = a.period
-g = a.grid_count(((-5, 37, step), (30, 46, step)), center=True)
-g.vars["count"] = g.vars["count"] / (step ** 2 * (t1 - t0))
-m = g.display(ax, name="count", vmin=0, vmax=2)
+g = a.grid_count(((-6, 37, step), (30, 46, step)), center=True)
+m = g.display(ax, name="count", vmin=0, vmax=2, factor=1 / (step ** 2 * (t1 - t0)))
 ax.grid()
 cb = plt.colorbar(m, cax=fig.add_axes([0.95, 0.05, 0.01, 0.9]))
 cb.set_label("Eddies by 1°^2 by day")
diff --git a/examples/10_tracking_diagnostics/pet_geographic_stats.py b/examples/10_tracking_diagnostics/pet_geographic_stats.py
@@ -0,0 +1,58 @@
+"""
+Geographical statistics
+=======================
+
+"""
+
+from matplotlib import pyplot as plt
+from py_eddy_tracker.observations.tracking import TrackEddiesObservations
+import py_eddy_tracker_sample
+
+
+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)
+
+fig = plt.figure(figsize=(15, 7))
+ax = fig.add_subplot(111)
+ax.set_xlim(-5, 37)
+ax.set_ylim(30, 46)
+ax.set_aspect("equal")
+step = 0.1
+g = a.grid_stat(((-5, 37, step), (30, 46, step)), "amplitude")
+m = g.display(ax, name="amplitude", vmin=0, vmax=0.2)
+ax.grid()
+cb = plt.colorbar(m, cax=fig.add_axes([0.92, 0.05, 0.01, 0.9]))
+cb.set_label("Amplitude (m)")
+ax.set_title("Amplitude mean by box of %s°" % step)
+
+fig = plt.figure(figsize=(15, 7))
+ax = fig.add_subplot(111)
+ax.set_xlim(-5, 37)
+ax.set_ylim(30, 46)
+ax.set_aspect("equal")
+step = 0.1
+g = a.grid_stat(((-5, 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=fig.add_axes([0.92, 0.05, 0.01, 0.9]))
+cb.set_label("Speed radius (km)")
+ax.set_title("Speed radius mean by box of %s°" % step)
+
+fig = plt.figure(figsize=(15, 7))
+ax = fig.add_subplot(111)
+ax.set_xlim(-5, 37)
+ax.set_ylim(30, 46)
+ax.set_aspect("equal")
+step = 0.1
+g = a.grid_stat(((-5, 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=fig.add_axes([0.92, 0.05, 0.01, 0.9]))
+cb.set_label("Percent of virtual (%)")
+ax.set_title("Percent of virtual by box of %s°" % step)
diff --git a/examples/10_tracking_diagnostics/pet_histo.py b/examples/10_tracking_diagnostics/pet_histo.py
diff --git a/examples/10_tracking_diagnostics/pet_lifetime.py b/examples/10_tracking_diagnostics/pet_lifetime.py
diff --git a/examples/10_tracking_diagnostics/pet_pixel_used.py b/examples/10_tracking_diagnostics/pet_pixel_used.py
@@ -15,6 +15,7 @@
 c = TrackEddiesObservations.load_file(
     py_eddy_tracker_sample.get_path("eddies_med_adt_allsat_dt2018/Cyclonic.zarr")
 )
+t0, t1 = a.period
 
 # Plot
 fig = plt.figure(figsize=(15, 20))
@@ -26,24 +27,21 @@
 ax_all.set_title("All eddies frequency")
 
 step = 0.1
+bins = ((-10, 37, step), (30, 46, step))
 kwargs_pcolormesh = dict(cmap="terrain_r", vmin=0, vmax=0.75)
-g_a = a.grid_count(((-5, 37, step), (30, 46, step)), intern=True)
-t0, t1 = a.period
-g_a.vars["count"] = g_a.vars["count"] / (t1 - t0)
-m = g_a.display(ax_a, name="count", **kwargs_pcolormesh)
+g_a = a.grid_count(bins, intern=True)
+m = g_a.display(ax_a, name="count", factor=1 / (t1 - t0), **kwargs_pcolormesh)
 
-g_c = c.grid_count(((-5, 37, step), (30, 46, step)), intern=True)
-t0, t1 = c.period
-g_c.vars["count"] = g_c.vars["count"] / (t1 - t0)
-m = g_c.display(ax_c, name="count", **kwargs_pcolormesh)
+g_c = c.grid_count(bins, intern=True)
+m = g_c.display(ax_c, name="count", factor=1 / (t1 - t0), **kwargs_pcolormesh)
 
 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, name="count", **kwargs_pcolormesh)
+m = g_c.display(ax_all, name="count", factor=1 / (t1 - t0), **kwargs_pcolormesh)
 
 for ax in (ax_a, ax_c, ax_all):
     ax.set_aspect("equal")
diff --git a/examples/10_tracking_diagnostics/pet_propagation.py b/examples/10_tracking_diagnostics/pet_propagation.py
diff --git a/examples/README.rst b/examples/README.rst
@@ -1,2 +1,7 @@
 Py eddy tracker toolbox
 =======================
+
+All figures in this gallery, used an experimental dataset, compute with this dataset : cmems_product_.
+
+
+.. _cmems_product: https://resources.marine.copernicus.eu/?option=com_csw&task=results?option=com_csw&view=details&product_id=SEALEVEL_MED_PHY_L4_REP_OBSERVATIONS_008_051
diff --git a/src/py_eddy_tracker/dataset/grid.py b/src/py_eddy_tracker/dataset/grid.py
@@ -1767,10 +1767,12 @@ def init_speed_coef(self, uname="u", vname="v"):
         """
         self._speed_ev = (self.grid(uname) ** 2 + self.grid(vname) ** 2) ** 0.5
 
-    def display(self, ax, name, **kwargs):
+    def display(self, ax, name, factor=1, **kwargs):
         if "cmap" not in kwargs:
             kwargs["cmap"] = "coolwarm"
-        return ax.pcolormesh(self.x_bounds, self.y_bounds, self.grid(name).T, **kwargs)
+        return ax.pcolormesh(
+            self.x_bounds, self.y_bounds, self.grid(name).T * factor, **kwargs
+        )
 
     def interp(self, grid_name, lons, lats):
         """
diff --git a/src/py_eddy_tracker/observations/observation.py b/src/py_eddy_tracker/observations/observation.py
@@ -1241,6 +1241,12 @@ def set_global_attr_netcdf(self, h_nc):
         for key, item in self.global_attr.items():
             h_nc.setncattr(key, item)
 
+    def scatter(self, ax, name, ref=None, ** kwargs):
+        x = self.longitude
+        if ref is not None:
+            x = (x - ref) % 360 + ref
+        return ax.scatter(x, self.latitude, c=self[name], **kwargs)
+
     def display(self, ax, ref=None, extern_only=False, **kwargs):
         if not extern_only:
             lon_s = flatten_line_matrix(self.obs["contour_lon_s"])
@@ -1264,8 +1270,8 @@ def display(self, ax, ref=None, extern_only=False, **kwargs):
 
     def grid_count(self, bins, intern=False, center=False):
         x_name, y_name = self.intern(intern)
-        x_bins = arange(*bins[0])
-        y_bins = arange(*bins[1])
+        x_bins, y_bins = arange(*bins[0]), arange(*bins[1])
+        x0 = bins[0][0]
         grid = ma.zeros((x_bins.shape[0] - 1, y_bins.shape[0] - 1), dtype="u4")
         from ..dataset.grid import RegularGridDataset
 
@@ -1274,17 +1280,35 @@ def grid_count(self, bins, intern=False, center=False):
             datas=dict(count=grid, lon=x_bins[:-1], lat=y_bins[:-1],),
         )
         if center:
-            x, y = self.longitude, self.latitude
+            x, y = (self.longitude - x0) % 360 + x0, self.latitude
             grid[:] = histogram2d(x, y, (x_bins, y_bins))[0]
             grid.mask = grid.data == 0
         else:
-            x, y = self[x_name], self[y_name]
+            x, y = (self[x_name] - x0) % 360 + x0, self[y_name]
             for x_, y_ in zip(x, y):
                 i, j = BasePath(custom_concat(x_, y_)).pixels_in(regular_grid)
                 grid_count_(grid, i, j)
             grid.mask = grid == 0
         return regular_grid
 
+    def grid_stat(self, bins, varname):
+        x_bins, y_bins = arange(*bins[0]), arange(*bins[1])
+        x0 = bins[0][0]
+        x, y = (self.longitude - x0) % 360 + x0, self.latitude
+        sum_obs = histogram2d(x, y, (x_bins, y_bins), weights=self[varname])[0]
+        nb_obs = histogram2d(x, y, (x_bins, y_bins))[0]
+        from ..dataset.grid import RegularGridDataset
+
+        regular_grid = RegularGridDataset.with_array(
+            coordinates=("x", "y"),
+            datas={
+                varname: ma.array(sum_obs / nb_obs, mask=nb_obs == 0),
+                "x": x_bins[:-1],
+                "y": y_bins[:-1],
+            },
+        )
+        return regular_grid
+
 
 @njit(cache=True)
 def grid_count_(grid, i, j):

Original file line number	Diff line number	Diff line change
`@@ -0,0 +1,2 @@`
	`1`	`+Tracking Manipulation`
	`2`	`+=======================`