Interpolating fields from the model llc grid to a regular lat lon grid

Objectives

  1. Learn how to interpolate scalar and vector fields from ECCOv4’s lat-lon-cap 90 (llc90) model grid to the more common regular latitude-longitude grid.
  2. Learn how to save these interpolated fields as netCDF for later analysis

Introduction

Recall the orientations of the 13 tiles of the ECCOv4 native llc90 model grid.

llc90 tile layout

Tiles 7-12 are rotated 90 degrees counter-clockwise relative to tiles 0-5.

In this tutorial we demonstrate two methods for mapping scalar and vector fields from the llc90 model grid to “regular” latitude-longitude grids of arbitrary resolution.

Note: There are many methods on can use to map between the grids (e.g., nearest neighbor, bilinear interpolation, bin-averaging, etc.), each with its own advantages.)

How to interpolate scalar ECCO fields to a lat-lon grid

Scalar fields are the most straightforward fields to interpolate.

Preliminaries: Load fields

First, let’s load the all 13 tiles for sea surface height and the model grid parameters.

[1]:
import numpy as np
import sys
import xarray as xr
import matplotlib.pyplot as plt
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')
from pprint import pprint
import importlib
[2]:
## Import the ecco_v4_py library into Python
## =========================================

## -- If ecco_v4_py is not installed in your local Python library,
##    tell Python where to find it.  For example, if your ecco_v4_py
##    files are in /Users/ifenty/ECCOv4-py/ecco_v4_py, then use:

from pathlib import Path
sys.path.append(str(Path('c:/Users/Ian/ECCOv4-py')))
import ecco_v4_py as ecco
[3]:
## Set top-level file directory for the ECCO NetCDF files
## =================================================================
# base_dir = homehome/username/'
ECCO_dir = Path('E:/inSync Share/Projects/ECCOv4/Release4/')
[4]:
## Load the model grid
grid_dir= ECCO_dir / 'nctiles_grid/'
[5]:
## Load the model grid
ecco_grid = ecco.load_ecco_grid_nc(grid_dir, 'ECCO-GRID.nc')
[6]:
## Load one year of 2D daily data, SSH, SST, and SSS
day_mean_dir= ECCO_dir / 'nctiles_monthly/'

ecco_vars = ecco.recursive_load_ecco_var_from_years_nc(day_mean_dir, \
                                           vars_to_load=['SSH'], \
                                           years_to_load=2000,less_output=True)

ecco_ds = []

## Merge the ecco_grid with the ecco_vars to make the ecco_ds
ecco_ds = xr.merge((ecco_grid , ecco_vars)).load()

pprint(ecco_ds.data_vars)
loading files of  SSH
Data variables:
    SSH      (time, tile, j, i) float32 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

Plotting the dataset

Plotting the ECCOv4 fields was covered in an earlier tutorial. Before demonstrating interpolation, we will first plot one of our SSH fields. Take a closer look at the arguments of plot_proj_to_latlon_grid, dx=2 and dy=2.

[7]:
plt.figure(figsize=(12,6), dpi= 90)

tmp_plt = ecco_ds.SSH.isel(time=1)
tmp_plt = tmp_plt.where(ecco_ds.hFacC.isel(k=0) !=0)

ecco.plot_proj_to_latlon_grid(ecco_ds.XC, \
                              ecco_ds.YC, \
                              tmp_plt, \
                              plot_type = 'pcolormesh', \
                              dx=2,\
                              dy=2, \
                              projection_type = 'robin',\
                              less_output = True);

_images/ECCO_v4_Interpolating_Fields_to_LatLon_Grid_9_0.png

These dx and dy arguments tell the plotting function to interpolate the native grid tiles onto a lat-lon grid with spacing of dx degrees longitude and dy degrees latitude. If we reduced dx and dy the resulting map would have finer features: Compare with the map when dx=dy=0.25 degrees:

[8]:
plt.figure(figsize=(12,6), dpi= 90)

tmp_plt = ecco_ds.SSH.isel(time=1)
tmp_plt = tmp_plt.where(ecco_ds.hFacC.isel(k=0) !=0)

ecco.plot_proj_to_latlon_grid(ecco_ds.XC, \
                              ecco_ds.YC, \
                              tmp_plt, \
                              plot_type = 'pcolormesh', \
                              dx=0.25,\
                              dy=0.25, \
                              projection_type = 'robin',\
                              less_output = True);

_images/ECCO_v4_Interpolating_Fields_to_LatLon_Grid_11_0.png

Of course you can interpolate to arbitrarily high resolution lat-lon grids, the model resolution will ultimately determine the smallest resolvable features.

Under the hood of plot_proj_to_lat_longrid is a call to the very useful routine resample_to_latlon which is the now described in more detail

resample_to_latlon

resample_to_latlon takes a field with a cooresponding set of lat lon coordinates (the source grid) and interpolates to a new lat-lon target grid. The arrangement of coordinates in the source grid is arbitrary. One also provides the bounds of the new lat lon grid. The example shown above uses -90..90N and -180..180E by default. In addition, one specifies which interpolation scheme to use (mapping method) and the radius of influence, the radius around the target grid cells within which to search for values from the source grid.

[9]:
help(ecco.resample_to_latlon)
Help on function resample_to_latlon in module ecco_v4_py.resample_to_latlon:

resample_to_latlon(orig_lons, orig_lats, orig_field, new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat, new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon, radius_of_influence=120000, fill_value=None, mapping_method='bin_average')
    Take a field from a source grid and interpolate to a target grid.

    Parameters
    ----------
    orig_lons, orig_lats, orig_field : xarray DataArray or numpy array  :
        the lons, lats, and field from the source grid

        new_grid_min_lat, new_grid_max_lat : float
                latitude limits of new lat-lon grid

    new_grid_delta_lat : float
        latitudinal extent of new lat-lon grid cells in degrees (-90..90)

    new_grid_min_lon, new_grid_max_lon : float
                longitude limits of new lat-lon grid (-180..180)

    new_grid_delta_llon : float
         longitudinal extent of new lat-lon grid cells in degrees

    radius_of_influence : float, optional.  Default 120000 m
        the radius of the circle within which the input data is search for
        when mapping to the new grid

    fill_value : float, optional. Default None
                value to use in the new lat-lon grid if there are no valid values
                from the source grid

        mapping_method : string, optional. Default 'bin_average'
        denote the type of interpolation method to use.
        options include
            'nearest_neighbor' - Take the nearest value from the source grid
                                                 to the target grid
            'bin_average'      - Use the average value from the source grid
                                                                 to the target grid

    RETURNS:
    new_grid_lon, new_grid_lat : ndarrays
        2D arrays with the lon and lat values of the new grid

    data_latlon_projection:
        the source field interpolated to the new grid

Demonstrations of resample_to_latlon

Global

First we will map to a global lat-lon grid at 1degree using nearest neighbor

[10]:
new_grid_delta_lat = 1
new_grid_delta_lon = 1

new_grid_min_lat = -90+new_grid_delta_lat/2
new_grid_max_lat = 90-new_grid_delta_lat/2

new_grid_min_lon = -180+new_grid_delta_lon/2
new_grid_max_lon = 180-new_grid_delta_lon/2

new_grid_lon, new_grid_lat, field_nearest_1deg =\
        ecco.resample_to_latlon(ecco_ds.XC, \
                                ecco_ds.YC, \
                                ecco_ds.SSH.isel(time=0),\
                                new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,\
                                new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,\
                                fill_value = np.NaN, \
                                mapping_method = 'nearest_neighbor',
                                radius_of_influence = 120000)
[11]:
# Dimensions of the new grid:
pprint(new_grid_lat.shape)
pprint(new_grid_lon.shape)
(180, 360)
(180, 360)
[12]:
pprint(new_grid_lon)

# the second dimension of new_grid_lon has the center longitude of the new grid cells
pprint(new_grid_lon[0,0:10])
array([[-179.5, -178.5, -177.5, ...,  177.5,  178.5,  179.5],
       [-179.5, -178.5, -177.5, ...,  177.5,  178.5,  179.5],
       [-179.5, -178.5, -177.5, ...,  177.5,  178.5,  179.5],
       ...,
       [-179.5, -178.5, -177.5, ...,  177.5,  178.5,  179.5],
       [-179.5, -178.5, -177.5, ...,  177.5,  178.5,  179.5],
       [-179.5, -178.5, -177.5, ...,  177.5,  178.5,  179.5]])
array([-179.5, -178.5, -177.5, -176.5, -175.5, -174.5, -173.5, -172.5,
       -171.5, -170.5])
[13]:
# The set of lat points
pprint(new_grid_lat)

# or as a 1D vector
# the first dimension of new_grid_lat has the center latitude of the new grid cells
pprint(new_grid_lat[0:10,0])
array([[-89.5, -89.5, -89.5, ..., -89.5, -89.5, -89.5],
       [-88.5, -88.5, -88.5, ..., -88.5, -88.5, -88.5],
       [-87.5, -87.5, -87.5, ..., -87.5, -87.5, -87.5],
       ...,
       [ 87.5,  87.5,  87.5, ...,  87.5,  87.5,  87.5],
       [ 88.5,  88.5,  88.5, ...,  88.5,  88.5,  88.5],
       [ 89.5,  89.5,  89.5, ...,  89.5,  89.5,  89.5]])
array([-89.5, -88.5, -87.5, -86.5, -85.5, -84.5, -83.5, -82.5, -81.5,
       -80.5])
[14]:
# plot the whole field
plt.figure(figsize=(12,6), dpi= 90)
plt.imshow(field_nearest_1deg,origin='lower')
[14]:
<matplotlib.image.AxesImage at 0x1ff70b5f9a0>
_images/ECCO_v4_Interpolating_Fields_to_LatLon_Grid_19_1.png

Notice that although we specified a fill_value of np.nan, the values on land are zero. This is because we did not first mask out original field with nans over dry points. Let’s nan out the land points first and interpolate again:

[15]:
original_field_with_land_mask= np.where(ecco_ds.maskC.isel(k=0)>0, ecco_ds.SSH.isel(time=0), np.nan)

new_grid_delta_lat = 1
new_grid_delta_lon = 1

new_grid_min_lat = -90+new_grid_delta_lat/2
new_grid_max_lat = 90-new_grid_delta_lat/2

new_grid_min_lon = -180+new_grid_delta_lon/2
new_grid_max_lon = 180-new_grid_delta_lon/2
new_grid_lon, new_grid_lat, field_nearest_1deg =\
        ecco.resample_to_latlon(ecco_ds.XC, \
                                ecco_ds.YC, \
                                original_field_with_land_mask,\
                                new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,\
                                new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,\
                                fill_value = np.NaN, \
                                mapping_method = 'nearest_neighbor',
                                radius_of_influence = 120000)

# plot the whole field, this time land values are nans.
plt.figure(figsize=(12,6), dpi= 90)
plt.imshow(field_nearest_1deg,origin='lower')
[15]:
<matplotlib.image.AxesImage at 0x1ff70989280>
_images/ECCO_v4_Interpolating_Fields_to_LatLon_Grid_21_1.png

Regional

1 degree, nearest neighbor

First we’ll interpolate only to a subset of the N. Atlantic at 1 degree.

[16]:
original_field_with_land_mask= np.where(ecco_ds.maskC.isel(k=0)>0, ecco_ds.SSH.isel(time=0), np.nan)

new_grid_delta_lat = 1
new_grid_delta_lon = 1

new_grid_min_lat = 30+new_grid_delta_lat/2
new_grid_max_lat = 82-new_grid_delta_lat/2

new_grid_min_lon = -90+new_grid_delta_lon/2
new_grid_max_lon = 10-new_grid_delta_lon/2

new_grid_lon, new_grid_lat, field_nearest_1deg =\
        ecco.resample_to_latlon(ecco_ds.XC, \
                                ecco_ds.YC, \
                                original_field_with_land_mask,\
                                new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,\
                                new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,\
                                fill_value = np.NaN, \
                                mapping_method = 'nearest_neighbor',
                                radius_of_influence = 120000)

# plot the whole field, this time land values are nans.
plt.figure(figsize=(6,6), dpi= 90)
plt.imshow(field_nearest_1deg,origin='lower',cmap='jet')
[16]:
<matplotlib.image.AxesImage at 0x1ff708d1ca0>
_images/ECCO_v4_Interpolating_Fields_to_LatLon_Grid_23_1.png

0.05 degree, nearest neighbor

Next we’ll interpolate the same domain at a much higher resolution, 0.05 degree:

[17]:
original_field_with_land_mask= np.where(ecco_ds.maskC.isel(k=0)>0, ecco_ds.SSH.isel(time=0), np.nan)

new_grid_delta_lat = 0.05
new_grid_delta_lon = 0.05

new_grid_min_lat = 30+new_grid_delta_lat/2
new_grid_max_lat = 82-new_grid_delta_lat/2

new_grid_min_lon = -90+new_grid_delta_lon/2
new_grid_max_lon = 10-new_grid_delta_lon/2


new_grid_lon, new_grid_lat, field_nearest_1deg =\
        ecco.resample_to_latlon(ecco_ds.XC, \
                                ecco_ds.YC, \
                                original_field_with_land_mask,\
                                new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,\
                                new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,\
                                fill_value = np.NaN, \
                                mapping_method = 'nearest_neighbor',
                                radius_of_influence = 120000)

# plot the whole field, this time land values are nans.
plt.figure(figsize=(6,6), dpi= 90)
plt.imshow(field_nearest_1deg,origin='lower',cmap='jet');
_images/ECCO_v4_Interpolating_Fields_to_LatLon_Grid_25_0.png

The new grid is much finer than the llc90 grid. Because we are using the nearest neighbor method, you can make out the approximate size and location of the llc90 model grid cells (think about it: nearest neighbor).

0.05 degree, bin average

With bin averaging many values from the source grid are ‘binned’ together and then averaged to determine the value in the new grid. The numer of grid cells from the source grid depends on the source grid resolution and the radius of influence. If we were to choose a radius of influence of 120000 m (a little longer than 1 degree longitude at the equator) the interpolated lat-lon map would be smoother than if we were to choose a smaller radius. To demonstrate:

bin average with 120000 m radius

[18]:
original_field_with_land_mask= np.where(ecco_ds.maskC.isel(k=0)>0, ecco_ds.SSH.isel(time=0), np.nan)

new_grid_lon, new_grid_lat, field_nearest_1deg =\
        ecco.resample_to_latlon(ecco_ds.XC, \
                                ecco_ds.YC, \
                                original_field_with_land_mask,\
                                new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,\
                                new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,\
                                fill_value = np.NaN, \
                                mapping_method = 'bin_average',
                                radius_of_influence = 120000)

# plot the whole field, this time land values are nans.
plt.figure(figsize=(6,6), dpi= 90)
plt.imshow(field_nearest_1deg,origin='lower',cmap='jet');
_images/ECCO_v4_Interpolating_Fields_to_LatLon_Grid_28_0.png

bin average with 40000m radius

[19]:
original_field_with_land_mask= np.where(ecco_ds.maskC.isel(k=0)>0, ecco_ds.SSH.isel(time=0), np.nan)

new_grid_lon, new_grid_lat, field_nearest_1deg =\
        ecco.resample_to_latlon(ecco_ds.XC, \
                                ecco_ds.YC, \
                                original_field_with_land_mask,\
                                new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,\
                                new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,\
                                fill_value = np.NaN, \
                                mapping_method = 'bin_average',
                                radius_of_influence = 40000)

# plot the whole field, this time land values are nans.
plt.figure(figsize=(6,6), dpi= 90)
plt.imshow(field_nearest_1deg,origin='lower',cmap='jet');
_images/ECCO_v4_Interpolating_Fields_to_LatLon_Grid_30_0.png

You may wonder why there are missing (white) values in the resulting map. It’s because there are some lat-lon grid cells whose center is more than 40km away from the center of any grid cell on the source grid. There are ways to get around this problem by specifying spatially-varying radii of influence, but that will have to wait for another tutorial.

If you want to explore on your own, explore some of the low-level routines of the pyresample package: https://pyresample.readthedocs.io/en/latest/

Interpolating ECCO vectors fields to lat-lon grids

Vector fields require a few more steps to interpolate to the lat-lon grid. At least if what you want are the zonal and meridional vector components. Why? Because in the llc90 grid, vector fields like UVEL and VVEL are defined with positive in the +x and +y directions of the local coordinate system, respectively. A term like oceTAUX correponds to ocean wind stress in the +x direction and does not correspond to zonal wind stress. To calculate zonal wind stress we need to rotate the vector from the llc90 x-y grid system to the lat-lon grid.

To demonstrate why rotation is necessary, consider the model field oceTAUX

[20]:
## Load vector fields
day_mean_dir= ECCO_dir / 'nctiles_monthly/'

ecco_vars = ecco.recursive_load_ecco_var_from_years_nc(day_mean_dir, \
                                           vars_to_load=['oceTAUX', 'oceTAUY'], \
                                           years_to_load=2000,less_output=True)

ecco_ds = []

## Merge the ecco_grid with the ecco_vars to make the ecco_ds
ecco_ds = xr.merge((ecco_grid , ecco_vars)).load()

pprint(ecco_ds.data_vars)
loading files of  oceTAUX
loading files of  oceTAUY
Data variables:
    oceTAUX  (time, tile, j, i_g) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    oceTAUY  (time, tile, j_g, i) float32 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
[21]:
tmp = ecco_ds.oceTAUX.isel(time=0)
tmp_masked = np.where(ecco_ds.maskC.isel(k=0)==1, tmp, np.nan)
ecco.plot_tiles(tmp_masked);
_images/ECCO_v4_Interpolating_Fields_to_LatLon_Grid_34_0.png

We can sees the expected positive zonal wind stress in tiles 0-5 because the x-y coordinates of tiles 0-5 happen to approximately line up with the meridians and parallels of the lat-lon grid. Importantly, for tiles 7-12 wind stress in the +x direction corresponds to mainly wind stress in the SOUTH direction. To see the zonal wind stress in tiles 7-12 one has to plot oceTAUY and recognize that for those tiles positive values correspond with wind stress in the tile’s +y direction, which is approximately east-west. To wit,

[22]:
tmp = ecco_ds.oceTAUY.isel(time=0)
tmp_masked = np.where(ecco_ds.maskC.isel(k=0)==1, tmp, np.nan)
ecco.plot_tiles(tmp_masked);
_images/ECCO_v4_Interpolating_Fields_to_LatLon_Grid_36_0.png

Great, but if you want zonal and meridional wind stress, we need to determine the zonal and meridional components of these vectors.

Vector rotation

For vector rotation we leverage the very useful XGCM package (https://xgcm.readthedocs.io/en/latest/).

Note: The XGCM documentation contains an MITgcm ECCOv4 Example Page: https://xgcm.readthedocs.io/en/latest/example_eccov4.html. In that example the dimension ‘tile’ is called ‘face’ and the fields were loaded from the binary output of the MITgcm, not the netCDF files that we produce for the official ECCOv4r4 product. Differences are largely cosmetic.)

We use XGCM to map the +x and +y vector components to the grid cell centers from the u and v grid points of the Arakawa C grid. The ecco-v4-py routine get_llc_grid creates the XGCM grid object using a DataSet object containing the following information about the model grid:

i,j,i_g,j_g,k,k_l,k_u,k_p1 .

Our ecco_ds DataSet does have these fields, they come from the ECCO-GRID.nc object:

[23]:
# dimensions of the ecco_ds DataSet
ecco_ds.dims
[23]:
Frozen(SortedKeysDict({'k_p1': 51, 'j_g': 90, 'i_g': 90, 'k': 50, 'j': 90, 'k_u': 50, 'i': 90, 'k_l': 50, 'tile': 13, 'time': 12, 'nv': 2}))
[24]:
# Make the XGCM object
XGCM_grid = ecco.get_llc_grid(ecco_ds)

# look at the XGCM object.
XGCM_grid

# Depending on how much you want to geek out, you can learn about this fancy XGCM_grid object here:
# https://xgcm.readthedocs.io/en/latest/grid_topology.html

[24]:
<xgcm.Grid>
T Axis (not periodic):
  * center   time
X Axis (not periodic):
  * center   i --> left
  * left     i_g --> center
Y Axis (not periodic):
  * center   j --> left
  * left     j_g --> center
Z Axis (not periodic):
  * center   k --> left
  * left     k_l --> center
  * outer    k_p1 --> center
  * right    k_u --> center

Once we have the XGCM_grid object, we can use built-in routines of XGCM to interpolate the x and y components of a vector field to the cell centers.

[25]:
import xgcm
xfld = ecco_ds.oceTAUX.isel(time=0)
yfld = ecco_ds.oceTAUY.isel(time=0)

velc = XGCM_grid.interp_2d_vector({'X': xfld, 'Y': yfld},boundary='fill')

velc is a dictionary of the x and y vector components taken to the model grid cell centers. At this point they are not yet rotated!

[26]:
velc.keys()
[26]:
dict_keys(['X', 'Y'])

The magic comes here, with the use of the grid cosine ‘cs’ and grid sine ‘cs’ values of the ECCO-GRID object:

[27]:
# Compute the zonal and meridional vector components of oceTAUX and oceTAUY
oceTAU_E  = velc['X']*ecco_ds['CS'] - velc['Y']*ecco_ds['SN']
oceTAU_N  = velc['X']*ecco_ds['SN'] + velc['Y']*ecco_ds['CS']

Now we have the zonal and meridional components of the vectors, albeit still on the llc90 grid.

[28]:
oceTAU_E_masked = np.where(ecco_ds.maskC.isel(k=0)>0, oceTAU_E, np.nan)
ecco.plot_tiles(oceTAU_E_masked);
_images/ECCO_v4_Interpolating_Fields_to_LatLon_Grid_48_0.png

So let’s resample oceTAU_E to a lat-lon grid (that’s why you’re here, right?) and plot

[29]:
new_grid_delta_lat = .5
new_grid_delta_lon = .5

new_grid_min_lat = -90+new_grid_delta_lat/2
new_grid_max_lat = 90-new_grid_delta_lat/2

new_grid_min_lon = -180+new_grid_delta_lon/2
new_grid_max_lon = 180-new_grid_delta_lon/2

new_grid_lon, new_grid_lat, oceTAU_E_masked_latlon =\
        ecco.resample_to_latlon(ecco_ds.XC, \
                                ecco_ds.YC, \
                                oceTAU_E_masked,\
                                new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,\
                                new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,\
                                fill_value = np.NaN, \
                                mapping_method = 'nearest_neighbor',
                                radius_of_influence = 120000)

# plot the whole field, this time land values are nans.
plt.figure(figsize=(12,8), dpi= 90);
plt.imshow(oceTAU_E_masked_latlon,origin='lower',vmin=-.2,vmax=.3,cmap='jet');
plt.title('ECCOv4 Jan 2000 Zonal Wind Stress (N/m^2)');
plt.colorbar(orientation='horizontal');
_images/ECCO_v4_Interpolating_Fields_to_LatLon_Grid_50_0.png

Or with extra fanciness:

[30]:
plt.figure(figsize=(12,5), dpi= 90)

ecco.plot_proj_to_latlon_grid(ecco_ds.XC, ecco_ds.YC, \
                              oceTAU_E_masked, \
                              user_lon_0=180,\
                              projection_type='PlateCaree',\
                              plot_type = 'pcolormesh', \
                              dx=1,dy=1,show_colorbar=True,cmin=-.2, cmax=0.3);
_images/ECCO_v4_Interpolating_Fields_to_LatLon_Grid_52_0.png

Compare vs “The Scatterometer Climatology of Ocean Winds (SCOW)”,

https://nanoos.ceoas.oregonstate.edu/scow/index.html

Risien, C.M., and D.B. Chelton, 2008: A Global Climatology of Surface Wind and Wind Stress Fields from Eight Years of QuikSCAT Scatterometer Data. J. Phys. Oceanogr., 38, 2379-2413.

image0

And of course we can’t forget about our meridional wind stress:

[31]:
oceTAU_N_masked = np.where(ecco_ds.maskC.isel(k=0)>0, oceTAU_N, np.nan)

plt.figure(figsize=(12,6), dpi= 90)
ecco.plot_proj_to_latlon_grid(ecco_ds.XC, ecco_ds.YC, \
                              oceTAU_N_masked, \
                              user_lon_0=180,\
                              projection_type='PlateCaree',\
                              plot_type = 'pcolormesh', \
                              dx=1,dy=1,cmin=-.2, cmax=0.3,show_colorbar=True);
_images/ECCO_v4_Interpolating_Fields_to_LatLon_Grid_56_0.png

Which we also compare against SCOW:

image0

UEVNfromUXVY

The ecco-v4-py library includes a routine, UEVNfromUXVY which does the interpolation to the grid cell centers and the rotation in one call:

[32]:
xfld = ecco_ds.oceTAUX.isel(time=0)
yfld = ecco_ds.oceTAUY.isel(time=0)

# Compute the zonal and meridional vector components of oceTAUX and oceTAUY
oceTAU_E, oceTAU_N  = ecco.vector_calc.UEVNfromUXVY(xfld, yfld, ecco_ds)

Before interpolating to the lat-lon grid it is convenient to mask out the land values

[33]:
# mask the rotated vectors
oceTAU_E=oceTAU_E.where(ecco_ds.maskC.isel(k=0))
oceTAU_N=oceTAU_N.where(ecco_ds.maskC.isel(k=0))

Now plot to verify

[34]:
# interpolate to lat-lon
new_grid_lon, new_grid_lat, oceTAU_E_masked_latlon =\
        ecco.resample_to_latlon(ecco_ds.XC, \
                                ecco_ds.YC, \
                                oceTAU_E,\
                                new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,\
                                new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,\
                                fill_value = np.NaN, \
                                mapping_method = 'nearest_neighbor',
                                radius_of_influence = 120000)

# plot the whole field, this time land values are nans.
plt.figure(figsize=(12,8), dpi= 90);
plt.imshow(oceTAU_E_masked_latlon,origin='lower',vmin=-.2,vmax=.3,cmap='jet');
plt.title('ECCOv4 Jan 2000 Zonal Wind Stress (N/m^2)');
plt.colorbar(orientation='horizontal');
_images/ECCO_v4_Interpolating_Fields_to_LatLon_Grid_64_0.png

Saving interpolated fields to netCDF

For this demonstration we will rotate the 12 monthly-mean records of oceTAUX and oceTAUY to their zonal and meridional components, interpolate to a lat-lon grids, and then save the output as netCDF format.

[35]:
pprint(ecco_ds.oceTAUX.time)
<xarray.DataArray 'time' (time: 12)>
array(['2000-01-16T12:00:00.000000000', '2000-02-15T12:00:00.000000000',
       '2000-03-16T12:00:00.000000000', '2000-04-16T00:00:00.000000000',
       '2000-05-16T12:00:00.000000000', '2000-06-16T00:00:00.000000000',
       '2000-07-16T12:00:00.000000000', '2000-08-16T12:00:00.000000000',
       '2000-09-16T00:00:00.000000000', '2000-10-16T12:00:00.000000000',
       '2000-11-16T00:00:00.000000000', '2000-12-16T12:00:00.000000000'],
      dtype='datetime64[ns]')
Coordinates:
    timestep  (time) int64 70860 71556 72300 73020 ... 76692 77436 78156 78900
  * time      (time) datetime64[ns] 2000-01-16T12:00:00 ... 2000-12-16T12:00:00
Attributes:
    long_name:  center time of averaging period
    bounds:     time_bnds
    axis:       T

We will loop through each month, use UEVNfromUXVY to determine the zonal and meridional components of the ocean wind stress vectors, and then interpolate to a 0.5 degree lat-lon grid.

[36]:
oceTAUE = np.zeros((12, 360,720))
oceTAUN = np.zeros((12, 360,720))

new_grid_delta_lat = .5
new_grid_delta_lon = .5

new_grid_min_lat = -90+new_grid_delta_lat/2
new_grid_max_lat = 90-new_grid_delta_lat/2

new_grid_min_lon = -180+new_grid_delta_lon/2
new_grid_max_lon = 180-new_grid_delta_lon/2


for m in range(12):
    cur_oceTAUX = ecco_ds.oceTAUX.isel(time=m)
    cur_oceTAUY = ecco_ds.oceTAUY.isel(time=m)

    # Compute the zonal and meridional vector components of oceTAUX and oceTAUY
    tmp_e, tmp_n  = ecco.vector_calc.UEVNfromUXVY(cur_oceTAUX, cur_oceTAUY, ecco_ds)

    # apply landmask
    tmp_e_masked = np.where(ecco_ds.maskC.isel(k=0)>0, tmp_e, np.nan)
    tmp_n_masked = np.where(ecco_ds.maskC.isel(k=0)>0, tmp_n, np.nan)

    # zonal component
    new_grid_lon, new_grid_lat, tmp_e_masked_latlon =\
        ecco.resample_to_latlon(ecco_ds.XC, \
                                ecco_ds.YC, \
                                tmp_e_masked,\
                                new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,\
                                new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,\
                                fill_value = np.NaN, \
                                mapping_method = 'nearest_neighbor',
                                radius_of_influence = 120000)

    # meridional component
    new_grid_lon, new_grid_lat, tmp_n_masked_latlon =\
        ecco.resample_to_latlon(ecco_ds.XC, \
                                ecco_ds.YC, \
                                tmp_n_masked,\
                                new_grid_min_lat, new_grid_max_lat, new_grid_delta_lat,\
                                new_grid_min_lon, new_grid_max_lon, new_grid_delta_lon,\
                                fill_value = np.NaN, \
                                mapping_method = 'nearest_neighbor',
                                radius_of_influence = 120000)

    oceTAUE[m,:] = tmp_e_masked_latlon
    oceTAUN[m,:] = tmp_n_masked_latlon
[37]:
# make the new data array structures for the zonal and meridional wind stress fields
oceTAUE_DA = xr.DataArray(oceTAUE,  name = 'oceTAUE',
                      dims = ['time','latitude','longitude'],
                      coords = {'latitude': new_grid_lat[:,0],
                                'longitude': new_grid_lon[0,:],
                                'time': ecco_ds.time})

# make the new data array structures for the zonal and meridional wind stress fields
oceTAUN_DA = xr.DataArray(oceTAUN,  name = 'oceTAUN',
                      dims = ['time','latitude','longitude'],
                      coords = {'latitude': new_grid_lat[:,0],
                                'longitude': new_grid_lon[0,:],
                                'time': ecco_ds.time})

Plot the zonal and meridional wind stresses in these new DataArray objects:

[38]:
oceTAUE_DA.dims
[38]:
('time', 'latitude', 'longitude')
[39]:
oceTAUE_DA.isel(time=0).plot()
[39]:
<matplotlib.collections.QuadMesh at 0x1ff001ebee0>
_images/ECCO_v4_Interpolating_Fields_to_LatLon_Grid_72_1.png
[40]:
oceTAUN_DA.isel(time=0).plot()
[40]:
<matplotlib.collections.QuadMesh at 0x1ff00c6d730>
_images/ECCO_v4_Interpolating_Fields_to_LatLon_Grid_73_1.png

Saving is easy with xarray

[41]:
# meridional component
output_path = Path('C:/Users/Ian/Downloads/oceTAUN_latlon_2000.nc')
oceTAUN_DA.to_netcdf(output_path)

# zonal component
output_path = Path('C:/Users/Ian/Downloads/oceTAUE_latlon_2000.nc')
oceTAUE_DA.to_netcdf(output_path)

… done!