Spatial Coordinate Reference Systems and Transformations for 3D GIS Applications

1. Introduction

The Zhengtong Crystal Twin Platform achieves multi‑scale digital simulation from global to city‑level, relying on spatial coordinate reference concepts covered in geodesy and GIS textbooks, as well as computer graphics and solid geometry.

2. Spatial Coordinate Reference

2.1. Geodetic Coordinate System

Understanding the Earth ellipsoid is essential; the geodetic latitude‑longitude‑height (BLH) system represents points on the ellipsoid surface.

Latitude is the angle between the normal at point P and the equatorial plane; longitude is the angle between the meridian plane of P and the prime meridian, which explains the ranges –90° to +90° for latitude and –180° to +180° for longitude.

2.2. Earth‑Centered Earth‑Fixed (ECEF) Coordinate System

ECEF is a right‑handed Cartesian system with its origin at the ellipsoid’s centre; X‑axis points to the intersection of the prime meridian and equator, Y‑axis points eastward along the equator, and Z‑axis aligns with the rotation axis.

Both BLH and ECEF describe the same point, but in different forms; conversion between them is possible.

2.3. Web Mercator Projection

Projecting curved surface coordinates onto a plane yields a projected coordinate system. Web Mercator is a simplified transverse Mercator projection widely used in Web GIS because of its simple conversion formulas and square tile layout.

Typical projected X range is [-20037508.3427892, 20037508.3427892] m; Y range is limited to latitudes between –85.05112878° and +85.05112878°.

2.4. Local (Station) Coordinate System (ENU)

Global coordinates have large numeric values, which reduce precision for local calculations. An ENU (East‑North‑Up) system centered at a chosen station point provides small, numerically stable coordinates for local spatial computation.

3. Spatial Coordinate Transformations

3.1. BLH ↔ ECEF Conversion

3.1.1. BLH → ECEF (XYZ)

Derivation starts by placing the point’s meridian ellipse on a plane, establishing a local 2‑D Cartesian system, then using ellipsoid geometry to compute the normal radius N and the Cartesian coordinates:

Key formulas include the first eccentricity e, the prime vertical radius N, and the final XYZ expressions:

3.1.2. ECEF → BLH

Iterative solution is required for latitude B because it appears in both trigonometric and algebraic terms. The iteration starts with an initial guess and converges to the correct latitude, after which height H is obtained directly.

3.1.3. Implementation

The following C++ code implements the BLH‑to‑Web‑Mercator and inverse conversions, using the WGS‑84 ellipsoid parameters (a = 6378137 m, f⁻¹ = 298.257223563).

#include
//#include
//#include
using namespace std;
const double epsilon = 0.000000000000001;
const double pi = 3.14159265358979323846;
const double d2r = pi / 180;
const double r2d = 180 / pi;
const double a = 6378137.0; // semi‑major axis
const double f_inverse = 298.257223563; // inverse flattening
const double b = a - a / f_inverse;
const double e = sqrt(a * a - b * b) / a;
static double mercatorAngleToGeodeticLatitude(double mercatorAngle) {
    return pi / 2.0 - (2.0 * atan(exp(-mercatorAngle)));
}
static double maximumLatitude = mercatorAngleToGeodeticLatitude(pi);
static double geodeticLatitudeToMercatorAngle(double latitude) {
    if (latitude > maximumLatitude) latitude = maximumLatitude;
    else if (latitude < -maximumLatitude) latitude = -maximumLatitude;
    double sinLatitude = sin(latitude);
    return 0.5 * log((1.0 + sinLatitude) / (1.0 - sinLatitude));
}
void Blh2Wmc(double &x, double &y, double &z) {
    x = x * d2r * a;
    y = geodeticLatitudeToMercatorAngle(y * d2r) * a;
}
void Wmc2Blh(double &x, double &y, double &z) {
    x = x / a * r2d;
    y = mercatorAngleToGeodeticLatitude(y / a) * r2d;
}
int main() {
    double x = 113.6;
    double y = 38.8;
    double z = 100;
    printf("%.10lf\n", maximumLatitude * r2d);
    printf("Original BLH: %.10lf\t%.10lf\t%.10lf\n", x, y, z);
    Blh2Wmc(x, y, z);
    printf("Web Mercator: %.10lf\t%.10lf\t%.10lf\n", x, y, z);
    Wmc2Blh(x, y, z);
    printf("Back to BLH: %.10lf\t%.10lf\t%.10lf\n", x, y, z);
}

Running the program prints the maximum latitude, the original BLH coordinates, the converted Web Mercator values, and the recovered BLH values, demonstrating round‑trip accuracy.

3.2. BLH ↔ Web Mercator Conversion

Web Mercator is a specialization of the transverse Mercator projection. Longitude maps linearly to X, while latitude is transformed using the Mercator formula, limiting usable latitudes to ±85.0511°.

3.3. ECEF ↔ ENU Conversion

3.3.1. Principle

Let the station point P have geodetic coordinates (B, L, H) and corresponding ECEF coordinates (X₀, Y₀, Z₀). ENU is defined with X → East, Y → North, Z → Up.

3.3.1.1. Translation

To convert ENU → ECEF, translate the origin from the station point to the Earth centre using the station’s ECEF coordinates as the translation vector.

3.3.1.2. Rotation

The rotation depends on the station’s latitude B and longitude L. For ENU → ECEF, rotate first around the X‑axis by –B, then around the Z‑axis by –L. The inverse order applies for ECEF → ENU.

The rotation matrices are:

Combining translation and rotation yields the full ENU→ECEF transformation matrix and its inverse for ECEF→ENU.

3.3.2. Implementation

A code snippet (similar to the one above) can be used to perform these conversions, applying the derived matrices to any point expressed in either coordinate system.

4. Extension

The presented concepts are essential for any 3D GIS developer working on the Zhengtong platform or similar systems. While the platform has its own proprietary spatial transformation pipeline, understanding these fundamental geodetic conversions greatly aids in debugging, extending, and optimizing 3‑D rendering pipelines.