Roms grid rotation and interpolation

[stef]

I'm trying to configure a model to practice data assimilation, but stumbled across a very basic question about interpolation (maybe embarrassing for a geophysicist , but I never had to do with that problem before.).

When I interpolate data on an equidistant lat/lon grid, do I need to do some spatial weighting to account for the sphericity? If yes, is this consistently done in the ROMS code, e.g. in the interpolation of wind stress data? I'm asking because I would think that in all the OA, OI and 3/4dvar routines, this is accounted for? Can you briefly comment on where one needs to worry about that, and where not?

Thanks a lot for your help!

[rduran]

It seems the key word is conservative interpolation or conservative remapping, I hope ROMS folks can chip in re what ROMS does when interpolating but in general this paper answers your question

http://journals.ametsoc.org/doi/pdf/10. ... 2.0.CO%3B2

The software implementing the algorithms in that paper is available in SCRIP
http://oceans11.lanl.gov/trac/SCRIP

I also hope others can chip in on how much / when one should worry about this.

[stef]

Very good, thanks! Meanwhile I also found some Python code, specific to ROMS:

https://github.com/ESMG/pyroms/blob/mas ... g/remap.py

with "conservative" and "bicubic" remapping. It seems they use the "SCRIP" library you linked to, so that starts to make sense to me Cheers, Stefan

EDIT: One more thing: I saw that ROMS does the bicubic interpolation of wind stress per default (subroutine 'cinterp2d', this has been discussed on other pages as well). Whether this is equivalent to the SCRIP version, I haven't found out yet. In particular, I'm not sure whether the ROMS version takes into account the sphericity. Maybe someone could confirm that. If it does, wouldn't it be handy to isolate that function (perhaps re-write it in Matlab/Python) and re-use it e.g. for bathymetry grid rotation?

[stef]

After spending some time today browsing through the SCRIP manual, I think I should better explain the motivation of my initial question. If I understand correctly none of the bilinear and bicubic methods take into account the sphericity, neither in SCRIP nor in ROMS. What they do is interpolate in a "logically rectangular" coordinate system (defined on a plane) onto points that don't coincide with the grid. While that's interesting and surely useful (I'd think for curvilinear coordinates), my question actually relates to whether the "initial" projection from the sphere to a plane has to be accounted for in the interpolation.

To make it clearer, consider e.g. the simplest type of "realistic" grids we use, i.e. lat-lon grids with equidistant spacing and dlat=dlon. These are not only "logically rectangular", they are practically rectangular too (after the projection from the sphere). So the procedure of determining the (i,j) coordinates of a point specified by a (lon,lat) pair is trivial (?), since both coordinate systems are equivalent. Lets say I want to interpolate onto a point that's exactly in the center of the (square) grid cell, I would do (I suppose this is bilinear):

P_c = (1/4)*(P_sw + P_se + P_ne + P_nw),

where P_c is the center point, and P_sw is the south-western corner, etc. However, in reality (on the sphere),

1) the distance between the two poleward grid points is less than the distance between the two equatorward grid points
2) Because of that, P_c is closer to the two poleward grid points

Is this true or does this not make sense? I found a post on the matlab forum
http://au.mathworks.com/matlabcentral/n ... read/58326
asking a related question, and somebody mentioned the term "spherical harmonics". Not sure what that means.

[kate]

It's true that we don't worry about this stuff, just let scrip do it's bilinear thing on a lat-lon grid. We're interpolating from fields with errors, our numerics have errors. I doubt that interpolation biases are going to be noticeable compared to all our other problems. You can make yourself crazy thinking about things like volume of the grid cell from all the bathymetry data - what's the best way to get a model bathymetry that's not too steep, conserves volume, doesn't destroy shelf-break currents, etc. Then do that for all your other data sources. I have first-order problems like mapping hydrology output to ROMS river sources while not losing the river mouths that don't line up with the ROMS coastline, trying to get the same volume flux. Or inventing BGC fluxes for all those rivers.

[stef]

Ok thanks a lot for your time Kate. Let me see if I understand correctly:

The real question to ask is what to use for pm and pn in the first place. This reflects how the projection is made from the sphere to the plane. It already involves a trade off: If one chooses Mercator, one doesn't conserve volume. If one chooses actual spherical distance to calculate pm and pn, then one doesn't get an orthogonal grid in two dimensional space.

Is this right so far, or not? (If yes, is there a rule of thumb which one to prefer? I guess not - see this post: viewtopic.php?f=14&t=801 )

If yes, then it may be inconsistent to *afterwards* worry about sphericity any longer anyways.

Does this sound reasonable?

[kate]

For pm, pn, I have an initial value based on the flat projection space. Then I later remap to the sphere so that I can assign lat,lon values to each point. At that time I recompute pm, pn based on the lat, lon values and distances on a sphere - with or without an ellipsoidal correction. As long as I use a conformal (angle-preserving) projection, it should still be orthogonal back on the sphere.

[stef]

Thanks for the help! Now I also found a presentation , where you explain this on page 4 (may be helpful for other readers).

However, strictly speaking, when we measure pm, pn on the sphere and use it in the horizontal plane, we don't get an orthogonal grid, right? So if we computed the ds quantities corresponding to dxi and deta (defined here), and tried to assembled them on a horizontal plane, the resulting grid would be non-orthogonal. (?)

To summarize, could one say that this is "the price we are willing to pay" to keep the distances (and in particular the grid-cell volumes) realistic? The latter is really the motivation for this approach, right?

[kate]

There are some map projections in which angles are preserved, called "conformal". Those are the ones we can use to have an orthogonal grid on a sphere.

[stef]

Ok, but ROMS doesn't care whether it's orthogonal on the sphere, right? What it wants is orthogonality in the plane. And this cannot be fulfilled if we measure pm and pn on the sphere and then use it on the plane. Is this right? Or is there a conformal projection that can do that?

The way I see it, is if I insist on an orthogonal grid on the sphere, I have the choice of either breaking orthogonality in the plane, or "breaking" pm and pn (or not using a conformal projection in the first place)?

In *theory*, couldn't one use an area-preserving projection, like Lambert cylindrical equal-area projection (here). Looking at the Wikipedia plot, that would obviously only work in the tropics and subtropics. The inverse map projection would yield a non-orthogonal grid on the sphere. But the algorithms to adjust "true" North etc. are there already anyways from the curvilinear grids? And the plots would look weird too. But the grids could be strictly orthogonal in the "ROMS-plane" and have the true pm an pn, i.e. volume conserving. For the plots maybe one could use a different projection. I guess that sounds weird, but is there a logical reason this wouldn't work?

[kate]

ROMS wants an orthogonal grid, either on the sphere or on a plane - it doesn't care which. If you want your grids to line up with lat,lon lines, the spacing doesn't matter - it will be orthogonal. The reason to map to flat space using a conformal map projection is if you want a curvilinear grid. The tools for creating curvilinear, orthogonal grids generally work on the complex plane, i.e. flat space.

[stef]

Very good, thanks for your patience and the information, it helped a lot. I will try to process all of this during the next couple of weeks

[stef]

figure_grid_rotation_tutorial_mercator.png (236.82 KiB) Viewed 23172 times

I'd like to show off a figure about one way to correctly (?) rotate the grid, where I'm trying to follow the steps from above post:

For pm, pn, I have an initial value based on the flat projection space. Then I later remap to the sphere so that I can assign lat,lon values to each point. At that time I recompute pm, pn based on the lat, lon values and distances on a sphere - with or without an ellipsoidal correction. As long as I use a conformal (angle-preserving) projection, it should still be orthogonal back on the sphere.

The objective is to obtain a ROMS grid aligned to the US west coast. The colours in the figure show depth (m) of a subsample of ETOPO2 in the North-East pacific. The bottom (left) axis shows the x (y) axis of the flat projection space (Mercator). The right (top) axis show latitude (longitude). Note the spacing between 5-degree increments of latitude increases polewards. This is characteristic of the Mercator projection.

The black grid is an equidistantly-spaced grid in the Mercator-plane. The red grid is a rotation thereof, again in Mercator-coordinates. Rotations preserve angles, so I guess they are trivial conformal maps (?). However, the rotation is an operation in the *plane*. Clearly, the circles that get rotated to the north encompass a smaller area on the sphere (i.e. after the inverse projection), but that doesn't matter; The only thing one cares about is that the red grid is still rectangular on the sphere. This way, ROMS' curvilinear coordinates xi and eta can be interpreted as *spherical* coorinates (in the non-rotated case that would be achieved by setting pm=1/(dlon*R*cos(lat)) and pn=1/(dlat*R)). That the grid is orthogonal on the sphere is guaranteed by the fact that one used

1) a conformal map to project into the plane
2) a rotation in the plane (which trivially preserves angles)
2) a conformal map to project back onto the sphere.

After the rotation, one has to interpolate the lat/lon data onto the red Mercator-plane grid. Then we get the great-circle distances for pm, pn for the rotated grid (e.g. from the Haversine formula).

Let me know in case it's still not right.

[kate]

Yes, that's all correct. Mercator is fine until you get to high latitudes. For mid-latitudes, a better projection is the Lambert conformal conic projection. We have a Northeast Pacific grid using the latter. Plus I'm sure you've heard of polar stereographic.

[stef]

Wow, Lambert conformal conic! On Wikipedia it says:

By scaling the resulting map, two parallels can be assigned unit scale, with scale decreasing between the two parallels and increasing outside them. This gives the map two standard parallels. In this way, deviation from unit scale can be minimized within a region of interest that lies largely between the two standard parallels.

I noticed that you even plotted it this way in the presentation (link above), i.e. the cone "cuts" through the sphere somewhere in the midlatitudes. So the advantage over Mercator is that (1) the lines of unit scale always coincide with parallels, and (2) can be chosen arbitrarily.

So with Mercator, if I wanted to have unit lines coincide with parallels, I'm restricted to choose them symmetrically around the equator, say 30deg North/South. Otherwise, if one tilted the "Mercator cylinder", the two standard lines wouldn't coincide with parallels any more? I think it makes sense, although I'm not sure if it's crucial to have the "unit lines" coincide with parallels in a modelling context?

But anyways, regarding modelling grids, the take-home message is that one can "tweak" the projection so that an equidistant grid in the projection plane yields approximately uniform grid-cell sizes on the sphere. To do that one can choose *two* "unit lines" (or circle in the stereographic projection) according to the extent/shape of the grid.

I think this relates closely to my initial question about interpolation. Perhaps not within a single grid cell, but across multiple cells within the domain.

Excellent, this is becoming much clearer now! Thanks again.

[stef]

It turns out that "oblique Mercator projection" is exactly what I need, I think. It lets you define two points in lat,lon coordinates through which the central line for the cylindrical projection goes (that's the line at which the projection cylinder touches the sphereoid), and this can be tilted with respect to parallels (hence the name "oblique"). That should be great for coastlines with a direction that is not predominantly zonally/meridionally aligned. Moreover, if I understand correctly the projection itself does the rotation for you, and the output is a coordinate system who's X-axis is aligned with the oblique central line.

Some links in case anyone is interested in further details:

basemap/proj4: Description (page 13), and the C code. It seems this is ellipsoidal only, but I'm not sure.

For the formulas for the spherical version, see

Snyder, John Parr. Map projections--A working manual. Vol. 1395. US Government Printing Office, 1987.
Link to Goolge books.

See particularly their page 69 for the spherical formulas. Maybe I'm going to write this down in Python to avoid a basemap dependency. It seems one can easily obtain the "secant" version by setting a scale factor to something other than 1.

This link has nice figures including oblique cylindrical projections.

[stef]

Just a quick warning in case anybody (like me) considers making oblique Mercator projections with Proj version (4.9.3): be careful