Add an entry to produce grid of frequency

Author: AntSimi

setup.py

@@ -40,6 +40,7 @@
             "EddyId = py_eddy_tracker.appli.grid:eddy_id",
             # eddies
             "MergeEddies = py_eddy_tracker.appli.eddies:merge_eddies",
+            "EddyFrequency = py_eddy_tracker.appli.eddies:get_frequency_grid",
             "EddyInfos = py_eddy_tracker.appli.eddies:display_infos",
             # network
             "EddyNetworkGroup = py_eddy_tracker.appli.network:build_network",

src/py_eddy_tracker/appli/eddies.py

@@ -57,6 +57,37 @@ def merge_eddies():
     obs.write_file(filename=args.out)
 
 
+def get_frequency_grid():
+    parser = EddyParser("Compute eddy frequency")
+    parser.add_argument("observations", help="Input observations to compute frequency")
+    parser.add_argument("out", help="Grid output file")
+    parser.add_argument(
+        "--intern",
+        help="Use speed contour instead of effective contour",
+        action="store_true",
+    )
+    parser.add_argument(
+        "--xrange", nargs="+", type=float, help="Horizontal range : START,STOP,STEP"
+    )
+    parser.add_argument(
+        "--yrange", nargs="+", type=float, help="Vertical range : START,STOP,STEP"
+    )
+    args = parser.parse_args()
+
+    if (args.xrange is None or len(args.xrange) not in (3,)) or (
+        args.yrange is None or len(args.yrange) not in (3,)
+    ):
+        raise Exception("Use START/STOP/STEP for --xrange and --yrange")
+
+    var_to_load = ["longitude"]
+    var_to_load.extend(EddiesObservations.intern(args.intern, public_label=True))
+    e = EddiesObservations.load_file(args.observations, include_vars=var_to_load)
+
+    bins = args.xrange, args.yrange
+    g = e.grid_count(bins, intern=args.intern)
+    g.write(args.out)
+
+
 def display_infos():
     parser = EddyParser("Display General inforamtion")
     parser.add_argument("observations", nargs='+', help="Input observations to compute frequency")

src/py_eddy_tracker/appli/gui.py

@@ -109,7 +109,7 @@ def prev(self):
         return self.draw_contour()
 
     def draw_contour(self):
-        t0, t1 = self.period
+        # t0, t1 = self.period
         # select contour for this time step
         m = self.t == self.now
         self.ax.figure.canvas.restore_region(self.bg_cache)

src/py_eddy_tracker/dataset/grid.py

@@ -584,7 +584,8 @@ def eddy_identification(
         :param float,int step: Height between two layers in m
         :param float,int shape_error: Maximal error allow for outter contour in %
         :param int sampling: Sampling of contour and speed profile
-        :param (int,int),None pixel_limit: Min and max of pixel which must be inside inner and outer contour to be considered like an eddy
+        :param (int,int),None pixel_limit:
+            Min and max of pixel which must be inside inner and outer contour to be considered like an eddy
         :param float,None precision: Truncate value at the defined precision in m
         :param str force_height_unit: Unit to used for height unit
         :param str force_speed_unit: Unit to used for speed unit

src/py_eddy_tracker/observations/observation.py

@@ -49,6 +49,7 @@
     histogram2d,
     linspace,
     sin,
+    histogram,
 )
 from netCDF4 import Dataset
 from datetime import datetime
@@ -1359,16 +1360,22 @@ def grid_count(self, bins, intern=False, center=False):
             ),
             centered=True,
         )
+        debug_active = logger.getEffectiveLevel() == logging.DEBUG
         if center:
             x, y = (self.longitude - x0) % 360 + x0, self.latitude
             grid[:] = histogram2d(x, y, (x_bins, y_bins))[0]
             grid.mask = grid.data == 0
         else:
             x_ref = ((self.longitude - x0) % 360 + x0 - 180).reshape(-1, 1)
             x, y = (self[x_name] - x_ref) % 360 + x_ref, self[y_name]
-            for x_, y_ in zip(x, y):
+            nb = x_ref.shape[0]
+            for i_, (x_, y_) in enumerate(zip(x, y)):
+                if debug_active and i_ % 10000 == 0:
+                    print(f"{i_}/{nb}", end="\r")
                 i, j = BasePath(create_vertice(x_, y_)).pixels_in(regular_grid)
                 grid_count_(grid, i, j)
+            if debug_active:
+                print()
             grid.mask = grid == 0
         return regular_grid