Request for help: No extrema found in contour of xxx pixels in level xxx

[cdmpbp123]

Dear Antoine, Thank you for providing the eddy detecting program. I am trying to detect eddies from a model result of NorthWest Pacific with grid of 2128*2129. I found a problem that several eddies were missing when running eddy identification module with many warnings of No extrema found in contour of xxx pixels in level xxx. The log file is: log_NWP.log

The version of python and other library: Python 3.10.13 pyEddyTracker # Version: 3.6.1 numpy# Version: 1.22.4 numba# Version: 0.55.2 matplotlib# Version: 3.7.0 xarray# Version: 2023.10.1 pandas# Version: 2.1.2 Here is the example data:NWP_20220101_sla.zip The test code is here: `

from py_eddy_tracker import start_logger   
start_logger().setLevel('DEBUG')   
from datetime import datetime  
from matplotlib import pyplot as plt  
from matplotlib.path import Path  
from py_eddy_tracker import data  

from py_eddy_tracker.dataset.grid import RegularGridDataset
from py_eddy_tracker.poly import create_vertice
from py_eddy_tracker.eddy_feature import Contours
from netCDF4 import Dataset

import xarray as xr  
import pandas as pd
import numpy as np

global lat_s,lat_n,lon_w,lon_e,domain_name
domain_name = 'NWP'
lat_s = -17
lat_n = 48
lon_w = 99
lon_e = 170

def start_axes(title):
    fig = plt.figure(figsize=(18, 18))
    ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])
    ax.set_xlim(lon_w, lon_e), ax.set_ylim(lat_s, lat_n)
    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]))
if __name__ == '__main__':
    fn = 'NWP_20220101_sla.nc'
    nf = Dataset(fn,'r')
    lat0 = nf.variables['latitude'][:].data
    lon0 = nf.variables['longitude'][:].data
    time = nf.variables['time'][:].data
    lat_idx = np.where((lat0>=lat_s) & (lat0<=lat_n))[0].tolist()
    lon_idx = np.where((lon0>=lon_w) & (lon0<=lon_e))[0].tolist()
    lat = lat0[lat_idx]
    lon = lon0[lon_idx]
    nf.close()
    da_time = xr.load_dataset(fn).time[0].values
    date_str = pd.to_datetime(str(da_time)).strftime('%Y%m%d')
    date = datetime.strptime(date_str, '%Y%m%d')
    margin = 0
    h = RegularGridDataset(fn, "longitude", "latitude",
        indexs=dict(
            longitude=slice(np.min(lon_idx).astype(int) - margin, np.max(lon_idx).astype(int) + margin),
            latitude=slice(np.min(lat_idx).astype(int) - margin, np.max(lat_idx).astype(int) + margin),
            ),   
    )

    h.add_uv("sla", "ugosa", "vgosa")
    a, c = h.eddy_identification(
    'sla', 'ugosa', 'vgosa', # Variables used for identification
    date, # Date of identification
    0.002, # step between two isolines of detection (m)
    pixel_limit=(20, 10000), # Min and max pixel count for valid contour
    )

    kwargs_a_adt = dict(lw=0.5, label="Anticyclonic ADT ({nb_obs} eddies)", ref=-10, color="b")
    kwargs_c_adt = dict(lw=0.5, label="Cyclonic ADT ({nb_obs} eddies)", ref=-10, color="r")
    ax = start_axes(f"SLA (m): {date_str}")
    m = h.display(ax, "sla", cmap="RdBu_r")
    a.display(ax, **kwargs_a_adt)
    c.display(ax, **kwargs_c_adt)
    ax.legend()
    update_axes(ax, m)
    plt.savefig('NWP_test_eddy_detect.png',facecolor ="w",dpi=600)
    plt.close()   

    ax = start_axes(f"SLA (m): {date_str}")
    m = h.display(ax, "sla", cmap="RdBu_r")
    ax.legend()
    great_current = Contours(h.x_c, h.y_c, h.grid("sla"), levels=(0.3,), keep_unclose=True)
    great_current.display(ax, color="k")
    update_axes(ax, m)
    plt.savefig('NWP_test_ssh_contour.png',facecolor ="w",dpi=600)
    plt.close()  

`

I add coutour of sla=0.3 to this map, it shows that many eddies are not identified. However, when I change the area to South China Sea while using same data and parameter. The eddies in the SCS were identified correctly. `

domain_name = 'SCS'
lat_s = 0
lat_n = 25
lon_w = 99
lon_e = 130

`

I have reviewed previous issues, but still can not solve this problem. Thank you for your time. Best wishes.

[AntSimi]

Several things come to me, you need to apply a large scale filter to remove large scale structure like in this example. Input grid have which resolution?

[cdmpbp123]

Thanks for your response and suggestion. I have tried different scale filters with 400, 700 and 1000km. The results with 1000km seems better than without filter, but most eddies detected in small region still can not be detected. And warning 'No extrema found in contour of xxx pixels in level xxx) still existed. Could you briefly explain why this warning happen? Is this because sla contour of my data so complicated that the algorithm can not find closed contour?
The input grid I used is 0.03°. How should high-pass filter value be determined, and could you provide a recommendation for my case? Thanks a lot!

[AntSimi]

My first guess. Contour detection algorithm could be high sensitive to high frequency data or noise. You could try to apply a low pass filter to smoth/simplify a little bit.

[cdmpbp123]

Thanks for your suggestion. I've tried several filters but still can't get satisfied results. I then divide the domain to several pieces and do the detection then concatenate them. The results seem OK now. Maybe the reason is the data quality, I will check the data later. Thanks for your help!

[AntSimi]

If you divide domain, result are different? i don't like this algotrithm behaviour., There is something to dig.

[cdmpbp123]

Yes. I test without any filter. Results vary significantly depending on the region dividing.