Ellipsoids

Ellipsoid is a three-dimensional shape in geometry formed by rotating an ellipse around its major or minor axis. The shape of the Earth is not a perfect sphere but is closer to a flattened sphere, hence referred to as the “Earth ellipsoid.” Specifically, the Earth is slightly bulging at the equator and slightly flattened at the poles, which is known as an oblate spheroid, a special type of ellipsoid.

When conducting Earth science research, calculating distances is unavoidable.

Here I have compiled an introduction to different ellipsoids and their parameters:

1. Airy (1830):

• Major axis radius: 6377.563 km, Minor axis radius: 6356.257 km
• This ellipsoid was proposed by George Biddell Airy in 1830 for the geographical measurements of the British Empire. It was the standard ellipsoid used in early geographical surveys in the UK.

Bessel:

• Major axis radius: 6377.397 km, Minor axis radius: 6356.079 km
• Proposed by German mathematician Friedrich Wilhelm Bessel in 1841, it was widely used in some surveying works in Europe and Asia. Its calculations were based on a large amount of geodetic measurement data.

Clarke (1880):

• Major axis radius: 6378.249145 km, Minor axis radius: 6356.51486955 km
• The Clarke 1880 ellipsoid was proposed by Alexander Ross Clarke, based on updated measurement data. It was widely used in Africa and Europe.

FAI sphere:

• Both major and minor axis radii: 6371 km
• This is an idealized spherical model that uses a uniform radius (6371 km) to represent the Earth. The International Aeronautical Federation (FAI) sometimes uses it for simple calculations.

GRS-67:

• Major axis radius: 6378.160 km, Minor axis radius: 6356.775 km
• The Geodetic Reference System 1967 (GRS-67) was initially proposed by the International Association of Geodesy (IAG) in 1967. This model provided a reference for many maps and measurements in the 1970s.

International:

• Major axis radius: 6378.388 km, Minor axis radius: 6356.912 km
• This is an ellipsoid widely used in 20th-century international geodetic measurements, especially in Europe. It is often referred to as the “International 1924” ellipsoid.

Krasovsky:

• Major axis radius: 6378.245 km, Minor axis radius: 6356.863 km
• Proposed by Soviet scientist Fyodor Krasovsky in 1940, it was used for geographical measurements in the former Soviet Union and its allies.

NAD27 (North American Datum 1927):

• Major axis radius: 6378.206 km, Minor axis radius: 6356.584 km
• NAD27 was developed based on the Clarke 1866 ellipsoid and is the standard for the North American geographical reference system, widely used in 20th-century maps and measurements.

WGS66 (World Geodetic System 1966):

• Major axis radius: 6378.145 km, Minor axis radius: 6356.758 km
• The World Geodetic System 1966 is an early global Earth reference system, primarily developed by the U.S. Department of Defense, initially used for navigation and satellite orbit calculations.

WGS72 (World Geodetic System 1972):

• Major axis radius: 6378.135 km, Minor axis radius: 6356.751 km
• The World Geodetic System 1972 is a geodetic reference system developed by the U.S. in the 1970s, applied to navigation for Earth and satellite observations.

WGS84 (World Geodetic System 1984):

• Major axis radius: 6378.1370 km, Minor axis radius: 6356.7523 km
• WGS84 is the most widely used Earth ellipsoid model, proposed by the U.S. Department of Defense in 1984 and continuously updated, serving as the standard ellipsoid for GPS (Global Positioning System).

IERS (2003):

• Major axis radius: 6378.1366 km, Minor axis radius: 6356.7519 km
• The International Earth Reference System (IERS) proposed in 2003 provides more accurate Earth ellipsoid parameters and is used for high-precision scientific calculations and navigation.

These different ellipsoid models were proposed based on data from various periods and different measurement needs. Each ellipsoid model reflects the best estimate of the Earth’s shape at the time, and the differences among these ellipsoids lie in their measurement accuracy and applicability in different regions and contexts. Modern global positioning and navigation systems primarily use the WGS84 ellipsoid, which provides a consistent and accurate reference across the globe.

Implementation in Python

import math
import numpy as np
def haversine(origin, destination,units='km',ellipsoid="WGS84"):
    """
    To calculate distance between two points in Earth giving their coordinates (lat,lon)
    Parameters
    ----------
    origin:array like (lat,lon)
    coordinates of origin point

    destination: array like (lat,lon)
    coordinates of destinations points

    units: str
    units to return distance
    available units are kilometers (km), meters (m) and miles (mi)
    
    ellipsoid: String: type of projection
    available: Airy (1830), Bessel, Clarke (1880), FAI sphere, GRS-67, International, Krasovsky, NAD27, WGS66, WGS72, IERS (2003), WGS84 - default WGS84
    return
    distance between points
    """
    if units=="km" or units=="kilometers":
        factor=1
    elif units=="m" or  units=="meters":
        factor=1000
    elif units=="miles" or units=="mi":
        factor=0.621371
    else:
        raise ValueError('available units are kilometers (km), meters (m) and miles (mi)')
    lat0, lon0 = origin
    
    ellipsoids = {
            "Airy (1830)": (6377.563, 6356.257), # Ordnance Survey default
            "Bessel": (6377.397, 6356.079),
            "Clarke (1880)": (6378.249145, 6356.51486955),
            "FAI sphere": (6371, 6371), # Idealised
            "GRS-67": (6378.160, 6356.775),
            "International": (6378.388, 6356.912),
            "Krasovsky": (6378.245, 6356.863),
            "NAD27": (6378.206, 6356.584),
            "WGS66": (6378.145, 6356.758),
            "WGS72": (6378.135, 6356.751),
            "WGS84": (6378.1370, 6356.7523), # GPS default
            "IERS (2003)": (6378.1366, 6356.7519),
        }
    
    r1, r2 = ellipsoids[ellipsoid]
    lat, lon = destination
    mean_latitude = (lat0 + lat)/2
    A = (r1*r1 * math.cos(mean_latitude))**2
    B = (r2*r2 * math.sin(mean_latitude))**2
    C = (r1 * math.cos(mean_latitude))**2
    D = (r2 * math.sin(mean_latitude))**2
    radius  = np.sqrt(( A + B )/( C + D)) #radius of the earth in km
    
    dlat = math.radians(lat-lat0)
    dlon = math.radians(lon-lon0)
    a = math.sin(dlat/2) * math.sin(dlat/2) + math.cos(math.radians(lat0)) \
        * math.cos(math.radians(lat)) * math.sin(dlon/2) * math.sin(dlon/2)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
    distance=radius * c*factor

    return distance

Example Usage:

lat1=22.5
lon1=-74.3
lat2=23.8
lon2=-83.2

haversine((lat1,lon1),(lat2,lon2),units='km')
Out[7]: 919.627811821197

Previous Articles

• Python | MJO | Phase Diagram
• Python | Longitude Conversion Between Hemispheres
• Python | Atmospheric Science | Partial Correlation
• Python | Meteorological Plotting | Typhoon Precipitation
• Python | Resolving Issues with the cmaps Library
• Python | Batch Downloading NCEP2 Reanalysis Data
• Python | North Atlantic Oscillation | NAO Index | EOF