Chapter 3: WGS84 Coordinate System and ENU

Introduction

WGS84, ECEF, and ENU are coordinate systems familiar to GIS professionals, but since GIS software handles the internal processing, there are few opportunities to implement the internals yourself. In this chapter, we will build the coordinate transformations that are typically used as a black box, step by step.

As a 3D engine-specific consideration, there is the axis mapping between Three.js's Y-up convention and ECEF's Z-up (North Pole direction). By clicking on city markers, you can inspect the values in each coordinate system and experience how the transformation pipeline works.

What You Will Learn

• WGS84 Ellipsoid definition and parameters — Visualizing the GIS foundation as a Three.js sphere
• ECEF Coordinate System: Earth-centered, Earth-fixed Cartesian coordinates — Experience GPS's internal coordinate system in 3D space
• Conversion from geodetic coordinates (latitude, longitude, altitude) to ECEF — The geometric meaning of the prime vertical radius of curvature N
• ECEF to Three.js scene coordinate axis remapping — Conversion rules between Y-up (Three.js) and Z-up (ECEF)
• ENU Coordinate System: A local coordinate system on the surface — Verify cardinal directions in 3D
• Rotation matrix for ECEF to ENU conversion — Derivation and implementation of the orthogonal transformation matrix

Demo Visual Elements

Why an Ellipsoid Is Needed

The Earth is not a perfect sphere. Due to centrifugal force from its rotation, it bulges slightly at the equator, forming a shape close to an oblate spheroid. The mathematical model of this shape is called a reference ellipsoid.

WGS84 (World Geodetic System 1984) defines the standard reference ellipsoid that best approximates the Earth's shape. When working with latitude and longitude in GIS, you are effectively working with coordinates on the WGS84 ellipsoid.

WGS84 Ellipsoid Parameters

The WGS84 ellipsoid is defined by two fundamental parameters.

The difference between the equatorial and polar radii is about 21 km. While this is only about 0.3% of the Earth's diameter (approx. 12,756 km), it cannot be ignored in precise coordinate transformations.

export const WGS84 = {
  /** Semi-major axis (equatorial radius) [m] */
  a: 6_378_137.0,
  /** Flattening */
  f: 1 / 298.257_223_563,
  /** Semi-minor axis (polar radius) [m] */
  get b(): number {
    return this.a * (1 - this.f);
  },
  /** First eccentricity squared */
  get e2(): number {
    return 2 * this.f - this.f * this.f;
  }
} as const;

The reason e² is derived directly from f rather than using a and b is to maintain floating-point precision.

In the demo on the left, the Y-axis direction is scaled by WGS84.b / WGS84.a to reproduce this oblate shape. The transparency is reduced so that the internal structure is visible.

What Is the ECEF Coordinate System

ECEF (Earth-Centered, Earth-Fixed) is a Cartesian coordinate system with its origin at the center of the Earth, rotating together with the Earth. Position information received from GPS satellites is internally processed in ECEF.

Geodetic to ECEF Conversion

The conversion from latitude, longitude, and altitude to ECEF uses the ellipsoid's prime vertical radius of curvature (N). N is the "east-west radius of curvature" at a given latitude on the ellipsoid surface, which is the distance from a point on the ellipsoid surface along the normal direction to the Z-axis (rotation axis).

Note the formula for the Z component. While the X and Y components use (N + h), the Z component has (1 - e²) multiplied. This is due to the oblateness of the ellipsoid. If e² = 0 (a perfect sphere), then (1-e²) = 1, and this difference vanishes.

export function geodeticToEcef(
  lat: number, lon: number, alt: number
): EcefCoord {
  const phi = degToRad(lat);
  const lam = degToRad(lon);
  const sinPhi = Math.sin(phi);
  const cosPhi = Math.cos(phi);
  const sinLam = Math.sin(lam);
  const cosLam = Math.cos(lam);

  const N = WGS84.a / Math.sqrt(
    1 - WGS84.e2 * sinPhi * sinPhi
  );

  return {
    x: (N + alt) * cosPhi * cosLam,
    y: (N + alt) * cosPhi * sinLam,
    z: (N * (1 - WGS84.e2) + alt) * sinPhi
  };
}

Inverse Conversion: ECEF to Geodetic

The inverse conversion (ECEF to geodetic) is not as straightforward as the forward conversion. In practice, Bowring's iterative method is widely used. In this implementation, the initial estimate using the auxiliary angle θ = atan2(z·a, p·b) is sufficiently accurate, so sub-millimeter precision is effectively achieved in a single calculation.

Processing steps:

• Longitude: Directly computed as atan2(Y, X) (the angle in the ECEF X-Y plane = longitude)
• Auxiliary angle: Geometrically estimated via atan2(z*a, p*b) (this initial value alone provides sufficient precision)
• Latitude calculation: Latitude φ is computed using the auxiliary angle θ
• Ellipsoidal height: Different formulas are used near the equator vs. near the poles

export function ecefToGeodetic(
  x: number, y: number, z: number
): GeodeticCoord {
  const p = Math.sqrt(x*x + y*y);
  const lam = Math.atan2(y, x);

  // Bowring initial estimate
  let theta = Math.atan2(z*a, p*b);
  let phi = Math.atan2(
    z + ep2 * b * Math.sin(theta)**3,
    p - e2 * a * Math.cos(theta)**3
  );

  // 3 iterations are sufficient for precision
  for (let i = 0; i < 3; i++) {
    theta = Math.atan2(z*a, p*b);
    phi = Math.atan2(/* ... */);
  }
  // ... height calculation
}

Coordinate Inspection

Click on a city marker in the demo on the left to inspect the values in each coordinate system. Experience how the transformation pipeline works.

Click a city marker to display its coordinates

Why Axis Remapping Is Needed

The ECEF coordinate system and Three.js scene coordinate system have different axis orientations.

To align these, an axis remapping from ECEF coordinates to Three.js scene coordinates is required.

Mapping to Scene Coordinates

The sign inversion (-Y) is necessary to maintain a right-handed coordinate system. When swapping two axes in the conversion from ECEF (right-handed) to Three.js (right-handed), one axis must be negated to avoid producing a left-handed system.

export function ecefToScenePosition(
  x: number, y: number, z: number
): EcefCoord {
  return {
    x: x,     // ECEF X → Scene X
    y: z,     // ECEF Z → Scene Y (up)
    z: -y     // ECEF Y → Scene Z (negated)
  };
}

In this book, 1 Three.js unit = 1 meter, with no scaling applied. CesiumJS also adopts a 1:1 scale. Introducing scaling would make tile texture resolution calculations scale-dependent and complicate coordinate alignment with external data such as 3D Tiles.

Axis Remapping in Matrix Form

The axis remapping can also be expressed as a 4x4 matrix. The 3D Tiles adapter uses this matrix form.

// ECEF → Scene coordinate axis remapping matrix
const ecefToScene = new THREE.Matrix4().set(
  1,  0,  0,  0,   // X → X
  0,  0,  1,  0,   // Z → Y
  0, -1,  0,  0,   // -Y → Z
  0,  0,  0,  1
);
this.group.applyMatrix4(ecefToScene);

The custom 3D Tiles adapter implemented in Chapter 9 places models in ECEF coordinates, so by applying the axis remapping matrix to the Group, all child objects are automatically transformed into scene coordinates.

ENU Coordinate System

ENU (East-North-Up) is a local Cartesian coordinate system centered on a reference point on the Earth's surface.

Since the ENU coordinate system changes when the reference point changes, it is called a "local" coordinate system. ENU is used in the following contexts:

• Converting local coordinates in 3D Tiles tilesets to ECEF
• Placing objects near the surface (buildings, roads, etc.)
• Survey data processing, GPS relative positioning

Rotation Matrix for ECEF to ENU Conversion

The ECEF-to-ENU conversion is achieved by multiplying the difference vector by a rotation matrix R.

Meaning of each row:

• Row 1 (East): Unit vector in the east direction. A horizontal vector perpendicular to the longitude direction
• Row 2 (North): Unit vector in the north direction. A horizontal vector along the latitude direction
• Row 3 (Up): Normal direction of the ellipsoid. Coincides with the normal vector in ECEF

The inverse conversion (ENU to ECEF) uses the transpose of the rotation matrix (since the matrix is orthogonal, the transpose equals the inverse).

export function enuToEcef(
  e: number, n: number, u: number,
  latRef: number, lonRef: number,
  altRef: number
): EcefCoord {
  const ref = geodeticToEcef(
    latRef, lonRef, altRef
  );
  const phi = degToRad(latRef);
  const lam = degToRad(lonRef);
  // Transpose of rotation matrix (inverse conversion)
  return {
    x: ref.x
      + (-sinLam*e - sinPhi*cosLam*n
         + cosPhi*cosLam*u),
    y: ref.y
      + (cosLam*e - sinPhi*sinLam*n
         + cosPhi*sinLam*u),
    z: ref.z
      + (cosPhi*n + sinPhi*u)
  };
}

Drawing the ENU Axes

In the demo, enuToEcef converts each axis's unit vector to ECEF, and the direction is determined from the difference with the origin to draw lines.

// Convert the east unit vector to ECEF
const tipEcef = enuToEcef(
  1, 0, 0,       // E=1, N=0, U=0
  tokyoLat, tokyoLon, 0
);
// Get direction from difference with origin
const dir = tip.sub(origin).normalize();
// Multiply direction by enuAxisLen for line endpoint
const end = tokyoOrigin.clone()
  .add(dir.multiplyScalar(enuAxisLen));

Round-Trip Conversion Accuracy

The correctness of the implementation can be verified with round-trip conversion tests. The tests perform round-trip conversions at various points including on the equator, at the poles, Tokyo, and the Southern Hemisphere, confirming that latitude and longitude errors are below 10^-10 degrees.

it('round-trip for Tokyo', () => {
  const ecef =
    geodeticToEcef(35.6762, 139.6503, 100);
  const result =
    ecefToGeodetic(ecef.x, ecef.y, ecef.z);
  expect(result.lat)
    .toBeCloseTo(35.6762, 9);
  expect(result.lon)
    .toBeCloseTo(139.6503, 9);
  expect(result.alt)
    .toBeCloseTo(100, 3); // within 1mm
});

Summary

In this chapter, we implemented the coordinate systems and transformations that form the foundation of GIS.

• WGS84 Parameters: Derived semi-minor axis b and eccentricity e² from semi-major axis a and flattening f
• Forward Conversion: Geodetic to ECEF conversion using the prime vertical radius of curvature N
• Inverse Conversion: ECEF to geodetic conversion using Bowring's iterative method (3 iterations)
• Axis Remapping: X→X, Z→Y, -Y→Z for ECEF to scene coordinates
• ENU Coordinate System: East-North-Up local coordinates centered on a reference point, with rotation matrix

In the next chapter, we will cover the Web Mercator projection and tile coordinate system, which form the basis for map tiles.