Remapping ocean vector quantities
DamienIrving opened this issue · 5 comments
I'm working with ocean heat flux (hf) and surface wind stress (tau) data from the CMIP5 project. These are vector quantities (i.e. there are hfx/tauuo and hfy/tauvo components) that are provided on the native model grid, which is often a curvilinear or rotated pole grid.
I'm particularly interested in the meridional ocean heat flux and the surface zonal wind stress. The complication in calculating these quantities is that the x and y directions on a curvilinear or rotated pole grid are not the same as the x (zonal) and y (meridional) directions on a geographic latitude/longitude grid. This means you can't simply apply existing software tools (see #2) that are used for remapping scalar quantities from an ocean model grid to geographic grid. Any remapping requires the use of both the x and y components.
How do people go about remapping vector quantities from an ocean model grid to a regular geographic grid?
I contacted the authors of ESMPy (the Python interface to the Earth System Modeling Framework regridding utility) and they said they don’t currently have the capability to interpolate vectors taking different orientations of the source and destination grids into account.
However, that being said, they are currently working on a way to do that for a project they are helping with. What they are doing is converting the vectors from the local direction to a consistent underlying basis (3D Cartesian), doing the regridding, and then converting the resulting vector back again to the local directions in the destination grid. They are talking about bringing that capability into ESMF, bu don’t have a timeframe for when or even if that might happen.
The details of the process are as follows:
To do the regridding we first project the local vectors to 3D Cartesian space. To do that we multiply the magnitude of the local vector by a 3D Cartesian unit vector pointing in the same direction and then sum the two Cartesian vectors representing the local vectors together to give one 3D Cartesian vector representing the direction. In our case we have a north and east vector, so this transformation looks like this:
C = N * Un + E * Ue
Where:
C is the 3D Cartesian vector
N, E are the north and east vector magnitudes
Un, Ue are the north and east Cartesian unit vectors
And:
Un = [-sin(lat)sin(lon), sin(lat)cos(lon), cos(lat)]
Ue = [cos(lon), sin(lon), 0]
Where:
lat,lon are the latitude and longitude of the location of the N and E vectors
(The above might be different for you if you have a different idea of what the local vectors are. If so, Un and Ue should be 3D Cartesian unit vectors in your particular local directions.)
Once you have a 3D Cartesian vector you regrid it as usual. If you put the 3 vector components first in memory in the Field, ESMF_FieldRegrid() will regrid all three at once. (e.g. for a 2D grid creating a Field like this would look like: cart_vec_field=ESMF_FieldCreate(grid,typekind=ESMF_TYPEKIND_R8, gridToFieldMap=(/2,3/), ungriddedLBound=(/1/), ungriddedUBound=(/3/), rc=rc))
This regridding yields C’ (the Cartesian vector in the new location). After regridding you project it back to the local vector representation:
N' = C' . Un'
E' = C' . Ue'
(note that the “.” indicates a dot product)
Where:
C' is the regridded 3D Cartesian vector
N’,E' are the north and east vector magnitudes in the new location
Un', Ue' are the north and east Cartesian unit vectors in the new location.
An alternative to writing a vector regridding package from scratch (using the advice from ESMPy above) could possibly be to use the iris.analysis.cartography.rotate_winds function. I started an issue to ask about this, but I suspect it won't work (you have to specify the details of the input grid in terms of an iris.coord_system, which isn't possible).
For the surface wind stress, it looks like in most cases people just use the atmospheric variable (tauu; which is on a lat/lon grid) instead of the ocean variable (tauuo).
http://journals.ametsoc.org/doi/full/10.1175/JCLI-D-12-00591.1
Ryan Abernathey commented on twitter that it's hard to ensure that the remapped vector has the same conservation properties as it does on its native grid. For this reason, many people prefer to always work on the native grid. It's better to interpolate frame-invariant scalar quantities, such as divergence, than vector ones.
This motivates his preferred way to interpolate 2D (u, v) vectors: use helmholtz decomposition to write vector in terms of streamfunction ψ and potential φ. Then interpolate ψ and φ to new grid, and take gradients on new grid to recover u and v.
From He et al (2019): "Since the CCSM4 ocean model is in curvilinear grids, we re-grid the rotated Eulerian velocity and temperature to high-resolution Cartesian grids..."
They use pyresample for this - their code can be found at https://github.com/Yefee/xcesm/blob/master/xcesm/core/xcesm.py, line 435 onwards.