Error: cannot compute length

[santoine2710]

I am downloading my data from here;
https://data.marine.copernicus.eu/product/SEALEVEL_EUR_PHY_L4_NRT_OBSERVATIONS_008_060/download
I am able to plot the grid and run the bessel_high_filter (only over 30km, not higher).

But everytime I try to run the eddy_identification function I get the following error;

levels = arange(z_min - z_min % step, z_max - z_max % step + 2 * step, step)
Traceback (most recent call last):
  File "C:\Users\Sarah\Documents\pyeddytracker\Ibma_csv.py", line 65, in <module>
    a, c = g.eddy_identification(
  File "C:\ProgramData\Anaconda3\lib\site-packages\py_eddy_tracker\dataset\grid.py", line 724, in eddy_identification
    levels = arange(z_min - z_min % step, z_max - z_max % step + 2 * step, step)
ValueError: arange: cannot compute length

I played around with some parameters but nothing changed.
Do you have an idea why this dataset is not working?

Thank you,
Sarah

[AntSimi]

Could you share script use to run eddy identification

[santoine2710]

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
from py_eddy_tracker.observations.observation import EddiesObservations
from py_eddy_tracker import data
from py_eddy_tracker.eddy_feature import Contours

##data source https://data.marine.copernicus.eu/product/SEALEVEL_EUR_PHY_L4_NRT_OBSERVATIONS_008_060/download
path = r"C:\Users\Sarah\OneDrive - Universidade do Algarve\Documentos\EMSO\2022_Ibma_csv_Paper\altimetry\1406.nc"

start_logger().setLevel("DEBUG")

##load grid
grid_name, lon_name, lat_name = (
    path,
    "longitude",
    "latitude",
)
h = RegularGridDataset(grid_name, lon_name, lat_name)

def start_axes(title):
    fig = plt.figure(figsize=(13, 5))
    ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])
    ax.set_aspect("equal")
    ax.set_title(title, weight="bold")
    return ax


def update_axes(ax, mappable=None):
    ax.grid()
    if mappable:
        plt.colorbar(mappable, cax=ax.figure.add_axes([0.94, 0.05, 0.01, 0.9]))

ax = start_axes("ADT (m)")
m = h.display(ax, "adt", cmap="RdBu_r")
update_axes(ax, m)
plt.show()

ax = start_axes("U/V deduce from ADT (m)")
m = h.display(ax, "adt", cmap="RdBu_r",vmin=0.02, vmax=0.23)
ax.set_xlim(-11, -7), ax.set_ylim(35,38)
u, v = h.grid("ugos").T, h.grid("vgos").T
ax.quiver(h.x_c, h.y_c, u, v, scale=10)
update_axes(ax, m)
plt.show()

h.bessel_high_filter("adt",30, order=3)
date = datetime(2022, 6, 17)
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
)

[AntSimi]

Could you share ncdump -h of 1406.nc?

[santoine2710]

dimensions:
        time = 1 ;
        latitude = 73 ;
        longitude = 91 ;
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 ;
        double adt(time, latitude, longitude) ;
                adt: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" ;
                adt:coordinates = "longitude latitude" ;
                adt:grid_mapping = "crs" ;
                adt:long_name = "Absolute dynamic topography" ;
                adt:standard_name = "sea_surface_height_above_geoid" ;
                adt:units = "m" ;
                adt:_ChunkSizes = 1, 50, 50 ;
        double vgos(time, latitude, longitude) ;
                vgos:coordinates = "longitude latitude" ;
                vgos:grid_mapping = "crs" ;
                vgos:long_name = "Absolute geostrophic velocity: meridian component" ;
                vgos:standard_name = "surface_geostrophic_northward_sea_water_velocity" ;
                vgos:units = "m/s" ;
                vgos:_ChunkSizes = 1, 50, 50 ;
        int crs ;
                crs:comment = "This is a container variable that describes the grid_mapping used by the data in this file. This variable does not contain any data; only information about the geographic coordinate system." ;
                crs:inverse_flattening = 298.257 ;
                crs:grid_mapping_name = "latitude_longitude" ;
                crs:semi_major_axis = 6378136.3 ;
                crs:_CoordinateTransformType = "Projection" ;
                crs:_CoordinateAxisTypes = "GeoX GeoY" ;
        float latitude(latitude) ;
                latitude:axis = "Y" ;
                latitude:bounds = "lat_bnds" ;
                latitude:long_name = "Latitude" ;
                latitude:standard_name = "latitude" ;
                latitude:units = "degrees_north" ;
                latitude:valid_max = 40.3125f ;
                latitude:valid_min = 31.3125f ;
                latitude:_ChunkSizes = 50 ;
                latitude:_CoordinateAxisType = "Lat" ;
        double sla(time, latitude, longitude) ;
                sla:ancillary_variables = "err_sla" ;
                sla:comment = "The sea level anomaly is the sea surface height above mean sea surface; it is referenced to the [1993, 2012] period; see the product user manual for details" ;
                sla:coordinates = "longitude latitude" ;
                sla:grid_mapping = "crs" ;
                sla:long_name = "Sea level anomaly" ;
                sla:standard_name = "sea_surface_height_above_sea_level" ;
                sla:units = "m" ;
                sla:_ChunkSizes = 1, 50, 50 ;
        float time(time) ;
                time:axis = "T" ;
                time:calendar = "gregorian" ;
                time:long_name = "Time" ;
                time:standard_name = "time" ;
                time:units = "days since 1950-01-01 00:00:00" ;
                time:_ChunkSizes = 1 ;
                time:_CoordinateAxisType = "Time" ;
                time:valid_min = 26462.f ;
                time:valid_max = 26462.f ;
        float longitude(longitude) ;
                longitude:axis = "X" ;
                longitude:bounds = "lon_bnds" ;
                longitude:long_name = "Longitude" ;
                longitude:standard_name = "longitude" ;
                longitude:units = "degrees_east" ;
                longitude:valid_max = -3.8125f ;
                longitude:valid_min = -15.0625f ;
                longitude:_ChunkSizes = 50 ;
                longitude:_CoordinateAxisType = "Lon" ;

// global attributes:
                :Conventions = "CF-1.6" ;
                :FROM_ORIGINAL_FILE__Metadata_Conventions = "Unidata Dataset Discovery v1.0" ;
                :cdm_data_type = "Grid" ;
                :comment = "Sea Surface Height measured by Altimetry and derived variables" ;
                :contact = "[email protected]" ;
                :creator_email = "[email protected]" ;
                :creator_name = "CMEMS - Sea Level Thematic Assembly Center" ;
                :creator_url = "http://marine.copernicus.eu" ;
                :date_created = "2023-06-20T02:03:15Z" ;
                :date_issued = "2023-06-20T02:03:15Z" ;
                :date_modified = "2023-06-20T02:03:15Z" ;
                :FROM_ORIGINAL_FILE__geospatial_lat_max = 66.0625 ;
                :FROM_ORIGINAL_FILE__geospatial_lat_min = 19.9375 ;
                :FROM_ORIGINAL_FILE__geospatial_lat_resolution = 0.125 ;
                :FROM_ORIGINAL_FILE__geospatial_lat_units = "degrees_north" ;
                :FROM_ORIGINAL_FILE__geospatial_lon_max = 42.0625 ;
                :FROM_ORIGINAL_FILE__geospatial_lon_min = -30.0625 ;
                :FROM_ORIGINAL_FILE__geospatial_lon_resolution = 0.125 ;
                :FROM_ORIGINAL_FILE__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 = "2023-06-20 02:03:15Z: 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" ;
                :FROM_ORIGINAL_FILE__platform = "Altika, Cryosat-2 New Orbit, Haiyang-2B, Jason-3 Interleaved, Sentinel-3A, Sentinel-3B, Sentinel-6A" ;
                :processing_level = "L4" ;
                :FROM_ORIGINAL_FILE__product_version = "vDec2021" ;
                :project = "COPERNICUS MARINE ENVIRONMENT MONITORING SERVICE (CMEMS)" ;
                :references = "http://marine.copernicus.eu" ;
                :FROM_ORIGINAL_FILE__software_version = "19.2.0_DUACS_DT2021_baseline" ;
                :source = "Altimetry measurements" ;
                :ssalto_duacs_comment = "Sentinel-6A 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 European Seas." ;
                :time_coverage_duration = "P1D" ;
                :time_coverage_end = "2023-06-20T12:00:00Z" ;
                :time_coverage_resolution = "P1D" ;
                :time_coverage_start = "2023-06-19T12:00:00Z" ;
                :title = "NRT merged all satellites European Seas Gridded SSALTO/DUACS Sea Surface Height L4 product and derived variables" ;
                :_CoordSysBuilder = "ucar.nc2.dataset.conv.CF1Convention" ;

[AntSimi]

When you run bessel_high_filter you can't set wavelength higher than 30 km? Could you share python traceback throw when your run with higher value?
This step is necessary to remove large scale pattern which could produce weak detection.
I think you do a subset on cmems grid, but to produce a good filtering it must be important to keep a larger area.
Maybe your issue is close to #173.