How to avoid anomalous strip when executing the difference in LLC90 data, like ECCOv4

When I try to calculate the difference of velocity of ECCOv4 (LLC90 data) using the following python codes:
div_uv = (grid_llc.diff(UVELMASS.mean(dim=‘time’) * data_uv.dyG, ‘X’, boundary=‘fill’) +
grid_llc.diff(VVELMASS.mean(dim=‘time’) * data_uv.dxG, ‘Y’, boundary=‘fill’))/(data_uv.rA)

I get the result with anomalous strip:

Does anyone know how to avoid anomalous strip when executing the difference in LLC90 data, like ECCOv4?

I do find a way to avoid it, like using grid.diff_2d_vector, how ever this function seems only can work with (i,j) grid, but now work for () grid, for example
zeta=grid_llc.diff_2d_vector({‘X’: (Vm.mean(dim=‘time’) * data.dyC),
‘Y’: (Um.mean(dim=‘time’) * data.dxC)},boundary=‘fill’)

which I get the bug reminding:-------NotImplementedError: Only vector interpolation to cell center is implemented, but vector X component is defined at center (dims: (‘face’, ‘j_g’, ‘i’))

1 Like

You might get more help if you give us a way to reproduce this? E.g. post a notebook here that has access to the output you use?

1 Like

I’ve created a topic with information on how to share notebook code

##load
import numpy as np
import xarray as xr
import pandas as pd
import os
from xmitgcm import open_mdsdataset
import xgcm
#ecco_v4_py
import sys
from os.path import join,expanduser
user_ECCO_dir = ('/g/data/p10/rh1363/ECCO/ECCOv4-py/')
sys.path.append(join(user_ECCO_dir,'ECCOv4-py'))
import ecco_v4_py as ecco
import warnings
warnings.simplefilter('ignore')
os.environ['PYTHONWARNOINGS'] = 'ignore'
#for plot
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
from cartopy.util import add_cyclic_point
##load data
data_uv = xr.open_dataset('/g/data/p10/rh1363/ECCO/uv_test.nc')
data_uv.attrs=[]
data_uv.load()
# define the connectivity between faces
face_connections = {'face':
                    {0: {'X':  ((12, 'Y', False), (3, 'X', False)),
                         'Y':  (None,             (1, 'Y', False))},
                     1: {'X':  ((11, 'Y', False), (4, 'X', False)),
                         'Y':  ((0, 'Y', False),  (2, 'Y', False))},
                     2: {'X':  ((10, 'Y', False), (5, 'X', False)),
                         'Y':  ((1, 'Y', False),  (6, 'X', False))},
                     3: {'X':  ((0, 'X', False),  (9, 'Y', False)),
                         'Y':  (None,             (4, 'Y', False))},
                     4: {'X':  ((1, 'X', False),  (8, 'Y', False)),
                         'Y':  ((3, 'Y', False),  (5, 'Y', False))},
                     5: {'X':  ((2, 'X', False),  (7, 'Y', False)),
                         'Y':  ((4, 'Y', False),  (6, 'Y', False))},
                     6: {'X':  ((2, 'Y', False),  (7, 'X', False)),
                         'Y':  ((5, 'Y', False),  (10, 'X', False))},
                     7: {'X':  ((6, 'X', False),  (8, 'X', False)),
                         'Y':  ((5, 'X', False),  (10, 'Y', False))},
                     8: {'X':  ((7, 'X', False),  (9, 'X', False)),
                         'Y':  ((4, 'X', False),  (11, 'Y', False))},
                     9: {'X':  ((8, 'X', False),  None),
                         'Y':  ((3, 'X', False),  (12, 'Y', False))},
                     10: {'X': ((6, 'Y', False),  (11, 'X', False)),
                          'Y': ((7, 'Y', False),  (2, 'X', False))},
                     11: {'X': ((10, 'X', False), (12, 'X', False)),
                          'Y': ((8, 'Y', False),  (1, 'X', False))},
                     12: {'X': ((11, 'X', False), None),
                          'Y': ((9, 'Y', False),  (0, 'X', False))}}}
# create the grid object
grid_llc = xgcm.Grid(data_uv, periodic=False, face_connections=face_connections)
zeta = (-grid_llc.diff(data_uv.UVELMASS * data_uv.dxC, 'Y',boundary='fill',fill_value='NaN') +
        grid_llc.diff(data_uv.VVELMASS * data_uv.dyC, 'X',boundary='fill',fill_value='NaN'))/data_uv.rAz

# plot the wrong result
plt.figure(figsize=(12,5), dpi= 90)
ecco.plot_proj_to_latlon_grid(data_uv.XG, data_uv.YG, \
                              zeta.isel(k=0), \
                              user_lon_0=180,\
                              projection_type='robin',\
                              plot_type = 'pcolormesh', \
                              dx=0.25,dy=0.25,show_colorbar=True,\
                              cmin=-1.e-6,cmax=1.e-6);
# the following do not work
zeta0=grid_llc.diff_2d_vector({'X': (data_uv.VVELMASS * data_uv.dyC), \
                               'Y': (data_uv.UVELMASS * data_uv.dxC)},\
                              boundary='fill')
zeta2 = (zeta0['X'] -zeta0['Y']) / (data_uv.rA)


Thank you so much. Now I add the codes sucecessfully.

Thank you for you advice. Now I add the codes here.

The data can be download from here calculate vorticity using diff_2d_vector in LLC90 grid · Issue #599 · xgcm/xgcm · GitHub

Seems like access to p10 is required, right?

No you don’t need to, just because I put the data in p10. You can download from here (calculate vorticity using diff_2d_vector in LLC90 grid · Issue #599 · xgcm/xgcm · GitHub)

Ruhui, when you’re asking for help it is best to make it as easy as possible for people to assist you. One way you can do that is to put the data you’re using somewhere others can access it. The /scratch/public/ area on gadi is for just this purpose. You can make a directory in there, make sure it is world readable and copy the data there, again taking care to make it world readable.

See the NCI docs for more detail.

Thank you so much. Now I put the test data on /scratch/public/rh1363/
Anyone can use the follow code to test it without any change.

##load
import numpy as np
import xarray as xr
import pandas as pd
import os
from xmitgcm import open_mdsdataset
import xgcm
#ecco_v4_py
import sys
from os.path import join,expanduser
user_ECCO_dir = ('/scratch/public/rh1363/ECCOv4-py/')
sys.path.append(join(user_ECCO_dir,'ECCOv4-py'))
import ecco_v4_py as ecco
import warnings
warnings.simplefilter('ignore')
os.environ['PYTHONWARNOINGS'] = 'ignore'
#for plot
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
from cartopy.util import add_cyclic_point
##load data
data_uv = xr.open_dataset('/scratch/public/rh1363/uv_test.nc')
data_uv.attrs=[]
data_uv.load()
# define the connectivity between faces
face_connections = {'face':
                    {0: {'X':  ((12, 'Y', False), (3, 'X', False)),
                         'Y':  (None,             (1, 'Y', False))},
                     1: {'X':  ((11, 'Y', False), (4, 'X', False)),
                         'Y':  ((0, 'Y', False),  (2, 'Y', False))},
                     2: {'X':  ((10, 'Y', False), (5, 'X', False)),
                         'Y':  ((1, 'Y', False),  (6, 'X', False))},
                     3: {'X':  ((0, 'X', False),  (9, 'Y', False)),
                         'Y':  (None,             (4, 'Y', False))},
                     4: {'X':  ((1, 'X', False),  (8, 'Y', False)),
                         'Y':  ((3, 'Y', False),  (5, 'Y', False))},
                     5: {'X':  ((2, 'X', False),  (7, 'Y', False)),
                         'Y':  ((4, 'Y', False),  (6, 'Y', False))},
                     6: {'X':  ((2, 'Y', False),  (7, 'X', False)),
                         'Y':  ((5, 'Y', False),  (10, 'X', False))},
                     7: {'X':  ((6, 'X', False),  (8, 'X', False)),
                         'Y':  ((5, 'X', False),  (10, 'Y', False))},
                     8: {'X':  ((7, 'X', False),  (9, 'X', False)),
                         'Y':  ((4, 'X', False),  (11, 'Y', False))},
                     9: {'X':  ((8, 'X', False),  None),
                         'Y':  ((3, 'X', False),  (12, 'Y', False))},
                     10: {'X': ((6, 'Y', False),  (11, 'X', False)),
                          'Y': ((7, 'Y', False),  (2, 'X', False))},
                     11: {'X': ((10, 'X', False), (12, 'X', False)),
                          'Y': ((8, 'Y', False),  (1, 'X', False))},
                     12: {'X': ((11, 'X', False), None),
                          'Y': ((9, 'Y', False),  (0, 'X', False))}}}
# create the grid object
grid_llc = xgcm.Grid(data_uv, periodic=False, face_connections=face_connections)
zeta = (-grid_llc.diff(data_uv.UVELMASS * data_uv.dxC, 'Y',boundary='fill',fill_value='NaN') +
        grid_llc.diff(data_uv.VVELMASS * data_uv.dyC, 'X',boundary='fill',fill_value='NaN'))/data_uv.rAz

# plot the wrong result
plt.figure(figsize=(12,5), dpi= 90)
ecco.plot_proj_to_latlon_grid(data_uv.XG, data_uv.YG, \
                              zeta.isel(k=0), \
                              user_lon_0=180,\
                              projection_type='robin',\
                              plot_type = 'pcolormesh', \
                              dx=0.25,dy=0.25,show_colorbar=True,\
                              cmin=-1.e-6,cmax=1.e-6);
# the following do not work
zeta0=grid_llc.diff_2d_vector({'X': (data_uv.VVELMASS * data_uv.dyC), \
                               'Y': (data_uv.UVELMASS * data_uv.dxC)},\
                              boundary='fill')
zeta2 = (zeta0['X'] -zeta0['Y']) / (data_uv.rA)
2 Likes

I started to take a look, and quickly realised I had no idea how this data is structured. I searched and found this reference

and checked out some of the ECCO Version 4 Python Tutorial. They have specialised functions for doing global plots

https://ecco-v4-python-tutorial.readthedocs.io/ECCO_v4_Plotting_Tiles.html#Plotting-all-13-tiles-with-plot_proj_to_latlon_grid

Have you tried replicating this for the data you’re interested in?

I also see the there is a worked example for ECCOv4 in the xgcm docs

https://xgcm.readthedocs.io/en/latest/xgcm-examples/01_eccov4.html

It would be a good idea to adapt the LLCMapper class provided in that example for your uses.