Eddy detection for POP3

[zmbvsky]

Hi,
I'm a Master's Student at the University of Hamburg and I'd like to apply the eddy detection and also the tracking to a dataset derived from the Parallel Ocean Program (POP).
The installation went well and the detection is working for the given example file (nrt_global_allsat_phy_l4_20190223_20190226.nc).
I tried to put all my data into a NetCDF with similar format as shown in the example.

 Source:
           20010101format.nc
Format:
           netcdf4_classic
Dimensions:
           time      = 1
           latitude  = 1200
           longitude = 3600
Variables:
    time     
           Size:       1x1
           Dimensions: time
           Datatype:   single
    latitude 
           Size:       1200x1
           Dimensions: latitude
           Datatype:   single
    longitude
           Size:       3600x1
           Dimensions: longitude
           Datatype:   single
    SSH      
           Size:       3600x1200x1
           Dimensions: longitude,latitude,time
           Datatype:   single
    UVEL     
           Size:       3600x1200x1
           Dimensions: longitude,latitude,time
           Datatype:   single
    VVEL     
            Size:       3600x1200x1
            Dimensions: longitude,latitude,time
            Datatype:   single

When I try to run the eddy detection on this dataset, I get the following error message:

  C:\Users\felix\anaconda3>EddyId share/20010101format.nc 20010101 SSH UVEL VVEL longitude latitude ./data_out -v INFO
 ←[30;47mWARNING←[0m 2022-02-21 15:38:37,567 ←[30;47mgrid←[0m.←[30;47m__init__←[0m :
        ←[30;47mWe assume pixel position of grid is centered for share/20010101format.nc←[0m
←[36mINFO←[0m 2022-02-21 15:38:37,571 ←[36mgrid←[0m.←[36msetup_coordinates←[0m :
        ←[36mGrid center←[0m
Traceback (most recent call last):
  File "c:\users\felix\anaconda3\lib\runpy.py", line 194, in _run_module_as_main
    return _run_code(code, main_globals, None,
  File "c:\users\felix\anaconda3\lib\runpy.py", line 87, in _run_code
    exec(code, run_globals)
  File "C:\Users\felix\anaconda3\Scripts\EddyId.exe\__main__.py", line 7, in <module>
  File "c:\users\felix\anaconda3\lib\site-packages\py_eddy_tracker\appli\grid.py", line 135, in eddy_id
    a, c = identification(
  File "c:\users\felix\anaconda3\lib\site-packages\py_eddy_tracker\appli\grid.py", line 185, in identification
    grid.bessel_high_filter(h, cut_wavelength, **kw_filter)
  File "c:\users\felix\anaconda3\lib\site-packages\py_eddy_tracker\dataset\grid.py", line 1496, in bessel_high_filter
    data_out = self.convolve_filter_with_dynamic_kernel(
  File "c:\users\felix\anaconda3\lib\site-packages\py_eddy_tracker\dataset\grid.py", line 1383, in convolve_filter_with_dynamic_kernel
    kernel = kernel_func(lat, **kwargs_func)
  File "c:\users\felix\anaconda3\lib\site-packages\py_eddy_tracker\dataset\grid.py", line 1342, in kernel_bessel
    kernel[0, half_y_pt * order] = pi
IndexError: index 0 is out of bounds for axis 0 with size 0

Since I'm not very familiar with Python, it's very hard for me to understand the problem in this case and I would be very happy about any help or suggestions.
The only difference I see between my file and the given example are the missing longitude and latitude bounds. Are they required for the eddy detection?

I hope I provided all necessary informations and would be very happy about help, solving this problem.
Cheers
Felix

[zmbvsky]

Hey,
thanks for editing my post, this looks much better.
I added the latitude and longitude bounds and tried it again, but I'm still getting the same error as before. I also changed the datatype for some variables.

Source:
           C:\Users\felix\anaconda3\share\20010101format.nc
Format:
           netcdf4_classic
Dimensions:
           time      = 1
           latitude  = 1200
           longitude = 3600
           nv        = 2
Variables:
    time     
           Size:       1x1
           Dimensions: time
           Datatype:   single
    latitude 
           Size:       1200x1
           Dimensions: latitude
           Datatype:   single
    longitude
           Size:       3600x1
           Dimensions: longitude
           Datatype:   single
    SSH      
           Size:       3600x1200x1
           Dimensions: longitude,latitude,time
           Datatype:   int32
    UVEL     
           Size:       3600x1200x1
           Dimensions: longitude,latitude,time
           Datatype:   int32
    VVEL     
           Size:       3600x1200x1
           Dimensions: longitude,latitude,time
           Datatype:   int32
    lat_bnds 
           Size:       2x1200
           Dimensions: nv,latitude
           Datatype:   single
    lon_bnds 
           Size:       2x3600
           Dimensions: nv,longitude
           Datatype:   single

Below you can find the working example file.

Source:
           C:\Users\felix\anaconda3\share\nrt_global_allsat_phy_l4_20190223_20190226.nc
Format:
           netcdf4
Global Attributes:
           Conventions                    = 'CF-1.6'
           Metadata_Conventions           = 'Unidata Dataset Discovery v1.0'
           cdm_data_type                  = 'Grid'
           comment                        = 'Sea Surface Height measured by Altimetry and derived variables'
           contact                        = 'servicedesk.cmems@mercator-ocean.eu'
           creator_email                  = 'servicedesk.cmems@mercator-ocean.eu'
           creator_name                   = 'CMEMS - Sea Level Thematic Assembly Center'
           creator_url                    = 'http://marine.copernicus.eu'
           date_created                   = '2019-02-25T23:27:16Z'
           date_issued                    = '2019-02-25T23:27:16Z'
           date_modified                  = '2019-02-25T23:27:16Z'
           geospatial_lat_max             = 89.875
           geospatial_lat_min             = -89.875
           geospatial_lat_resolution      = 0.25
           geospatial_lat_units           = 'degrees_north'
           geospatial_lon_max             = 359.875
           geospatial_lon_min             = 0.125
           geospatial_lon_resolution      = 0.25
           geospatial_lon_units           = 'degrees_east'
           geospatial_vertical_max        = 0
           geospatial_vertical_min        = 0
           geospatial_vertical_positive   = 'down'
           geospatial_vertical_resolution = 'point'
           geospatial_vertical_units      = 'm'
           history                        = 'Thu Apr 11 17:31:18 2019: ncks -v longitude,latitude,adt,ugos,vgos -4 -L4 nrt_global_allsat_phy_l4_20190223_20190226.nc nrt_global_allsat_phy_l4_20190223_20190226_.nc
                                            2019-02-25 23:27:18Z: Creation'
           institution                    = 'CLS, CNES'
           keywords                       = 'Oceans > Ocean Topography > Sea Surface Height'
           keywords_vocabulary            = 'NetCDF COARDS Climate and Forecast Standard Names'
           license                        = 'http://marine.copernicus.eu/web/27-service-commitments-and-licence.php'
           platform                       = 'Altika, Cryosat-2, Jason-3, OSTM/Jason-2 Long Repeat Orbit, Sentinel-3A'
           processing_level               = 'L4'
           product_version                = '4.0.0'
           project                        = 'COPERNICUS MARINE ENVIRONMENT MONITORING SERVICE (CMEMS)'
           references                     = 'http://marine.copernicus.eu'
           software_version               = '18.2.0'
           source                         = 'Altimetry measurements'
           ssalto_duacs_comment           = 'Jason-3 is the reference mission used for the altimeter inter-calibration processing'
           standard_name_vocabulary       = 'NetCDF Climate and Forecast (CF) Metadata Convention Standard Name Table v37'
           summary                        = 'SSALTO/DUACS Near-Real-Time Level-4 sea surface height and derived variables measured by multi-satellite altimetry observations over Global Ocean.'
           time_coverage_duration         = 'P1D'
           time_coverage_end              = '2019-02-23T00:00:00Z'
           time_coverage_resolution       = 'P1D'
           time_coverage_start            = '2019-02-23T00:00:00Z'
           title                          = 'NRT merged all satellites Global Ocean Gridded SSALTO/DUACS Sea Surface Height L4 product and derived variables'
           NCO                            = '"4.5.4"'
Dimensions:
           time      = 1
           latitude  = 720
           longitude = 1440
           nv        = 2
Variables:
    adt      
           Size:       1440x720x1
           Dimensions: longitude,latitude,time
           Datatype:   int32
           Attributes:
                       _FillValue    = -2147483647
                       comment       = 'The absolute dynamic topography is the sea surface height above geoid; the adt is obtained as follows: adt=sla+mdt where mdt is the mean dynamic topography; see the product user manual for details'
                       coordinates   = 'longitude latitude'
                       grid_mapping  = 'crs'
                       long_name     = 'Absolute dynamic topography'
                       scale_factor  = 0.0001
                       standard_name = 'sea_surface_height_above_geoid'
                       units         = 'm'
    lat_bnds 
           Size:       2x720
           Dimensions: nv,latitude
           Datatype:   single
           Attributes:
                       comment = 'latitude values at the north and south bounds of each pixel.'
                       units   = 'degrees_north'
    latitude 
           Size:       720x1
           Dimensions: latitude
           Datatype:   single
           Attributes:
                       axis          = 'Y'
                       bounds        = 'lat_bnds'
                       long_name     = 'Latitude'
                       standard_name = 'latitude'
                       units         = 'degrees_north'
                       valid_max     = 89.875
                       valid_min     = -89.875
    lon_bnds 
           Size:       2x1440
           Dimensions: nv,longitude
           Datatype:   single
           Attributes:
                       comment = 'longitude values at the west and east bounds of each pixel.'
                       units   = 'degrees_east'
    longitude
           Size:       1440x1
           Dimensions: longitude
           Datatype:   single
           Attributes:
                       axis          = 'X'
                       bounds        = 'lon_bnds'
                       long_name     = 'Longitude'
                       standard_name = 'longitude'
                       units         = 'degrees_east'
                       valid_max     = 359.875
                       valid_min     = 0.125
    time     
           Size:       1x1
           Dimensions: time
           Datatype:   single
           Attributes:
                       axis          = 'T'
                       calendar      = 'gregorian'
                       long_name     = 'Time'
                       standard_name = 'time'
                       units         = 'days since 1950-01-01 00:00:00'
    ugos     
           Size:       1440x720x1
           Dimensions: longitude,latitude,time
           Datatype:   int32
           Attributes:
                       _FillValue    = -2147483647
                       coordinates   = 'longitude latitude'
                       grid_mapping  = 'crs'
                       long_name     = 'Absolute geostrophic velocity: zonal component'
                       scale_factor  = 0.0001
                       standard_name = 'surface_geostrophic_eastward_sea_water_velocity'
                       units         = 'm/s'
    vgos     
           Size:       1440x720x1
           Dimensions: longitude,latitude,time
           Datatype:   int32
           Attributes:
                       _FillValue    = -2147483647
                       coordinates   = 'longitude latitude'
                       grid_mapping  = 'crs'
                       long_name     = 'Absolute geostrophic velocity: meridian component'
                       scale_factor  = 0.0001
                       standard_name = 'surface_geostrophic_northward_sea_water_velocity'
                       units         = 'm/s'

Except for the "Global Attributes" and "Attributes" specified for each variable the two files share the same format. Could this be the problem?
I will try to get deeper into the code to find the problem but I'm still happy about any suggestions that could lead to a solution.
Right now, I don't understand what axis 0 is ment to be and why it is empty.
Cheers
Felix

[AntSimi]

Several things:

• You must specify unit for variable
• Could you give min and max of latitude
• How was the original format of your data
• SSH, UVEL, VVEL use int32 without scale_factor is it correct?

[zmbvsky]

Oh, I thought the attributes are only for documentation and not used by the algorithm.

• Min(latitude)=-78.472862, Max(latitude)=1.7497509 since POP is using a tripole grid I want to focus on the southern hemisphere where the grid is regular.
• The original data comes from several NetCDFs with the variables given in double prec. *e: right now Im extracting the data from the original NetCDFs and creating a new file with Matlab to match the given format.
• You are right, I'm missing the scale_factor here

I will add units and scale_factor to the attributes and try again. *e: I guess specifying the axis, bounds, valid_max, valid_min and _FillValue is also necessary and I will add those attributes aswell.
Thank you a lot.

[AntSimi]

only following attribut is mandatory : units
scale_factor/add_offset and _FillValue must be necessary in specific case

[zmbvsky]

Ok, thank you a lot for you time and help.

[zmbvsky]

Hi,
I added units, scale_factor and the according _FillValues to my data but I'm still getting the same error.

Source:
           C:\Users\felix\anaconda3\share\20010101format.nc
Format:
           netcdf4
Dimensions:
           time      = 1
           latitude  = 1200
           longitude = 3600
           nv        = 2
Variables:
    time     
           Size:       1x1
           Dimensions: time
           Datatype:   single
           Attributes:
                       units = 'days since 1950-01-01 00:00:00'
    latitude 
           Size:       1200x1
           Dimensions: latitude
           Datatype:   single
           Attributes:
                       units  = 'degrees_north'
                       bounds = 'lat_bnds'
    longitude
           Size:       3600x1
           Dimensions: longitude
           Datatype:   single
           Attributes:
                       units  = 'degrees_east'
                       bounds = 'lon_bnds'
    SSH      
           Size:       3600x1200x1
           Dimensions: longitude,latitude,time
           Datatype:   int32
           Attributes:
                       _FillValue   = -9999
                       units        = 'm'
                       scale_factor = 0.0001
    UVEL     
           Size:       3600x1200x1
           Dimensions: longitude,latitude,time
           Datatype:   int32
           Attributes:
                       _FillValue   = 0
                       units        = 'm/s'
                       scale_factor = 0.0001
    VVEL     
           Size:       3600x1200x1
           Dimensions: longitude,latitude,time
           Datatype:   int32
           Attributes:
                       _FillValue   = 0
                       units        = 'm/s'
                       scale_factor = 0.0001
    lat_bnds 
           Size:       2x1200
           Dimensions: nv,latitude
           Datatype:   single
           Attributes:
                       units = 'degrees_north'
    lon_bnds 
           Size:       2x3600
           Dimensions: nv,longitude
           Datatype:   single
           Attributes:
                       units = 'degrees_east'

Do you have any suggestion what I can do to rule out what's going wrong?
I thought about filling up the "working example" with my own data, to see if it has something to do with my data itself.

[AntSimi]

Could you share netcdf file use in input of EddyDetection?

[zmbvsky]

Hey,
I had to reduce the dataset from 1200x3600 (latitude x longitude) to 1200x1200 grid points to make it smaller than 25mb.
The given error is still the same.
20010101format.zip

[AntSimi]

Error message was not obvious, maybe i will add some code for warning.
But coordinates must be strictly increasing. Your longitude variable is incorrect.
Your SSH variable have an incorrect fillvalue, this code must return matrix without -9999....:

from netCDF4 import Dataset
with Dataset('20010101format.nc') as h:
    print(h.variables['SSH'][:])

output

[[[-9999.0 -9999.0 -9999.0 ... -9999.0 -9999.0 -9999.0]
  [-9999.0 -9999.0 -9999.0 ... -9999.0 -9999.0 -9999.0]
  [-9999.0 -9999.0 -9999.0 ... -9999.0 -9999.0 -9999.0]
  ...
  [0.12240000000000001 0.1228 0.1221 ... -9999.0 -9999.0 -9999.0]
  [0.1297 0.1288 0.1291 ... -9999.0 -9999.0 -9999.0]
  [0.1365 0.1364 0.1353 ... -9999.0 -9999.0 -9999.0]]]

Last advice for filtering step, 700 km is better than 500km and maybe use only order 1.

[zmbvsky]

Thank you very much for your help. I really apreciate it.
I will fix the issues and I'm really looking forward to study the results.