GELE Geodesy — Geodetic and Cartesian CoordinatesStudy Notes

A coordinate system in geodesy defines the spatial location of points on or above the Earth's surface. Two primary systems exist: **Geodetic Coordinates (φ, λ, h):** - Latitude (φ): angular distance from the equator, ranging from –90° (South Pole) to +90° (North Pole). In the Philippines, latitudes range approximately 5°N to 21°N. - Longitude (λ): angular distance from the Prime Meridian (Greenwich), ranging from –180° to +180° (or 0° to 360°). The Philippines spans roughly 117°E to 127°E. - Ellipsoidal Height (h): vertical distance from the reference ellipsoid surface to the point, measured along the ellipsoid normal. This differs from orthometric height (which references the geoid). **Cartesian Coordinates (X, Y, Z):** - A 3D right-handed system with origin at Earth's center of mass. - X-axis: points toward the intersection of the Prime Meridian and the equator (Greenwich meridian, equatorial plane). - Y-axis: points 90° east in the equatorial plane (toward 90°E longitude). - Z-axis: points toward the North Pole along the rotational axis. - Units: meters (m), referenced to WGS84 (a = 6,378,137 m; semi-major axis) or other reference ellipsoids. **Why Two Systems?** GNSS receivers (GPS, GNSS) intrinsically compute Cartesian positions because the satellite constellation orbits in a geocentric frame. However, land surveyors, mapmakers, and property owners need geodetic coordinates and projected coordinates (UTM, PPCS) to work with maps and cadastral records under Philippine surveying standards (PD 1529, CA 141).

A reference ellipsoid is a mathematical model of Earth's shape, defined by two primary parameters: **Semi-major axis (a):** The equatorial radius of the ellipsoid. **Semi-minor axis (b):** The polar radius. From these, we derive: - **Flattening:** f = (a − b) / a - **Eccentricity (first):** e² = (a² − b²) / a² = 2f − f² ≈ 2f for small f - **Second eccentricity:** e'² = (a² − b²) / b² **WGS84 (World Geodetic System 1984):** - a = 6,378,137.0 m - f = 1/298.257223563 - e² = 0.00669438002290 - Defined by the International Earth Rotation and Reference Systems Service (IERS). - Used globally and is the de facto standard for GNSS. **PRS92 (Philippine Reference System):** - Based on WGS84 datum (since 1994). - Legal standard for surveying and mapping in the Philippines under RA 4374. - Requires conversion of historical data (e.g., from Clarke 1880 ellipsoid used in older surveys). - PPCS (Philippine Plane Coordinate System) uses PRS92 geodetic datum with UTM-based projection. **Prime Vertical Radius (N):** The radius of curvature in the prime vertical plane (the plane containing the meridian and the Z-axis) is critical for forward and inverse conversions: $$N = \frac{a}{\sqrt{1 - e^2 \sin^2 \phi}}$$ Note that N varies with latitude: at the equator, N ≈ a; at the poles, N ≈ b. This variation is why iterative methods are needed for the inverse transformation.

The forward transformation takes a geodetic position (latitude, longitude, ellipsoidal height) and converts it to geocentric Cartesian coordinates. This is the primary conversion needed when receiving geodetic data and requiring Cartesian output. **Mathematical Formulas:** Given: - Latitude: φ - Longitude: λ - Ellipsoidal height: h - Reference ellipsoid parameters: a (semi-major axis), e² (first eccentricity squared) Prime vertical radius: $$N = \frac{a}{\sqrt{1 - e^2 \sin^2 \phi}}$$ Cartesian coordinates: $$X = (N + h) \cos \phi \cos \lambda$$ $$Y = (N + h) \cos \phi \sin \lambda$$ $$Z = [N(1 - e^2) + h] \sin \phi$$ **Physical Interpretation:** - The term (N + h) is the distance from the Earth's rotational axis to the point in the equatorial plane projection. - The factor cos(φ) in X and Y reflects the decreasing distance from the rotational axis as latitude increases. - The (1 − e²) factor in the Z formula accounts for the ellipsoid's flattening; it converts the vertical component to account for the ellipsoid's non-spherical shape. - Without (1 − e²), you would effectively model a sphere, introducing errors up to ±21 m in Z at high latitudes. **Step-by-Step Procedure:** 1. Compute the prime vertical radius N using the given latitude. 2. Calculate X = (N + h) cos(φ) cos(λ). 3. Calculate Y = (N + h) cos(φ) sin(λ). 4. Calculate Z = [N(1 − e²) + h] sin(φ). 5. Round to appropriate precision (typically ±0.001 m for survey-grade GNSS).

The inverse transformation is more complex because latitude (φ) and height (h) are coupled through the prime vertical radius N. Closed-form solutions are available (e.g., Bowring's method, Heikkinen's method) but iteration is conceptually simpler and widely taught in surveying programs. **Longitude (Direct Calculation):** Longitude is decoupled and can be found directly: $$\lambda = \arctan\left(\frac{Y}{X}\right)$$ **Caution:** The two-argument arctangent function atan2(Y, X) must be used to resolve the quadrant correctly: - If X > 0, Y > 0: λ is in the first quadrant (0° to 90°). - If X < 0, Y > 0: λ is in the second quadrant (90° to 180°). - If X < 0, Y < 0: λ is in the third quadrant (180° to 270° or −180° to −90°). - If X > 0, Y < 0: λ is in the fourth quadrant (270° to 360° or 0° to −90°). Always express λ in the range [−180°, +180°] or [0°, 360°] as convention dictates. **Latitude and Height (Iterative Method):** 1. **Compute the cylindrical radius:** $$p = \sqrt{X^2 + Y^2}$$ 2. **Initialize latitude (first guess):** $$\phi_0 = \arctan\left(\frac{Z}{p(1 - e^2)}\right)$$ This initial guess accounts for the ellipsoid's flattening and works well for initial approximation. 3. **Iteration loop (repeat until convergence):** a. Compute N(φ): $$N_i = \frac{a}{\sqrt{1 - e^2 \sin^2 \phi_i}}$$ b. Compute height: $$h_i = \frac{p}{\cos \phi_i} - N_i$$ (Note: If |cos φ| is very small, use h = Z/sin φ − N(1 − e²) instead.) c. Update latitude: $$\phi_{i+1} = \arctan\left(\frac{Z}{p} \left(1 - \frac{e^2 N_i}{N_i + h_i}\right)^{-1}\right)$$ Or equivalently: $$\phi_{i+1} = \arctan\left(\frac{Z + e^2 N_i \sin \phi_i}{p}\right)$$ 4. **Convergence criterion:** Stop when |φ_{i+1} − φ_i| < ε (typically 10⁻¹² radians or 10⁻⁸ seconds of arc). **Closed-Form Alternatives:** Bowring's method and other closed-form solutions avoid iteration but are more algebraically complex. They are useful for embedded systems or high-throughput calculations. For licensure exams, understanding the iterative method is usually sufficient. **Physical Interpretation:** - p is the distance from the Z-axis (the perpendicular distance to Earth's rotational axis). - The initial φ₀ guess uses the spheroid approximation; successive iterations refine it. - Convergence is usually achieved in 2–3 iterations for terrestrial survey points.

**Pitfall 1: Forgetting the (1 − e²) Factor in Z** A frequent mistake, especially under exam pressure, is to write: $$Z_\text{wrong} = (N + h) \sin \phi$$ Instead of: $$Z_\text{correct} = [N(1 - e^2) + h] \sin \phi$$ Error magnitude: At φ = 45°, this omission causes an error of approximately: $$\Delta Z \approx N e^2 \sin \phi \approx 6,379,000 \times 0.00669438 \times 0.707 \approx 30,200 \ \text{m}$$ This is catastrophic and will fail any licensure exam problem. **Pitfall 2: Incorrect Longitude Quadrant** Using the single-argument arctan (or misinterpreting its output) leads to ambiguous results: $$\lambda_\text{wrong} = \arctan(Y/X)$$ Example: For X = −3M, Y = 5M: - arctan(5/(−3)) = arctan(−1.667) ≈ −59.04°, which is wrong. - Correct: atan2(5M, −3M) = 180° − 59.04° = 120.96°. ✓ **Solution:** Always use the two-argument arctangent (atan2 in most programming languages and modern calculators). **Pitfall 3: Confusing Ellipsoidal Height (h) with Orthometric Height (H)** - **Ellipsoidal height (h):** Geometric distance from the WGS84 ellipsoid surface to the point, measured along the ellipsoid normal. This is what GNSS provides directly. - **Orthometric height (H):** Height above mean sea level (the geoid), determined by leveling and gravity measurements. Used for engineering, construction, and legal purposes in the Philippines. - **Geoid undulation (N):** h = H + N (approximately). Under RA 4374 and PD 1529, cadastral maps in the Philippines use orthometric heights. Many engineering projects require conversion from GNSS-derived ellipsoidal heights to orthometric heights using a geoid model (e.g., PRS92 geoid model). **Pitfall 4: Using the Wrong Reference Ellipsoid or Datum** - **WGS84:** Global standard, used by GNSS satellites (GPS, GLONASS, Galileo, BeiDou). Parameters: a = 6,378,137 m, e² = 0.00669438. - **PRS92:** Philippines national datum, legally required for surveying under RA 4374. Defined as WGS84 as of 1994 (no datum transformation needed if both use WGS84). - **Clarke 1880:** Old British survey ellipsoid, used in some historic Philippine surveys. Datum transformation required (typically via 7-parameter Helmert transformation or Molodensky). Always verify the datum of your input data. Mixing datums without transformation introduces errors of 100s of meters. **Pitfall 5: Numerical Precision and Rounding** - Use at least **double-precision floating-point** (64-bit) for all calculations. Single precision (32-bit) introduces ±1–10 m errors. - For survey-grade GNSS (±0.01 m), compute N to at least 0.01 m precision. - In iterative methods, iterate until |Δφ| < 10⁻¹² radians (approximately 10⁻⁸ seconds of arc). - Round final coordinates only after all computations are complete; do not round intermediate values. **Pitfall 6: Singularities Near the Poles** - At the North Pole (φ = 90°), cos(φ) = 0, making h = p / cos(φ) → ∞ (numerically unstable). - **Solution:** Near the poles, use h = Z / sin(φ) − N(1 − e²) instead. - For latitudes |φ| > 84°, use specialized polar-coordinate algorithms or the Bowring closed-form method. **Pitfall 7: Sign Conventions for Longitude** - Positive λ for east (Eastern Hemisphere); negative λ for west (Western Hemisphere). - The Philippines uses positive longitudes (117°E to 127°E). - Some legacy systems use 0°–360° notation; convert to ±180° convention to avoid confusion. **Pitfall 8: Iteration Not Converging** Rare, but can occur if: - Initial guess φ₀ is very poor (e.g., initial guess φ₀ = 0° for a polar point). - **Solution:** Use a more robust initial guess, e.g., φ₀ = atan(Z/p) without the (1−e²) factor for initial rapid convergence. - Use Bowring's or other closed-form methods if iteration is unreliable.

**Application 1: GNSS Survey Integration with Cadastral Maps** Under RA 4374 (Geodetic Engineers Law) and CA 141 (Public Land Act), property boundaries and lot surveys must be referenced to an approved national datum. The legal procedure in the Philippines: 1. Conduct GNSS survey on site, receiving WGS84 Cartesian (X, Y, Z) from the receiver. 2. Convert to WGS84 geodetic (φ, λ, h) using the inverse transformation. 3. Project to PPCS (Philippine Plane Coordinate System) using UTM conversion formulas. PPCS coordinates are easting (E) and northing (N) in meters, aligned with cadastral maps. 4. Adjust for local datum if historical surveys used Clark 1880 (rare now, but still found in some old titles). 5. Report results in PPCS coordinates on the cadastral sketch (plan), as required by PD 1529 (Decree on Boundary Markers). **Application 2: Baseline Correction in Trilateration and GPS Networks** When establishing control networks (RA 4374 mandates precise control networks for large surveys), distances measured on the ground must be reduced to the ellipsoid and then to the UTM projection plane: 1. Measure baseline distance **on ground** (slope distance, corrected for refraction). 2. Convert slope distance to ellipsoidal chord using ellipsoidal height from GNSS. 3. Convert ellipsoidal chord to **grid distance** (UTM) using scale factors (typically 0.99954–1.00046 depending on zone location). 4. Use grid distances in adjustment (least-squares network adjustment). Cartesian to geodetic conversion is the bridge step for determining ellipsoidal heights. **Application 3: Road and Infrastructure Design** A road project spanning multiple municipalities requires: 1. GNSS control points (Cartesian from receiver) → geodetic (φ, λ, h) → PPCS grid (E, N). 2. Align design centerline to PPCS coordinates. 3. Compute road chainage and cross-sections in grid coordinates. 4. Report tie-ins to existing cadastral corners (also in PPCS) for legal and construction purposes. **Application 4: Monitoring Ground Deformation and Structural Movement** Repeat GNSS campaigns (e.g., annual measurements at bridge pylons, dams, or landslide sites): 1. Time 1: GNSS → (X₁, Y₁, Z₁) in Cartesian → (φ₁, λ₁, h₁) geodetic. 2. Time 2: GNSS → (X₂, Y₂, Z₂) in Cartesian → (φ₂, λ₂, h₂) geodetic. 3. Displacement vector: ΔX = X₂ − X₁, etc. 4. Convert displacement to local topocentric frame (East, North, Up) for easier interpretation. The inverse transformation enables precise height changes (Δh) to be detected—critical for dam safety, tunnel convergence, and seismic monitoring. **Application 5: Establishing Vertical Datums for Leveling Networks** Philippine leveling networks refer to mean sea level (MSL) or orthometric height (H). However, GNSS provides ellipsoidal height (h). The relationship: $$h = H + N_\text{geoid}$$ where N_geoid is the geoid undulation (geoidal height). To tie GNSS benchmarks to the leveling network: 1. Establish GNSS position at known leveling benchmark (e.g., tide gauge station). 2. Compute ellipsoidal height h from inverse transformation. 3. Compute geoid undulation: N_geoid = h − H. 4. Apply this correction to all GNSS-derived heights in the region. **Application 6: Datum Transition in Legacy Surveys** Many Philippine cadastral records (pre-1994) reference **Clarke 1880** or **Everest ellipsoid**. To update them to PRS92/WGS84: 1. Extract old geodetic coordinates (φ_old, λ_old) from historical records. 2. Convert to Cartesian using Clarke 1880 parameters. 3. Apply 7-parameter Helmert transformation (translation, rotation, scale) to convert to WGS84 Cartesian. 4. Convert WGS84 Cartesian back to WGS84 geodetic. 5. Update survey records under RA 4374 requirements. This process, required for legal compliance, relies entirely on coordinate conversions.

Quick reference for coordinate system parameters and conversion methods: | Aspect | Geodetic (φ, λ, h) | Cartesian (X, Y, Z) | |--------|-------------------|---------------------| | **Definition** | Angular latitude/longitude + height | Geocentric 3D Cartesian | | **Units** | Degrees (°) or radians; meters (m) | Meters (m) | | **Typical Range (Philippines)** | 5°N–21°N, 117°E–127°E, 0–3000 m | X ≈ ±4M, Y ≈ ±6M, Z ≈ 0–2M | | **Reference Surface** | WGS84 ellipsoid (or PRS92) | WGS84 ellipsoid (or PRS92) | | **Forward Conversion** | φ, λ, h → X, Y, Z | Formula: X = (N+h)cosφcosλ, etc. | | **Inverse Conversion** | X, Y, Z → φ, λ, h | Longitude direct; φ, h iterative | | **Computational Complexity** | Low (direct formulas) | Medium (iteration required) | | **Primary Use** | Maps, land records, engineering design | GNSS output, coordinate transformation | | **Legal Context (PH)** | PPCS (UTM-projected) for cadastre | Intermediate for conversion | **Convergence Behavior for Inverse Iteration:** | Latitude φ | Initial Guess Error | Iterations to Convergence (ε = 10⁻¹²) | Final Accuracy (h) | |-----------|---------------------|----------------------------------------|---------| | 0° (Equator) | ±10° | 2–3 | ±0.001 m | | 15° (Manila) | ±10° | 2–4 | ±0.001 m | | 45° (Mid-latitude) | ±10° | 2–4 | ±0.001 m | | 80° (Near pole) | ±5° | 3–5 | ±0.01 m | Note: Convergence is very fast; even a poor initial guess converges within a few iterations.