43
43
from cv2 import filter2D
44
44
from numba import njit , types as numba_types
45
45
from matplotlib .path import Path as BasePath
46
- from pyproj import Proj
47
46
from pint import UnitRegistry
48
47
from ..observations .observation import EddiesObservations
49
48
from ..eddy_feature import Amplitude , Contours
@@ -701,28 +700,29 @@ def eddy_identification(
701
700
for contour in contour_paths :
702
701
if contour .used :
703
702
continue
704
- centlon_e , centlat_e , eddy_radius_e , aerr = contour .fit_circle ()
703
+ # FIXME : center could be not in contour and fit on raw sampling
704
+ _ , _ , _ , aerr = contour .fit_circle ()
705
705
706
706
# Filter for shape
707
707
if aerr < 0 or aerr > shape_error or isnan (aerr ):
708
708
continue
709
- # Get indices of centroid
710
- # Give only 1D array of lon and lat not 2D data
711
- i_x , i_y = self .nearest_grd_indice (centlon_e , centlat_e )
712
- i_x = self .normalize_x_indice (i_x )
713
-
714
- # Check if centroid is on define value
715
- if data .mask [i_x , i_y ]:
716
- continue
717
- # Test to know cyclone or anticyclone
718
- acyc_not_cyc = data [i_x , i_y ] >= cvalues
719
- if anticyclonic_search != acyc_not_cyc :
720
- continue
721
709
722
710
# Find all pixels in the contour
723
711
i_x_in , i_y_in = contour .pixels_in (self )
724
712
725
- # Maybe limit max must be replace with a maximum of surface
713
+ # Check if pixels in contour are masked
714
+ if data .mask [i_x_in , i_y_in ].any ():
715
+ continue
716
+
717
+ # Test to know cyclone or anticyclone
718
+ if anticyclonic_search :
719
+ if (data [i_x_in , i_y_in ] < cvalues ).any ():
720
+ continue
721
+ else :
722
+ if (data [i_x_in , i_y_in ] > cvalues ).any ():
723
+ continue
724
+
725
+ # FIXME : Maybe limit max must be replace with a maximum of surface
726
726
if (
727
727
contour .nb_pixel < pixel_limit [0 ]
728
728
or contour .nb_pixel > pixel_limit [1 ]
@@ -775,30 +775,7 @@ def eddy_identification(
775
775
pixel_min = pixel_limit [0 ],
776
776
)
777
777
778
- # Use azimuth equal projection for radius
779
- proj = Proj (
780
- "+proj=aeqd +ellps=WGS84 +lat_0={1} +lon_0={0}" .format (
781
- * inner_contour .mean_coordinates
782
- )
783
- )
784
- # First, get position based on innermost
785
- # contour
786
- centx_i , centy_i , _ , _ = fit_circle (
787
- * proj (inner_contour .lon , inner_contour .lat )
788
- )
789
- centlon_i , centlat_i = proj (centx_i , centy_i , inverse = True )
790
- # Second, get speed-based radius based on
791
- # contour of max uavg
792
- centx_s , centy_s , eddy_radius_s , aerr_s = fit_circle (
793
- * proj (speed_contour .lon , speed_contour .lat )
794
- )
795
- # Computed again to be coherent with speed_radius, we will be compute in same reference
796
- _ , _ , eddy_radius_e , aerr_e = fit_circle (
797
- * proj (contour .lon , contour .lat )
798
- )
799
- centlon_s , centlat_s = proj (centx_s , centy_s , inverse = True )
800
-
801
- # Instantiate new EddyObservation object (high cost need to be review)
778
+ # FIXME : Instantiate new EddyObservation object (high cost need to be review)
802
779
obs = EddiesObservations (
803
780
size = 1 ,
804
781
track_extra_variables = track_extra_variables ,
@@ -818,23 +795,38 @@ def eddy_identification(
818
795
else :
819
796
obs .obs ["uavg_profile" ] = raw_resample (speed_array , sampling )
820
797
obs .obs ["amplitude" ] = amp .amplitude
821
- obs .obs ["radius_s" ] = eddy_radius_s
822
798
obs .obs ["speed_average" ] = max_average_speed
799
+ obs .obs ["num_point_e" ] = contour .lon .shape [0 ]
800
+ xy_e = uniform_resample (contour .lon , contour .lat , ** out_sampling )
801
+ obs .obs ["contour_lon_e" ], obs .obs ["contour_lat_e" ] = xy_e
802
+ obs .obs ["num_point_s" ] = speed_contour .lon .shape [0 ]
803
+ xy_s = uniform_resample (
804
+ speed_contour .lon , speed_contour .lat , ** out_sampling
805
+ )
806
+ obs .obs ["contour_lon_s" ], obs .obs ["contour_lat_s" ] = xy_s
807
+
808
+ # FIXME : we use a contour without resampling
809
+ # First, get position based on innermost contour
810
+ centlon_i , centlat_i , _ , _ = _fit_circle_path (
811
+ create_vertice (inner_contour .lon , inner_contour .lat )
812
+ )
813
+ # Second, get speed-based radius based on contour of max uavg
814
+ centlon_s , centlat_s , eddy_radius_s , aerr_s = _fit_circle_path (
815
+ create_vertice (* xy_s )
816
+ )
817
+ # Computed again to use resample contour
818
+ _ , _ , eddy_radius_e , aerr_e = _fit_circle_path (
819
+ create_vertice (* xy_e )
820
+ )
821
+
822
+ obs .obs ["radius_s" ] = eddy_radius_s
823
823
obs .obs ["radius_e" ] = eddy_radius_e
824
824
obs .obs ["shape_error_e" ] = aerr_e
825
825
obs .obs ["shape_error_s" ] = aerr_s
826
826
obs .obs ["lon" ] = centlon_s
827
827
obs .obs ["lat" ] = centlat_s
828
828
obs .obs ["lon_max" ] = centlon_i
829
829
obs .obs ["lat_max" ] = centlat_i
830
- obs .obs ["num_point_e" ] = contour .lon .shape [0 ]
831
- xy = uniform_resample (contour .lon , contour .lat , ** out_sampling )
832
- obs .obs ["contour_lon_e" ], obs .obs ["contour_lat_e" ] = xy
833
- obs .obs ["num_point_s" ] = speed_contour .lon .shape [0 ]
834
- xy = uniform_resample (
835
- speed_contour .lon , speed_contour .lat , ** out_sampling
836
- )
837
- obs .obs ["contour_lon_s" ], obs .obs ["contour_lat_s" ] = xy
838
830
if aerr > 99.9 or aerr_s > 99.9 :
839
831
logger .warning (
840
832
"Strange shape at this step! shape_error : %f, %f" ,
0 commit comments