Coordinate Systems & Geodesy

Overview

Accurate position calculations are fundamental to balloon trajectory prediction. The simulator implements sophisticated geodetic calculations to handle the non-spherical shape of Earth, variations in gravity, and transformations between different coordinate systems. This section covers the mathematical foundations and practical implementations of these systems.

WGS84 Reference Ellipsoid

The World Geodetic System 1984 (WGS84) defines Earth's shape and gravitational field. It serves as the reference for GPS and most modern mapping systems.

Fundamental Parameters

Derived Relationships

Coordinate System Types

1. Geodetic Coordinates (LLA)

Geodetic coordinates specify position using latitude, longitude, and altitude above the ellipsoid:

2. Earth-Centered Earth-Fixed (ECEF)

ECEF is a Cartesian coordinate system that rotates with Earth:

3. Local Tangent Plane (ENU)

East-North-Up coordinates are useful for local calculations:

Coordinate Transformations

Geodetic to ECEF Transformation

Converting from latitude/longitude/height to ECEF coordinates requires the prime vertical radius of curvature:

The transformation equations are:

Implementation

class WGS84:
    def __init__(self):
        self.a = 6378137.0  # Semi-major axis [m]
        self.f = 1/298.257223563  # Flattening
        self.b = self.a * (1 - self.f)  # Semi-minor axis
        self.e2 = 2*self.f - self.f**2  # First eccentricity squared
        self.ep2 = self.e2 / (1 - self.e2)  # Second eccentricity squared

    def geodetic_to_ecef(self, lat_deg, lon_deg, height_m):
        """Convert geodetic coordinates to ECEF"""
        # Convert to radians
        lat = np.radians(lat_deg)
        lon = np.radians(lon_deg)

        # Prime vertical radius of curvature
        N = self.a / np.sqrt(1 - self.e2 * np.sin(lat)**2)

        # ECEF coordinates
        X = (N + height_m) * np.cos(lat) * np.cos(lon)
        Y = (N + height_m) * np.cos(lat) * np.sin(lon)
        Z = (N * (1 - self.e2) + height_m) * np.sin(lat)

        return X, Y, Z

ECEF to Geodetic Transformation

The inverse transformation is more complex and requires iteration. We use Bowring's method:

def ecef_to_geodetic(self, X, Y, Z):
    """Convert ECEF to geodetic coordinates using Bowring's method"""
    # Distance from Z-axis
    p = np.sqrt(X**2 + Y**2)

    # Longitude is straightforward
    lon = np.arctan2(Y, X)

    # Initial values for iteration
    lat = np.arctan2(Z, p * (1 - self.e2))

    # Iterate to refine latitude
    for _ in range(5):  # Usually converges in 2-3 iterations
        N = self.a / np.sqrt(1 - self.e2 * np.sin(lat)**2)
        h = p / np.cos(lat) - N
        lat = np.arctan2(Z, p * (1 - self.e2 * N / (N + h)))

    # Final height calculation
    N = self.a / np.sqrt(1 - self.e2 * np.sin(lat)**2)
    h = p / np.cos(lat) - N

    return np.degrees(lat), np.degrees(lon), h

Local Tangent Plane Transformation

To convert ECEF to local East-North-Up coordinates:

where the rotation matrix $R$ is:

def ecef_to_enu(self, X, Y, Z, lat0_deg, lon0_deg, h0_m):
    """Convert ECEF to local ENU coordinates"""
    # Reference point in ECEF
    X0, Y0, Z0 = self.geodetic_to_ecef(lat0_deg, lon0_deg, h0_m)

    # Difference vector
    dX = X - X0
    dY = Y - Y0
    dZ = Z - Z0

    # Convert reference point to radians
    lat0 = np.radians(lat0_deg)
    lon0 = np.radians(lon0_deg)

    # Rotation matrix elements
    sin_lat = np.sin(lat0)
    cos_lat = np.cos(lat0)
    sin_lon = np.sin(lon0)
    cos_lon = np.cos(lon0)

    # Apply rotation
    E = -sin_lon * dX + cos_lon * dY
    N = -sin_lat * cos_lon * dX - sin_lat * sin_lon * dY + cos_lat * dZ
    U = cos_lat * cos_lon * dX + cos_lat * sin_lon * dY + sin_lat * dZ

    return E, N, U

Gravitational Model

Earth's gravity varies with latitude and altitude due to rotation and the non-spherical shape. The WGS84 gravity model includes these effects.

Normal Gravity Formula

The Somigliana formula gives theoretical gravity on the ellipsoid:

Gravity at Altitude

For points above the ellipsoid, we apply altitude and free-air corrections:

J₂ Gravitational Correction

For high-precision calculations, we include the J₂ term (Earth's oblateness):

Implementation

def calculate_gravity(self, lat_deg, height_m):
    """Calculate gravitational acceleration using WGS84 model"""
    lat = np.radians(lat_deg)

    # Constants
    gamma_e = 9.7803253359  # Equatorial gravity [m/s²]
    k = 0.00193185265241

    # Normal gravity on ellipsoid (Somigliana formula)
    gamma_0 = gamma_e * (1 + k * np.sin(lat)**2) / np.sqrt(1 - self.e2 * np.sin(lat)**2)

    # Altitude correction (to second order)
    f1 = 5.3024e-3
    f2 = -5.9e-6

    # Combined formula
    g = gamma_0 * (1 + f1 * np.sin(lat)**2 + f2 * np.sin(2*lat)**2) * (1 - 2*height_m/self.a)

    return g

Local Tangent Plane Coordinates

Azimuth and Elevation

For tracking and communication, we often need azimuth and elevation angles from a ground station to the balloon:

Great Circle Distance

The shortest distance between two points on the ellipsoid uses Vincenty's formulae:

def vincenty_distance(self, lat1_deg, lon1_deg, lat2_deg, lon2_deg):
    """Calculate geodesic distance using Vincenty's formula"""
    lat1, lon1 = np.radians([lat1_deg, lon1_deg])
    lat2, lon2 = np.radians([lat2_deg, lon2_deg])

    # Reduced latitudes
    U1 = np.arctan((1 - self.f) * np.tan(lat1))
    U2 = np.arctan((1 - self.f) * np.tan(lat2))

    L = lon2 - lon1
    lambda_prev = L

    # Iterate to convergence
    for _ in range(100):
        sin_sigma = np.sqrt(
            (np.cos(U2) * np.sin(lambda_prev))**2 +
            (np.cos(U1) * np.sin(U2) - np.sin(U1) * np.cos(U2) * np.cos(lambda_prev))**2
        )
        cos_sigma = np.sin(U1) * np.sin(U2) + np.cos(U1) * np.cos(U2) * np.cos(lambda_prev)
        sigma = np.arctan2(sin_sigma, cos_sigma)

        sin_alpha = np.cos(U1) * np.cos(U2) * np.sin(lambda_prev) / sin_sigma
        cos2_alpha = 1 - sin_alpha**2
        cos_2sigma_m = cos_sigma - 2 * np.sin(U1) * np.sin(U2) / cos2_alpha

        C = self.f/16 * cos2_alpha * (4 + self.f * (4 - 3 * cos2_alpha))
        lambda_new = L + (1 - C) * self.f * sin_alpha * (
            sigma + C * sin_sigma * (
                cos_2sigma_m + C * cos_sigma * (-1 + 2 * cos_2sigma_m**2)
            )
        )

        if abs(lambda_new - lambda_prev) < 1e-12:
            break
        lambda_prev = lambda_new

    # Calculate distance
    u2 = cos2_alpha * self.ep2
    A = 1 + u2/16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)))
    B = u2/1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2)))
    delta_sigma = B * sin_sigma * (
        cos_2sigma_m + B/4 * (
            cos_sigma * (-1 + 2 * cos_2sigma_m**2) -
            B/6 * cos_2sigma_m * (-3 + 4 * sin_sigma**2) * (-3 + 4 * cos_2sigma_m**2)
        )
    )

    distance = self.b * A * (sigma - delta_sigma)
    return distance

Implementation Examples

Complete Coordinate System Class

class CoordinateSystem:
    """Complete implementation of WGS84 coordinate transformations"""

    def __init__(self):
        # WGS84 parameters
        self.a = 6378137.0
        self.f = 1/298.257223563
        self.b = self.a * (1 - self.f)
        self.e2 = 2*self.f - self.f**2
        self.ep2 = self.e2 / (1 - self.e2)

        # Derived parameters
        self.e = np.sqrt(self.e2)
        self.ep = np.sqrt(self.ep2)

        # Earth rotation rate
        self.omega = 7.292115e-5  # rad/s

    def get_earth_radius(self, lat_deg):
        """Get Earth radius at given latitude"""
        lat = np.radians(lat_deg)

        # Radius of curvature in meridian
        M = self.a * (1 - self.e2) / (1 - self.e2 * np.sin(lat)**2)**1.5

        # Radius of curvature in prime vertical
        N = self.a / np.sqrt(1 - self.e2 * np.sin(lat)**2)

        # Mean radius
        R = np.sqrt(M * N)

        return R

    def get_local_vertical(self, lat_deg, lon_deg):
        """Get local vertical unit vector in ECEF"""
        lat = np.radians(lat_deg)
        lon = np.radians(lon_deg)

        # Unit vector along ellipsoid normal
        nx = np.cos(lat) * np.cos(lon)
        ny = np.cos(lat) * np.sin(lon)
        nz = np.sin(lat)

        return np.array([nx, ny, nz])

    def get_rotation_matrix_ecef_to_enu(self, lat_deg, lon_deg):
        """Get rotation matrix from ECEF to ENU"""
        lat = np.radians(lat_deg)
        lon = np.radians(lon_deg)

        sin_lat = np.sin(lat)
        cos_lat = np.cos(lat)
        sin_lon = np.sin(lon)
        cos_lon = np.cos(lon)

        R = np.array([
            [-sin_lon, cos_lon, 0],
            [-sin_lat*cos_lon, -sin_lat*sin_lon, cos_lat],
            [cos_lat*cos_lon, cos_lat*sin_lon, sin_lat]
        ])

        return R

Trajectory Calculations

class TrajectoryCoordinates:
    """Handle coordinate transformations for balloon trajectories"""

    def __init__(self):
        self.cs = CoordinateSystem()

    def calculate_downrange_distance(self, launch_lat, launch_lon, current_lat, current_lon):
        """Calculate great circle distance from launch"""
        return self.cs.vincenty_distance(launch_lat, launch_lon, current_lat, current_lon)

    def calculate_bearing(self, lat1, lon1, lat2, lon2):
        """Calculate initial bearing from point 1 to point 2"""
        lat1, lon1, lat2, lon2 = np.radians([lat1, lon1, lat2, lon2])

        dlon = lon2 - lon1

        y = np.sin(dlon) * np.cos(lat2)
        x = np.cos(lat1) * np.sin(lat2) - np.sin(lat1) * np.cos(lat2) * np.cos(dlon)

        bearing_rad = np.arctan2(y, x)
        bearing_deg = np.degrees(bearing_rad)

        # Normalize to 0-360
        return (bearing_deg + 360) % 360

    def project_position(self, lat, lon, bearing_deg, distance_m):
        """Project a position given bearing and distance"""
        lat1 = np.radians(lat)
        lon1 = np.radians(lon)
        bearing = np.radians(bearing_deg)

        # Angular distance
        delta = distance_m / self.cs.get_earth_radius(lat)

        # Calculate new position
        lat2 = np.arcsin(
            np.sin(lat1) * np.cos(delta) +
            np.cos(lat1) * np.sin(delta) * np.cos(bearing)
        )

        lon2 = lon1 + np.arctan2(
            np.sin(bearing) * np.sin(delta) * np.cos(lat1),
            np.cos(delta) - np.sin(lat1) * np.sin(lat2)
        )

        return np.degrees(lat2), np.degrees(lon2)

Accuracy Considerations

Sources of Error

Best Practices

# Good practice: Keep calculations in consistent units
class SafeCoordinateCalculator:
    def __init__(self):
        self.cs = CoordinateSystem()

    def safe_distance_calculation(self, lat1_deg, lon1_deg, alt1_m,
                                 lat2_deg, lon2_deg, alt2_m):
        """Calculate 3D distance with proper altitude handling"""

        # Convert to ECEF for accurate 3D distance
        x1, y1, z1 = self.cs.geodetic_to_ecef(lat1_deg, lon1_deg, alt1_m)
        x2, y2, z2 = self.cs.geodetic_to_ecef(lat2_deg, lon2_deg, alt2_m)

        # Euclidean distance in ECEF
        distance_3d = np.sqrt((x2-x1)**2 + (y2-y1)**2 + (z2-z1)**2)

        # Surface distance for comparison
        surface_distance = self.cs.vincenty_distance(lat1_deg, lon1_deg,
                                                    lat2_deg, lon2_deg)

        return {
            'distance_3d': distance_3d,
            'surface_distance': surface_distance,
            'altitude_difference': alt2_m - alt1_m
        }