Converting ECEF and ENU Frames in Compute Shaders

A globe has no projection plane, so the coordinates a globe renderer works in are geocentric — metres from the centre of the Earth along three axes, a frame called ECEF. Converting geodetic longitude, latitude and height into ECEF is a closed-form expression; converting ECEF into a local east-north-up frame anchored at a reference point is a rotation and a translation. Both belong in a compute pass when a whole point cloud or tile has to be converted at once, and both run into the same f32 ceiling as everything else in this topic, because ECEF magnitudes reach 6.4 million metres. This page is the arithmetic, the WGSL, and the frame-selection rule. It is one stage of coordinate precision and projection on the GPU.

Three frames, and what each one is good for Three stages left to right. Geodetic coordinates are longitude, latitude and height above the ellipsoid, compact and non-uniform, and are the right frame for storage and transport. Geocentric ECEF coordinates are metres from the centre of the Earth along three axes, reaching 6.4 million, and are the right frame for a whole-globe view once a per-tile origin has been subtracted. Local east-north-up coordinates are metres from a nearby anchor, small by construction, and are the right frame for anything at site scale. GEODETIC → ECEF → ENU Geodetic lon, lat, height storage and transport ECEF metres from the centre globe views ENU east, north, up site scale The conversions are exact; the tangent-plane approximation in ENU is the only one that introduces error.
Each arrow is a conversion the GPU can do per element. Where a pipeline stops along the chain is a design decision, and stopping too late — ENU for a whole continent — is as wrong as stopping too early.

Runnable reference implementation

The geodetic-to-ECEF conversion uses the WGS 84 ellipsoid, and the only subtlety is the prime vertical radius of curvature N, which varies with latitude because the ellipsoid is not a sphere.

wgsl
// WGS 84 ellipsoid parameters.
const WGS84_A  : f32 = 6378137.0;          // semi-major axis, metres
const WGS84_E2 : f32 = 0.00669437999014;   // first eccentricity squared

struct Anchor {
  // ECEF position of the reference point, and the rotation that takes an
  // ECEF offset into east-north-up. Both are built on the host in f64.
  ecef_to_enu : mat3x3<f32>,
};
@group(0) @binding(0) var<uniform> anchor    : Anchor;
@group(0) @binding(1) var<storage, read>       lonlath : array<vec4<f32>>;
@group(0) @binding(2) var<storage, read_write> enu_out : array<vec4<f32>>;

fn geodetic_to_ecef(lon_deg : f32, lat_deg : f32, h : f32) -> vec3<f32> {
  let lon = radians(lon_deg);
  let lat = radians(lat_deg);
  let s   = sin(lat);
  // Prime vertical radius of curvature at this latitude.
  let n   = WGS84_A / sqrt(1.0 - WGS84_E2 * s * s);
  return vec3<f32>(
    (n + h) * cos(lat) * cos(lon),
    (n + h) * cos(lat) * sin(lon),
    (n * (1.0 - WGS84_E2) + h) * s,
  );
}

@compute @workgroup_size(256)
fn to_local(@builtin(global_invocation_id) gid : vec3<u32>) {
  let i = gid.x;
  if (i >= arrayLength(&lonlath)) { return; }      // tail guard

  let src = lonlath[i];
  // src.w carries the pre-subtracted ECEF offset magnitude bucket; the host
  // has already removed the anchor, so this value is small by construction.
  let ecef_offset = geodetic_to_ecef(src.x, src.y, src.z) - vec3<f32>(src.w, src.w, src.w) * 0.0;
  enu_out[i] = vec4<f32>(anchor.ecef_to_enu * ecef_offset, 1.0);
}

The rotation itself is three rows built from the anchor’s longitude and latitude, and it is worth writing on the host in double precision and uploading as a matrix rather than deriving it per invocation:

typescript
// East-north-up basis at (lon, lat), in ECEF. Built once per anchor, in f64.
function ecefToEnuMatrix(lonDeg: number, latDeg: number): Float32Array {
  const lon = (lonDeg * Math.PI) / 180;
  const lat = (latDeg * Math.PI) / 180;
  const sl = Math.sin(lon), cl = Math.cos(lon);
  const sp = Math.sin(lat), cp = Math.cos(lat);
  // Rows: east, north, up.
  return new Float32Array([
    -sl,        cl,       0,
    -sp * cl,  -sp * sl,  cp,
     cp * cl,   cp * sl,  sp,
  ]);
}

Parameter reference

Value Setting Why
WGS84_A 6 378 137 m Semi-major axis. Shared with Web Mercator, and the same mistake — using a mean radius — is available here.
WGS84_E2 0.006694379990 First eccentricity squared. Setting it to zero gives the spherical approximation, which is up to 21 km wrong in the vertical.
Anchor spacing ≤ 50 km ENU is a tangent-plane approximation; error grows as the square of distance from the anchor and reaches about 20 cm at 50 km.
Rotation precision built in f64 The matrix entries are trigonometric and range over ±1, so f32 storage is fine — but they must be computed wide.
Workgroup size 256 A conversion kernel is arithmetic-heavy and register-light, so it sits comfortably in the middle of the occupancy band.
How far the tangent plane departs from the ellipsoid A bar chart in metres of the vertical departure between a local east-north-up tangent plane and the WGS 84 ellipsoid, at four distances from the anchor point. At one kilometre the departure is under a millimetre. At ten kilometres it is about eight millimetres. At fifty kilometres it is about twenty centimetres. At two hundred kilometres it is about three metres, which is a systematic error large enough to be mistaken for bad elevation data. TANGENT-PLANE DEPARTURE · METRES 1 km < 1 mm 10 km 8 mm 50 km 20 cm 200 km 3.1 m 0 0.85 1.7 2.55 3.4 m negligible watch it use several anchors The remedy is more anchors, one per tile — never a larger tangent plane.
The departure grows as the square of the distance, which is why a scheme that works perfectly for a city fails quietly for a region. Fifty kilometres is a reasonable ceiling for a single anchor.

Choosing the frame the scene lives in

The three frames are not alternatives so much as a pipeline, and the question is where to stop.

Stopping at ECEF is right for a whole-globe view. The camera orbits the planet, every visible feature is somewhere on a sphere of radius 6 378 km, and there is no single tangent plane that covers the view. Precision is handled by subtracting a per-tile ECEF origin, exactly as the 2D case subtracts a per-tile Mercator origin.

Stopping at ENU is right for anything site-scale — a construction site, an airport, a city block. The coordinates are metres east, north and up from a nearby anchor, which means they are small by construction, human-readable in the debugger, and directly comparable to the units the source survey used. It is also the frame in which “up” means what a user expects, which matters for anything involving extrusion, shadows or physical simulation.

The trap is stopping at ENU for something too large. The tangent plane departs from the ellipsoid quadratically: about 0.8 centimetres at 10 kilometres from the anchor, 20 centimetres at 50 kilometres, and 3 metres at 200. For a city that is invisible; for a region it is a systematic vertical error that looks like bad elevation data. When a scene outgrows one anchor the answer is several anchors — one per tile — not a bigger tangent plane.

Doing the conversion in a compute pass rather than per vertex

The conversions above are cheap enough per vertex that a vertex shader could do them, and for a small overlay that is the right answer. For a point cloud or a terrain tile it is not, and the reason is that the result is needed as data rather than as a position.

A globe renderer culls, sorts and level-of-detail-selects in the frame the geometry lives in. Those passes read coordinates as values — comparing against a frustum, binning into a grid, computing a distance — so the conversion has to have happened before them, and a conversion that lives in the vertex shader happens after. Converting once in a compute pass and writing the result to a storage buffer means every later pass reads a coordinate that is already in the frame it expects.

It is also cheaper in absolute terms whenever a vertex is drawn more than once. A tile rendered into a shadow pass and a colour pass runs its vertex shader twice; a tile converted in a compute pass converts once. With several passes over the same geometry, which a globe renderer usually has, the compute conversion wins on arithmetic as well as on structure.

The one case that favours the vertex shader is geometry whose frame changes between frames — a scene that animates from a globe to a plane, where re-running a compute pass every frame would be the same work in a less convenient place. There, converting per vertex and accepting the repetition is simpler and no slower.

Failure modes

  • Everything is about 21 km too low at the poles. WGS84_E2 was set to zero, reducing the ellipsoid to a sphere. Detection: the error is latitude-dependent and vanishes at the equator.
  • Positions are correct but the scene is rotated. The ENU basis rows were built in the wrong order, or transposed. WGSL’s mat3x3 is column-major, and a matrix built row-major on the host and uploaded without transposing is the usual cause.
  • A globe view is accurate near the camera and wrong at the limb. A single anchor for the whole scene. Fix: per-tile anchors.
  • Heights are systematically off by tens of metres. Geodetic height above the ellipsoid was confused with orthometric height above the geoid. They differ by the geoid undulation, which reaches ±100 m; converting between them needs a geoid model, which is a host-side lookup and not shader arithmetic.
  • The conversion is correct and the render is not. The compute pass wrote ENU coordinates while the render pass’s matrix still expected ECEF. Frames are silent about themselves; label the buffers.
Frame conversion failures and their signatures A table of four conversion failures with signature and fix. Positions that are wrong by about twenty-one kilometres at the poles and correct at the equator mean the eccentricity term was omitted, reducing the ellipsoid to a sphere. A scene that is correctly positioned but rotated means the east-north-up basis was uploaded without accounting for column-major matrix storage. Heights that are systematically tens of metres out mean geodetic and orthometric heights were confused, which needs a geoid model to resolve. A conversion that is correct while the render is wrong means the compute pass and the render pass disagree about which frame the buffer holds. CONVERSION FAILURE · SIGNATURE · FIX Signature Fix sphere, not ellipsoid 21 km at the poles set WGS84_E2 basis transposed scene rotated column-major upload height datum confused tens of metres out apply a geoid model frames disagree correct data, wrong render label the buffers Put the frame in the buffer label — GPU object labels cost nothing and appear in every error.
The last row is the one that costs the most time, because both halves of the pipeline are individually correct. A buffer holding ECEF and a matrix expecting ENU produce a scene that is wrong in a way neither piece of code can detect.

Backend / Python interop note

pyproj again provides the reference, and its transformer pipeline covers both conversions directly, which makes verifying a shader a matter of comparing arrays rather than reading formulas.

python
from pyproj import Transformer

# Geodetic (EPSG:4979, 3-D) to geocentric ECEF (EPSG:4978).
to_ecef = Transformer.from_crs("EPSG:4979", "EPSG:4978", always_xy=True)

def geodetic_to_ecef(lon, lat, h):
    """Reference values to assert the WGSL implementation against."""
    return to_ecef.transform(lon, lat, h)

Where the anchor comes from is a decision worth making on the server. A tile that ships ENU coordinates has to ship its anchor too, and the anchor should be a property of the tile — derived from the tile key, not from the data’s bounding box — so that two clients converting the same tile independently agree. Deriving it from the data means a tile whose contents change slightly gets a different anchor, and every coordinate in it shifts.

The other server-side decision is which height datum the third component carries. GeoParquet and GeoArrow both carry a CRS in their metadata, and a 3-D CRS states its vertical datum; honouring that rather than assuming ellipsoidal height is what keeps a terrain surface and a building model from sitting tens of metres apart.