- Add method to colorize contour with a field

- Add option to force align on to return all step for reference dataset
- Add method and property to network to easily select segment and network
- Add method to found same track/segment/network in  dataset
- Rewrite particle candidate to be easily parallelize

Author: AntSimi

setup.cfg

@@ -1,8 +1,22 @@
 [flake8]
+max-line-length = 140
 ignore =
     E203, # whitespace before ':'
     W503,  # line break before binary operator
 
+[isort]
+combine_as_imports=True
+force_grid_wrap=0
+force_sort_within_sections=True
+force_to_top=typing
+include_trailing_comma=True
+line_length=140
+multi_line_output=3
+skip=
+    build
+    docs/source/conf.py
+
+
 [versioneer]
 VCS = git
 style = pep440

setup.py

@@ -48,6 +48,7 @@
             "EddyNetworkGroup = py_eddy_tracker.appli.network:build_network",
             "EddyNetworkBuildPath = py_eddy_tracker.appli.network:divide_network",
             "EddyNetworkSubSetter = py_eddy_tracker.appli.network:subset_network",
+            "EddyNetworkQuickCompare = py_eddy_tracker.appli.network:quick_compare",
             # anim/gui
             "EddyAnim = py_eddy_tracker.appli.gui:anim",
             "GUIEddy = py_eddy_tracker.appli.gui:guieddy",

src/py_eddy_tracker/__init__.py

@@ -422,14 +422,20 @@ def identify_time(str_date):
         nc_name="previous_cost",
         nc_type="float32",
         nc_dims=("obs",),
-        nc_attr=dict(long_name="Previous cost for previous observation", comment="",),
+        nc_attr=dict(
+            long_name="Previous cost for previous observation",
+            comment="",
+        ),
     ),
     next_cost=dict(
         attr_name=None,
         nc_name="next_cost",
         nc_type="float32",
         nc_dims=("obs",),
-        nc_attr=dict(long_name="Next cost for next observation", comment="",),
+        nc_attr=dict(
+            long_name="Next cost for next observation",
+            comment="",
+        ),
     ),
     n=dict(
         attr_name=None,
@@ -640,7 +646,8 @@ def identify_time(str_date):
         nc_type="f4",
         nc_dims=("obs",),
         nc_attr=dict(
-            long_name="Log base 10 background chlorophyll", units="Log(Chl/[mg/m^3])",
+            long_name="Log base 10 background chlorophyll",
+            units="Log(Chl/[mg/m^3])",
         ),
     ),
     year=dict(

src/py_eddy_tracker/appli/eddies.py

@@ -243,7 +243,8 @@ def browse_dataset_in(
         filenames = bytes_(glob(full_path))
 
     dataset_list = empty(
-        len(filenames), dtype=[("filename", "S500"), ("date", "datetime64[s]")],
+        len(filenames),
+        dtype=[("filename", "S500"), ("date", "datetime64[s]")],
     )
     dataset_list["filename"] = filenames
 
@@ -371,7 +372,8 @@ def track(
 
     logger.info("Longer track saved have %d obs", c.nb_obs_by_tracks.max())
     logger.info(
-        "The mean length is %d observations for long track", c.nb_obs_by_tracks.mean(),
+        "The mean length is %d observations for long track",
+        c.nb_obs_by_tracks.mean(),
     )
 
     long_track.write_file(**kw_write)
@@ -381,7 +383,14 @@ def track(
 
 
 def get_group(
-    dataset1, dataset2, index1, index2, score, invalid=2, low=10, high=60,
+    dataset1,
+    dataset2,
+    index1,
+    index2,
+    score,
+    invalid=2,
+    low=10,
+    high=60,
 ):
     group1, group2 = dict(), dict()
     m_valid = (score * 100) >= invalid
@@ -490,7 +499,8 @@ def get_values(v, dataset):
             ]
 
     labels = dict(
-        high=f"{high:0.0f} <= high", low=f"{invalid:0.0f} <= low < {low:0.0f}",
+        high=f"{high:0.0f} <= high",
+        low=f"{invalid:0.0f} <= low < {low:0.0f}",
     )
 
     keys = [labels.get(key, key) for key in list(gr_ref.values())[0].keys()]

src/py_eddy_tracker/appli/network.py

@@ -8,6 +8,7 @@
 from .. import EddyParser
 from ..observations.network import Network, NetworkObservations
 from ..observations.tracking import TrackEddiesObservations
+from numpy import in1d, zeros
 
 logger = logging.getLogger("pet")
 
@@ -128,3 +129,135 @@ def subset_network():
     if args.period is not None:
         n = n.extract_with_period(args.period)
     n.write_file(filename=args.out)
+
+
+def quick_compare():
+    parser = EddyParser(
+        """Tool to have a quick comparison between several network:
+        - N : network
+        - S : segment
+        - Obs : observations
+        """
+        
+    )
+    parser.add_argument("ref", help="Identification file of reference")
+    parser.add_argument("others", nargs="+", help="Identifications files to compare")
+    parser.add_argument(
+        "--path_out", default=None, help="Save each group in separate file"
+    )
+    args = parser.parse_args()
+
+    kw = dict(
+        include_vars=['longitude', 'latitude', 'time', 'track', 'segment', 'next_obs', 'previous_obs']
+    )
+
+    if args.path_out is not None:
+        kw = dict()
+
+    ref = NetworkObservations.load_file(args.ref, **kw)
+    print(
+        f"[ref] {args.ref} -> {ref.nb_network} network / {ref.nb_segment} segment / {len(ref)} obs "
+        f"-> {ref.network_size(0)} trash obs, "
+        f"{len(ref.merging_event())} merging, {len(ref.splitting_event())} spliting"
+        )
+    others = {other: NetworkObservations.load_file(other, **kw) for other in args.others}
+
+    if args.path_out is not None:
+        groups_ref, groups_other = run_compare(ref, others, **kwargs)
+        if not exists(args.path_out):
+            mkdir(args.path_out)
+        for i, other_ in enumerate(args.others):
+            dirname_ = f"{args.path_out}/{other_.replace('/', '_')}/"
+            if not exists(dirname_):
+                mkdir(dirname_)
+            for k, v in groups_other[other_].items():
+                basename_ = f"other_{k}.nc"
+                others[other_].index(v).write_file(filename=f"{dirname_}/{basename_}")
+            for k, v in groups_ref[other_].items():
+                basename_ = f"ref_{k}.nc"
+                ref.index(v).write_file(filename=f"{dirname_}/{basename_}")
+        return
+    display_compare(ref, others)
+
+
+def run_compare(ref, others):
+    outs = dict()
+    for i, (k, other) in enumerate(others.items()):
+        out = dict()
+        print(
+            f"[{i}]   {k} -> {other.nb_network} network / {other.nb_segment} segment / {len(other)} obs "
+            f"-> {other.network_size(0)} trash obs, "
+            f"{len(other.merging_event())} merging, {len(other.splitting_event())} spliting"
+            )
+        ref_id, other_id = ref.identify_in(other, size_min=2)
+        m = other_id != -1
+        ref_id, other_id = ref_id[m], other_id[m]
+        out['same N(N)'] = m.sum()
+        out['same N(Obs)'] = ref.network_size(ref_id).sum()
+
+        # For network which have same obs
+        ref_, other_ = ref.networks(ref_id), other.networks(other_id)
+        ref_segu, other_segu = ref_.identify_in(other_, segment=True)
+        m = other_segu==-1
+        ref_track_no_match, _ = ref_.unique_segment_to_id(ref_segu[m])
+        ref_segu, other_segu = ref_segu[~m], other_segu[~m]
+        m = ~in1d(ref_id, ref_track_no_match)
+        out['same NS(N)'] = m.sum()
+        out['same NS(Obs)'] = ref.network_size(ref_id[m]).sum()
+
+        # Check merge/split
+        def follow_obs(d, i_follow):
+            m = i_follow != -1
+            i_follow = i_follow[m]
+            t, x, y = zeros(m.size, d.time.dtype), zeros(m.size, d.longitude.dtype), zeros(m.size, d.latitude.dtype)
+            t[m], x[m], y[m] = d.time[i_follow], d.longitude[i_follow], d.latitude[i_follow]
+            return t, x, y
+        def next_obs(d, i_seg):
+            last_i = d.index_segment_track[1][i_seg] - 1
+            return follow_obs(d, d.next_obs[last_i])
+        def previous_obs(d, i_seg):
+            first_i = d.index_segment_track[0][i_seg]
+            return follow_obs(d, d.previous_obs[first_i])
+
+        tref, xref, yref = next_obs(ref_, ref_segu)
+        tother, xother, yother = next_obs(other_, other_segu)
+
+        m = (tref == tother) & (xref == xother) & (yref == yother) 
+        print(m.sum(), m.size, ref_segu.size, ref_track_no_match.size)
+
+        tref, xref, yref = previous_obs(ref_, ref_segu)
+        tother, xother, yother = previous_obs(other_, other_segu)
+
+        m = (tref == tother) & (xref == xother) & (yref == yother) 
+        print(m.sum(), m.size, ref_segu.size, ref_track_no_match.size)
+
+
+
+        ref_segu, other_segu = ref.identify_in(other, segment=True)
+        m = other_segu != -1
+        out['same S(S)'] = m.sum()
+        out['same S(Obs)'] = ref.segment_size()[ref_segu[m]].sum()
+
+        outs[k] = out
+    return outs
+
+def display_compare(ref, others):
+    def display(value, ref=None):
+        if ref:
+            outs = [f"{v/ref[k] * 100:.1f}% ({v})" for k, v in value.items()]
+        else:
+            outs = value
+        return "".join([f"{v:^18}" for v in outs])
+    
+    datas = run_compare(ref, others)
+    ref_ = {
+        'same N(N)' : ref.nb_network,
+        "same N(Obs)": len(ref),
+        'same NS(N)' : ref.nb_network,
+        'same NS(Obs)' : len(ref),
+        'same S(S)' : ref.nb_segment,
+        'same S(Obs)' : len(ref),
+        }
+    print("     ", display(ref_.keys()))
+    for i, (_, v) in enumerate(datas.items()):
+        print(f"[{i:2}] ", display(v, ref=ref_))

src/py_eddy_tracker/dataset/grid.py

@@ -310,6 +310,11 @@ def __init__(
             self.load_general_features()
             self.load()
 
+    def populate(self):
+        if self.dimensions is None:
+            self.load_general_features()
+            self.load()
+
     @property
     def is_centered(self):
         """Give True if pixel is described with its center's position or
@@ -539,7 +544,8 @@ def grid(self, varname, indexs=None):
                         self.vars[varname] = self.vars[varname].T
             if self.nan_mask:
                 self.vars[varname] = ma.array(
-                    self.vars[varname], mask=isnan(self.vars[varname]),
+                    self.vars[varname],
+                    mask=isnan(self.vars[varname]),
                 )
             if not hasattr(self.vars[varname], "mask"):
                 self.vars[varname] = ma.array(
@@ -869,7 +875,9 @@ def eddy_identification(
                         num_fac=presampling_multiplier,
                     )
                     xy_e = uniform_resample(
-                        contour.lon, contour.lat, num_fac=presampling_multiplier,
+                        contour.lon,
+                        contour.lat,
+                        num_fac=presampling_multiplier,
                     )
                     xy_s = uniform_resample(
                         speed_contour.lon,
@@ -1204,7 +1212,9 @@ def setup_coordinates(self):
         dx = self.x_bounds[1:] - self.x_bounds[:-1]
         dy = self.y_bounds[1:] - self.y_bounds[:-1]
         if (dx < 0).any() or (dy < 0).any():
-            raise Exception("Coordinates in RegularGridDataset must be strictly increasing")
+            raise Exception(
+                "Coordinates in RegularGridDataset must be strictly increasing"
+            )
         self._x_step = (self.x_c[1:] - self.x_c[:-1]).mean()
         self._y_step = (self.y_c[1:] - self.y_c[:-1]).mean()
 
@@ -1736,7 +1746,7 @@ def compute_stencil(
             self.x_c,
             self.y_c,
             data.data,
-            data.mask,
+            self.get_mask(data),
             self.EARTH_RADIUS,
             vertical=vertical,
             stencil_halfwidth=stencil_halfwidth,
@@ -2285,23 +2295,23 @@ def __init__(self):
         self.datasets = list()
 
     @classmethod
-    def from_netcdf_cube(cls, filename, x_name, y_name, t_name, heigth=None):
+    def from_netcdf_cube(cls, filename, x_name, y_name, t_name, heigth=None, **kwargs):
         new = cls()
         with Dataset(filename) as h:
             for i, t in enumerate(h.variables[t_name][:]):
-                d = RegularGridDataset(filename, x_name, y_name, indexs={t_name: i})
+                d = RegularGridDataset(filename, x_name, y_name, indexs={t_name: i}, **kwargs)
                 if heigth is not None:
                     d.add_uv(heigth)
                 new.datasets.append((t, d))
         return new
 
     @classmethod
-    def from_netcdf_list(cls, filenames, t, x_name, y_name, indexs=None, heigth=None):
+    def from_netcdf_list(cls, filenames, t, x_name, y_name, indexs=None, heigth=None, **kwargs):
         new = cls()
         for i, _t in enumerate(t):
             filename = filenames[i]
             logger.debug(f"load file {i:02d}/{len(t)} t={_t} : {filename}")
-            d = RegularGridDataset(filename, x_name, y_name, indexs=indexs)
+            d = RegularGridDataset(filename, x_name, y_name, indexs=indexs, **kwargs)
             if heigth is not None:
                 d.add_uv(heigth)
             new.datasets.append((_t, d))
@@ -2349,6 +2359,7 @@ def __iter__(self):
     def __getitem__(self, item):
         for t, d in self.datasets:
             if t == item:
+                d.populate()
                 return d
         raise KeyError(item)
 
@@ -2448,10 +2459,13 @@ def advect(
         :param array y: Latitude of obs to move
         :param str,array u_name: U field to advect obs
         :param str,array v_name: V field to advect obs
+        :param float t_init: time to start advection
+        :param array,None mask_particule: advect only i mask is True
         :param int nb_step: Number of iteration before to release data
         :param int time_step: Number of second for each advection
+        :param bool rk4: Use rk4 algorithm instead of finite difference
 
-        :return: x,y position
+        :return: t,x,y position
 
         .. minigallery:: py_eddy_tracker.GridCollection.advect
         """
@@ -2477,7 +2491,7 @@ def advect(
         else:
             mask_particule += isnan(x) + isnan(y)
         while True:
-            logger.debug(f"advect : t={t}")
+            logger.debug(f"advect : t={t/86400}")
             if (backward and t <= t1) or (not backward and t >= t1):
                 t0, u0, v0, m0 = t1, u1, v1, m1
                 t1, d1 = generator.__next__()
@@ -2507,7 +2521,7 @@ def get_next_time_step(self, t_init):
         for i, (t, dataset) in enumerate(self.datasets):
             if t < t_init:
                 continue
-
+            dataset.populate()
             logger.debug(f"i={i}, t={t}, dataset={dataset}")
             yield t, dataset
 
@@ -2517,7 +2531,7 @@ def get_previous_time_step(self, t_init):
             i -= 1
             if t > t_init:
                 continue
-
+            dataset.populate()
             logger.debug(f"i={i}, t={t}, dataset={dataset}")
             yield t, dataset