Give special treatment to phantom patterns

This is case 2 correction in the Matlab program of Rémi Laxenaire.

Trajectories/trajectories.py

@@ -16,6 +16,8 @@ import graph_tool
 from graph_tool import topology
 import json
 import sys
+import report_graph
+import numpy as np
 
 def new_traj(ind_traj, traj_prop, n, traj_vert_ind):
     """Assign new trajectory to vertex index n."""
@@ -61,6 +63,55 @@ traj_vert_ind = []
 # of vertex indices.
 
 for n in topology.topological_sort(g):
+    # First, we check whether we have a phantom pattern. Phantom
+    # eddies may appear in altimetric data because of
+    # interpolation. In a phantom pattern, we have splitting
+    # immediately followed by merging. We do not want to let phantom
+    # patterns interrupt trajectories.
+
+    if g.vertex(n).in_degree() == 2:
+        pred = g.get_in_neighbors(n)
+
+        if np.all(g.get_in_degrees(pred) == 1) \
+           and np.all(g.get_out_degrees(pred) == 1):
+            n1 = g.get_in_neighbors(pred[0])[0]
+            n2 = g.get_in_neighbors(pred[1])[0]
+
+            if n1 == n2 and g.vertex(n1).out_degree() == 2:
+                # We have a splitting event at n1 immediately followed
+                # by a merging at n. Check the time interval between
+                # splitting and merging. The date of merging is the
+                # date of the first eddy in segment n and the date of
+                # splitting is the date of the last eddy in segment
+                # n1.
+                date_merging = report_graph.node_to_date_eddy(
+                    g.vp.name[n], g.gp.e_overestim, only_date = True)
+                date_splitting = report_graph.node_to_date_eddy(
+                    g.vp.inst_eddies[n1][- 1], g.gp.e_overestim,
+                    only_date = True)
+
+                if date_merging - date_splitting <= 6:
+                    # We have a phantom pattern. Get the vertex on the
+                    # side of the shortest path:
+                    v_short = topology.shortest_path(
+                        g, n1, n, weights = g.ep.cost_function,
+                        dag = True)[0][1]
+
+                    if traj_prop[v_short] != traj_prop[n1]:
+                        # Correct the trajectories that were assigned
+                        # according to closest successor, that is:
+                        # swap the trajectories of pred[0] and
+                        # pred[1]:
+                        traj_prop[pred[0]], traj_prop[pred[1]] \
+                            = traj_prop[pred[1]], traj_prop[pred[0]]
+                        traj_vert_ind[traj_prop[pred[0]]][- 1], \
+                            traj_vert_ind[traj_prop[pred[1]]][- 1] \
+                            = traj_vert_ind[traj_prop[pred[1]]][- 1], \
+                            traj_vert_ind[traj_prop[pred[0]]][- 1]
+
+                    traj_prop[n] = traj_prop[n1]
+                    traj_vert_ind[traj_prop[n1]].append(n)
+
     if traj_prop[n] == - 1:
         if g.vertex(n).in_degree() >= 1:
             # Find the index of the closest direct predecessor of a node: