add merging and splitting triplets (AntSimi#51)

add triplet in merging and splitting

Author: CoriPegliasco

examples/16_network/pet_atlas.py

@@ -17,7 +17,7 @@
 step = 1 / 10.0
 bins = ((-10, 37, step), (30, 46, step))
 kw_time = dict(cmap="terrain_r", factor=100.0 / n.nb_days, name="count")
-kw_ratio = dict(cmap=plt.get_cmap("coolwarm_r", 10))
+kw_ratio = dict(cmap=plt.get_cmap("magma_r", 10))
 
 
 # %%
@@ -58,8 +58,10 @@ def update_axes(ax, mappable=None):
 m = g_10.display(ax, **kw_time, vmin=0, vmax=75)
 update_axes(ax, m).set_label("Pixel used in % of time")
 
+# %%
 # Display the ratio between the short and total presence.
-# Red = mostly short networks
+#
+# Light = mostly short networks
 ax = start_axes("")
 m = g_10.display(
     ax,
@@ -73,7 +75,7 @@ def update_axes(ax, mappable=None):
 # Network longer than 20 days
 # ---------------------------
 # Display the % of time each pixel (1/10°) is within an anticyclonic network
-# which total lifetime in longer than 20 days
+# which total lifetime is longer than 20 days
 ax = start_axes("")
 n20 = n.longer_than(20)
 g_20 = n20.grid_count(bins)
@@ -82,7 +84,8 @@ def update_axes(ax, mappable=None):
 
 # %%
 # Display the ratio between the short and total presence.
-# Red = mostly short networks
+#
+# Light = mostly short networks
 ax = start_axes("")
 m = g_20.display(
     ax,
@@ -95,7 +98,9 @@ def update_axes(ax, mappable=None):
 
 # %%
 # Display the ratio between the short and total presence.
-# Red = mostly short networks
+#
+# Light = mostly short networks
+#
 # Networks shorter than 365 days are masked
 ax = start_axes("")
 m = g_20.display(
@@ -110,9 +115,10 @@ def update_axes(ax, mappable=None):
 update_axes(ax, m).set_label("Pixel used in % all atlas")
 # %%
 # Display the ratio between the short and total presence.
-# Red = mostly short networks
+#
 # Networks longer than 365 days are masked
-# -> Coastal areas are mostly populated by short networks
+#
+# # -> Coastal areas are mostly populated by short networks
 ax = start_axes("")
 m = g_20.display(
     ax,
@@ -133,8 +139,9 @@ def update_axes(ax, mappable=None):
 # Display the occurence of merging events
 ax = start_axes("")
 g_all_merging = n.merging_event().grid_count(bins)
-m = g_all_merging.display(ax, **kw_time, vmin=0, vmax=0.2)
+m = g_all_merging.display(ax, **kw_time, vmin=0, vmax=1)
 update_axes(ax, m).set_label("Pixel used in % of time")
+
 # %%
 # Ratio merging events / eddy presence
 ax = start_axes("")
@@ -152,15 +159,15 @@ def update_axes(ax, mappable=None):
 # ---------------------------------------
 ax = start_axes("")
 g_10_merging = n10.merging_event().grid_count(bins)
-m = g_10_merging.display(ax, **kw_time, vmin=0, vmax=0.2)
+m = g_10_merging.display(ax, **kw_time, vmin=0, vmax=1)
 update_axes(ax, m).set_label("Pixel used in % of time")
 # %%
 ax = start_axes("")
 m = g_10_merging.display(
     ax,
     **kw_ratio,
     vmin=0,
-    vmax=0.5,
+    vmax=2,
     name=ma.array(
         g_10_merging.vars["count"] * 100.0 / g_10.vars["count"],
         mask=g_10.vars["count"] < 365,

examples/16_network/pet_ioannou_2017_case.py

@@ -69,6 +69,7 @@ def update_axes(ax, mappable=None):
 
 # %%
 # We know the position and the time of a specific eddy
+#
 # `n.extract_with_mask` give us the corresponding network
 n = NetworkObservations.load_file(
     "med/Anticyclonic_seg.nc"

examples/16_network/pet_relative.py

@@ -32,6 +32,7 @@
 # %%
 # Display timeline with events
 # A segment generated by a splitting is marked with a star
+#
 # A segment merging in another is marked with an exagon
 fig = plt.figure(figsize=(15, 5))
 ax = fig.add_axes([0.04, 0.04, 0.92, 0.92])
@@ -82,6 +83,7 @@
 # Parameters timeline
 # -------------------
 # Scatter is usefull to display the parameters' temporal evolution
+#
 # Effective Radius and Amplitude
 kw = dict(s=25, cmap="Spectral_r", zorder=10)
 fig = plt.figure(figsize=(15, 8))
@@ -125,7 +127,7 @@
 # Remove all tiny segments with less than N obs which didn't join two segments
 #
 # .. warning::
-#     Must be explore, no solution to solve all the case
+#     Must be explore, no solution to solve all cases
 
 n_clean = n.remove_dead_end(nobs=10)
 fig = plt.figure(figsize=(15, 8))
@@ -140,6 +142,7 @@
 # Keep close relative
 # -------------------
 # When you want to investigate one particular observation and select only the closest segments
+#
 # First choose an observation in the network
 fig = plt.figure(figsize=(15, 5))
 ax = fig.add_axes([0.04, 0.06, 0.90, 0.88])
@@ -201,21 +204,36 @@
 # Display the position of the eddies after a merging
 fig = plt.figure(figsize=(15, 8))
 ax = fig.add_axes([0.04, 0.06, 0.90, 0.88], projection="full_axes")
-merging = close_to_i2.merging_event()
-merging.display(ax)
+close_to_i2.plot(ax, color_cycle=close_to_i2.COLORS)
+m1, m0, m0_stop = close_to_i2.merging_event(triplet=True)
+m1.display(ax, color="violet", lw=2, label="Eddies after merging")
+m0.display(ax, color="blueviolet", lw=2, label="Eddies before merging")
+m0_stop.display(ax, color="black", lw=2, label="Eddies stopped by merging")
+ax.plot(m1.lon, m1.lat, marker=".", color="purple", ls="")
+ax.plot(m0.lon, m0.lat, marker=".", color="blueviolet", ls="")
+ax.plot(m0_stop.lon, m0_stop.lat, marker=".", color="black", ls="")
+ax.legend()
 ax.set_xlim(-10, 20), ax.set_ylim(-37, -21), ax.grid()
-merging
+m1
+
 
 # %%
 # Get spliting event
 # ------------------
 # Display the position of the eddies before a splitting
 fig = plt.figure(figsize=(15, 8))
 ax = fig.add_axes([0.04, 0.06, 0.90, 0.88], projection="full_axes")
-spliting = close_to_i2.spliting_event()
-spliting.display(ax)
+close_to_i2.plot(ax, color_cycle=close_to_i2.COLORS)
+s0, s1, s1_start = close_to_i2.spliting_event(triplet=True)
+s0.display(ax, color="violet", lw=2, label="Eddies before splitting")
+s1.display(ax, color="blueviolet", lw=2, label="Eddies after splitting")
+s1_start.display(ax, color="black", lw=2, label="Eddies starting by splitting")
+ax.plot(s0.lon, s0.lat, marker=".", color="purple", ls="")
+ax.plot(s1.lon, s1.lat, marker=".", color="blueviolet", ls="")
+ax.plot(s1_start.lon, s1_start.lat, marker=".", color="black", ls="")
+ax.legend()
 ax.set_xlim(-10, 20), ax.set_ylim(-37, -21), ax.grid()
-spliting
+s1
 
 # %%
 # Get birth event

src/py_eddy_tracker/observations/network.py

@@ -113,8 +113,8 @@ def longer_than(self, nb_day_min=-1, nb_day_max=-1):
         """
         Select network on time duration
 
-        :param int nb_day_min: Minimal number of day which must be covered by one network, if negative -> not used
-        :param int nb_day_max: Maximal number of day which must be covered by one network, if negative -> not used
+        :param int nb_day_min: Minimal number of day covered by one network, if negative -> not used
+        :param int nb_day_max: Maximal number of day covered by one network, if negative -> not used
         """
         if nb_day_max < 0:
             nb_day_max = 1000000000000
@@ -132,7 +132,7 @@ def longer_than(self, nb_day_min=-1, nb_day_max=-1):
     @classmethod
     def from_split_network(cls, group_dataset, indexs, **kwargs):
         """
-        Build a NetworkObservations object with Group dataset and indexs
+        Build a NetworkObservations object with Group dataset and indexes
 
         :param TrackEddiesObservations group_dataset: Group dataset
         :param indexs: result from split_network
@@ -204,6 +204,9 @@ def __close_segment(cls, father, shift, connexions, distance):
             cls.__close_segment(son, shift, connexions, distance)
 
     def segment_relative_order(self, seg_origine):
+        """
+        Compute the relative order of each segment to the chosen segment
+        """
         i_s, i_e, i_ref = build_index(self.segment)
         segment_connexions = self.connexions()
         relative_tr = -ones(i_s.shape, dtype="i4")
@@ -216,6 +219,9 @@ def segment_relative_order(self, seg_origine):
         return d
 
     def relative(self, i_obs, order=2, direct=True, only_past=False, only_future=False):
+        """
+        Extract the segments at a certain order.
+        """
         d = self.segment_relative_order(self.segment[i_obs])
         m = (d <= order) * (d != -1)
         return self.extract_with_mask(m)
@@ -234,7 +240,7 @@ def only_one_network(self):
         """
         _, i_start, _ = self.index_network
         if len(i_start) > 1:
-            raise Exception("Several network")
+            raise Exception("Several networks")
 
     def position_filter(self, median_half_window, loess_half_window):
         self.median_filter(median_half_window, "time", "lon").loess_filter(
@@ -266,7 +272,15 @@ def display_timeline(
         self, ax, event=True, field=None, method=None, factor=1, **kwargs
     ):
         """
-        Must be call on only one network
+        Plot a timeline of a network.
+        Must be called on only one network.
+
+        :param matplotlib.axes.Axes ax: matplotlib axe used to draw
+        :param bool event: if True, draw the splitting and merging events
+        :param str,array field: yaxis values, if None, segments are used
+        :param str method: if None, mean values are used
+        :param float factor: to multiply field
+        :return: plot mappable
         """
         self.only_one_network()
         j = 0
@@ -516,6 +530,7 @@ def extract_event(self, indices):
 
     @property
     def segment_track_array(self):
+        """Return a unique segment id when multiple networks are considered"""
         return build_unique_array(self.segment, self.track)
 
     def birth_event(self):
@@ -542,27 +557,65 @@ def death_event(self):
                 indices.append(i.stop - 1)
         return self.extract_event(list(set(indices)))
 
-    def merging_event(self):
-        indices = list()
+    def merging_event(self, triplet=False):
+        """Return observation after a merging event.
+
+        If `triplet=True` return the eddy after a merging event, the eddy before the merging event,
+        and the eddy stopped due to merging.
+        """
+        idx_m1 = list()
+        if triplet:
+            idx_m0_stop = list()
+            idx_m0 = list()
+
         for i, _, _ in self.iter_on(self.segment_track_array):
             nb = i.stop - i.start
             if nb == 0:
                 continue
             i_n = self.next_obs[i.stop - 1]
             if i_n != -1:
-                indices.append(i.stop - 1)
-        return self.extract_event(list(set(indices)))
+                if triplet:
+                    idx_m0_stop.append(i.stop - 1)
+                    idx_m0.append(self.previous_obs[i_n])
+                idx_m1.append(i_n)
+
+        if triplet:
+            return (
+                self.extract_event(list(idx_m1)),
+                self.extract_event(list(idx_m0)),
+                self.extract_event(list(idx_m0_stop)),
+            )
+        else:
+            return self.extract_event(list(set(idx_m1)))
 
-    def spliting_event(self):
-        indices = list()
+    def spliting_event(self, triplet=False):
+        """Return observation before a splitting event.
+
+        If `triplet=True` return the eddy before a splitting event, the eddy after the splitting event,
+        and the eddy starting due to splitting.
+        """
+        idx_s0 = list()
+        if triplet:
+            idx_s1_start = list()
+            idx_s1 = list()
         for i, _, _ in self.iter_on(self.segment_track_array):
             nb = i.stop - i.start
             if nb == 0:
                 continue
             i_p = self.previous_obs[i.start]
             if i_p != -1:
-                indices.append(i.start)
-        return self.extract_event(list(set(indices)))
+                if triplet:
+                    idx_s1_start.append(i.start)
+                    idx_s1.append(self.next_obs[i_p])
+                idx_s0.append(i_p)
+        if triplet:
+            return (
+                self.extract_event(list(idx_s0)),
+                self.extract_event(list(idx_s1)),
+                self.extract_event(list(idx_s1_start)),
+            )
+        else:
+            return self.extract_event(list(set(idx_s0)))
 
     def dissociate_network(self):
         """
@@ -655,7 +708,7 @@ def remove_dead_end(self, nobs=3, recursive=0, mask=None):
         self.only_one_network()
         segments_keep = list()
         connexions = self.connexions()
-        for i, b0, b1 in self.iter_on("segment"):
+        for i, b0, _ in self.iter_on("segment"):
             nb = i.stop - i.start
             if mask and mask[i].any():
                 segments_keep.append(b0)
@@ -852,6 +905,7 @@ def apply_replace(x, x0, x1):
 
 @njit(cache=True)
 def build_unique_array(id1, id2):
+    """Give a unique id for each (id1, id2) with id1 and id2 increasing monotonically"""
     k = 0
     new_id = empty(id1.shape, dtype=id1.dtype)
     id1_previous = id1[0]