Local Coordinates and the Ellipsoid

library(geographiclib)

Example Locations
Geocentric (ECEF) Coordinates
Local Cartesian (ENU) Coordinates
WGS84 Ellipsoid Properties
Combining Coordinate Systems
Coordinate System Summary
See Also

Example Locations

# Locations at various latitudes
world_pts <- cbind(
  lon = c(151.21, 0, -74.01, 0, 0),
  lat = c(-33.87, 51.51, 40.71, 0, -90)
)
rownames(world_pts) <- c("Sydney", "London", "New York", "Equator/Prime", "South Pole")

# Antarctic stations
antarctic <- cbind(
  lon = c(166.67, 77.85, -68.13, 39.58, 0),
  lat = c(-77.85, -67.60, -67.57, -69.41, -90)
)
rownames(antarctic) <- c("McMurdo", "Davis", "Palmer", "Progress", "South Pole")

Geocentric (ECEF) Coordinates

Earth-Centered Earth-Fixed (ECEF) coordinates express positions as X, Y, Z relative to the Earth’s center:

X-axis: Points toward 0°N, 0°E (intersection of equator and prime meridian)
Y-axis: Points toward 0°N, 90°E (intersection of equator and 90°E)
Z-axis: Points toward the North Pole

Basic Conversion

geocentric_fwd(world_pts)
#>          X        Y        Z    lon    lat h
#> 1 -4646017  2553114 -3534483 151.21 -33.87 0
#> 2  3977778        0  4969055   0.00  51.51 0
#> 3  1333729 -4654333  4138065 -74.01  40.71 0
#> 4  6378137        0        0   0.00   0.00 0
#> 5        0        0 -6356752   0.00 -90.00 0

Understanding ECEF

# Points on the equator have Z ≈ 0
equator_pts <- cbind(
  lon = c(0, 90, 180, -90),
  lat = c(0, 0, 0, 0)
)
geocentric_fwd(equator_pts)
#>          X        Y Z lon lat h
#> 1  6378137        0 0   0   0 0
#> 2        0  6378137 0  90   0 0
#> 3 -6378137        0 0 180   0 0
#> 4        0 -6378137 0 -90   0 0

# Points at poles have X ≈ 0, Y ≈ 0
pole_pts <- cbind(lon = c(0, 0), lat = c(90, -90))
geocentric_fwd(pole_pts)
#>   X Y        Z lon lat h
#> 1 0 0  6356752   0  90 0
#> 2 0 0 -6356752   0 -90 0

Height Above Ellipsoid

ECEF can include height above the WGS84 ellipsoid:

# Ground level vs airplane altitude (10km) vs satellite (400km)
pt <- c(151.21, -33.87)

data.frame(
  location = c("Ground", "Aircraft (10km)", "ISS (~400km)", "GPS satellite (~20200km)"),
  h = c(0, 10000, 400000, 20200000),
  geocentric_fwd(cbind(rep(pt[1], 4), rep(pt[2], 4)), 
                 h = c(0, 10000, 400000, 20200000))
)
#>                   location        h         X        Y         Z    lon    lat
#> 1                   Ground        0  -4646017  2553114  -3534483 151.21 -33.87
#> 2          Aircraft (10km)    10000  -4653294  2557113  -3540056 151.21 -33.87
#> 3             ISS (~400km)   400000  -4937086  2713064  -3757407 151.21 -33.87
#> 4 GPS satellite (~20200km) 20200000 -19344970 10630591 -14792154 151.21 -33.87
#>        h.1
#> 1        0
#> 2    10000
#> 3   400000
#> 4 20200000

Antarctic Stations in ECEF

# Antarctic stations are all near the "bottom" of the coordinate system
geocentric_fwd(antarctic)
#>            X          Y        Z    lon    lat h
#> 1 -1310449.4   310501.7 -6213433 166.67 -77.85 0
#> 2   513025.6  2382902.9 -5874236  77.85 -67.60 0
#> 3   909126.8 -2264950.0 -5872961 -68.13 -67.57 0
#> 4  1733893.8  1433382.5 -5948210  39.58 -69.41 0
#> 5        0.0        0.0 -6356752   0.00 -90.00 0

# Note the large negative Z values (Southern Hemisphere)
# and relatively small X, Y values (near the axis)

Round-trip Conversion

fwd <- geocentric_fwd(world_pts)
rev <- geocentric_rev(fwd$X, fwd$Y, fwd$Z)

# Verify accuracy
max(abs(rev$lon - world_pts[,1]))
#> [1] 0
max(abs(rev$lat - world_pts[,2]))
#> [1] 1.421085e-14
max(abs(rev$h))  # Height should be ~0
#> [1] 1.409644e-09

GPS Applications

ECEF is the native coordinate system for GPS satellites:

# Convert GPS receiver position to geodetic
# (Example: receiver at Sydney at ~100m altitude)
X <- 4648241   # meters
Y <- -2560342
Z <- -3526276

geocentric_rev(X, Y, Z)
#>         lon       lat         h       X        Y        Z
#> 1 -28.84683 -33.78127 -55.64177 4648241 -2560342 -3526276

Local Cartesian (ENU) Coordinates

Local Cartesian, also called ENU (East-North-Up), creates a local coordinate system with:

East: Positive X direction
North: Positive Y direction
Up: Positive Z direction

This is ideal for local surveys, robotics, and navigation.

Basic Conversion

# Set up local system centered on Sydney
sydney <- c(151.21, -33.87)

nearby_pts <- cbind(
  lon = c(151.21, 151.31, 151.11, 151.21, 151.21),
  lat = c(-33.87, -33.87, -33.87, -33.77, -33.97)
)
rownames(nearby_pts) <- c("Origin", "East 10km", "West 10km", "North 10km", "South 10km")

localcartesian_fwd(nearby_pts, lon0 = sydney[1], lat0 = sydney[2])
#>               x             y         z    lon    lat h
#> 1  0.000000e+00      0.000000  0.000000 151.21 -33.87 0
#> 2  9.252524e+03     -4.499921 -6.704168 151.31 -33.87 0
#> 3 -9.252524e+03     -4.499921 -6.704168 151.11 -33.87 0
#> 4  5.275069e-11  11091.908279 -9.679492 151.21 -33.77 0
#> 5 -2.150955e-10 -11092.088563 -9.679702 151.21 -33.97 0

Local Survey Application

# Survey points around McMurdo Station
mcmurdo <- c(166.67, -77.85)

survey_pts <- cbind(
  lon = c(166.67, 166.70, 166.64, 166.67, 166.73),
  lat = c(-77.85, -77.85, -77.85, -77.84, -77.86)
)
rownames(survey_pts) <- c("Base", "Point A", "Point B", "Point C", "Point D")

result <- localcartesian_fwd(survey_pts, lon0 = mcmurdo[1], lat0 = mcmurdo[2])
result
#>               x            y          z    lon    lat h
#> 1  0.000000e+00     0.000000  0.0000000 166.67 -77.85 0
#> 2  7.051476e+02    -0.180472 -0.0388546 166.70 -77.85 0
#> 3 -7.051476e+02    -0.180472 -0.0388546 166.64 -77.85 0
#> 4 -1.907097e-11  1116.439380 -0.0974277 166.67 -77.84 0
#> 5  1.409152e+03 -1117.161493 -0.2527202 166.73 -77.86 0

# Distances from base in meters
data.frame(
  point = rownames(survey_pts),
  east_m = round(result$x),
  north_m = round(result$y),
  distance_m = round(sqrt(result$x^2 + result$y^2))
)
#>     point east_m north_m distance_m
#> 1    Base      0       0          0
#> 2 Point A    705       0        705
#> 3 Point B   -705       0        705
#> 4 Point C      0    1116       1116
#> 5 Point D   1409   -1117       1798

Including Height

# Local system with height differences
# Simulating a hill near Sydney
pts_with_height <- cbind(
  lon = c(151.21, 151.22, 151.20),
  lat = c(-33.87, -33.87, -33.86)
)
heights <- c(0, 50, 100)  # meters above ellipsoid

result <- localcartesian_fwd(pts_with_height, 
                              lon0 = 151.21, lat0 = -33.87, h0 = 0,
                              h = heights)
result
#>           x             y        z    lon    lat   h
#> 1    0.0000    0.00000000  0.00000 151.21 -33.87   0
#> 2  925.2601   -0.04499957 49.93296 151.22 -33.87  50
#> 3 -925.3752 1109.17194197 99.83615 151.20 -33.86 100

Round-trip Conversion

pts <- cbind(
  lon = c(166.67, 166.70, 166.64),
  lat = c(-77.85, -77.84, -77.86)
)

fwd <- localcartesian_fwd(pts, lon0 = 166.67, lat0 = -77.85)
rev <- localcartesian_rev(fwd$x, fwd$y, fwd$z, lon0 = 166.67, lat0 = -77.85)

max(abs(rev$lon - pts[,1]))
#> [1] 0
max(abs(rev$lat - pts[,2]))
#> [1] 0

Robotics/Navigation Application

# Robot path planning at Davis Station
davis <- c(77.85, -67.60)

# Waypoints for a robot traverse
waypoints_enu <- data.frame(
  name = c("Start", "WP1", "WP2", "WP3", "End"),
  x = c(0, 100, 200, 250, 300),    # East (meters)
  y = c(0, 50, 100, 100, 150),     # North (meters)
  z = c(0, 0, 0, 0, 0)             # Up (meters)
)

# Convert to geographic coordinates for GPS navigation
result <- localcartesian_rev(
  waypoints_enu$x, waypoints_enu$y, waypoints_enu$z,
  lon0 = davis[1], lat0 = davis[2]
)

data.frame(
  waypoint = waypoints_enu$name,
  lon = round(result$lon, 6),
  lat = round(result$lat, 6)
)
#>   waypoint      lon       lat
#> 1    Start 77.85000 -67.60000
#> 2      WP1 77.85235 -67.59955
#> 3      WP2 77.85470 -67.59910
#> 4      WP3 77.85588 -67.59910
#> 5      End 77.85705 -67.59865

WGS84 Ellipsoid Properties

The WGS84 ellipsoid is the reference surface for GPS and most modern mapping.

Basic Parameters

ellipsoid_params()
#> $a
#> [1] 6378137
#> 
#> $f
#> [1] 0.003352811
#> 
#> $b
#> [1] 6356752
#> 
#> $e2
#> [1] 0.00669438
#> 
#> $ep2
#> [1] 0.006739497
#> 
#> $n
#> [1] 0.00167922
#> 
#> $area
#> [1] 5.100656e+14
#> 
#> $volume
#> [1] 1.083207e+21

Key parameters: - a: Semi-major axis (equatorial radius) ≈ 6,378,137 m - b: Semi-minor axis (polar radius) ≈ 6,356,752 m - f: Flattening ≈ 1/298.257 - e2: First eccentricity squared - area: Total surface area - volume: Total volume

Earth’s Shape

params <- ellipsoid_params()

# Equatorial vs polar radius difference
equatorial_radius <- params$a
polar_radius <- params$b

cat("Equatorial radius:", equatorial_radius, "m\n")
#> Equatorial radius: 6378137 m
cat("Polar radius:", polar_radius, "m\n")
#> Polar radius: 6356752 m
cat("Difference:", equatorial_radius - polar_radius, "m\n")
#> Difference: 21384.69 m
cat("Flattening:", 1/params$f, "(1/f)\n")
#> Flattening: 298.2572 (1/f)

Radii of Curvature

The Earth’s curvature varies with latitude:

# Curvature at various latitudes
lats <- c(0, -33.87, -42.88, -67.60, -77.85, -90)
names_lat <- c("Equator", "Sydney", "Hobart", "Davis", "McMurdo", "South Pole")

result <- ellipsoid_curvature(lats)

data.frame(
  location = names_lat,
  latitude = lats,
  meridional_km = round(result$meridional / 1000, 2),
  transverse_km = round(result$transverse / 1000, 2)
)
#>     location latitude meridional_km transverse_km
#> 1    Equator     0.00       6335.44       6378.14
#> 2     Sydney   -33.87       6355.25       6384.78
#> 3     Hobart   -42.88       6365.01       6388.05
#> 4      Davis   -67.60       6390.21       6396.46
#> 5    McMurdo   -77.85       6396.73       6398.64
#> 6 South Pole   -90.00       6399.59       6399.59

Circle of Latitude Properties

# Properties of circles at different latitudes
lats <- c(0, -30, -45, -60, -75, -90)

result <- ellipsoid_circle(lats)

data.frame(
  latitude = lats,
  circle_radius_km = round(result$radius / 1000, 2),
  meridian_dist_km = round(result$meridian_distance / 1000, 2)
)
#>   latitude circle_radius_km meridian_dist_km
#> 1        0          6378.14             0.00
#> 2      -30          5528.26         -3320.11
#> 3      -45          4517.59         -4984.94
#> 4      -60          3197.10         -6654.07
#> 5      -75          1655.96         -8326.94
#> 6      -90             0.00        -10001.97

Auxiliary Latitudes

For advanced geodetic calculations, various auxiliary latitudes are used:

# Compare different latitude types at 45°S
lats <- c(0, -30, -45, -60, -90)

result <- ellipsoid_latitudes(lats)
result
#>   lat parametric geocentric rectifying  authalic conformal isometric
#> 1   0    0.00000    0.00000    0.00000   0.00000   0.00000   0.00000
#> 2 -30  -29.91675  -29.83364  -29.87515 -29.88900 -29.83368 -31.28104
#> 3 -45  -44.90379  -44.80758  -44.85568 -44.87170 -44.80768 -50.22747
#> 4 -60  -59.91661  -59.83308  -59.87489 -59.88879 -59.83322 -75.12340
#> 5 -90  -90.00000  -90.00000  -90.00000 -90.00000 -90.00000      -Inf

# The differences are small but matter for precise calculations

Combining Coordinate Systems

GNSS Data Processing

# Typical GNSS workflow:
# 1. Receive ECEF coordinates from GPS
# 2. Convert to geodetic (lat/lon/height)
# 3. Convert to local for navigation

# Example GPS receiver data (ECEF, meters)
gps_ecef <- data.frame(
  time = 1:5,
  X = c(4648241, 4648242, 4648240, 4648243, 4648241),
  Y = c(-2560342, -2560340, -2560343, -2560341, -2560342),
  Z = c(-3526276, -3526275, -3526277, -3526274, -3526276)
)

# Convert to geodetic
geodetic <- geocentric_rev(gps_ecef$X, gps_ecef$Y, gps_ecef$Z)

# Convert to local (relative to first point)
local <- localcartesian_fwd(
  cbind(geodetic$lon, geodetic$lat),
  lon0 = geodetic$lon[1], lat0 = geodetic$lat[1], h0 = geodetic$h[1],
  h = geodetic$h
)

data.frame(
  time = gps_ecef$time,
  east_m = round(local$x, 2),
  north_m = round(local$y, 2),
  up_m = round(local$z, 2)
)
#>   time east_m north_m  up_m
#> 1    1   0.00    0.00  0.00
#> 2    2   2.23    0.78 -0.63
#> 3    3  -1.36   -1.05  0.23
#> 4    4   1.84    2.37 -0.06
#> 5    5   0.00    0.00  0.00

Antarctic Field Survey

# Simulated survey data at McMurdo
mcmurdo <- c(166.67, -77.85, 10)  # lon, lat, height

# Survey points in local coordinates
survey_local <- data.frame(
  point = c("Control", "P1", "P2", "P3", "P4"),
  east = c(0, 500, -300, 200, -100),
  north = c(0, 200, 400, -150, -300),
  up = c(0, 5, -2, 8, -5)
)

# Convert to geodetic for GPS upload
geodetic <- localcartesian_rev(
  survey_local$east, survey_local$north, survey_local$up,
  lon0 = mcmurdo[1], lat0 = mcmurdo[2], h0 = mcmurdo[3]
)

# Convert to ECEF for satellite positioning
ecef <- geocentric_fwd(
  cbind(geodetic$lon, geodetic$lat),
  h = geodetic$h
)

data.frame(
  point = survey_local$point,
  lon = round(geodetic$lon, 5),
  lat = round(geodetic$lat, 5),
  X = round(ecef$X),
  Y = round(ecef$Y),
  Z = round(ecef$Z)
)
#>     point      lon       lat        X      Y        Z
#> 1 Control 166.6700 -77.85000 -1310451 310502 -6213443
#> 2      P1 166.6913 -77.84821 -1310758 310061 -6213406
#> 3      P2 166.6572 -77.84642 -1310762 310884 -6213357
#> 4      P3 166.6785 -77.85134 -1310357 310274 -6213482
#> 5      P4 166.6657 -77.85269 -1310142 310532 -6213501

Coordinate System Summary

System	Description	Use Case
Geodetic (lon/lat/h)	Geographic coordinates	Maps, GIS, human communication
ECEF (X/Y/Z)	Earth-centered Cartesian	GPS satellites, orbit calculations
ENU (E/N/U)	Local tangent plane	Robotics, local surveys, navigation

Contents

Example Locations

Geocentric (ECEF) Coordinates

Basic Conversion

Understanding ECEF

Height Above Ellipsoid

Antarctic Stations in ECEF

Round-trip Conversion

GPS Applications

Local Cartesian (ENU) Coordinates

Basic Conversion

Local Survey Application

Including Height

Round-trip Conversion

Robotics/Navigation Application

WGS84 Ellipsoid Properties

Basic Parameters

Earth’s Shape

Radii of Curvature

Circle of Latitude Properties

Auxiliary Latitudes

Combining Coordinate Systems

GNSS Data Processing

Antarctic Field Survey

Coordinate System Summary

See Also