Coordinate Systems

Cocoa supports multiple cylindrical map projections for representing spherical Earth coordinates in the computational domain. This document describes the available projections and their mathematical formulations.

Overview

Coastal ocean models operate on a spherical Earth but require planar coordinates for finite element computation. Cocoa implements five cylindrical projections, each with different trade-offs between preserving 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. 18 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

• Simple 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 relationships across the domain, making it suitable for analysis where accurate area representation is critical.

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), historically used for navigation. It provides correct local angles but 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, reducing pole distortion while maintaining reasonable shape preservation.

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 = \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

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 = \frac{g \bar{H}}{4 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.

GWCE Mass Matrix

Nodal contributions to the mass matrix are weighted by \(\cos\phi\) to account for the spherical area element:

\[M_{ij}^e = \frac{A_e}{12} \cos\phi_i \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\]

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, this should be near the center of the computational domain.

Global Mesh Support

For global ocean modeling, Cocoa provides two features to handle challenges that arise when meshes span the entire Earth:

1. Prime meridian crossing detection (automatic)

2. Coordinate rotation (opt-in via configuration)

Prime Meridian Crossing

Elements that straddle the \(\pm 180\)° longitude line (antimeridian) can cause incorrect Jacobian calculations due to the discontinuity in longitude values. Cocoa automatically detects and corrects these elements during mesh initialization.

Detection Algorithm:

1. Normalize all element vertex longitudes to [0°, 360°)

2. Compute the element Jacobian and total edge length

3. If Jacobian < 0 OR edge length > 360°, the element likely crosses the meridian

4. Apply correction:

   • If exactly one vertex has longitude > 180°: subtract 360° from that vertex

   • If exactly one vertex has longitude < 180°: add 360° to that vertex

5. Verify the correction improved the geometry (positive Jacobian)

6. If not improved, revert to original coordinates

This handling is automatic and requires no user configuration. It is implemented in GlobalMeshHandler.cpp and invoked during element projection from BasisFunctions.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 supports rotating the computational coordinate system so that the computational poles are located 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. Physics computation: Coriolis parameter and spherical metric corrections use the rotated latitude \(\phi'\) (not geographic latitude)

3. Tidal potential: Uses geographic latitude \(\phi\) (physical requirement)

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

5. Output: Velocities can be written in geographic frame for visualization

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.

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

Implementation Files

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

• geometry/ProjectionTransformer.cpp: Projection implementations

• geometry/ProjectionScaleFactor.hpp: Scale factor data structure

• geometry/BasisFunctions.cpp: Scale factor computations

• geometry/GlobalMeshHandler.hpp: Prime meridian crossing detection

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

• geometry/VelocityTransform.hpp: Velocity transformation matrices

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