Use last obs of network, manage wrapping of network around 0

Author: AntSimi

examples/12_external_data/pet_SST_collocation.py

@@ -27,7 +27,7 @@
 
 # %%
 # 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"
@@ -58,14 +58,14 @@ def update_axes(ax, mappable=None, unit=""):
 
 # %%
 # 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]")
 
 # %%
 # SST first display
-# -----------------------------
+# -----------------
 
 # %%
 # We can now plot SST from `sst`

src/py_eddy_tracker/appli/network.py

@@ -5,15 +5,9 @@
 
 import logging
 
-from netCDF4 import Dataset
-from numpy import arange, empty, zeros
-from Polygon import Polygon
-
 from .. import EddyParser
-from ..generic import build_index
 from ..observations.network import Network
 from ..observations.tracking import TrackEddiesObservations
-from ..poly import create_vertice_from_2darray, polygon_overlap
 
 logger = logging.getLogger("pet")
 
@@ -57,246 +51,6 @@ def divide_network():
         include_vars=("time", "track", "latitude", "longitude", *contour_name),
     )
     ids = e.split_network(intern=args.intern, window=args.window)
-
-
-def split_network(input, output):
-    """Divide each group in track"""
-    sl = slice(None)
-    with Dataset(input) as h:
-        group = h.variables["track"][sl]
-    track_s, track_e, track_ref = build_index(group)
-    # nb = track_e - track_s
-    # m = nb > 1500
-    # print(group[track_s[m]])
-
-    track_id = 12003
-    sls = [slice(track_s[track_id - track_ref], track_e[track_id - track_ref], None)]
-    for sl in sls:
-
-        print(sl)
-        with Dataset(input) as h:
-            time = h.variables["time"][sl]
-            group = h.variables["track"][sl]
-            x = h.variables["effective_contour_longitude"][sl]
-            y = h.variables["effective_contour_latitude"][sl]
-        print(group[0])
-        ids = empty(
-            time.shape,
-            dtype=[
-                ("group", group.dtype),
-                ("time", time.dtype),
-                ("track", "u2"),
-                ("previous_cost", "f4"),
-                ("next_cost", "f4"),
-                ("previous_observation", "i4"),
-                ("next_observation", "i4"),
-            ],
-        )
-        ids["group"] = group
-        ids["time"] = time
-        # To store id track
-        ids["track"] = 0
-        ids["previous_cost"] = 0
-        ids["next_cost"] = 0
-        ids["previous_observation"] = -1
-        ids["next_observation"] = -1
-        # Cost with previous
-        track_start, track_end, track_ref = build_index(group)
-        for i0, i1 in zip(track_start, track_end):
-            if (i1 - i0) == 0 or group[i0] == Network.NOGROUP:
-                continue
-            sl_group = slice(i0, i1)
-            set_tracks(
-                x[sl_group],
-                y[sl_group],
-                time[sl_group],
-                i0,
-                ids["track"][sl_group],
-                ids["previous_cost"][sl_group],
-                ids["next_cost"][sl_group],
-                ids["previous_observation"][sl_group],
-                ids["next_observation"][sl_group],
-                window=5,
-            )
-
-        new_i = ids.argsort(order=("group", "track", "time"))
-        ids_sort = ids[new_i]
-        # To be able to follow indices sorting
-        reverse_sort = empty(new_i.shape[0], dtype="u4")
-        reverse_sort[new_i] = arange(new_i.shape[0])
-        # Redirect indices
-        m = ids_sort["next_observation"] != -1
-        ids_sort["next_observation"][m] = reverse_sort[ids_sort["next_observation"][m]]
-        m = ids_sort["previous_observation"] != -1
-        ids_sort["previous_observation"][m] = reverse_sort[
-            ids_sort["previous_observation"][m]
-        ]
-        # print(ids_sort)
-        display_network(
-            x[new_i],
-            y[new_i],
-            ids_sort["track"],
-            ids_sort["time"],
-            ids_sort["next_cost"],
-        )
-
-
-def next_obs(
-    i_current, next_cost, previous_cost, polygons, t, t_start, t_end, t_ref, window
-):
-    t_max = t_end.shape[0] - 1
-    t_cur = t[i_current]
-    t0, t1 = t_cur + 1 - t_ref, t_cur + window - t_ref
-    if t0 > t_max:
-        return -1
-    t1 = min(t1, t_max)
-    for t_step in range(t0, t1 + 1):
-        i0, i1 = t_start[t_step], t_end[t_step]
-        # No observation at the time step !
-        if i0 == i1:
-            continue
-        sl = slice(i0, i1)
-        # Intersection / union, to be able to separte in case of multiple inside
-        c = polygon_overlap(polygons[i_current], polygons[sl])
-        # We remove low overlap
-        if (c > 0.1).sum() > 1:
-            print(c)
-        c[c < 0.1] = 0
-        # We get index of maximal overlap
-        i = c.argmax()
-        c_i = c[i]
-        # No overlap found
-        if c_i == 0:
-            continue
-        target = i0 + i
-        # Check if candidate is already used
-        c_target = previous_cost[target]
-        if (c_target != 0 and c_target < c_i) or c_target == 0:
-            previous_cost[target] = c_i
-        next_cost[i_current] = c_i
-        return target
-    return -1
-
-
-def set_tracks(
-    x,
-    y,
-    t,
-    ref_index,
-    track,
-    previous_cost,
-    next_cost,
-    previous_observation,
-    next_observation,
-    window,
-):
-    # Will split one group in tracks
-    t_start, t_end, t_ref = build_index(t)
-    nb = x.shape[0]
-    used = zeros(nb, dtype="bool")
-    current_track = 1
-    # build all polygon (need to check if wrap is needed)
-    polygons = list()
-    for i in range(nb):
-        polygons.append(Polygon(create_vertice_from_2darray(x, y, i)))
-
-    for i in range(nb):
-        # If observation already in one track, we go to the next one
-        if used[i]:
-            continue
-        build_track(
-            i,
-            current_track,
-            used,
-            track,
-            previous_observation,
-            next_observation,
-            ref_index,
-            next_cost,
-            previous_cost,
-            polygons,
-            t,
-            t_start,
-            t_end,
-            t_ref,
-            window,
-        )
-        current_track += 1
-
-
-def build_track(
-    first_index,
-    track_id,
-    used,
-    track,
-    previous_observation,
-    next_observation,
-    ref_index,
-    next_cost,
-    previous_cost,
-    *args,
-):
-    i_next = first_index
-    while i_next != -1:
-        # Flag
-        used[i_next] = True
-        # Assign id
-        track[i_next] = track_id
-        # Search next
-        i_next_ = next_obs(i_next, next_cost, previous_cost, *args)
-        if i_next_ == -1:
-            break
-        next_observation[i_next] = i_next_ + ref_index
-        if not used[i_next_]:
-            previous_observation[i_next_] = i_next + ref_index
-        # Target was previously used
-        if used[i_next_]:
-            if next_cost[i_next] == previous_cost[i_next_]:
-                m = track[i_next_:] == track[i_next_]
-                track[i_next_:][m] = track_id
-                previous_observation[i_next_] = i_next + ref_index
-            i_next_ = -1
-        i_next = i_next_
-
-
-def display_network(x, y, tr, t, c):
-    tr0, tr1, t_ref = build_index(tr)
-    import matplotlib.pyplot as plt
-
-    cmap = plt.get_cmap("jet")
-    from ..generic import flatten_line_matrix
-
-    fig = plt.figure(figsize=(20, 10))
-    ax = fig.add_subplot(121, aspect="equal")
-    ax.grid()
-    ax_time = fig.add_subplot(122)
-    ax_time.grid()
-    i = 0
-    for s, e in zip(tr0, tr1):
-        if s == e:
-            continue
-        sl = slice(s, e)
-        color = cmap((tr[s] - tr[tr0[0]]) / (tr[tr0[-1]] - tr[tr0[0]]))
-        ax.plot(
-            flatten_line_matrix(x[sl]),
-            flatten_line_matrix(y[sl]),
-            color=color,
-            label=f"{tr[s]} - {e-s} obs from {t[s]} to {t[e-1]}",
-        )
-        i += 1
-        ax_time.plot(
-            t[sl],
-            tr[s].repeat(e - s) + c[sl],
-            color=color,
-            label=f"{tr[s]} - {e-s} obs",
-            lw=0.5,
-        )
-        ax_time.plot(t[sl], tr[s].repeat(e - s), color=color, lw=1, marker="+")
-        ax_time.text(t[s], tr[s] + 0.15, f"{x[s].mean():.2f}, {y[s].mean():.2f}")
-        ax_time.axvline(t[s], color=".75", lw=0.5, ls="--", zorder=-10)
-        ax_time.text(
-            t[e - 1], tr[e - 1] - 0.25, f"{x[e-1].mean():.2f}, {y[e-1].mean():.2f}"
-        )
-    ax.legend()
-    ax_time.legend()
-    plt.show()
+    e = e.add_fields(("sub_track",))
+    e.sub_track[:] = ids["track"]
+    e.write_file(filename=args.out)

src/py_eddy_tracker/generic.py

@@ -83,7 +83,8 @@ def build_index(groups):
             first_index[group - i0 + 1 : next_group - i0 + 1] = i + 1
     last_index = zeros(amplitude, dtype=numba_types.int_)
     last_index[:-1] = first_index[1:]
-    last_index[-1] = i + 1
+    # + 2 because we iterate only until -2 and we want upper bound ( 1 + 1)
+    last_index[-1] = i + 2
     return first_index, last_index, i0
 
 

src/py_eddy_tracker/observations/tracking.py

@@ -26,11 +26,10 @@
     where,
     zeros,
 )
-from Polygon import Polygon
 
 from .. import VAR_DESCR_inv, __version__
 from ..generic import build_index, cumsum_by_track, distance, split_line, wrap_longitude
-from ..poly import create_vertice_from_2darray, merge, polygon_overlap
+from ..poly import bbox_intersection, merge, vertice_overlap
 from .observation import EddiesObservations
 
 logger = logging.getLogger("pet")
@@ -642,15 +641,14 @@ def set_tracks(self, x, y, ids, window):
         used = zeros(nb, dtype="bool")
         track_id = 1
         # build all polygon (need to check if wrap is needed)
-        polygons = [Polygon(create_vertice_from_2darray(x, y, i)) for i in range(nb)]
         for i in range(nb):
             # If observation already in one track, we go to the next one
             if used[i]:
                 continue
-            self.follow_obs(i, track_id, used, ids, polygons, *time_index, window)
+            self.follow_obs(i, track_id, used, ids, x, y, *time_index, window)
             track_id += 1
             # Search a possible ancestor
-            self.previous_obs(i, ids, polygons, *time_index, window)
+            self.previous_obs(i, ids, x, y, *time_index, window)
 
     @classmethod
     def follow_obs(cls, i_next, track_id, used, ids, *args):
@@ -676,7 +674,7 @@ def follow_obs(cls, i_next, track_id, used, ids, *args):
             i_next = i_next_
 
     @staticmethod
-    def previous_obs(i_current, ids, polygons, time_s, time_e, time_ref, window):
+    def previous_obs(i_current, ids, x, y, time_s, time_e, time_ref, window):
         time_cur = ids["time"][i_current]
         t0, t1 = time_cur - 1 - time_ref, max(time_cur - window - time_ref, 0)
         for t_step in range(t0, t1 - 1, -1):
@@ -685,7 +683,12 @@ def previous_obs(i_current, ids, polygons, time_s, time_e, time_ref, window):
             if i0 == i1:
                 continue
             # Intersection / union, to be able to separte in case of multiple inside
-            c = polygon_overlap(polygons[i_current], polygons[i0:i1], minimal_area=True)
+            xi, yi, xj, yj = x[[i_current]], y[[i_current]], x[i0:i1], y[i0:i1]
+            ii, ij = bbox_intersection(xi, yi, xj, yj)
+            if len(ii) == 0:
+                continue
+            c = zeros(len(xj))
+            c[ij] = vertice_overlap(xi[ii], yi[ii], xj[ij], yj[ij], minimal_area=True)
             # We remove low overlap
             c[c < 0.1] = 0
             # We get index of maximal overlap
@@ -699,7 +702,7 @@ def previous_obs(i_current, ids, polygons, time_s, time_e, time_ref, window):
             break
 
     @staticmethod
-    def next_obs(i_current, ids, polygons, time_s, time_e, time_ref, window):
+    def next_obs(i_current, ids, x, y, time_s, time_e, time_ref, window):
         time_max = time_e.shape[0] - 1
         time_cur = ids["time"][i_current]
         t0, t1 = time_cur + 1 - time_ref, min(time_cur + window - time_ref, time_max)
@@ -711,7 +714,12 @@ def next_obs(i_current, ids, polygons, time_s, time_e, time_ref, window):
             if i0 == i1:
                 continue
             # Intersection / union, to be able to separte in case of multiple inside
-            c = polygon_overlap(polygons[i_current], polygons[i0:i1], minimal_area=True)
+            xi, yi, xj, yj = x[[i_current]], y[[i_current]], x[i0:i1], y[i0:i1]
+            ii, ij = bbox_intersection(xi, yi, xj, yj)
+            if len(ii) == 0:
+                continue
+            c = zeros(len(xj))
+            c[ij] = vertice_overlap(xi[ii], yi[ii], xj[ij], yj[ij], minimal_area=True)
             # We remove low overlap
             c[c < 0.1] = 0
             # We get index of maximal overlap