Eddy identification fail

[Manuel-Fundora]

Bug report

Bug summary

Code for reproduction

# Paste your code here
#
from py_eddy_tracker import start_logger
from matplotlib import pyplot as plt
from py_eddy_tracker.dataset.grid import RegularGridDataset
from datetime import datetime

start_logger().setLevel("DEBUG")  # Available options: ERROR, WARNING, INFO, DEBUG

grid_name, lon_name, lat_name = (
    "cmems_obs-sl_glo_phy-ssh_my_allsat-l4-duacs-2014-11-22.nc",
    "longitude",
    "latitude",
)
# %%
h = RegularGridDataset(grid_name, lon_name, lat_name)
h.bessel_high_filter("adt", 500, order=3)
date = datetime(2014, 11, 22)
a, c = h.eddy_identification(
    "adt",
    "ugos",
    "vgos",  # Variables used for identification
    date,  # Date of identification
    0.002,  # step between two isolines of detection (m)
    pixel_limit=(5, 2000),  # Min and max pixel count for valid contour
    shape_error=55,  # Error max (%) between ratio of circle fit and contour
)
# %%
fig = plt.figure(figsize=(15, 7))
ax = fig.add_axes([0.03, 0.03, 0.94, 0.94])
ax.set_title("Eddies detected -- Cyclonic(red) and Anticyclonic(blue)")
ax.set_ylim(15, 25)
ax.set_xlim(-115, -105)
ax.set_aspect("equal")
a.display(ax, color="b", linewidth=0.5)
c.display(ax, color="r", linewidth=0.5)
ax.grid()
fig.savefig("eddies.png")

#

Actual outcome

# If applicable, paste the console output here
#runfile('C:/Users/LENOVO/Downloads/untitled1.py', wdir='C:/Users/LENOVO/Downloads')
WARNING 2022-11-03 21:13:54,904 grid.__init__ :
    We assume pixel position of grid is centered for cmems_obs-sl_glo_phy-ssh_my_allsat-l4-duacs-2014-11-22.nc
    DEBUG 2022-11-03 21:13:54,906   grid.   load_general_features :
        Load general feature from cmems_obs-sl_glo_phy-ssh_my_allsat-l4-duacs-2014-11-22.nc
    DEBUG 2022-11-03 21:13:54,927   grid.   bessel_high_filter :
        Run filtering with wavelength of 500 km and order of 3 ...
    DEBUG 2022-11-03 21:13:54,927   grid.   grid :
        Load adt from cmems_obs-sl_glo_phy-ssh_my_allsat-l4-duacs-2014-11-22.nc
    DEBUG 2022-11-03 21:13:55,111   grid.   bessel_high_filter :
        Filtering done
    DEBUG 2022-11-03 21:13:55,112   grid.   grid :
        Load ugos from cmems_obs-sl_glo_phy-ssh_my_allsat-l4-duacs-2014-11-22.nc
    DEBUG 2022-11-03 21:13:55,114   grid.   grid :
        Load vgos from cmems_obs-sl_glo_phy-ssh_my_allsat-l4-duacs-2014-11-22.nc
Remain  0:00:00.002639 ETA  2022-11-03 21:13:55.112125 current kernel size : (119, 107) Step : 40/41    
INFO 2022-11-03 21:13:55,267 grid.eddy_identification :
    We will apply on step a factor to be coherent with grid : 1.000000
    DEBUG 2022-11-03 21:13:55,269   grid.   eddy_identification :
        Levels from nan to nan
C:\Users\LENOVO\anaconda3\lib\site-packages\py_eddy_tracker\dataset\grid.py:724: RuntimeWarning: invalid value encountered in double_scalars
  levels = arange(z_min - z_min % step, z_max - z_max % step + 2 * step, step)
Traceback (most recent call last):

  File ~\Downloads\untitled1.py:25 in <module>
    a, c = h.eddy_identification(

  File ~\anaconda3\lib\site-packages\py_eddy_tracker\dataset\grid.py:724 in eddy_identification
    levels = arange(z_min - z_min % step, z_max - z_max % step + 2 * step, step)

ValueError: arange: cannot compute length
#

Expected outcome this code should provide an image with the cyclonic and anticyclonic eddies in the region

Copernicus_Data.zip

[AntSimi]

Ok if you want apply filter with your parameter on grid you need a much larger grid.

• For global detection, i use this parameter for ADT filtering 700km with order 1
• Grid must be mask and not have nan, original grid of CMEMS don't have nan, did you get grid with MOTU?

[Manuel-Fundora]

I am working on a local scale, only the Mexican Pacific region. download my data according to this region. im not using MOTU, do you advise me to use it?

[AntSimi]

I don't advise to use MOTU. When you do a ncdump on your file there are several thing which are modified from original file: Your file :

ncdump -h cmems_obs-sl_glo_phy-ssh_my_allsat-l4-duacs-2014-11-22.nc
{
dimensions:
        time = 1 ;
        latitude = 41 ;
        longitude = 41 ;
        nv = 2 ;
variables:
        ...
        double ugos(time, latitude, longitude) ;
                ugos:coordinates = "longitude latitude" ;
                ugos:grid_mapping = "crs" ;
                ugos:long_name = "Absolute geostrophic velocity: zonal component" ;
                ugos:standard_name = "surface_geostrophic_eastward_sea_water_velocity" ;
                ugos:units = "m/s" ;
                ugos:_ChunkSizes = 1, 50, 50 ;

// global attributes:
                ...
                :time_coverage_duration = "P1D" ;
                :time_coverage_end = "2022-02-09T12:00:00Z" ;
                :time_coverage_resolution = "P1D" ;
                :time_coverage_start = "2022-02-08T12:00:00Z" ;
                :title = "DT merged all satellites Global Ocean Gridded SSALTO/DUACS Sea Surface Height L4 product and derived variables" ;
                :_CoordSysBuilder = "ucar.nc2.dataset.conv.CF1Convention" ;
}

Original one :

ncdump -h  dt_global_allsat_phy_l4_20141122_20210726.nc 
{
dimensions:
        time = 1 ;
        latitude = 720 ;
        longitude = 1440 ;
        nv = 2 ;
variables:
        ....
        int ugos(time, latitude, longitude) ;
                ugos:_FillValue = -2147483647 ;
                ugos:coordinates = "longitude latitude" ;
                ugos:grid_mapping = "crs" ;
                ugos:long_name = "Absolute geostrophic velocity: zonal component" ;
                ugos:scale_factor = 0.0001 ;
                ugos:standard_name = "surface_geostrophic_eastward_sea_water_velocity" ;
                ugos:units = "m/s" ;

// global attributes:
               ...
                :time_coverage_duration = "P1D" ;
                :time_coverage_end = "2014-11-22T12:00:00Z" ;
                :time_coverage_resolution = "P1D" ;
                :time_coverage_start = "2014-11-21T12:00:00Z" ;
                :title = "DT merged all satellites Global Ocean Gridded SSALTO/DUACS Sea Surface Height L4 product and derived variables" ;
}

[Manuel-Fundora]

Thank you very much. Can you specify where I can download the data (the one you indicate as original) or if you can send it to me. To start work on the script with it. Excuse the inconvenience, but it would be of great help to my work to obtain the results that you expose in your algorithms.

[AntSimi]

you could get data on CMEMS

[Manuel-Fundora]

Thanks, I already downloaded the data from the ftp. Now when I run the script it generates these images. Apparently there is something with the coordinates, with the latitude. I tried to transform it in the file but I couldn't, I tried to add 360 to it but it didn't work either. What could I do? I send you the resulting images.

[AntSimi]

Quick solution: if you want detetion of eddies on CMEMS maps you could use AVISO atlas which are computed with CMEMS maps. With py eddy tracker tools you could do a subset on your area.

[Manuel-Fundora]

Thanks, I already downloaded the data for the eddies from A.VISO (META3.2_DT_allsat_Cyclonic_long_19930101_20220209.nc) for both cyclonic and anticyclonic. This data with what section of the eddie tracker can I process it.

[AntSimi]

# To select on area from [200W:250W] and [20N:40N]
EddySubSetter my_atlas_in.nc atlas_out.nc --area 200 20 250 40
# For documentation
EddySubSetter -h

[Manuel-Fundora]

hello, I have already managed to make the spatial cut, but the temporary one is not done, this is the line that I am putting

EddySubSetter "META3.2_DT_allsat_Anticyclonic_untracked_19930101_20220209.nc" "cutted3.nc" \
     --area 245 15 255 25 -p 26306 26309

I already tried with the days in julians or starting from scratch but when I add this part the cut is not made, how could I make the cut by time too.

[AntSimi]

Time is stored with a scale/offset options so use this option --no_raw_mode to filter time

[Manuel-Fundora]

Thank you very much for your help and time

[Manuel-Fundora]

hello, how can i create the .zarr archives?, to makes the trajectories and the lifetimes.

[AntSimi]

When you save use .zarr extension instead of .nc. But netcdf could also store trajetories.

[Manuel-Fundora]

thanks, So in the eddy_tracker the scripts that call .zarr files can be replaced by the same .nc that I cut from the AVI.SO Atlas? Or when I cut with the EddySubsetter should I save the file as .zarr?

[AntSimi]

zarr or nc is only extension for storage, those 2 formats could be use like you want at each processings step

[Manuel-Fundora]

ok, how can i know the id for a specific eddy? for example: in the track animation script, how can i change the directory to my .nc file and the id number for extract? And knowing the id how can I know its independent properties, radius, area, EKE, OW parameter, for a single eddy i mean.

[AntSimi]

• Just give filename to load_file to load another dataset
• You need to know which ID you want to display
• my_atlas.extract_ids([my_ids, ....]).speed_radius will give speed radius for all given ID

[Manuel-Fundora]

Hi, I'm trying to track a particular eddy for my study region, how can I get the ids of those eddies involved. On the other hand, in the scripts corresponding to the tracking (pet_one_track, pet_track_anim, pet_track_anim_matplotlib_animation) when I use my database and change the id it does not make the video. Could you help me. Thank you very much, I am attaching my cropped area from A.VISO atlas. Anticyclonic_PTFM.zip Scripts.zip thanks for your time.

[AntSimi]

Please for a new topic open a new issue or discussion thread :).