surface detection

[loic23physics]

Hi,

I'm running the code for surface detections using zeta variable, and I see that for a detected eddy, in this case a NBC ring, the contour size fluctuates a lot between the time steps. I checked this using a code to map SSH contours, and overlaying the detected contours. The parameters used here were : z_max = 0.5 m, z_min = -0.5 m and step = 0.0001. If I try a bigger step it will not detected the ring at all for a certain percentage of the timesteps, whereas if I make it smaller it doesn't improve the detection any further. I have tried with a range of zeta values from (0.1,-0.1) to (3,3), without success, as well as changing the error shape bounds.

Would you have an idea how this could be solved ?

Thanks,

Loïc

[AntSimi]

Could you explain how you run detection?

[loic23physics]

detection_gigatl6 (1).txt The detection is run using this code. It works with a 1h time step simulation, and uses the parameters set in the identification section

from datetime import datetime, timedelta
from py_eddy_tracker import start_logger
from py_eddy_tracker.dataset.grid import UnRegularGridDataset
from numpy import zeros, arange, float
from netCDF4 import Dataset
from matplotlib import pyplot as plt
import sys,os,shutil

class RomsDataset(UnRegularGridDataset):

    __slots__ = list()

    def init_speed_coef(self, uname, vname):
        # xi_u and eta_v must be specified because this dimension are not use in lon/lat
        u = self.grid(uname, indexs=self.indexs)
        v = self.grid(vname, indexs=self.indexs)
        u = self.rho_2d(u.T).T
        v = self.rho_2d(v)
        self._speed_norm = (v ** 2 + u ** 2) ** .5

    @staticmethod
    def rho_2d(x):
        """Transformation to have u or v on same grid than h
        """
        M, Lp = x.shape
        new_x = zeros((M + 1, Lp))
        new_x[1:-1] = .5 * (x[:-1] + x[1:])
        new_x[0] = new_x[1]
        new_x[-1] = new_x[-2]
        return new_x

    @staticmethod
    def psi2rho(var_psi):
        """Transformation to have vrt on same grid than h
        """
        M, L = var_psi.shape
        Mp=M+1; Lp=L+1
        Mm=M-1; Lm=L-1
        var_rho = zeros((Mp, Lp))
        var_rho[1:M,1:L]=0.25*(var_psi[0:Mm,0:Lm]+var_psi[0:Mm,1:L]+var_psi[1:M,0:Lm]+var_psi[1:M,1:L])
        var_rho[0,:]=var_rho[1,:]
        var_rho[Mp-1,:]=var_rho[M-1,:]
        var_rho[:,0]=var_rho[:,1]
        var_rho[:,Lp-1]=var_rho[:,L-1]
        return var_rho

varname = 'zeta'

if __name__ == '__main__':
    start_logger().setLevel('DEBUG')
    ####
    # create case-specific folder
    folder = varname #+ '/' + sys.argv[2] + '/'

    try:
        os.mkdir(folder)
    except OSError:
        print ("Directory %s already exists" % folder)

    # copy script in folder
    shutil.copy('detection_gigatl6.py',folder)

    # Pick a depth/isopycnal
    if varname == 'zeta':
        s_rho = -1
        depth = 0
    elif  varname == 'ow':
        #isopycnal
        #grid_name = './iso/gigatl6_1h_isopycnal_section.01440.nc'
        grid_name = '/home/datawork-lops-megatl/EDDY/gigatl6/iso/gigatl6_1h_isopycnal_section.01440.nc'
        nc = Dataset(grid_name,'r')
        isopycnals = nc.variables['isopycnal'][:]
        nc.close()

        s_rho =  5 #
        depth = isopycnals[s_rho]

    # Time loop
    for time in range(31056+120, 31176+120):
        # hourly data
        dtfile = 1
        # 12-hourly
        #dtfile = 12

        tfile = 5*24//dtfile

        ###########
        realyear_origin = datetime(2004,1,15)
        date = realyear_origin + timedelta(days=float(time)*dtfile/24.)

        infiletime = time%tfile
        filetime = time - time%tfile

        # Using times:
        #grid_name = './GIGATL6/gigatl6_1h_horizontal_section.' + '{0:05}'.format(filetime) + '.nc'

        if varname =='ow':
            grid_name = '/home/datawork-lops-megatl/EDDY/gigatl6/iso/gigatl6_1h_isopycnal_section.' + '{0:05}'.format(filetime) + '.nc'
        elif varname == 'zeta':
            # or using dates
            date1 = realyear_origin + timedelta(days=float(filetime)*dtfile/24.)
            date2 = date1 + timedelta(days=float(4.5))

            filedate =  '{0:04}'.format(date1.year)+'-'+\
                        '{0:02}'.format(date1.month) + '-'+\
                        '{0:02}'.format(date1.day)+ '-'+\
                        '{0:04}'.format(date2.year)+'-'+\
                        '{0:02}'.format(date2.month) + '-' +\
                        '{0:02}'.format(date2.day)

            #grid_name = './HIS/GIGATL6_1h_inst_surf_' + filedate + '.nc'
            grid_name = '/home/datawork-lops-megatl/GIGATL6/GIGATL6_1h/HIS/GIGATL6_1h_inst_surf_' + filedate + '.nc'

            print(f"file used is {grid_name}")
            print(f"infiletime is {infiletime}")

        # Identification
        if varname=='zeta':
            lon_name, lat_name = 'nav_lon_rho', 'nav_lat_rho'
        elif varname=='ow':
            lon_name, lat_name = 'lon', 'lat'

        # domain grid points: x1, x2, y1, y2
        #x1, x2 = 0, 1500
        #y1, y2 = 0, 1500

        x1, x2 = 400, 600
        y1, y2 = 900, 1200

        h = RomsDataset(grid_name, lon_name, lat_name,
                indexs=dict(time=infiletime,time_counter=infiletime,
                eta_rho=slice(y1, y2),
                xi_rho=slice(x1, x2),
                eta_v=slice(y1, y2-1),
                xi_u=slice(x1, x2-1),
                y_rho=slice(y1, y2),
                x_rho=slice(x1, x2),
                y_v=slice(y1, y2-1),
                y_u=slice(y1, y2),
                x_u=slice(x1, x2-1),
                x_v=slice(x1, x2),
                s_rho=s_rho)
        )

        # Identification
        if varname=='zeta':
            z_min = -0.5 ; z_max = 0.5; step = 0.0001        
        elif varname=='ow':
            # ow is multiplied by 1e10
            z_min = -1; z_max = -0.01; step = 0.01

        #continue

        a, c = h.eddy_identification(varname, 'u', 'v', date, z_min =  z_min, z_max = z_max, step = step, pixel_limit=(5, 5000), shape_error=55, force_height_unit='m',force_speed_unit='m/s',vorticity_name='vrt')

        #a, c = h.eddy_identification('zeta', 'u', 'v', date, step = 0.002, pixel_limit=(5, 2000), force_height_unit='m',force_speed_unit='m/s', shape_error=55)
        ####

        filename = 'gigatl6_1h_' + varname + '_'+  '{0:04}'.format(depth) + '_'

        with Dataset(date.strftime('/home2/datawork/leisenri/Desktop2/py-eddy-tracker-master/zetatest7/' + folder + 'Anticyclonic_' + filename + '%Y%m%d%H.nc'), 'w') as h:
            a.to_netcdf(h)
        with Dataset(date.strftime('/home2/datawork/leisenri/Desktop2/py-eddy-tracker-master/zetatest7/' + folder + 'Cyclonic_' + filename + '%Y%m%d%H.nc'), 'w') as h:
            c.to_netcdf(h)

        # PLOT
        fig = plt.figure(figsize=(15,7))
        ax = fig.add_axes([.03,.03,.94,.94])
        ax.set_title('Eddies detected -- Cyclonic(red) and Anticyclonic(blue)')
        #ax.set_ylim(-75,75)
        ax.set_xlim(250,360)
        ax.set_aspect('equal')
        a.display(ax, color='b', linewidth=.5)
        c.display(ax, color='r', linewidth=.5)
        ax.grid()
        fig.savefig('/home2/datawork/leisenri/Desktop2/py-eddy-tracker-master/zetatest7/' + folder + 'eddies_' + date.strftime( filename + '%Y%m%d%H') +'.png')

[AntSimi]

Firstly, i never work with other field than SLA(Sea Level Anomalie) or ADT(Absolute Dynamic Topography) for eddy detection. Also, i think also you used a modified version of py-eddy-tracker because in the orginal one you can't set zmin and zmax for eddy detection.

But i imagine your problem come from some noisy values which stop contour collecting:

• First solution could be to apply a filtering step before eddy detection, to remove this high frequency signal which stop the contour collecting
• Second solution, add some code to allow very small local extremum in eddies

[evanmason]

As I see it the code is working, but there is a large disparity in the detected areas of the ring between your two subplots. This will mean you have a very small speed_radius in the bottom subplot, and this is unlikely given that NBC rings are large.

The easy fix is to do some post-processing of your tracks using the Loess filter (available from py_eddy_tracker.observations.tracking import TrackEddiesObservations) on your speed_radius and other variables as necessary.

See https://py-eddy-tracker.readthedocs.io/en/stable/_autosummary/py_eddy_tracker.observations.tracking.track_loess_filter.html?highlight=loess

[Mesharou]

Hi guys,

Thanks for the very quick answer! Loic is working with us at LOPS and using a code we slightly modified (https://github.com/Mesharou/py-eddy-tracker). What we call 'zeta' is in fact ADT.

Here I believe the problem is due to multiple extrema inside the ring. So just using a higher value for 'mle' seems to solve the issue. Note that filtering small scales does the trick too:

Cheers, Jonathan

[AntSimi]

Bottom right subplot have some contours smaller than top left subplot. Did you do more changes than 'mle'?

[Mesharou]

Indeed, the shape_error has also been reduced for the experiment in the bottom right panel.

Here is what you get if you change only mle: