Effective radius will be compute in same geo reference than speed radius

Author: AntSimi

src/py_eddy_tracker/dataset/grid.py

@@ -588,13 +588,13 @@ def eddy_identification(self, grid_height, uname, vname, date, step=0.005, shape
                     proj = Proj('+proj=aeqd +ellps=WGS84 +lat_0={1} +lon_0={0}'.format(*inner_contour.mean_coordinates))
                     # First, get position based on innermost
                     # contour
-                    c_x, c_y = proj(inner_contour.lon, inner_contour.lat)
-                    centx_i, centy_i, _, _ = fit_circle(c_x, c_y)
+                    centx_i, centy_i, _, _ = fit_circle(*proj(inner_contour.lon, inner_contour.lat))
                     centlon_i, centlat_i = proj(centx_i, centy_i, inverse=True)
                     # Second, get speed-based radius based on
                     # contour of max uavg
-                    c_x, c_y = proj(speed_contour.lon, speed_contour.lat)
-                    centx_s, centy_s, eddy_radius_s, aerr_s = fit_circle(c_x, c_y)
+                    centx_s, centy_s, eddy_radius_s, aerr_s = fit_circle(*proj(speed_contour.lon, speed_contour.lat))
+                    # Computed again to be coherent with speed_radius, we will be compute in same reference
+                    _, _, eddy_radius_e, aerr_e = fit_circle(*proj(current_contour.lon, current_contour.lat))
                     centlon_s, centlat_s = proj(centx_s, centy_s, inverse=True)
 
                     # Instantiate new EddyObservation object (high cost need to be review)
@@ -615,7 +615,7 @@ def eddy_identification(self, grid_height, uname, vname, date, step=0.005, shape
                     properties.obs['radius_s'] = eddy_radius_s
                     properties.obs['speed_average'] = max_average_speed
                     properties.obs['radius_e'] = eddy_radius_e
-                    properties.obs['shape_error_e'] = aerr
+                    properties.obs['shape_error_e'] = aerr_e
                     properties.obs['shape_error_s'] = aerr_s
                     properties.obs['lon'] = centlon_s
                     properties.obs['lat'] = centlat_s

src/scripts/EddyId

@@ -7,7 +7,7 @@ from datetime import datetime
 from netCDF4 import Dataset
 from py_eddy_tracker import EddyParser
 from py_eddy_tracker.dataset.grid import RegularGridDataset
-
+import zarr
 
 def id_parser():
     parser = EddyParser('Eddy Identification')
@@ -30,6 +30,9 @@ def id_parser():
                         help='Force height unit')
     parser.add_argument('--speed_unit', default=None, type=str,
                         help='Force speed unit')
+    parser.add_argument('--zarr',
+                        action='store_true',
+                        help='Output will be wrote in zarr')
     return parser
 
 
@@ -48,7 +51,13 @@ if __name__ == '__main__':
     a, c = h.eddy_identification(args.h, u, v, date, args.isoline_step, pixel_limit=(5, 2000),
                                  shape_error=args.fir_errmax, force_height_unit=args.height_unit,
                                  force_speed_unit=args.speed_unit)
-    with Dataset(args.path_out + date.strftime('/Anticyclonic_%Y%m%d.nc'), 'w') as h:
-        a.to_netcdf(h)
-    with Dataset(args.path_out + date.strftime('/Cyclonic_%Y%m%d.nc'), 'w') as h:
-        c.to_netcdf(h)
+    if args.zarr:
+        h = zarr.open(args.path_out + date.strftime('/Anticyclonic_%Y%m%d.zarr'), 'w')
+        a.to_zarr(h)
+        h = zarr.open(args.path_out + date.strftime('/Cyclonic_%Y%m%d.zarr'), 'w')
+        c.to_zarr(h)
+    else:
+        with Dataset(args.path_out + date.strftime('/Anticyclonic_%Y%m%d.nc'), 'w') as h:
+            a.to_netcdf(h)
+        with Dataset(args.path_out + date.strftime('/Cyclonic_%Y%m%d.nc'), 'w') as h:
+            c.to_netcdf(h)