How to do an identification on ORCA025 irregular grid

[annesophiefortin]

The file attached contains the sea level anomalies ('sla') and the geostrophic velocities ('u','v'). The coordinates lat/lon have two dimensions. Note also that the grid is tripolar.

Because of the tri-polarity of the grid, we add the following line of code: if centlat_e > 60: continue to the eddy_identification function (on line 582) so that it skip identifying eddies too close to the north poles.

However, we get an out-of-bound index error message on line 80 of eddy_feature, i.e. for the init of the class Amplituder. IndexError: index 1442 is out of bounds for axis 1 with size 1442

Note that for small domain (where less than 1442 eddies are identified) the code successfully identifies eddies.

GIOPS_NativeGrid_20200424.zip

[AntSimi]

I am not sure, but i think "nav_lon" content produce problem with matplotlib contour:

When we apply matplotlib.pyplot.contour on all the ORCA data, we have that, i think contour need no jump in lon. i will try to find a solution...

from netCDF4 import Dataset
from matplotlib import pyplot as plt
import numpy as np
import pylook
if __name__ == '__main__':
    grid_name = 'GIOPS_NativeGrid_20200424.nc'
    lon_name, lat_name = 'nav_lon', 'nav_lat'
    ax = plt.subplot(111, projection='plat_carre')
    with Dataset(grid_name) as h:
        lon = h.variables[lon_name][:]
        lat = h.variables[lat_name][:]
        sla = h.variables['sla'][:]
        ax.contour(lon,lat, sla, levels=np.arange(-2,2,.25))
        plt.show()

[AntSimi]

First, i am not familiar with ORCAGrid.I don't know if all ORCAGrid have jump in nav_lon. Solution given below isn't added for now in py-eddy-tracker(PET), but anyone could use with current version of PET.

Explanations

Problem are due to nav_lon have jump around indice 450. Matplolib contour function doesn't work with this nav_lon. This function provide big contour which are not usable for detection. shift of nav_lon seems to be enough to have correct identification.

Warnings for ORCAGrid

• UnregularGrid doesn't try to close around sphere, so identification at grid border could be not accurate
• I didn't check identification result at high latitude near of ORCAGrid closing, above Russia and Canada

Solution used

I rewrote load function to add two lines to shift a part of nav_lon

import pylook
from netCDF4 import Dataset
from matplotlib import pyplot as plt
from datetime import datetime
from py_eddy_tracker import start_logger
from py_eddy_tracker.dataset.grid import UnRegularGridDataset

class ORCAGrid(UnRegularGridDataset):

    def load(self):
        """Load variable (data)
        """
        x_name, y_name = self.coordinates
        with Dataset(self.filename) as h:
            self.x_dim = h.variables[x_name].dimensions
            self.y_dim = h.variables[y_name].dimensions

            sl_x = [self.indexs.get(dim, slice(None)) for dim in self.x_dim]
            sl_y = [self.indexs.get(dim, slice(None)) for dim in self.y_dim]
            self.vars[x_name] = h.variables[x_name][sl_x]
            self.vars[y_name] = h.variables[y_name][sl_y]

            self.x_c = self.vars[x_name]
            ###
            # To avoid jump of longitude, needed for file in example
            x_v = self.x_c[:,:500]
            x_v[x_v > 0] -= 360
            ###
            self.y_c = self.vars[y_name]

            self.init_pos_interpolator()

if __name__ == '__main__':
    start_logger().setLevel('DEBUG')
    grid_name = 'GIOPS_NativeGrid_20200424.nc'
    lon_name, lat_name = 'nav_lon', 'nav_lat'

    # Must be set with time of grid
    date = datetime(2020, 4, 24)  

    # Identification every 2 mm
    args = ('sla', 'u', 'v', date, .002)
    kwargs = dict(pixel_limit=(5, 2000), shape_error=55)
    # without filter
    h = ORCAGrid(grid_name, lon_name, lat_name, indexs=dict(time=0))
    a_unfiltered, c_unfiltered = h.eddy_identification(*args, **kwargs)
    # with filter
    h = ORCAGrid(grid_name, lon_name, lat_name, indexs=dict(time=0))
    # Ugly filter for unregular
    h.high_filter('sla', 500)
    a, c = h.eddy_identification(*args, **kwargs)
    # plot
    ax = plt.subplot(111, projection='plat_carre')
    ax.grid()
    kwargs_display = dict(lw=.5)
    a_unfiltered.display(ax, label='Anticyclonic unfiltered', color='r', **kwargs_display)
    c_unfiltered.display(ax, label='Cyclonic unfiltered', color='b', **kwargs_display)
    a.display(ax, label='Anticyclonic', color='green', **kwargs_display)
    c.display(ax, label='Cyclonic', color='purple', **kwargs_display)
    ax.legend()
    ax.set_title('Identification on ORCA grid')
    plt.show()

Quantities seems to be coherent with classic AVISO ADT grid

Filter

Filter step is important to remove long scale in grid, in order to highlight mesoscale structure which could hide by large scale. I used in this example for unregular grid an ugly solution to filter sea level grid, which also could add some other problems ...