"""
This Python 3.11 code implements ellipsoids of revolution.
- Alexander Raichev (AR), 2012-01-26: Refactored code from release 0.3.
NOTE:
All lengths are measured in meters and all angles are measured in radians
unless indicated otherwise.
By 'ellipsoid' throughout, I mean an ellipsoid of revolution and *not* a general (triaxial) ellipsoid.
Points lying on an ellipsoid are given in geodetic (longitude, latitude) coordinates.
"""
# *****************************************************************************
# Copyright (C) 2013 Alexander Raichev <alex.raichev@gmail.com>
#
# Distributed under the terms of the GNU Lesser General Public License (LGPL)
# http://www.gnu.org/licenses/
# *****************************************************************************
from random import uniform
import pyproj
from numpy import arcsin, arctanh, cos, deg2rad, pi, rad2deg, sin, sqrt
from rhealpixdggs.utils import auth_lat, auth_rad, my_round
# Parameters of some common ellipsoids.
WGS84_A = pyproj.get_ellps_map()["WGS84"]["a"] # 6378137.0
WGS84_F = 1 / pyproj.get_ellps_map()["WGS84"]["rf"] # 298.257223563
WGS84_B = WGS84_A * (1 - WGS84_F)
WGS84_E = sqrt(WGS84_F * (2 - WGS84_F))
WGS84_R_A = sqrt(WGS84_A**2 / 2 + WGS84_B**2 / 2 * (arctanh(WGS84_E) / WGS84_E))
R_EM = pyproj.get_ellps_map()["sphere"]["a"] # 6371000 (Earth's mean radius)
[docs]
class Ellipsoid:
"""
Represents an ellipsoid of revolution (possibly a sphere) with a
geodetic longitude-latitude coordinate frame.
INSTANCE ATTRIBUTES:
- `sphere` - True if the ellipsoid is a sphere, and False otherwise.
- `R` - The radius of the ellipsoid in meters, implying that it is a
sphere.
- `a` - Major radius of the ellipsoid in meters.
- `b` - Minor radius of the ellipsoid in meters.
- `e` - Eccentricity of the ellipsoid.
- `f` - Flattening of the ellipsoid.
- `R_A` - Authalic radius of the ellipsoid in meters.
- `lon_0` - Central meridian. Projections and grids built on the
ellipsoid are rotated about the polar axis so that this meridian
becomes the planar x = 0 line.
- `lat_0` - Latitude of origin. Must be 0. Shifting latitude is not a
rotation of an ellipsoid of revolution: a grid built on a shifted
latitude is neither equal-area nor geographically coherent near the
poles (issue #93). The attribute is kept so that the names match
PROJ's; reorient a grid with `lon_0` and the DGGS's `north_square`
and `south_square` instead.
- `radians` - If True, use angles measured in radians for all calculations.
Use degrees otherwise.
- `phi_0` - The latitude separating the equatorial region and
the north polar region in the context of the (r)HEALPix projection.
Except for phi_0, these attribute names match the names of the
`PROJ ellipsoid parameters <https://proj.org/en/9.4/usage/ellipsoids.html#ellipsoid-spherification-parameters>`
"""
def __init__(
self,
R: float | None = None,
a: float = WGS84_A,
b: float | None = None,
e: float | None = None,
f: float = WGS84_F,
lon_0: float = 0,
lat_0: float = 0,
radians: bool = False,
) -> None:
if lat_0 != 0:
raise ValueError(
f"lat_0 must be 0, got {lat_0!r}: shifting latitude is not a "
"rotation of the ellipsoid, so a grid built on it is neither "
"equal-area nor geographically coherent near the poles. "
"Reorient the grid with lon_0 and the DGGS's "
"north_square/south_square instead."
)
self.lon_0 = lon_0
self.lat_0 = lat_0
self.radians = radians
if R is not None:
# The ellipsoid is a sphere.
# Override the other geometric parameters.
self.sphere = True
self.R = R
self.a = R
self.b = R
self.e: float = 0
self.f: float = 0
self.R_A = R
else:
self.sphere = False
self.a = a
if b is not None:
# Derive the other geometric parameters from a and b.
self.b = b
self.e = sqrt(1 - (b / a) ** 2)
self.f = (a - b) / a
elif e is not None:
# Derive the other geometric parameters from a and e.
self.e = e
self.b = a * sqrt(1 - e**2)
self.f = 1 - sqrt(1 - e**2)
else:
self.f = f
self.b = self.a * (1 - f)
self.e = sqrt(f * (2 - f))
self.R_A = auth_rad(self.a, self.e)
self.phi_0 = auth_lat(arcsin(2.0 / 3), e=self.e, radians=True, inverse=True)
if not self.radians:
# Convert to degrees.
self.phi_0 = rad2deg(self.phi_0)
def __str__(self) -> str:
result = ["ellipsoid:"]
# result.append('lengths measured in meters')
for k, v in sorted(self.__dict__.items()):
if k == "phi_0":
continue
if k in {"sphere", "radians"}:
result.append(" " + k + " = " + str(v))
else:
result.append(" " + k + " = " + str(my_round(v, 15)))
return "\n".join(result)
def __eq__(self, other: object) -> bool:
if not isinstance(other, Ellipsoid):
return NotImplemented
return self.a == other.a and self.b == other.b
def __ne__(self, other: object) -> bool:
"""
The inequality relation on cells.
Since Python 3.3 doesn't automatically create reverse relations
from given ones, i must define this seemingly redundant relation.
"""
return not self.__eq__(other)
[docs]
def pi(self) -> float:
"""
Return pi if `self.radians` = True and 180 otherwise.
"""
if self.radians:
return pi
else:
return 180.0
[docs]
def random_point(
self,
lam_min: float | None = None,
lam_max: float | None = None,
phi_min: float | None = None,
phi_max: float | None = None,
) -> tuple[float, float]:
"""
Return a point (given in geodetic coordinates) sampled uniformly at
random from the section of this ellipsoid with longitude in the range
`lam_min <= lam < lam_max` and latitude in the range
`phi_min <= phi < phi_max`.
But avoid the poles.
EXAMPLES::
>>> E = UNIT_SPHERE
>>> print(E.random_point()) # doctest: +SKIP
(-1.0999574573422948, 0.21029104897701129)
"""
PI = self.pi()
if lam_min is None:
lam_min = -PI
if lam_max is None:
lam_max = PI
if phi_min is None:
phi_min = -PI / 2
if phi_max is None:
phi_max = PI / 2
if not self.radians:
# Convert to radians.
lam_min, lam_max, phi_min, phi_max = deg2rad(
[lam_min, lam_max, phi_min, phi_max]
)
# Pick a longitude.
while True:
u = uniform(0, 1)
lam = (lam_max - lam_min) * u + lam_min
# Don't include lam_max.
if lam < lam_max:
# Success.
break
# Pick a latitude.
delta = pi / 360
while True:
v = uniform(0, 1)
if self.sphere:
phi = arcsin((sin(phi_max) - sin(phi_min)) * v + sin(phi_min))
else:
# Sample from the authalic sphere.
# The map from the ellipsoid to the authalic sphere is
# an equiareal diffeomorphism.
# So a uniform distribution on the authalic sphere gives
# rise to a uniform distribution on the ellipsoid.
beta0 = auth_lat(phi_min, e=self.e, radians=True)
beta1 = auth_lat(phi_max, e=self.e, radians=True)
beta = arcsin((sin(beta1) - sin(beta0)) * v + sin(beta0))
phi = auth_lat(beta, e=self.e, radians=True, inverse=True)
# Avoid the poles.
if abs(phi) <= pi / 2 - delta:
# Success.
break
if not self.radians:
# Convert back to degrees.
lam, phi = rad2deg([lam, phi])
return lam, phi
[docs]
def xyz(self, lam: float, phi: float) -> tuple[float, float, float]:
"""
Given a point on this ellipsoid with longitude-latitude coordinates
`(lam, phi)`, return the point's 3D rectangular coordinates.
EXAMPLES::
>>> E = UNIT_SPHERE
>>> print(tuple(x.tolist() for x in my_round(E.xyz(0, 45), 15)))
(0.707106781186548, 0.0, 0.707106781186548)
"""
a = self.a
e = self.e
if not self.radians:
lam, phi = deg2rad([lam, phi])
# Equals a iff e = 0 (sphere):
N = a / sqrt(1 - e**2 * sin(phi) ** 2)
return (
N * cos(lam) * cos(phi),
N * sin(lam) * cos(phi),
N * (1 - e**2) * sin(phi),
)
# Define some common ellipsoids.
WGS84_ELLIPSOID = Ellipsoid(a=WGS84_A, f=WGS84_F)
WGS84_ELLIPSOID_RADIANS = Ellipsoid(a=WGS84_A, f=WGS84_F, radians=True)
WGS84_ASPHERE = Ellipsoid(R=WGS84_R_A)
WGS84_ASPHERE_RADIANS = Ellipsoid(R=WGS84_R_A, radians=True)
EMR_SPHERE = Ellipsoid(R=R_EM)
EMR_SPHERE_RADIANS = Ellipsoid(R=R_EM, radians=True)
UNIT_SPHERE = Ellipsoid(R=1)
UNIT_SPHERE_RADIANS = Ellipsoid(R=1, radians=True)