================== 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 :math:`(\lambda_0, \phi_0)` as the origin and Earth's radius :math:`R = 6378206.4` m (Clarke 1866 equatorial radius). Available Projections --------------------- .. list-table:: :header-rows: 1 :widths: 25 20 20 20 15 * - 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 .. figure:: ../_static/images/diagrams/projection_comparison.svg :alt: Side-by-side comparison of equidistant cylindrical, equal area, and Mercator projections :width: 650px :align: center Comparison of three cylindrical projections supported by Cocoa. Green ellipses (Tissot indicatrices) show local distortion at different latitudes. .. _equidistant-cylindrical: Equidistant Cylindrical (Default) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The equidistant cylindrical projection, also known as the Plate Carr\ |eacute|\ e or Carte Parallelo-grammatique (CPP) projection, is Cocoa's default. It corresponds to ADCIRC's ICS=21 setting. **Projection equations:** .. math:: 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 :math:`\phi_0` is near the domain center **Scale factors:** The governing equations require derivatives in physical space. Because this projection uses :math:`\cos\phi_0` rather than the local :math:`\cos\phi`, x-derivatives must be corrected: .. math:: S_x = \frac{\cos\phi_0}{\cos\phi} For nodes near the reference latitude, :math:`S_x \approx 1`. At higher latitudes where :math:`\cos\phi < \cos\phi_0`, the scale factor exceeds unity. .. _equal-area: Equal Area ^^^^^^^^^^ The cylindrical equal-area projection preserves area across the domain. **Projection equations:** .. math:: 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: .. math:: S_x^{cont} = \cos\phi, \quad S_y^{cont} = \frac{1}{\cos\phi_0} .. _mercator: Mercator ^^^^^^^^ The Mercator projection is conformal (angle-preserving) and severely distorts areas at high latitudes. **Projection equations:** .. math:: \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} **Properties:** - Conformal (preserves local angles and shapes) - Severe area distortion approaching poles - Singularity at :math:`\phi = \pm 90^\circ` - Best for equatorial and mid-latitude regions **Scale factors:** The projection exponent is 2, yielding: .. math:: S_x = S_y = \frac{\cos\phi_0}{\cos\phi} The equal x and y scale factors reflect the conformal property. .. _miller: Miller Cylindrical ^^^^^^^^^^^^^^^^^^ The Miller projection is a compromise between Mercator and equal-area that reduces distortion near the poles. **Projection equations:** .. math:: \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} **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: Gall Stereographic ^^^^^^^^^^^^^^^^^^ The Gall stereographic projection uses stereographic latitude scaling with standard parallels at approximately :math:`\pm 45^\circ`. It does not depend on a reference origin. **Projection equations:** .. math:: 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 :math:`\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 :math:`\frac{\sqrt{2}}{2}\lambda`, without :math:`R`, while its own Gall scale factors (for example :math:`\text{SFMX} = \frac{\sqrt{2}/2}{\cos\phi}`) are the stretch of an easting that carries :math:`R`, so the two disagree by a factor of :math:`R` in :math:`x`. Cocoa keeps ADCIRC's scale factors and gives the easting its :math:`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: .. list-table:: :header-rows: 1 :widths: 30 70 * - Scale Factor - Purpose * - ``scale_factor`` - General coordinate distortion correction * - ``scale_factor_continuity_x`` - Corrects :math:`\partial u/\partial x` in continuity equation * - ``scale_factor_continuity_y`` - Corrects :math:`\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 (:math:`\cos\phi` for spherical) * - ``scale_factor_tan_phi`` - Spherical geometry correction: :math:`\tan\phi / R` General Scale Factor Formulas ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ For cylindrical projections with exponent :math:`p` (where :math:`p=0` for equal-area, :math:`p=1` for equidistant, :math:`p=2` for Mercator): .. math:: 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} The :math:`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 :math:`\nabla \cdot (gH\nabla\zeta)` term uses scale factors to correct gradient products: .. math:: 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 :math:`C_x = 1/S_x^2` and :math:`C_y = 1` for the equidistant cylindrical projection (Equation :eq:`gwce-stiffness-matrix` gives the same term in the code's :math:`FDX_i = 2A_e\,\partial N_i/\partial x` form). GWCE Mass Matrix ^^^^^^^^^^^^^^^^ Nodal contributions to the mass matrix are weighted by :math:`\cos\phi` to account for the spherical area element. The weight is that of the column node :math:`j`: .. math:: 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: .. math:: \frac{\partial \zeta}{\partial x} \rightarrow S_x^{mom} \frac{\partial \zeta}{\partial x} The spherical metric correction in the effective Coriolis parameter uses :math:`S_{\tan\phi}`: .. math:: 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 :math:`\phi` and the metric term takes the rotated :math:`\phi'`. See `Coordinate Rotation`_. Configuration ------------- Specify the projection in the YAML configuration file: .. code-block:: yaml 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 :math:`(\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. .. _prime-meridian-handling: Longitude Branch Cut ^^^^^^^^^^^^^^^^^^^^ A global mesh stores longitudes on a 360\ |deg| 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\ |deg|, 360\ |deg|) 2. If they span at most 180\ |deg|, the element does not cross the cut and is left alone 3. Otherwise shift by 360\ |deg| every vertex on the side of 180\ |deg| 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\ |deg|) 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\ |deg| one. Node coordinates are projected individually, so a difference of projected :math:`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 (:doc:`wetdry`) 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: 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 :math:`\mathbf{R}` transforms geographic coordinates to a rotated system: .. math:: \begin{pmatrix} x' \\ y' \\ z' \end{pmatrix} = \mathbf{R} \begin{pmatrix} x \\ y \\ z \end{pmatrix} where :math:`(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 :math:`(\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 :math:`(\lambda, \phi)` to rotated :math:`(\lambda', \phi')` coordinates 2. Metric terms: spherical metric corrections and projection scale factors use the rotated latitude :math:`\phi'` (not geographic latitude), since they are properties of the frame the equations are discretized in 3. Coriolis parameter: uses geographic latitude :math:`\phi`. The planetary rotation term is frame-independent, since the orthonormal rotation carries the pole direction and the local vertical together and leaves :math:`\hat{z}\cdot\hat{r} = \sin\phi` unchanged 4. Tidal potential: uses geographic latitude :math:`\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:** .. code-block:: yaml 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 :ref:`rotation-configuration` for full configuration options and :doc:`../user_guide/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