Coordinate Systems

Overview

The finite element computation needs planar coordinates, so Cocoa maps the sphere onto one of five cylindrical projections. Each preserves a different combination of angles, areas and distances.

All projections use a reference point \((\lambda_0, \phi_0)\) as the origin and Earth’s radius \(R = 6378206.4\) m (Clarke 1866 equatorial radius).

Available Projections

Projection	Config Name	Angles	Areas	Pole Distortion
Equidistant Cylindrical	EquidistantCylindrical	No	No	Moderate
Equal Area	EqualArea	No	Yes	None
Mercator	Mercator	Yes	No	Extreme
Miller	Miller	Partial	No	Reduced
Gall Stereographic	GallStereographic	Partial	Partial	Moderate

Fig. 37 Comparison of three cylindrical projections supported by Cocoa. Green ellipses (Tissot indicatrices) show local distortion at different latitudes.

Equidistant Cylindrical (Default)

The equidistant cylindrical projection, also known as the Plate Carrée or Carte Parallelo-grammatique (CPP) projection, is Cocoa’s default. It corresponds to ADCIRC’s ICS=21 setting.

Projection equations:

\[x = R \cos\phi_0 (\lambda - \lambda_0), \quad y = R \phi\]

Properties:

• Preserves distances along meridians (constant latitude spacing)

• Distorts east-west distances away from the reference latitude

• Linear relationship between grid and geographic coordinates

• Suitable for regional models where the reference latitude \(\phi_0\) is near the domain center

Scale factors:

The governing equations require derivatives in physical space. Because this projection uses \(\cos\phi_0\) rather than the local \(\cos\phi\), x-derivatives must be corrected:

\[S_x = \frac{\cos\phi_0}{\cos\phi}\]

For nodes near the reference latitude, \(S_x \approx 1\). At higher latitudes where \(\cos\phi < \cos\phi_0\), the scale factor exceeds unity.

Equal Area

The cylindrical equal-area projection preserves area across the domain.

Projection equations:

\[x = R \cos\phi_0 (\lambda - \lambda_0), \quad y = R \frac{\sin\phi}{\cos\phi_0}\]

Properties:

• Preserves area (no distortion in integrated quantities)

• Distorts shapes and angles, especially at high latitudes

• No singularity at poles

• Useful for mass conservation verification

Scale factors:

The projection exponent is 0, yielding scale factors:

\[S_x^{cont} = \cos\phi, \quad S_y^{cont} = \frac{1}{\cos\phi_0}\]

Mercator

The Mercator projection is conformal (angle-preserving) and severely distorts areas at high latitudes.

Projection equations:

\[\begin{split}\begin{aligned} x &= R \cos\phi_0 (\lambda - \lambda_0) \\[4pt] y &= R \cos\phi_0 \ln\left[\tan\left(\frac{\pi}{4} + \frac{\phi}{2}\right)\right] \end{aligned}\end{split}\]

Properties:

• Conformal (preserves local angles and shapes)

• Severe area distortion approaching poles

• Singularity at \(\phi = \pm 90^\circ\)

• Best for equatorial and mid-latitude regions

Scale factors:

The projection exponent is 2, yielding:

\[S_x = S_y = \frac{\cos\phi_0}{\cos\phi}\]

The equal x and y scale factors reflect the conformal property.

Miller Cylindrical

The Miller projection is a compromise between Mercator and equal-area that reduces distortion near the poles.

Projection equations:

\[\begin{split}\begin{aligned} x &= R \cos\phi_0 (\lambda - \lambda_0) \\[4pt] y &= 1.25\, R \cos\phi_0 \ln\left[\tan\left(\frac{\pi}{4} + 0.4\phi\right)\right] \end{aligned}\end{split}\]

Properties:

• Reduces extreme distortion near poles compared to Mercator

• Neither conformal nor equal-area (a compromise)

• Useful for whole-Earth visualizations

• Less common for hydrodynamic modeling

Scale factors:

Uses modified latitude scaling with factor 0.4 instead of 1.0, reducing the rate of scale change with latitude.

Gall Stereographic

The Gall stereographic projection uses stereographic latitude scaling with standard parallels at approximately \(\pm 45^\circ\). It does not depend on a reference origin.

Projection equations:

\[x = R\,\frac{\sqrt{2}}{2}\,\lambda, \quad y = R \left(1 + \frac{\sqrt{2}}{2}\right) \tan\left(\frac{\phi}{2}\right)\]

Properties:

• Standard parallels at \(\pm 45^\circ\) where distortion is minimized

• Compromise between area and shape preservation

• Independent of reference latitude (uses fixed standard parallels)

• Moderate distortion throughout

Note

This departs from ADCIRC on purpose. ADCIRC’s CYLINDERMAP (ICS = 24) writes the easting as \(\frac{\sqrt{2}}{2}\lambda\), without \(R\), while its own Gall scale factors (for example \(\text{SFMX} = \frac{\sqrt{2}/2}{\cos\phi}\)) are the stretch of an easting that carries \(R\), so the two disagree by a factor of \(R\) in \(x\). Cocoa keeps ADCIRC’s scale factors and gives the easting its \(R\). Before 0.6.0 Cocoa copied the radius-less easting, and a Gall run diverged within its first steps.

Scale Factor System

Cocoa precomputes seven scale factors per element to correct the shallow water equations for projection distortions:

Scale Factor	Purpose
scale_factor	General coordinate distortion correction
scale_factor_continuity_x	Corrects \(\partial u/\partial x\) in continuity equation
scale_factor_continuity_y	Corrects \(\partial v/\partial y\) in continuity equation
scale_factor_momentum_x	Corrects x-velocity derivatives in momentum equations
scale_factor_momentum_y	Corrects y-velocity derivatives in momentum equations
scale_factor_y_coordinate	Y-coordinate correction (\(\cos\phi\) for spherical)
scale_factor_tan_phi	Spherical geometry correction: \(\tan\phi / R\)

General Scale Factor Formulas

For cylindrical projections with exponent \(p\) (where \(p=0\) for equal-area, \(p=1\) for equidistant, \(p=2\) for Mercator):

\[\begin{split}S_x^{cont} &= \cos\phi_0 \cdot \cos^{p-1}\phi \\ S_y^{cont} &= \cos^{p-1}\phi_0 \\ S_x^{mom} &= \frac{\cos\phi_0}{\cos\phi} \\ S_y^{mom} &= \left(\frac{\cos\phi_0}{\cos\phi}\right)^{p-1} \\ S_{y,coord} &= \cos\phi \\ S_{\tan\phi} &= \frac{\tan\phi}{R}\end{split}\]

The \(S_{\tan\phi}\) factor appears in the spherical metric correction for Coriolis and lateral stress terms.

Application in Equations

GWCE Stiffness Matrix

The stiffness matrix for the \(\nabla \cdot (gH\nabla\zeta)\) term uses scale factors to correct gradient products:

\[K_{ij}^e = g \bar{H} A_e \left( C_x \frac{\partial N_i}{\partial x} \frac{\partial N_j}{\partial x} + C_y \frac{\partial N_i}{\partial y} \frac{\partial N_j}{\partial y} \right)\]

where \(C_x = 1/S_x^2\) and \(C_y = 1\) for the equidistant cylindrical projection (Equation (13) gives the same term in the code’s \(FDX_i = 2A_e\,\partial N_i/\partial x\) form).

GWCE Mass Matrix

Nodal contributions to the mass matrix are weighted by \(\cos\phi\) to account for the spherical area element. The weight is that of the column node \(j\):

\[M_{ij}^e = \frac{A_e}{12} \cos\phi_j \cdot (\text{mass coefficients})\]

Momentum Equations

Velocity derivatives in the momentum equations are corrected by the momentum scale factors. For the pressure gradient:

\[\frac{\partial \zeta}{\partial x} \rightarrow S_x^{mom} \frac{\partial \zeta}{\partial x}\]

The spherical metric correction in the effective Coriolis parameter uses \(S_{\tan\phi}\):

\[f_{eff} = 2\Omega\sin\phi + U \cdot S_{\tan\phi} \cdot R = 2\Omega\sin\phi + U\tan\phi\]

On a rotated-pole mesh the two terms take different latitudes: the planetary term keeps the geographic \(\phi\) and the metric term takes the rotated \(\phi'\). See Coordinate Rotation.

Configuration

Specify the projection in the YAML configuration file:

mesh:
  filename: domain.nc
  projection:
    type: EquidistantCylindrical       # Default
    center: [-90.0, 29.0]              # [longitude, latitude]

Available projection type values:

• EquidistantCylindrical (default)

• EqualArea

• Mercator

• Miller

• GallStereographic

The center specifies \((\lambda_0, \phi_0)\) in decimal degrees. For most coastal applications, put it near the center of the computational domain.

Global Mesh Support

For meshes that span the entire Earth, Cocoa handles the longitude branch cut automatically and offers coordinate rotation as a configuration option.

Longitude Branch Cut

A global mesh stores longitudes on a 360° branch, so an element that straddles the branch cut has vertices at opposite ends of the range although they are neighbors on the sphere. Projected as stored, such an element would be a full revolution wide. Cocoa puts each element’s vertices on one branch before it computes the element’s geometry.

Unwrapping Algorithm (triangles and quadrilaterals alike):

1. Normalize the element’s vertex longitudes to [0°, 360°)

2. If they span at most 180°, the element does not cross the cut and is left alone

3. Otherwise shift by 360° every vertex on the side of 180° that holds no more vertices than the other (a straddling triangle’s lone minority vertex; on a quadrilateral split 2-2, the two vertices above 180°)

4. Keep the shift only if it narrowed the span

Only an actual span over half a revolution triggers the unwrap. The element’s Jacobian sign plays no part: a clockwise element is a connectivity property, fixed by reordering its nodes when the mesh is read, and treating a negative Jacobian as a crossing would widen a compact element into a 360° one.

Node coordinates are projected individually, so a difference of projected \(x\) between two nodes on opposite sides of the cut is reduced to its shortest equivalent modulo the map’s easting period (Projection::shortest_easting_delta) wherever such a difference is formed: the wet/dry head gradient (Wetting and Drying) and land-boundary chords and edge lengths.

The unwrapping is geometry/MeshProjector.hpp, called when elements are oriented (HostMesh.cpp) and when element geometry is evaluated (ElementGeometry.cpp).

Coordinate Rotation

For global models that include polar regions, the geographic poles create singularities in the shallow water equations (infinite Coriolis at poles, converging meridians). Cocoa can rotate the computational coordinate system to move the computational poles away from the geographic poles.

Mathematical Basis:

A 3x3 rotation matrix \(\mathbf{R}\) transforms geographic coordinates to a rotated system:

\[\begin{split}\begin{pmatrix} x' \\ y' \\ z' \end{pmatrix} = \mathbf{R} \begin{pmatrix} x \\ y \\ z \end{pmatrix}\end{split}\]

where \((x, y, z)\) are Cartesian coordinates on the unit sphere.

The rotation is specified by providing the location of the new North Pole in geographic coordinates \((\lambda_{pole}, \phi_{pole})\). The rotation matrix is constructed using ADCIRC’s GET_ROTMAT_ZNVEC algorithm, which creates a rotation that moves the geographic North Pole to the specified location.

Workflow:

1. Coordinate transformation: all node coordinates are rotated from geographic \((\lambda, \phi)\) to rotated \((\lambda', \phi')\) coordinates

2. Metric terms: spherical metric corrections and projection scale factors use the rotated latitude \(\phi'\) (not geographic latitude), since they are properties of the frame the equations are discretized in

3. Coriolis parameter: uses geographic latitude \(\phi\). The planetary rotation term is frame-independent, since the orthonormal rotation carries the pole direction and the local vertical together and leaves \(\hat{z}\cdot\hat{r} = \sin\phi\) unchanged

4. Tidal potential: uses geographic latitude \(\phi\) (physical requirement)

5. Velocity transformation: per-node 2x2 matrices transform velocities between geographic (east, north) and rotated coordinate frames

6. Output: velocities can be written in the geographic frame

Configuration:

mesh:
  rotation:
    enabled: true
    pole_location: [0.0, 0.0]  # Move North Pole to equator at prime meridian

The pole_location specifies where the geographic North Pole will appear in the rotated coordinate system, using the same format as ADCIRC’s znorth_in_spherical_coors.

On a rotated mesh the projection is chosen as ADCIRC’s CYLINDERMAP chooses it, by a scale-factor exponent: EqualArea, EquidistantCylindrical and Mercator keep their own projection, while Miller and GallStereographic, which have no exponent of their own, are projected and scaled as Mercator.

See Coordinate Rotation (Global Meshes) for full configuration options and Global Simulations for a walkthrough.

Implementation Files

• geometry/ProjectionTransformer.hpp: Projection type enum and conversions

• geometry/ProjectionTransformer.cpp: Projection dispatch

• geometry/ProjectionFormulas.hpp: Projection formulas and shortest_easting_delta

• geometry/ProjectionScaleFactor.hpp: Scale factor data structure

• geometry/ProjectionScaleFactors.cpp: Scale factor computations

• geometry/MeshProjector.hpp: Longitude branch-cut unwrapping

• geometry/CoordinateRotation.hpp: Rotation matrix and spherical rotation

• geometry/VelocityTransform.hpp: Velocity transformation matrices

• geometry/RotationData.hpp: Container for rotation-related data