Last correction of identification code to solve x bounds problem, contour was mis-weigthed and test of poly which contain poly was inssuficient around x bounds

Author: AntSimi

README.md

@@ -40,7 +40,8 @@ for tracking.
 Loading grid
 ```python
 from py_eddy_tracker.dataset.grid import RegularGridDataset
-h = RegularGridDataset('share/nrt_global_allsat_phy_l4_20190223_20190226.nc', 'longitude', 'latitude')
+grid_name, lon_name, lat_name = 'share/nrt_global_allsat_phy_l4_20190223_20190226.nc', 'longitude', 'latitude'
+h = RegularGridDataset(grid_name, lon_name, lat_name)
 ```
 
 Plotting grid
@@ -51,13 +52,14 @@ ax = fig.add_axes([.02, .51, .9, .45])
 ax.set_title('ADT (m)')
 ax.set_ylim(-75, 75)
 ax.set_aspect('equal')
-m=ax.pcolormesh(h.x_bounds, h.y_bounds, h.grid('adt').T.copy(), vmin=-1, vmax=1, cmap='coolwarm')
+m = h.display(ax, name='adt', vmin=-1, vmax=1)
 ax.grid(True)
 plt.colorbar(m, cax=fig.add_axes([.94, .51, .01, .45]))
 ```
 
 Filtering
 ```python
+h = RegularGridDataset(grid_name, lon_name, lat_name)
 h.bessel_high_filter('adt', 500, order=3)
 ```
 
@@ -72,7 +74,7 @@ ax = fig.add_axes([.02, .02, .9, .45])
 ax.set_title('ADT Filtered (m)')
 ax.set_aspect('equal')
 ax.set_ylim(-75, 75)
-m=ax.pcolormesh(h.x_bounds, h.y_bounds, h.grid('adt').T, vmin=-.1, vmax=.1, cmap='coolwarm')
+m = h.display(ax, name='adt', vmin=-.1, vmax=.1)
 ax.grid(True)
 plt.colorbar(m, cax=fig.add_axes([.94, .02, .01, .45]))
 fig.savefig('share/png/filter.png')
@@ -83,7 +85,6 @@ fig.savefig('share/png/filter.png')
 #### Compute spectrum and spectrum ratio on some area ####
 Load data
 ```python
-grid_name, lon_name, lat_name = 'share/nrt_global_allsat_phy_l4_20190223_20190226.nc', 'longitude', 'latitude'
 raw = RegularGridDataset(grid_name, lon_name, lat_name)
 filtered = RegularGridDataset(grid_name, lon_name, lat_name)
 filtered.bessel_low_filter('adt', 150, order=3)
@@ -169,8 +170,8 @@ ax.set_title('Eddies detected -- Cyclonic(red) and Anticyclonic(blue)')
 ax.set_ylim(-75,75)
 ax.set_xlim(0,360)
 ax.set_aspect('equal')
-ax.plot(a.obs['contour_lon_s'].T, a.obs['contour_lat_s'].T, 'b', linewidth=.5)
-ax.plot(c.obs['contour_lon_s'].T, c.obs['contour_lat_s'].T, 'r', linewidth=.5)
+a.display(ax, color='b', linewidth=.5)
+c.display(ax, color='r', linewidth=.5)
 ax.grid()
 fig.savefig('share/png/eddies.png')
 ```

share/nrt_global_allsat_phy_l4_20190223_20190226.nc

@@ -34,25 +34,22 @@ def raw_resample(datas, fixed_size):
 
 @property
 def mean_coordinates(self):
-    return self.vertices.mean(axis=0)
-
-
-BasePath.mean_coordinates = mean_coordinates
+    # last coordinates == first
+    return self.vertices[1:].mean(axis=0)
 
 
 @property
 def lon(self):
     return self.vertices[:, 0]
 
 
-BasePath.lon = lon
-
-
 @property
 def lat(self):
     return self.vertices[:, 1]
 
 
+BasePath.mean_coordinates = mean_coordinates
+BasePath.lon = lon
 BasePath.lat = lat
 
 
@@ -89,7 +86,8 @@ def fit_circle_path(self):
 @njit(cache=True, fastmath=True)
 def _fit_circle_path(vertice):
     lons, lats = vertice[:, 0], vertice[:, 1]
-    lon0, lat0 = lons.mean(), lats.mean()
+    # last coordinates == first
+    lon0, lat0 = lons[1:].mean(), lats[1:].mean()
     c_x, c_y = coordinates_to_local(lons, lats, lon0, lat0)
     # Some time, edge is only a dot of few coordinates
     d_lon = lons.max() - lons.min()

share/png/eddies.png

@@ -309,7 +309,7 @@ def find_wrapcut_path_and_join(self, x0, x1):
                                 path_left.vertices[-1, 1] - path_right.vertices[0, 1]) < self.DELTA_PREC:
                             polys_to_pop_left.append(i_left)
                             path_right.vertices[:, 0] -= 360
-                            path_left.vertices = concatenate((path_left.vertices, path_right.vertices))
+                            path_left.vertices = concatenate((path_left.vertices, path_right.vertices[1:]))
                             path_left.vertices[-1] = path_left.vertices[0]
                             paths_solve.append(path_left)
                             paths_right.pop(i_right)

share/png/filter.png

@@ -128,18 +128,19 @@ def fit_circle(x_vec, y_vec):
     p_inon_x = empty(nb_elt)
     p_inon_y = empty(nb_elt)
 
-    x_mean = x_vec.mean()
-    y_mean = y_vec.mean()
+    # last coordinates == first
+    x_mean = x_vec[1:].mean()
+    y_mean = y_vec[1:].mean()
 
-    norme = (x_vec - x_mean) ** 2 + (y_vec - y_mean) ** 2
+    norme = (x_vec[1:] - x_mean) ** 2 + (y_vec[1:] - y_mean) ** 2
     norme_max = norme.max()
     scale = norme_max ** .5
 
     # Form matrix equation and solve it
     # Maybe put f4
-    datas = ones((nb_elt, 3))
-    datas[:, 0] = 2. * (x_vec - x_mean) / scale
-    datas[:, 1] = 2. * (y_vec - y_mean) / scale
+    datas = ones((nb_elt - 1, 3))
+    datas[:, 0] = 2. * (x_vec[1:] - x_mean) / scale
+    datas[:, 1] = 2. * (y_vec[1:] - y_mean) / scale
 
     (center_x, center_y, radius), residuals, rank, s = lstsq(datas, norme / norme_max)
 

src/py_eddy_tracker/dataset/grid.py

@@ -728,9 +728,13 @@ def write_netcdf(self, path='./', filename='%(path)s/%(sign_type)s.nc'):
     def set_global_attr_netcdf(self, h_nc):
         pass
 
-    def display(self, ax, ref=0, **kwargs):
-        ax.plot((self.obs['contour_lon_s'].T - ref) % 360 + ref, self.obs['contour_lat_s'].T, ** kwargs)
-        ax.plot((self.obs['contour_lon_e'].T - ref) % 360 + ref, self.obs['contour_lat_e'].T, linestyle='-.', ** kwargs)
+    def display(self, ax, ref=None, **kwargs):
+        if ref is None:
+            ax.plot(self.obs['contour_lon_s'].T, self.obs['contour_lat_s'].T, ** kwargs)
+            ax.plot(self.obs['contour_lon_e'].T, self.obs['contour_lat_e'].T, linestyle='-.', ** kwargs)
+        else:
+            ax.plot((self.obs['contour_lon_s'].T - ref) % 360 + ref, self.obs['contour_lat_s'].T, ** kwargs)
+            ax.plot((self.obs['contour_lon_e'].T - ref) % 360 + ref, self.obs['contour_lat_e'].T, linestyle='-.', ** kwargs)
 
 
 @njit(cache=True, fastmath=True)

src/py_eddy_tracker/eddy_feature.py

@@ -28,7 +28,7 @@ def poly_contain_poly(xy_poly_out, xy_poly_in):
     x_ref = xy_poly_out[0,0]
     d_x = abs(x[0] - x_ref)
     if d_x > 180:
-        x = (x - x_ref + 180) % 360 - 180
+        x = (x - x_ref + 180) % 360 + x_ref - 180
     for i_elt in prange(nb_elt):
         wn = winding_number_poly(x[i_elt], xy_poly_in[i_elt, 1], xy_poly_out)
         if wn == 0: