diff --git a/src/py_eddy_tracker/dataset/grid.py b/src/py_eddy_tracker/dataset/grid.py index 2bb5b70d..5b884b68 100644 --- a/src/py_eddy_tracker/dataset/grid.py +++ b/src/py_eddy_tracker/dataset/grid.py @@ -607,6 +607,7 @@ def eddy_identification( date, step=0.005, shape_error=55, + presampling_multiplier=10, sampling=50, sampling_method="visvalingam", pixel_limit=None, @@ -624,8 +625,10 @@ def eddy_identification( :param datetime.datetime date: Date to be stored in object to date data :param float,int step: Height between two layers in m :param float,int shape_error: Maximal error allowed for outermost contour in % + :param int presampling_multiplier: + Evenly oversample the initial number of points in the contour by nb_pts x presampling_multiplier to fit circles :param int sampling: Number of points to store contours and speed profile - :param str sampling_method: Method to resample, 'uniform' or 'visvalingam' + :param str sampling_method: Method to resample the stored contours, 'uniform' or 'visvalingam' :param (int,int),None pixel_limit: Min and max number of pixels inside the inner and the outermost contour to be considered as an eddy :param float,None precision: Truncate values at the defined precision in m @@ -849,42 +852,59 @@ def eddy_identification( obs.amplitude[:] = amp.amplitude obs.speed_average[:] = max_average_speed obs.num_point_e[:] = contour.lon.shape[0] - xy_e = resample(contour.lon, contour.lat, **out_sampling) - obs.contour_lon_e[:], obs.contour_lat_e[:] = xy_e obs.num_point_s[:] = speed_contour.lon.shape[0] - xy_s = resample( - speed_contour.lon, speed_contour.lat, **out_sampling - ) - obs.contour_lon_s[:], obs.contour_lat_s[:] = xy_s - # FIXME : we use a contour without resampling - # First, get position based on innermost contour - centlon_i, centlat_i, _, _ = _fit_circle_path( - create_vertice(inner_contour.lon, inner_contour.lat) + # Evenly resample contours with nb_pts = nb_pts_original x presampling_multiplier + xy_i = uniform_resample( + inner_contour.lon, + inner_contour.lat, + num_fac=presampling_multiplier + ) + xy_e = uniform_resample( + contour.lon, + contour.lat, + num_fac=presampling_multiplier, ) - # Second, get speed-based radius based on contour of max uavg + xy_s = uniform_resample( + speed_contour.lon, + speed_contour.lat, + num_fac=presampling_multiplier, + ) + + # First, get position of max SSH based on best fit circle with resampled innermost contour + centlon_i, centlat_i, _, _ = _fit_circle_path(create_vertice(*xy_i)) + obs.lon_max[:] = centlon_i + obs.lat_max[:] = centlat_i + + # Second, get speed-based radius, shape error, eddy center, area based on resampled contour of max uavg centlon_s, centlat_s, eddy_radius_s, aerr_s = _fit_circle_path( create_vertice(*xy_s) ) - # Compute again to use resampled contour - _, _, eddy_radius_e, aerr_e = _fit_circle_path( - create_vertice(*xy_e) - ) - obs.radius_s[:] = eddy_radius_s - obs.radius_e[:] = eddy_radius_e - obs.shape_error_e[:] = aerr_e obs.shape_error_s[:] = aerr_s obs.speed_area[:] = poly_area( *coordinates_to_local(*xy_s, lon0=centlon_s, lat0=centlat_s) ) + obs.lon[:] = centlon_s + obs.lat[:] = centlat_s + + # Third, compute effective radius, shape error, area from resampled effective contour + _, _, eddy_radius_e, aerr_e = _fit_circle_path( + create_vertice(*xy_e) + ) + obs.radius_e[:] = eddy_radius_e + obs.shape_error_e[:] = aerr_e obs.effective_area[:] = poly_area( *coordinates_to_local(*xy_e, lon0=centlon_s, lat0=centlat_s) ) - obs.lon[:] = centlon_s - obs.lat[:] = centlat_s - obs.lon_max[:] = centlon_i - obs.lat_max[:] = centlat_i + + # Finally, resample contours with output parameters + xy_e_f = resample(*xy_e, **out_sampling) + xy_s_f = resample(*xy_s, **out_sampling) + + obs.contour_lon_s[:], obs.contour_lat_s[:] = xy_s_f + obs.contour_lon_e[:], obs.contour_lat_e[:] = xy_e_f + if aerr > 99.9 or aerr_s > 99.9: logger.warning( "Strange shape at this step! shape_error : %f, %f", diff --git a/src/py_eddy_tracker/poly.py b/src/py_eddy_tracker/poly.py index 0f0271ee..56fb55e7 100644 --- a/src/py_eddy_tracker/poly.py +++ b/src/py_eddy_tracker/poly.py @@ -86,7 +86,7 @@ def poly_area_vertice(v): @njit(cache=True) def poly_area(x, y): """ - Must be call with local coordinates (in m, to get an area in m²). + Must be called with local coordinates (in m, to get an area in m²). :param array x: :param array y: