- Add example to illustrate grouping step:

	- add feature to group with in-memory object
Simplification of transfer method when 2 network are joined

Author: AntSimi

examples/16_network/pet_group_anim.py

@@ -0,0 +1,150 @@
+"""
+Network group process
+=====================
+"""
+
+import re
+from datetime import datetime
+
+from matplotlib import pyplot as plt
+from matplotlib.animation import FuncAnimation
+from matplotlib.colors import ListedColormap
+from numba import njit
+from numpy import arange, array, empty, ones
+
+from py_eddy_tracker import data
+from py_eddy_tracker.generic import flatten_line_matrix
+from py_eddy_tracker.observations.network import Network
+from py_eddy_tracker.observations.observation import EddiesObservations
+
+
+# %%
+class VideoAnimation(FuncAnimation):
+    def _repr_html_(self, *args, **kwargs):
+        """To get video in html and have a player"""
+        content = self.to_html5_video()
+        return re.sub(
+            'width="[0-9]*"\sheight="[0-9]*"', 'width="100%" height="100%"', content
+        )
+
+    def save(self, *args, **kwargs):
+        if args[0].endswith("gif"):
+            # In this case gif is use to create thumbnail which are not use but consume same time than video
+            # So we create an empty file, to save time
+            with open(args[0], "w") as _:
+                pass
+            return
+        return super().save(*args, **kwargs)
+
+
+# %%
+NETWORK_GROUPS = list()
+
+
+@njit(cache=True)
+def apply_replace(x, x0, x1):
+    nb = x.shape[0]
+    for i in range(nb):
+        if x[i] == x0:
+            x[i] = x1
+
+
+# %%
+# Modified class to catch group process at each step in order to illustrate processing
+class MyNetwork(Network):
+    def get_group_array(self, results, nb_obs):
+        """With a loop on all pair of index, we will label each obs with a group
+        number
+        """
+        nb_obs = array(nb_obs, dtype="u4")
+        day_start = nb_obs.cumsum() - nb_obs
+        gr = empty(nb_obs.sum(), dtype="u4")
+        gr[:] = self.NOGROUP
+
+        id_free = 1
+        for i, j, ii, ij in results:
+            gr_i = gr[slice(day_start[i], day_start[i] + nb_obs[i])]
+            gr_j = gr[slice(day_start[j], day_start[j] + nb_obs[j])]
+            # obs with no groups
+            m = (gr_i[ii] == self.NOGROUP) * (gr_j[ij] == self.NOGROUP)
+            nb_new = m.sum()
+            gr_i[ii[m]] = gr_j[ij[m]] = arange(id_free, id_free + nb_new)
+            id_free += nb_new
+            # associate obs with no group with obs with group
+            m = (gr_i[ii] != self.NOGROUP) * (gr_j[ij] == self.NOGROUP)
+            gr_j[ij[m]] = gr_i[ii[m]]
+            m = (gr_i[ii] == self.NOGROUP) * (gr_j[ij] != self.NOGROUP)
+            gr_i[ii[m]] = gr_j[ij[m]]
+            # case where 2 obs have a different group
+            m = gr_i[ii] != gr_j[ij]
+            if m.any():
+                # Merge of group, ref over etu
+                for i_, j_ in zip(ii[m], ij[m]):
+                    g0, g1 = gr_i[i_], gr_j[j_]
+                    apply_replace(gr, g0, g1)
+            NETWORK_GROUPS.append((i, j, gr.copy()))
+        return gr
+
+
+# %%
+# Movie period
+t0 = (datetime(2005, 5, 1) - datetime(1950, 1, 1)).days
+t1 = (datetime(2005, 6, 1) - datetime(1950, 1, 1)).days
+
+# %%
+# Get data from period and area
+e = EddiesObservations.load_file(data.get_path("Anticyclonic_seg.nc"))
+e = e.extract_with_mask((e.time >= t0) * (e.time < t1)).extract_with_area(
+    dict(llcrnrlon=25, urcrnrlon=35, llcrnrlat=31, urcrnrlat=37.5)
+)
+# %%
+# Reproduce individual daily identification(for demonstration)
+EDDIES_BY_DAYS = list()
+for i, b0, b1 in e.iter_on("time"):
+    EDDIES_BY_DAYS.append(e.index(i))
+# need for display
+e = EddiesObservations.concatenate(EDDIES_BY_DAYS)
+
+# %%
+# Run network building group to intercept every step
+n = MyNetwork.from_eddiesobservations(EDDIES_BY_DAYS, window=7)
+_ = n.group_observations(minimal_area=True)
+
+
+# %%
+def update(frame):
+    i_current, i_match, gr = NETWORK_GROUPS[frame]
+    current = EDDIES_BY_DAYS[i_current]
+    x = flatten_line_matrix(current.contour_lon_e)
+    y = flatten_line_matrix(current.contour_lat_e)
+    current_contour.set_data(x, y)
+    match = EDDIES_BY_DAYS[i_match]
+    x = flatten_line_matrix(match.contour_lon_e)
+    y = flatten_line_matrix(match.contour_lat_e)
+    matched_contour.set_data(x, y)
+    groups.set_array(gr)
+    txt.set_text(f"Day {i_current} match with day {i_match}")
+    s = 80 * ones(gr.shape)
+    s[gr == 0] = 4
+    groups.set_sizes(s)
+
+
+# %%
+# Anim
+# ----
+fig = plt.figure(figsize=(16, 9), dpi=50)
+ax = fig.add_axes([0, 0, 1, 1])
+ax.set_aspect("equal"), ax.grid(), ax.set_xlim(26, 34), ax.set_ylim(32, 36.5)
+groups = ax.scatter(
+    e.lon,
+    e.lat,
+    c=NETWORK_GROUPS[0][2],
+    cmap=ListedColormap(["gray", *e.COLORS[:-1]], name="from_list", N=30),
+    vmin=0,
+    vmax=30,
+)
+current_contour = ax.plot([], [], "k", lw=2, label="Current contour")[0]
+matched_contour = ax.plot([], [], "r", lw=1, ls="--", label="Candidate contour")[0]
+txt = ax.text(31.5, 35.5, "", fontsize=25)
+ax.legend(fontsize=25)
+ani = VideoAnimation(fig, update, frames=len(NETWORK_GROUPS), interval=220)

src/py_eddy_tracker/observations/network.py

@@ -43,6 +43,8 @@ def __init__(self, buffersize, intern=False, memory=False):
         self.memory = memory
 
     def load_contour(self, filename):
+        if isinstance(filename, EddiesObservations):
+            return filename[self.xname], filename[self.yname]
         if filename not in self.DATA:
             if len(self.FLIST) > self.buffersize:
                 self.DATA.pop(self.FLIST.pop(0))
@@ -815,10 +817,18 @@ def __init__(self, input_regex, window=5, intern=False, memory=False):
         """
         self.window = window
         self.buffer = Buffer(window, intern, memory)
+        self.memory = memory
+
         self.filenames = glob(input_regex)
         self.filenames.sort()
         self.nb_input = len(self.filenames)
-        self.memory = memory
+
+    @classmethod
+    def from_eddiesobservations(cls, observations, *args, **kwargs):
+        new = cls("", *args, **kwargs)
+        new.filenames = observations
+        new.nb_input = len(new.filenames)
+        return new
 
     def get_group_array(self, results, nb_obs):
         """With a loop on all pair of index, we will label each obs with a group
@@ -849,10 +859,13 @@ def get_group_array(self, results, nb_obs):
             if m.any():
                 # Merge of group, ref over etu
                 for i_, j_ in zip(ii[m], ij[m]):
-                    merge_id.append((gr_i[i_], gr_j[j_]))
+                    g0, g1 = gr_i[i_], gr_j[j_]
+                    if g0 > g1:
+                        g0, g1 = g1, g0
+                    merge_id.append((g0, g1))
 
         gr_transfer = arange(id_free, dtype="u4")
-        for i, j in merge_id:
+        for i, j in set(merge_id):
             gr_i, gr_j = gr_transfer[i], gr_transfer[j]
             if gr_i != gr_j:
                 apply_replace(gr_transfer, gr_i, gr_j)

src/py_eddy_tracker/observations/observation.py

@@ -490,6 +490,8 @@ def iter_on(self, xname, bins=None):
         d = x[1:] - x[:-1]
         if bins is None:
             bins = arange(x.min(), x.max() + 2)
+        elif not isinstance(bins, ndarray):
+            bins = array(bins)
         nb_bins = len(bins) - 1
         i = numba_digitize(x, bins) - 1
         # Not monotonous
@@ -1609,7 +1611,9 @@ def extract_with_area(self, area, **kwargs):
 
         .. minigallery:: py_eddy_tracker.EddiesObservations.extract_with_area
         """
-        mask = (self.latitude > area["llcrnrlat"]) * (self.latitude < area["urcrnrlat"])
+        lat0 = area.get("llcrnrlat", -90)
+        lat1 = area.get("urcrnrlat", 90)
+        mask = (self.latitude > lat0) * (self.latitude < lat1)
         lon0 = area["llcrnrlon"]
         lon = (self.longitude - lon0) % 360 + lon0
         mask *= (lon > lon0) * (lon < area["urcrnrlon"])