Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Suggestion: correct coordinates after curl/norm #34

Open
simon3122 opened this issue Jun 15, 2016 · 4 comments
Open

Suggestion: correct coordinates after curl/norm #34

simon3122 opened this issue Jun 15, 2016 · 4 comments

Comments

@simon3122
Copy link

simon3122 commented Jun 15, 2016 •
edited
Loading

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
Copy link
Collaborator

apatlpo commented Jun 16, 2016

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
Copy link
Author

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
Copy link
Collaborator

apatlpo commented Jul 4, 2016 •
edited
Loading

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
Copy link
Collaborator

apatlpo commented Jul 4, 2016 •
edited
Loading

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}

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@apatlpo @simon3122