#!/usr/bin/env python
"""
datum.py
Written by Tyler Sutterley (03/2023)
Gravitational and ellipsoidal parameters
CALLING SEQUENCE
wgs84 = datum('WGS84', units='MKS')
INPUT:
ellipsoid - reference ellipsoid name
CLK66 = Clarke 1866
GRS67 = Geodetic Reference System 1967 (IAU ellipsoid)
GRS80 = Geodetic Reference System 1980
WGS72 = World Geodetic System 1972
WGS84 = World Geodetic System 1984
ATS77 = Quasi-earth centred ellipsoid for ATS77
NAD27 = North American Datum 1927
NAD83 = North American Datum 1983
INTER = International
KRASS = Krassovsky (USSR)
MAIRY = Modified Airy (Ireland 1965/1975)
HGH80 = Hughes 1980 Ellipsoid used in some NSIDC data
TOPEX = TOPEX/POSEIDON ellipsoid
EGM96 = EGM 1996 gravity model
IERS = IERS Numerical Standards (2010)
PYTHON DEPENDENCIES:
numpy: Scientific Computing Tools For Python
https://numpy.org
https://numpy.org/doc/stable/user/numpy-for-matlab-users.html
REFERENCE:
Hofmann-Wellenhof and Moritz, "Physical Geodesy" (2006)
IERS Numerical Standards
https://iers-conventions.obspm.fr/content/tn36.pdf
UPDATE HISTORY:
Updated 03/2023: include typing of variables
set ellipsoid name as constants attribute
Updated 01/2023: include main ellipsoid attributes in docstring
Forked 12/2022 from ref_ellipsoid.py
"""
from __future__ import annotations
import numpy as np
_ellipsoids = [
'CLK66',
'GRS67',
'GRS80',
'WGS72',
'WGS84',
'ATS77',
'NAD27',
'NAD83',
'INTER',
'KRASS',
'MAIRY',
'HGH80',
'TOPEX',
'EGM96',
'IERS',
]
_units = ['MKS', 'CGS']
[docs]
class datum(object):
"""
Class for gravitational and ellipsoidal parameters
Parameters
----------
ellipsoid: str, default 'WGS84'
Reference ellipsoid name
- ``'CLK66'``: Clarke 1866
- ``'GRS67'``: Geodetic Reference System 1967
- ``'GRS80'``: Geodetic Reference System 1980
- ``'HGH80'``: Hughes 1980 Ellipsoid
- ``'WGS72'``: World Geodetic System 1972
- ``'WGS84'``: World Geodetic System 1984
- ``'ATS77'``: Quasi-earth centred ellipsoid for ATS77
- ``'NAD27'``: North American Datum 1927
- ``'NAD83'``: North American Datum 1983
- ``'INTER'``: International
- ``'KRASS'``: Krassovsky (USSR)
- ``'MAIRY'``: Modified Airy (Ireland 1965/1975)
- ``'TOPEX'``: TOPEX/POSEIDON ellipsoid
- ``'EGM96'``: EGM 1996 gravity model
- ``'IERS'``: IERS Numerical Standards (2010)
units: str, default `MKS`
Output units
- ``'MKS'``: meters, kilograms, seconds
- ``'CGS'``: centimeters, grams, seconds
Attributes
----------
a_axis: float
Semi-major axis of the ellipsoid
flat: float
Flattening of the ellipsoid
omega: float
Angular velocity of the Earth
GM: float
Geocentric gravitational constant
"""
np.seterr(invalid='ignore')
def __init__(self, ellipsoid: str = 'WGS84', units: str = 'MKS'):
# set ellipsoid name and units
self.name = ellipsoid.upper()
self.units = units.upper()
# validate ellipsoid and units
assert self.name in _ellipsoids
assert self.units in _units
# set parameters for ellipsoid
if self.name in ('CLK66', 'NAD27'):
# Clarke 1866
self.a_axis = 6378206.4 # [m] semimajor axis of the ellipsoid
self.flat = 1.0 / 294.9786982 # flattening of the ellipsoid
elif self.name in ('GRS80', 'NAD83'):
# Geodetic Reference System 1980
# North American Datum 1983
self.a_axis = 6378135.0 # [m] semimajor axis of the ellipsoid
self.flat = 1.0 / 298.26 # flattening of the ellipsoid
self.GM = 3.986005e14 # [m^3/s^2] Geocentric Gravitational Constant
elif self.name == 'GRS67':
# Geodetic Reference System 1967
# International Astronomical Union (IAU ellipsoid)
self.a_axis = 6378160.0 # [m] semimajor axis of the ellipsoid
self.flat = 1.0 / 298.247167427 # flattening of the ellipsoid
self.GM = 3.98603e14 # [m^3/s^2] Geocentric Gravitational Constant
self.omega = (
7292115.1467e-11 # angular velocity of the Earth [rad/s]
)
elif self.name == 'WGS72':
# World Geodetic System 1972
self.a_axis = 6378135.0 # [m] semimajor axis of the ellipsoid
self.flat = 1.0 / 298.26 # flattening of the ellipsoid
elif self.name == 'WGS84':
# World Geodetic System 1984
self.a_axis = 6378137.0 # [m] semimajor axis of the ellipsoid
self.flat = 1.0 / 298.257223563 # flattening of the ellipsoid
elif self.name == 'ATS77':
# Quasi-earth centred ellipsoid for ATS77
self.a_axis = 6378135.0 # [m] semimajor axis of the ellipsoid
self.flat = 1.0 / 298.257 # flattening of the ellipsoid
elif self.name == 'KRASS':
# Krassovsky (USSR)
self.a_axis = 6378245.0 # [m] semimajor axis of the ellipsoid
self.flat = 1.0 / 298.3 # flattening of the ellipsoid
elif self.name == 'INTER':
# International
self.a_axis = 6378388.0 # [m] semimajor axis of the ellipsoid
self.flat = 1 / 297.0 # flattening of the ellipsoid
elif self.name == 'MAIRY':
# Modified Airy (Ireland 1965/1975)
self.a_axis = 6377340.189 # [m] semimajor axis of the ellipsoid
self.flat = 1 / 299.3249646 # flattening of the ellipsoid
elif self.name == 'HGH80':
# Hughes 1980 Ellipsoid used in some NSIDC data
self.a_axis = 6378273.0 # [m] semimajor axis of the ellipsoid
self.flat = 1.0 / 298.279411123064 # flattening of the ellipsoid
elif self.name == 'TOPEX':
# TOPEX/POSEIDON ellipsoid
self.a_axis = 6378136.3 # [m] semimajor axis of the ellipsoid
self.flat = 1.0 / 298.257 # flattening of the ellipsoid
self.GM = 3.986004415e14 # [m^3/s^2]
elif self.name == 'EGM96':
# EGM 1996 gravity model
self.a_axis = 6378136.3 # [m] semimajor axis of the ellipsoid
self.flat = 1.0 / 298.256415099 # flattening of the ellipsoid
self.GM = 3.986004415e14 # [m^3/s^2]
elif self.name == 'IERS':
# IERS Numerical Standards
self.a_axis = 6378136.6 # [m] semimajor axis of the ellipsoid
self.flat = 1.0 / 298.25642 # flattening of the ellipsoid
# set default parameters if not listed as part of ellipsoid
if self.name not in ('GRS80', 'GRS67', 'NAD83', 'TOPEX', 'EGM96'):
# for ellipsoids not listing the Geocentric Gravitational Constant
self.GM = 3.986004418e14 # [m^3/s^2]
if self.name not in ('GRS67'):
# for ellipsoids not listing the angular velocity of the Earth
self.omega = 7292115e-11 # [rad/s]
# universal gravitational constant [N*m^2/kg^2]
self.G = 6.67430e-11
# convert units to CGS
if self.units == 'CGS':
self.a_axis *= 100.0
self.GM *= 1e6
self.G *= 1000.0 # [dyn*cm^2/g^2]
# mean radius of the Earth having the same volume
# (4pi/3)R^3 = (4pi/3)(a^2)b = (4pi/3)(a^3)(1 - f)
@property
def rad_e(self) -> float:
"""Average radius of the Earth with same volume as ellipsoid"""
return self.a_axis * (1.0 - self.flat) ** (1.0 / 3.0)
# semiminor axis of the ellipsoid
@property
def b_axis(self) -> float:
"""Semi-minor axis of the ellipsoid"""
return (1.0 - self.flat) * self.a_axis
# Ratio between ellipsoidal axes
@property
def ratio(self) -> float:
"""Ratio between ellipsoidal axes"""
return 1.0 - self.flat
# Polar radius of curvature
@property
def rad_p(self) -> float:
"""Polar radius of curvature"""
return self.a_axis / (1.0 - self.flat)
# Linear eccentricity
@property
def ecc(self) -> float:
"""Linear eccentricity"""
return np.sqrt((2.0 * self.flat - self.flat**2) * self.a_axis**2)
# first numerical eccentricity
@property
def ecc1(self) -> float:
"""First numerical eccentricity"""
return self.ecc / self.a_axis
# second numerical eccentricity
@property
def ecc2(self) -> float:
"""Second numerical eccentricity"""
return self.ecc / self.b_axis
# m parameter [omega^2*a^2*b/(GM)]
# p. 70, Eqn.(2-137)
@property
def m(self) -> float:
"""m Parameter"""
return self.omega**2 * ((1 - self.flat) * self.a_axis**3) / self.GM
# q
# p. 67, Eqn.(2-113)
@property
def q(self) -> float:
"""q Parameter"""
return 0.5 * (
(1.0 + 3.0 / (self.ecc2**2)) * np.arctan(self.ecc2)
- 3.0 / self.ecc2
)
# q_0
# p. 67, Eqn.(2-113)
@property
def q0(self) -> float:
r"""q\ :sub:`0`\ Parameter"""
return (
3
* (1.0 + 1.0 / (self.ecc2**2))
* (1.0 - 1.0 / self.ecc2 * np.arctan(self.ecc2))
- 1.0
)
# J_2 p. 75 Eqn.(2-167), p. 76 Eqn.(2-172)
@property
def J2(self) -> float:
"""Oblateness coefficient"""
return (
(self.ecc1**2)
* (1.0 - 2.0 * self.m * self.ecc2 / (15.0 * self.q))
/ 3.0
)
# Normalized C20 harmonic
# p. 60, Eqn.(2-80)
@property
def C20(self) -> float:
r"""Normalized C\ :sub:`20`\ harmonic"""
return -self.J2 / np.sqrt(5.0)
# Normal gravity at the equator
# p. 71, Eqn.(2-141)
@property
def ga(self) -> float:
"""Normal gravity at the equator"""
return (
self.GM
/ (self.a_axis * self.b_axis)
* (1.0 - self.m - self.m * self.ecc2 * self.q0 / (6.0 * self.q))
)
# Normal gravity at the pole
# p. 71, Eqn.(2-142)
@property
def gb(self) -> float:
"""Normal gravity at the pole"""
return (
self.GM
/ (self.a_axis**2.0)
* (1.0 + self.m * self.ecc2 * self.q0 / (3.0 * self.q))
)
# ratio between gravity at pole versus gravity at equator
@property
def dk(self) -> float:
"""Ratio between gravity at pole versus gravity at equator"""
return self.b_axis * self.gb / (self.a_axis * self.ga) - 1.0
# Normal potential at the ellipsoid
# p. 68, Eqn.(2-123)
@property
def U0(self) -> float:
"""Normal potential at the ellipsoid"""
return (
self.GM / self.ecc * np.arctan(self.ecc2)
+ (1.0 / 3.0) * self.omega**2 * self.a_axis**2
)
# Surface area of the reference ellipsoid
@property
def area(self) -> float:
"""Surface area of the ellipsoid"""
return (
np.pi
* self.a_axis**2.0
* (
2.0
+ ((1.0 - self.ecc1**2) / self.ecc1)
* np.log((1.0 + self.ecc1) / (1.0 - self.ecc1))
)
)
# Volume of the reference ellipsoid
@property
def volume(self) -> float:
"""Volume of the ellipsoid"""
return (
(4.0 * np.pi / 3.0)
* (self.a_axis**3.0)
* (1.0 - self.ecc1**2.0) ** 0.5
)
# Average density
@property
def rho_e(self) -> float:
"""Average density"""
return self.GM / (self.G * self.volume)