Suggestion: correct coordinates after curl/norm

[simon3122]

Hello,

Here is my example: I am reading a u,v field, respectively on U-grid and V-grid but also with respective vertical coordinates: depthu and depthv. I want to compute curl, which will be centered on f-grid. I point out that the issue is similar with vector norm (centered on T-grid)

            vec_xr = gridsc.VectorField2d(vi_xr, vj_xr,
                                   x_component_grid_location='u',
                                   y_component_grid_location='v')

            #EE mode
            vcurl_xr=grd.vertical_component_of_curl(vec_xr)
            print vi_xr
            print vj_xr
            print vec_xr
            print vcurl_xr

Here is the output:

<xarray.DataArray (t: 5, depthu: 19, y: 3454, x: 5422)>
dask.array<elemwis..., shape=(5, 19, 3454, 5422), dtype=float64, chunksize=(1, 1, 1727, 2711)>
Coordinates:
    nav_lat        (y, x) float32 26.5648 26.5648 26.5648 26.5648 26.5648 ...
    nav_lon        (y, x) float32 -81.4512 -81.4346 -81.4179 -81.4012 ...
  * depthu         (depthu) float32 0.480455 1.55879 2.79421 4.18731 5.73867 ...
  * x              (x) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ...
  * y              (y) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ...
  * t              (t) datetime64[ns] 2008-01-07T12:00:00 ...
    time_centered  (t) datetime64[ns] 2008-01-07T12:00:00 ...
Attributes:
    grid_location: u

<xarray.DataArray (t: 5, depthv: 19, y: 3454, x: 5422)>
dask.array<elemwis..., shape=(5, 19, 3454, 5422), dtype=float64, chunksize=(1, 1, 1727, 2711)>
Coordinates:
    nav_lat        (y, x) float32 26.5648 26.5648 26.5648 26.5648 26.5648 ...
    nav_lon        (y, x) float32 -81.4512 -81.4346 -81.4179 -81.4012 ...
  * depthv         (depthv) float32 0.480455 1.55879 2.79421 4.18731 5.73867 ...
  * x              (x) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ...
  * y              (y) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ...
  * t              (t) datetime64[ns] 2008-01-07T12:00:00 ...
    time_centered  (t) datetime64[ns] 2008-01-07T12:00:00 ...
Attributes:
    grid_location: v

VectorField2d(x_component=<xarray.DataArray (t: 5, depthu: 19, y: 3454, x: 5422)>
dask.array<elemwis..., shape=(5, 19, 3454, 5422), dtype=float64, chunksize=(1, 1, 1727, 2711)>
Coordinates:
    nav_lat        (y, x) float32 26.5648 26.5648 26.5648 26.5648 26.5648 ...
    nav_lon        (y, x) float32 -81.4512 -81.4346 -81.4179 -81.4012 ...
  * depthu         (depthu) float32 0.480455 1.55879 2.79421 4.18731 5.73867 ...
  * x              (x) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ...
  * y              (y) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ...
  * t              (t) datetime64[ns] 2008-01-07T12:00:00 ...
    time_centered  (t) datetime64[ns] 2008-01-07T12:00:00 ...
Attributes:
    grid_location: u, y_component=<xarray.DataArray (t: 5, depthv: 19, y: 3454, x: 5422)>
dask.array<elemwis..., shape=(5, 19, 3454, 5422), dtype=float64, chunksize=(1, 1, 1727, 2711)>
Coordinates:
    nav_lat        (y, x) float32 26.5648 26.5648 26.5648 26.5648 26.5648 ...
    nav_lon        (y, x) float32 -81.4512 -81.4346 -81.4179 -81.4012 ...
  * depthv         (depthv) float32 0.480455 1.55879 2.79421 4.18731 5.73867 ...
  * x              (x) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ...
  * y              (y) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ...
  * t              (t) datetime64[ns] 2008-01-07T12:00:00 ...
    time_centered  (t) datetime64[ns] 2008-01-07T12:00:00 ...
Attributes:
    grid_location: v, x_component_grid_location='u', y_component_grid_location='v')

<xarray.DataArray 'vertical component of the curl' (t: 5, depthv: 19, y: 3454, x: 5422, depthu: 19)>
dask.array<elemwis..., shape=(5, 19, 3454, 5422, 19), dtype=float64, chunksize=(1, 1, 1727, 2711, 1)>
Coordinates:
  * depthv         (depthv) float32 0.480455 1.55879 2.79421 4.18731 5.73867 ...
  * x              (x) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ...
  * y              (y) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ...
  * t              (t) datetime64[ns] 2008-01-07T12:00:00 ...
    time_centered  (t) datetime64[ns] 2008-01-07T12:00:00 ...
  * depthu         (depthu) float32 0.480455 1.55879 2.79421 4.18731 5.73867 ...
Attributes:
    grid_location: f

Original coordinates of the u,v components are broadcast through the operation (here, curl) and the result is a 4-coordinate variable.

It can be very consequent in memory and writing operations, because it creates a 4-D variable which is N times bigger than the expected 3-D variable, with N as the vertical size.

Best regards,
Simon

[apatlpo]

Hi Simon,
The operation grd.vertical_component_of_curl does not handle depth dependent arrays at the moment.
This will come once the treatment of the vertical dimension will be added to oocgcm.
Meanwhile you can manually loop around the vertical coordinate and compute the curl on horizontal slices of the arrays.
Let us know if this works and please post the solution once it's works,
cheers
aurelien

[simon3122]

Hi Aurelien,
Here is my solution:

I switched off chunking along vertical coordinate and passed an explicit loop. It is possible to fill back a 3D-array with 'dict' and the .values attribute:

            # initialises vcurl
            print vi_xr
            vcurl_xr=vi_xr

            # loop over depth levels
            for k in range(lev1,lev2):

                 # constructs a VectorField2d object: NY*NX
                 vec_xr = gridsc.VectorField2d(vi_xr.isel(depthu=k), 
                                               vj_xr.isel(depthv=k),
                                   x_component_grid_location='u',
                                   y_component_grid_location='v')

                 # compute curl
                 vcurl_xr_1lev=grd.vertical_component_of_curl(vec_xr)
                 # stores in NY*NX*NZ array
                 vcurl_xr[dict(depthu=k)].values=vcurl_xr_1lev.values

            # renaming
            vcurl_xr.name=vcurl_xr_1lev.name
            print vcurl_xr

with the output:

<xarray.DataArray 'vozocrtx' (t: 5, depthu: 19, y: 3454, x: 5422)>
dask.array<getitem..., shape=(5, 19, 3454, 5422), dtype=float64, chunksize=(1, 19, 1727, 2711)>
Coordinates:
    nav_lat        (y, x) float32 26.5648 26.5648 26.5648 26.5648 26.5648 ...
    nav_lon        (y, x) float32 -81.4512 -81.4346 -81.4179 -81.4012 ...
  * depthu         (depthu) float32 0.480455 1.55879 2.79421 4.18731 5.73867 ...
  * x              (x) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ...
  * y              (y) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ...
  * t              (t) datetime64[ns] 2008-01-07T12:00:00 ...
    time_centered  (t) datetime64[ns] 2008-01-07T12:00:00 ...
Attributes:
    long_name: ocean current along i-axis
    units: m/s
    online_operation: average
    interval_operation: 40s
    interval_write: 5d

<xarray.DataArray 'vertical component of the curl' (t: 5, depthu: 19, y: 3454, x: 5422)>
dask.array<getitem..., shape=(5, 19, 3454, 5422), dtype=float64, chunksize=(1, 19, 1727, 2711)>
Coordinates:
    nav_lat        (y, x) float32 26.5648 26.5648 26.5648 26.5648 26.5648 ...
    nav_lon        (y, x) float32 -81.4512 -81.4346 -81.4179 -81.4012 ...
  * depthu         (depthu) float32 0.480455 1.55879 2.79421 4.18731 5.73867 ...
  * x              (x) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ...
  * y              (y) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ...
  * t              (t) datetime64[ns] 2008-01-07T12:00:00 ...
    time_centered  (t) datetime64[ns] 2008-01-07T12:00:00 ...
Attributes:
    long_name: ocean current along i-axis
    units: m/s
    online_operation: average
    interval_operation: 40s
    interval_write: 5d

[apatlpo]

Hi all,
The following commit proposes a fix where all vertical coordinates variables are referred to as 'depth':
commit

The following code tests the implementation:

""" Illustrates a bug with oogcm when loading u/v variables which are described
by coordinates depthu and depthv and are not handled properly when curl/divergence/norm
are computed
"""


# - Modules
#
import sys, os
import xarray as xr
import time, calendar
import matplotlib.pyplot as plt
from matplotlib import animation
import cartopy.crs as ccrs
import math #SF
import numpy as np #SF

#sys.path.append("../../oocgcm/")
from oocgcm.oceanmodels.nemo import grids
from oocgcm.oceanmodels.nemo import io
from oocgcm.core import grids as gridsc

chunks = (1727, 2711)
xr_chunks = {'x': chunks[-1], 'y': chunks[-2]}
xr_chunks_t = {'x': chunks[-1], 'y': chunks[-2], 'time_counter':1}
xr_chunks_tz = {'x': chunks[-1], 'y': chunks[-2], 'time_counter':1, 'deptht': 1}
xr_chunks_tzu = {'x': chunks[-1], 'y': chunks[-2], 'time_counter':1, 'depthu': 1}
xr_chunks_tzv = {'x': chunks[-1], 'y': chunks[-2], 'time_counter':1, 'depthv': 1}

# - Variables names and plotting parameters
#
# dictionary contains: netcdf variable name, colormap to use, vmin, vmax
vdict={}
vdict['current']={'vnameu': 'vozocrtx', 'vnamev': 'vomecrty',
                  'chunksu': xr_chunks_tzu, 'chunksv': xr_chunks_tzv,
                  'flag3D': True, 'flagVector':True,
                  'cmap': plt.cm.YlOrBr_r, 'vmin': 0., 'vmax': 1.}


#- picks which variable is considered
vkey='current'

# - Parameter
natl60_path = '/home7/pharos/othr/NATL60/'
coordfile = natl60_path + 'NATL60-I/NATL60_coordinates_v4.nc'
maskfile = natl60_path + 'NATL60-I/NATL60_v4.1_cdf_byte_mask.nc'


filenatl60u = natl60_path+'NATL60-MJM155-S/5d/2008/NATL60-MJM155_y2008m01d*gridU.nc'
filenatl60v = natl60_path+'NATL60-MJM155-S/5d/2008/NATL60-MJM155_y2008m01d*gridV.nc'

# - creating the grid object
grd = grids.nemo_2d_grid(nemo_coordinate_file=coordfile, nemo_byte_mask_file=maskfile, chunks=xr_chunks)

# - defining a 2D xarray
vi_xr = io.return_xarray_mfdataset(filenatl60u, chunks=vdict[vkey]['chunksu'])[vdict[vkey]['vnameu']][:]
vj_xr = io.return_xarray_mfdataset(filenatl60v, chunks=vdict[vkey]['chunksv'])[vdict[vkey]['vnamev']][:]    #vi_xr = xr.open_mfdataset(filenatl60u, chunks=vdict[vkey]['chunks'], lock=False)[vdict[vkey]['vnameu']][:]

vi_xr=vi_xr.isel(t=0)
vj_xr=vj_xr.isel(t=0)

vec_xr = gridsc.VectorField2d(vi_xr, vj_xr,
                              x_component_grid_location='u',
                              y_component_grid_location='v')
vcurl_xr = grd.vertical_component_of_curl(vec_xr)
v_xr = grd.norm_of_vectorfield(vec_xr)


print "----"
print vec_xr
print vcurl_xr

[apatlpo]

An update on the previous fix proposal is found at: ce67ffd

The chunks that needs to be provided need now to be consistent with final names, e.g.:

xr_chunks_tz = {'x': chunks[-1], 'y': chunks[-2], 'time_counter':1, 'deptht': 1}

should now be:

xr_chunks_tz = {'x': chunks[-1], 'y': chunks[-2], 't':1, 'depth': 1}