Converts geographic coordinates from WGS84 to NZGD2000.

transform_wgs84_to_nzgd2000.py

"""Converts geographic coordinates from WGS84 to NZGD2000."""
import math

AXIS_SEMI_MAJOR = 6378137.0  # Semi-major axis for GRS 80
AXIS_SEMI_MINOR = 6356752.314245179  # Semi-minor axis for GRS 80
FLATTENING = (AXIS_SEMI_MAJOR -
              AXIS_SEMI_MINOR) / AXIS_SEMI_MAJOR  # Flattening
# FLATTENING = 1 / 298.257222101  # or use the standard value for GRS 80
E2 = 2 * FLATTENING - FLATTENING**2  # Eccentricity squared
FALSE_EASTING = 1600000
FALSE_NORTHING = 10000000
SCALE_FACTOR = 0.9996
CENTRAL_MERIDIAN = math.radians(173)  # Central meridian in radians


def transform_wgs84_to_nzgd2000(latitude: float, longitude: float) -> tuple:
    """Convert geographic coordinates from WGS84 to NZGD2000.

    This function calculates the Easting and Northing coordinates for the
    NZTM2000 projection using the Transverse Mercator projection method. The
    calculation involves several steps, including calculating the meridional
    arc, which is the distance on the ellipsoid from the equator to the
    latitude phi.

    Args:
        latitude: The latitude in degrees.
        longitude: The longitude in degrees.

    Returns:
        A tuple containing the Easting and Northing coordinates in meters.

    Key Terms:

    - latitude in radians (phi): geographic latitude in radians.

    - longitude in radians (lambda): geographic longitude in radians.

    - longitude of the central meridian (lambda0): for NZTM2000, the central
      meridian is 173 degrees east.

    - semi-major axis (alpha): the semi-major axis of the ellipsoid, the
      longest radius of the ellipsoid.

    - flattening (f): a measure of the compression of a sphere (in this
      context, a planet or moon) along its axis, resulting in an ellipsoid with
      a shorter axis from pole to pole than at the equator. It quantifies how
      much an ellipsoid deviates from being a perfect sphere.
      - formula: (alpha - beta) / alpha
        - alpha: semi-major axis
        - beta: semi-minor axis
      - purpose
        - Shape Approximation: Flattening is critical in accurately describing
          the Earth's shape. The Earth is not a perfect sphere but an oblate
          spheroid, meaning it is flattened at the poles due to its rotation.
          This flattening must be accurately represented in models used for
          navigation, surveying, and mapping.
        - Distortion Minimization in Mapping: In map projections, understanding
          the flattening of the Earth allows cartographers to minimize
          distortions when converting the three-dimensional surface of the
          Earth to two-dimensional maps. This is particularly important for
          ensuring accuracy in areas of critical navigation, land ownership
          mapping, and other applications where spatial precision is mandatory.

    - eccentricity squared (e^2): the square of the eccentricity of the
      ellipsoid.
      - formula: 2 * f - f^2 = (alpha^2 - beta^2) / alpha^2
      - purpose
        - Ellipsoidal Shape Modeling: Eccentricity squared is critical for
          accurately modeling the ellipsoidal shape of the Earth. It directly
          influences calculations related to the Earth's curvature in different
          geographic coordinate systems and geodetic calculations.
        - Map Projections: In map projections, `e^2` is crucial for adjusting
          the scaling and distortion that occur when the Earth's surface is
          represented on a flat plane, ensuring that these representations are
          as accurate as possible.

    - meridional arc (M):

      alpha * (A0 * phi - A2 * sin(2 phi) + A4 * sin(4 phi) - A6 * sin(6 phi))

      where coefficients A0, A2, A4, and A6 are derived from the eccentricity
      of the ellipsoid:

      - A0 = 1 - (1 / 4 * e^2) - (3 / 64 * e^4) - (5 / 256 * e^6)
      - A2 = (3 / 8) * (e^2 + (1 / 4 * e^4) + (15 / 128 * e^6))
      - A4 = (15 / 256) * (e^4 + (3 / 4 * e^6))
      - A6 = (35 / 3072) * e^6

    - radius of curvature in the prime vertical (N): the prime vertical is the
      east-west oriented plane that intersects the ellipsoid at a given
      latitude and is perpendicular to the meridian plane.

      - Formula: alpha / sqrt(1 - e^2 * sin^2(phi))
      - Purpose
        - measures how much the ellipsoid flattens at the poles and bulges at
          the equator by providing the curvature radius at any given latitude.
        - serves as a scaling factor for the latitude to adjust the distance
          calculations to the actual surface curve of the ellipsoid.
        - In practical terms, N adjusts the easting and northing calculations
          to be more accurate according to the true distances along the surface
          of the Earth.

    - tangent of latitude (T):

      - Formula: tan(phi)
      - Purpose
        - directly impacts the calculation of the northing coordinate in the
          projection, particularly influencing the quadratic and higher-order
          terms.
        - The tangent function increases dramatically as latitude approaches
          +/-90 degrees (the poles), which is critical because the projection
          needs to account for greater distortions as one moves away from the
          equator.
        - also represents the slope of the tangent line to the meridian at
          latitude phi, relating to how steeply the Earth curves away at that
          latitude.

    - second eccentricity squared (C):

      - Formula: (e^2 / (1 - e^2)) * cos^2(phi)
      - Purpose
        - accounts for the ellipsoidal shape's effect on meridian convergence,
          modifying how distances calculated along meridians and parallels vary
          with latitude.
        - adjusts the distance calculations to factor in the flattening effect
          of the ellipsoid as one moves north or south from the equator. The
          Earth is not a perfect sphere, so the distances along the parallels
          (east-west direction) and meridians (north-south direction) do not
          uniformly scale; `C` helps correct for these variances.
        - primarily affects the higher-order corrections in the easting and
          northing formulas, ensuring that these adjustments become more
          significant as one moves further from the equator and as the effects
          of the Earth's flattening become more pronounced.

    - adjusted longitude difference (A): adjusted longitudinal difference
      between a given point and the central meridian, factoring in the
      latitude. It essentially adjusts for the convergence of meridians toward
      the poles.

      - formula: (lambda - lambda0) * cos(phi)
      - purpose
        - Eastward Displacement: `A` quantifies the eastward displacement from
          the central meridian but scales this displacement by the cosine of
          the latitude. This scaling is crucial because as you move away from
          the equator toward the poles, the physical distance represented by a
          degree of longitude decreases. `A` compensates for this by reducing
          the eastward displacement as the cosine of the latitude decreases
          (i.e., as one moves closer to the poles).
        - Projection Scaling: By multiplying the longitudinal difference by the
          cosine of the latitude, `A` helps in maintaining accurate scale
          across the map. In a cylindrical projection like the Transverse
          Mercator, if this scaling were not applied, areas farther from the
          equator would appear disproportionately stretched in the east-west
          direction relative to their actual ground size.
        - Reduces Distortion: This term is critical in reducing distortion in
          the projected coordinates. By adjusting for the latitude's influence
          on the distance covered by a degree of longitude, `A` helps ensure
          that the projection retains more accurate geographical proportions,
          which is particularly important in high-accuracy mapping applications
          like cadastral and navigation systems.
        - Simplifies Higher-Order Corrections: `A` is used in the calculation
          of higher-order terms for both easting and northing, which further
          correct for the curvature and the elliptical shape of the earth.
          These corrections are essential for minimizing errors in the
          projection, especially over larger areas or regions with significant
          relief.

    - Easting (x):

      - Formula: N * (Linear Term + Cubic Correction + Quintic Correction)
        - Linear Term: Basic Eastward Displacement
          - Formula: A
          - Purpose: represents the basic eastward displacement from the
            central meridian, adjusted by the cosine of latitude to account for
            convergence of meridians towards the poles.
        - Cubic Correction: Primary Ellipsoidal Correction
          - Formula: (1 - T^2 + C) * A^3 / 6
          - Purpose: Corrects for the ellipsoidal shape of the earth, taking
            into account the squared tangent of latitude `T^2` and the square
            of the first eccentricity `C`, which adjusts for the flattening
            effects near the poles.
        - Quintic Correction: Higher-Order Ellipsoidal Correction
          - Formula: (5 - 18 * T^2 + T^4 + 72 * C - 58 * e^2) * A^5 / 120
          - Purpose: Further refinement based on higher-order terms of latitude
            and eccentricity, ensuring accuracy over larger areas and extreme
            latitudinal positions.

    - Northing (y): Incorporates the meridional arc and additional terms that
      correct for north-south displacement and the curvature effects.

      - Formula:
        Meridional Arc + N * T * (Quadratic + Quartic + Sextic Correction)
        - Meridional Arc (M): Meridional Distance
          - Purpose: Represents the northward distance from the equator to the
            latitude `phi` accounting for the curvature of the ellipsoidal
            earth.
        - Quadratic Correction: Basic Curvature Correction
          - Formula: A^2 / 2
          - Purpose: Adjusts for basic curvature effects due to the projection
            of the spherical coordinates onto a plane, particularly noticeable
            at higher latitudes.
        - Quartic Correction: Moderate Ellipsoidal Correction
          - Formula: (5 - T^2 + 9 * C + 4 * C^2) * A^4 / 24
          - Purpose: Corrects for moderate distortions due to the ellipsoid
            shape, incorporating terms for tangent and eccentricity squared.
        - Sextic Correction: Advanced Ellipsoidal Correction
          - Formula: (61 - 58 * T^2 + T^4 + 600 * C - 330 * e^2) * A^6 / 720
          - Purpose: A higher-order correction for more significant
            distortions, integrating complex interactions of latitude,
            flattening, and eccentricity for precise mapping across extensive
            geographical areas.
    """
    # Convert latitude to radians
    phi = math.radians(latitude)

    # Calculate coefficients for the meridional arc
    a0 = 1 - (E2 / 4) - (3 * E2**2 / 64) - (5 * E2**3 / 256)
    a2 = (3 / 8) * (E2 + (E2**2 / 4) + (15 * E2**3 / 128))
    a4 = (15 / 256) * (E2**2 + (3 * E2**3 / 4))
    a6 = (35 / 3072) * E2**3

    # Calculate the meridional arc
    m = AXIS_SEMI_MAJOR * (a0 * phi - a2 * math.sin(2 * phi) +
                           a4 * math.sin(4 * phi) - a6 * math.sin(6 * phi))

    # Recalculate phi and lambda for longitude conversion
    lambda_ = math.radians(longitude)

    # Calculate the radius of curvature in the prime vertical
    n = AXIS_SEMI_MAJOR / math.sqrt(1 - E2 * math.sin(phi)**2)
    t = math.tan(phi)
    c = (E2 / (1 - E2)) * math.cos(phi)**2
    a = (lambda_ - CENTRAL_MERIDIAN) * math.cos(phi)

    # Calculate easting and northing
    x = n * (a + (1 - t**2 + c) * a**3 / 6 +
             (5 - 18 * t**2 + t**4 + 72 * c - 58 * E2) * a**5 / 120)
    y = m + n * math.tan(phi) * (
        a**2 / 2 + (5 - t**2 + 9 * c + 4 * c**2) * a**4 / 24 +
        (61 - 58 * t**2 + t**4 + 600 * c - 330 * E2) * a**6 / 720)

    # Apply the scale factor and the false easting/northing
    easting = x * SCALE_FACTOR + FALSE_EASTING
    northing = y * SCALE_FACTOR + FALSE_NORTHING

    return easting, northing


if __name__ == "__main__":
    # Example usage:
    easting, northing = transform_wgs84_to_nzgd2000(-41.2865, 174.7762)
    print(f"Easting: {easting:.2f} m, Northing: {northing:.2f} m")
    easting, northing = transform_wgs84_to_nzgd2000(-43.5321, 172.6362)
    print(f"Easting: {easting:.2f} m, Northing: {northing:.2f} m")
    easting, northing = transform_wgs84_to_nzgd2000(-36.8485, 176.2920)
    print(f"Easting: {easting:.2f} m, Northing: {northing:.2f} m")