Stefan's website
\( \renewcommand{\b}[1]{\symbf{#1}} \)

Grid generation for numerical ocean models

Riha, S.

1 Spherical, Curvilinear or Non-curvilinear?

ROMS is formulated in curvilinear coordinates. Curvilinear coordinates can be used to create boundary-following coordinate systems. However, the most basic and common application of curvilinear coordinates, is to formulate the equations of motion on the surface of a sphere. The use of curvilinear coordinates is mandatory in this case, independent of whether the projected planar geometry follows any coastlines or not. For example, Haidvogel et al. (2000) state (in their section 2.1) that

Note that [ROMS’] general formulation of curvilinear coordinates includes Cartesian coordinates [...] as well as spherical coordinates [...].

It is interesting to note that some experienced regional and coastal ocean modellers seem to avoid boundary-following grids altogether (see e.g. Arango 2006). For a realistic model with a grid that is not boundary-following, one has the choice between:

  1. Interpreting the modeling domain as (part of) a sphere by using curvilinear coordinates.
  2. Interpreting the modeling domain as a flat plane and use conventional Cartesian coordinates. This implies that one is willing to accept any errors (or, perhaps more appropriately, lack of realism), which are unavoidably caused by conformal (and thus area-distorting) projections.

In both cases, the final grids need to be orthogonal. Suppose one uses a conformal projection, for example Mercator, and sets up a perfectly rectangular grid in the projection plane, i.e. the spacing between two adjacent grid points in x- and y- direction is constant. Then the orthogonality for option 1 is obtained by measuring the grid spacing that enters the mathematical model along great circles of the sphere. To this end, one needs the geographical coordinates (latitude and longitude) of the grid points. On the other hand, for option 2, the spacing which enters the mathematical model is set to a constant, that is the grid spacing of the flat projection space. Note that with option 1, the actual cell area of the mathematical model will vary with latitude, since the Mercator projection (as most, if not all conformal projections) is not area preserving. With option 2, the cell area is constant in the entire domain, yet there is necessarily a lack of realism. But how large are the differences?

2 Mercator projection

To quantify the differences, Fig. 1 shows a Mercator projection of ETOPO2 bathymetry (U.S. National Geophysical Data Center, 2016) in the Southern California Bight. The coastline of the bight is oriented neither along a meridian nor a parallel, so it makes sense to use a rotated grid here, whose extent is outlined by the black rectangle. Note that the actual grid boundaries shown in the figure may not make sense from a modeling viewpoint (but see e.g. Kumar et al. 2015), but are chosen arbitrarily for this example.

Mercator projection
Figure 1: Mercator projection of bathymetry in the Southern California Bight. Coast lines are shown in red. The black rectangle shows the outline of a rotated grid, which could be used as the computational domain of a numerical ocean model. Elevations larger than sea level have been clipped.

To set up a curvilinear grid, we follow kate’s suggestion in Hedstrom et al. (2017), and first choose a grid spacing in the projection plane. Let’s say dx=dy=3km. With this constant spacing, we set up a grid which is initially aligned with the coordinate axis. To rotate the grid, the grid point coordinates are multiplied with a rotation matrix. To get geographical coordinates, we use proj ( Evenden 1990; Evenden 2005) to convert the projection-plane coordinates of the grid-points to latitude-longitude coordinates with an inverse Mercator projection. After that, we use geod (Evenden, 1990) to get the geodetic distances between the latitude-longitude coordinates of the grid points. Fig. 2 shows the grid spacing in ξ direction. proj let’s you tweak the (inverse) Mercator projection to have an arbitrary latitude of true scale. Here it is located roughly in the center of the domain. Along this latitude, the spherical grid spacing corresponds exactly to the constant spacing of dx=3km in the planar projection plane.

Cell width of a Mercator projection
Figure 2: Grid spacing in ξ direction. Projections and geodetic calculations were obtained from proj and geod using the WGS84 spheroid.

To make this more obvious, Fig. 3 shows the relative deviation of the magnitude of the grid spacing from dx=3km. According to these figures, the variation of the spacing is between two orders of magnitude lower than the actual values.

Deviation of the cell width from the mean (Mercator)
Figure 3: Relative deviation (percent) of the magnitude of grid spacing in ξ direction from dx=3km.

3 Oblique Mercator projection

The Oblique Mercator projection is a cylindrical and conformal projection just like Mercator, but with a tilted (oblique) projection cylinder. How does it differ from an ordinary Mercator projection?

  • Ordinary Mercator: The central line (that’s the great circle line along which the cylinder is exactly tangent to the sphere) of an ordinary Mercator projection is always at the equator. This means that the region of low variation of the scaling factor (that’s the region around the central line) is also around the equator. Although we can adjust the radius of the cylinder to obtain a true scale at any latitude, there is no way to avoid the increasingly strong variation of the scaling away from the equator.
  • Oblique Mercator: For the Oblique Mercator projection, one can choose an arbitrary great circle line to be the central line. For the above example grid, we choose the great circle to have an angle of about 60° in anti-clockwise direction from North (i.e. an azimuth of about -60°), resulting in the projection shown in Fig. 4
Oblique Mercator projection
Figure 4: Oblique Mercator projection of bathymetry in the Southern California Bight. The central line of the projection intersects with the origin of the coordinate system, at an angle of roughly 60° in anti-clockwise direction from North. Elevations larger than sea level have been clipped.

Given the relatively small region, Fig. 4 looks very similar to Fig. 1, but note the slight differences along the edges of the figures. The plotted data is extracted from ETOPO2, and is gridded on a perfect rectangle in latitude-longitude space. While the data boundaries (meridians and parallels) exactly align with the planar coordinate system in Fig. 1, this is obviously not the case in Fig. 4. In fact, all meridians and parallels are generally curved lines in an Oblique Mercator projection. This is important because the angle between grid-cell boundaries and true North varies in a grid generated with the Oblique Mercator projection. This has to be accounted for when further configuring ROMS.

Let's take a look of the variation of the grid spacing. Fig. 5 and Fig. 6 show that the variation is much smaller than in the Mercator projection. According to Fig. 6, the variation of the spacing is now three to four orders of magnitude smaller than the actual spacing.

Cell width of an oblique Mercator projection
Figure 5: Grid spacing in ξ direction for the Oblique Mercator projection.
Deviation of the cell width from the mean (oblique Mercator)
Figure 6: Relative deviation (percent) of the magnitude of grid spacing in ξ direction from dx=3km.

It is possible to decrease the variation even further by using a secant projection in which the radius of the projection cylinder is slightly decreased, such that it cuts through the spheroid. Fig. 7 shows a secant projection with proj's parameter k set to k=0.99993. Note that there are now two swaths where the grid spacing is equal to dx=3km. In between these swaths the spacing is slightly larger, outside it is slightly lower (the smaller/larger regions are not discernable in the figure since absolute values of deviation are shown).

Deviation of cell width with modified secant parameter
Figure 7: Relative deviation (percent) of the magnitude of grid spacing in ξ direction from dx=3km for a secant Oblique Mercator projection with k=0.99993.

As noted above, latitude/longitude grids appear as curved lines on an Oblique Mercator projection. Fig. 8 shows the angle between the ξ-axis and true East. The variation is roughly an order of magnitude lower than the average value.

Variation of cell boundary orientation with respect to true directions for an oblique Mercator projection
Figure 8: Angle (degrees) between the XI-axis and true East.

4 Discussion

Are there any advantages in tuning the projection? I don’t know yet. Some questions and points for consideration:

  • Arguably the most important thing to keep in mind is that the simulated flow should be independent of grid topology (see e.g. Arango, 2006).
  • If the modeling domain is interpreted as part of a sphere (or spheroid), which implies that the curvilinear terms of the numerical model are used, is there any advantage in having a grid with relatively uniform cell sizes? The answer is not obvious to me, but I would think that it certainly isn't a drawback. If it is indeed an advantage, then this means that in extratropical regions, a (conformal) projection other than standard Mercator must be used, with minimal scale variation in the modeling domain. The only additional ”burden” relative to standard Mercator seems to be the varying angle of the grid alignment to true East. However, with programs like proj/geod this is straightforward to calculate. What one gets in return is a closer correspondence between projected and unprojected grid cell sizes.
  • kate suggests in Hedstrom et al. (2017) that the Lambert conformal projection is often a good choice for mid-latitudes. The Lambert conformal projection can be tweaked to have the lines of true scale along arbitrary parallels, which means that the region of low scale variation can be chosen to be around any desired latitude (except perhaps the poles). I'm not sure about any possible (dis-) advantages of the Lambert conformal projection compared to Oblique Mercator. A difference is certainly that Oblique Mercator can have true scale along lines other than parallels, which may be an advantage if the coast (and hence current systems) are not aligned with parallels.
  • Is it meaningful or advantageous to interpret the modeling domain as a flat space instead of as a part of a spheroid? Given the distortions of area (and possibly direction, in the case of the oblique Mercator projection), it is clear that the conservation laws governing "real" geophysical flow are violated when modeling on the plane. For example, the conservation of momentum (mass times velocity) is violated due to distortions of area (mass) and direction. For small domains, are these errors potentially small enough to neglect them?
    • The only benefit I can think of is the lower computational cost after switching off the curvilinear terms in the model, since those aren't used anyways. A downside are possible complications regarding the processing of the forcing fields. I suppose these would need to be projected onto the plane as well, which may increase the computational cost again? On the other hand, given that e.g. linear momentum is not conserved by the projection anyways, it may not make sense to modify the forcing. What I do know is that renting a computer from a cloud provider is terribly expensive, so it’s an attractive prospect to save some money for those of us who don’t have a supercomputer facility available for free in the basement of their institution.
    • Unless anybody can answer these questions, it’s probably prudent to set up a realistic (spherical curvilinear) simulation and then start to simplify things and see if it makes a difference.
    • In another post, we use a plane model to do some simple storm surge modeling. This allows using relatively realistic (but slightly distorted) topography, while not having to formulate/implement the model in curvilinear coordinates. Pretty sloppy, but still useful.

If you have any comments or suggestions, please let me know.

References