# -*- coding: utf-8 -*-
"""Functions for handling geographical coordinate systems
and reference ellipsoids.
Unless otherwise stated functions are ported from atmlab-2-3-181.
"""
import numpy as np
from typhon import constants
__all__ = [
'ellipsoidmodels',
'ellipsoid_r_geodetic',
'ellipsoid_r_geocentric',
'ellipsoid2d',
'ellipsoidcurvradius',
'sind',
'cosd',
'tand',
'asind',
'inrange',
'cart2geocentric',
'geocentric2cart',
'cart2geodetic',
'geodetic2cart',
'geodetic2geocentric',
'geocentric2geodetic',
'great_circle_distance',
'geographic_mean',
'cartposlos2geocentric',
'geocentricposlos2cart',
'sphere_plane_intersection',
'tunnel_distance',
]
[docs]def sind(x):
"""Sine of argument in degrees."""
return np.sin(np.deg2rad(x))
[docs]def cosd(x):
"""Cosine of argument in degrees."""
return np.cos(np.deg2rad(x))
[docs]def tand(x):
"""Tangent of argument in degrees."""
return np.tan(np.deg2rad(x))
[docs]def asind(x):
"""Inverse sine in degrees."""
return np.arcsin(np.deg2rad(x))
[docs]def inrange(x, minx, maxx, exclude='none', text=None):
"""Test if x is within given bounds.
Parameters:
x: Variable to test.
minx: Lower boundary.
maxx: Upper boundary.
exclude (str): Exclude boundaries. Possible values are:
'none', 'lower', 'upper' and 'both'
text (str): Addiitional warning text.
Raises:
Exception: If value is out of bounds.
"""
compare = {'none': (np.greater_equal, np.less_equal),
'lower': (np.greater, np.less_equal),
'upper': (np.greater_equal, np.less),
'both': (np.greater, np.less),
}
greater, less = compare[exclude]
if less(x, minx) or greater(x, maxx):
if text is None:
raise Exception('Range out of bound [{}, {}]'.format(minx, maxx))
else:
raise Exception(
'Range out of bound [{}, {}]: {}'.format(minx, maxx, text)
)
[docs]class ellipsoidmodels():
"""Provide data for different reference ellipsoids.
The following models are covered:
* SphericalEarth (radius set as constants.earth_radius)
* WGS84
* SphericalVenus (radius same as used in ARTS)
* SphericalMars (radius same as used in ARTS)
* EllipsoidMars
* SphericalJupiter (radius same as used in ARTS)
Examples:
>>> e = ellipsoidmodels()
.. Ported from atmlab. Original author: Patrick Eriksson
"""
[docs] def __init__(self):
self._data = {
"SphericalEarth": (constants.earth_radius, 0),
"WGS84": (6378137, 0.0818191908426),
"SphericalVenus": (6051800.0, 0),
"SphericalMars": (3389500.0, 0),
"EllipsoidMars": (3396190.0, 0.1083),
"SphericalJupiter": (69911000.0, 0),
}
def __getitem__(self, model):
return self.get(model)
[docs] def get(self, model='WGS84'):
"""Return data for different reference ellipsoids.
Parameters:
model (str): Model ellipsoid.
Returns:
tuple: Equatorial radius (r), eccentricity (e)
Examples:
>>> e['WGS84']
(6378137, 0.0818191908426)
>>> e.get('WGS84')
(6378137, 0.0818191908426)
>>> ellipsoidmodel()['WGS84']
(6378137, 0.0818191908426)
"""
if model in self._data:
return self._data.__getitem__(model)
else:
raise Exception('Unknown ellipsoid model "{}".'.format(model))
@property
def models(self):
"""List of available models.
Examples:
>>> e.models
['SphericalVenus',
'SphericalMars',
'WGS84',
'SphericalJupiter',
'EllipsoidMars',
'SphericalEarth']
"""
return list(self._data.keys())
[docs]def ellipsoid_r_geocentric(ellipsoid, lat):
"""Geocentric radius of a reference ellipsoid.
Gives the distance from the Earth's centre and the reference ellipsoid
as a function of geo\ **centric** latitude.
Note:
To obtain the radii for **geodetic** latitude,
use :func:`ellipsoid_r_geodetic`.
Parameters:
ellipsoid (tuple): Model ellipsoid as returned
by :class:`ellipsoidmodels`.
lat (float or ndarray): Geocentric latitudes.
Returns:
float or ndarray: Radii.
.. Ported from atmlab. Original author: Patrick Eriksson
"""
errtext = 'Invalid excentricity value in ellipsoid model.'
inrange(ellipsoid[1], 0, 1, exclude='upper', text=errtext)
if ellipsoid[1] == 0:
r = np.ones(np.shape(lat)) * ellipsoid[0]
else:
c = 1 - ellipsoid[1]**2
b = ellipsoid[0] * np.sqrt(c)
r = b / np.sqrt(c * cosd(lat)**2 + sind(lat)**2)
return r
[docs]def ellipsoid_r_geodetic(ellipsoid, lat):
"""Geodetic radius of a reference ellipsoid.
Gives the distance from the Earth's centre and the reference ellipsoid
as a function of geo\ **detic** latitude.
Note:
To obtain the radii for **geocentric** latitude,
use :func:`ellipsoid_r_geocentric`.
Parameters:
ellipsoid (tuple): Model ellipsoid as returned
by :class:`ellipsoidmodels`.
lat (float or ndarray): Geodetic latitudes.
Returns:
float or ndarray: Radii.
.. Ported from atmlab. Original author: Patrick Eriksson
"""
errtext = 'Invalid excentricity value in ellipsoid model.'
inrange(ellipsoid[1], 0, 1, exclude='upper', text=errtext)
if ellipsoid[1] == 0:
r = np.ones(np.shape(lat)) * ellipsoid[0]
else:
e2 = ellipsoid[1]**2
sin2 = sind(lat)**2
r = (ellipsoid[0] * np.sqrt((1 - e2)**2 * sin2 +
cosd(lat) ** 2) / np.sqrt(1 - e2 * sin2))
return r
[docs]def ellipsoid2d(ellipsoid, orbitinc):
"""Approximate ellipsoid for 2D calculations.
Determines an approximate reference ellipsoid following an orbit track.
The new ellipsoid is determined simply, by determining the radius at the
maximum latitude and from this value calculate a new eccentricity.
The orbit is specified by giving the orbit inclination, that is
normally a value around 100 deg for polar sun-synchronous orbits.
Parameters:
ellipsoid (tuple): Model ellipsoid as returned
by :class:`ellipsoidmodels`.
orbitinc (float): Orbit inclination.
Returns:
tuple: Modified ellipsoid vector.
.. Ported from atmlab. Original author: Patrick Eriksson
"""
errtext = 'Invalid excentricity value in ellipsoid model.'
inrange(ellipsoid[1], 0, 1, exclude='upper', text=errtext)
inrange(orbitinc, 0, 180,
exclude='both',
text='Invalid orbit inclination.')
rp = ellipsoid_r_geocentric(ellipsoid, orbitinc)
return ellipsoid[0], np.sqrt(1 - (rp / ellipsoid[0])**2)
[docs]def ellipsoidcurvradius(ellipsoid, lat_gd, azimuth):
"""Sets ellispoid to local curvature radius
Calculates the curvature radius for the given latitude and azimuth
angle, and uses this to set a spherical reference ellipsoid
suitable for 1D calculations. The curvature radius is a better
local approximation than using the local ellipsoid radius.
For exact result the *geodetic* latitude shall be used.
Parameters:
lat_gd: Geodetic latitude.
azimuth: Azimuthal angle (angle from NS plane).
If given curvature radii are returned, see above.
Returns:
tuple: Modified ellipsoid.
.. Ported from atmlab. Original author: Patrick Eriksson
"""
errtext = 'Invalid excentricity value in ellipsoid model.'
inrange(ellipsoid[1], 0, 1, exclude='upper', text=errtext)
aterm = 1 - ellipsoid[1]**2 * sind(lat_gd)**2
rn = 1 / np.sqrt(aterm)
rm = (1 - ellipsoid[1]**2) * (rn / aterm)
e0 = (ellipsoid[0] / (cosd(azimuth)**2.0 / rm + sind(azimuth)**2.0 / rn))
e1 = 0
return e0, e1
[docs]def cart2geocentric(x, y, z, lat0=None, lon0=None, za0=None, aa0=None):
"""Convert cartesian position to spherical coordinates.
The geocentric Cartesian coordinate system is fixed with respect to the
Earth, with its origin at the center of the ellipsoid and its X-, Y-,
and Z-axes intersecting the surface at the following points:
* X-axis: Equator at the Prime Meridian (0°, 0°)
* Y-axis: Equator at 90-degrees East (0°, 90°)
* Z-axis: North Pole (90°, 0°)
A common synonym is Earth-Centered, Earth-Fixed coordinates, or ECEF.
If the optional arguments are given, it is ensured that latitude and
longitude are kept constant for zenith or nadir cases, and the longitude
for N-S cases. The optional input shall be interpreted as the [x,y,z]
is obtained by moving from [lat0,lon0] in the direction of [za0,aa0].
Parameters:
x: Coordinate in x dimension.
y: Coordinate in y dimension.
z: Coordinate in z dimension.
lat0: Original latitude.
lon0: Original longitude.
za0: Orignal zenith angle.
aa0: Orignal azimuth angle.
Returns:
tuple: Radius, Latitude, Longitude
.. Ported from atmlab. Original author: Bengt Rydberg
"""
r = np.sqrt(x**2 + y**2 + z**2)
if np.any(r == 0):
raise Exception("This set of functions does not handle r = 0.")
lat = np.rad2deg(np.arcsin(z / r))
lon = np.rad2deg(np.arctan2(y, x))
if all(x is not None for x in [lat0, lon0, za0, aa0]):
for i in range(np.size(r)):
if za0[i] < 1e-06 or za0[i] > 180 - 1e-06:
lat[i] = lat0[i]
lon[i] = lon0[i]
if (abs(lat0[i]) < 90 - 1e-08 and
(abs(aa0[i]) < 1e-06 or abs(aa0[i] - 180) < 1e-06)):
if abs(lon[i] - lon0[i]) < 1:
lon[i] = lon0[i]
else:
if lon0[i] > 0:
lon[i] = lon0[i] - 180
else:
lon[i] = lon0[i] + 180
return r, lat, lon
[docs]def geocentric2cart(r, lat, lon):
"""Convert from spherical coordinate to a cartesian position.
See :func:`cart2geocentric` for a defintion of the geocentric
coordinate system.
Parameters:
r: Radius.
lat: Latitude in degree.
lon Longitude in degree.
Returns:
tuple: Coordinate in x, y, z dimension.
.. Ported from atmlab. Original author: Bengt Rydberg
"""
if np.any(r == 0):
raise Exception("This set of functions does not handle r = 0.")
latrad = np.deg2rad(lat)
lonrad = np.deg2rad(lon)
x = r * np.cos(latrad)
y = x * np.sin(lonrad)
x = x * np.cos(lonrad)
z = r * np.sin(latrad)
return x, y, z
[docs]def cart2geodetic(x, y, z, ellipsoid=None):
"""Convert from cartesian to geodetic coordinates.
The geodetic coordinates refer to the reference ellipsoid
specified by input ellipsoid.
See module docstring for a defintion of the geocentric coordinate system.
Parameters:
x: Coordinates in x dimension.
y: Coordinates in y dimension.
z: Coordinates in z dimension.
ellipsoid: A tuple with the form (semimajor axis, eccentricity).
Default is 'WGS84' from :class:`ellipsoidmodels`.
Returns:
tuple: Geodetic height, latitude and longitude
.. Ported from atmlab. Original author: Bengt Rydberg
"""
if ellipsoid is None:
ellipsoid = ellipsoidmodels()['WGS84']
errtext = 'Invalid excentricity value in ellipsoid model.'
inrange(ellipsoid[1], 0, 1, exclude='upper', text=errtext)
lon = np.rad2deg(np.arctan2(y, x))
B0 = np.arctan2(z, np.hypot(x, y))
B = np.ones(B0.shape)
e2 = ellipsoid[1]**2
if e2 == 0.0:
h, lat, lon = cart2geocentric(x, y, z)
h -= ellipsoid[0]
else:
while (np.any(np.abs(B - B0) > 1e-10)):
N = ellipsoid[0] / np.sqrt(1 - e2 * np.sin(B0)**2)
h = np.hypot(x, y) / np.cos(B0) - N
B = B0.copy()
B0 = np.arctan(z/np.hypot(x, y) * ((1-e2*N/(N+h))**(-1)))
lat = np.rad2deg(B)
return h, lat, lon
[docs]def geodetic2cart(h, lat, lon, ellipsoid=None):
"""Convert from geodetic to geocentric cartesian coordinates.
The geodetic coordinates refer to the reference ellipsoid
specified by input ellipsoid.
See module docstring for a defintion of the geocentric coordinate system.
Parameters:
h: Geodetic height (height above the reference ellipsoid).
lat: Geodetic latitude.
lon: Geodetic longitude.
ellipsoid: A tuple with the form (semimajor axis, eccentricity).
Default is 'WGS84' from :class:`ellipsoidmodels`.
Returns:
tuple: x, y, z coordinates.
.. Ported from atmlab. Original author: Bengt Rydberg
"""
if ellipsoid is None:
ellipsoid = ellipsoidmodels()['WGS84']
errtext = 'Invalid excentricity value in ellipsoid model.'
inrange(ellipsoid[1], 0, 1, exclude='upper', text=errtext)
a = ellipsoid[0]
e2 = ellipsoid[1] ** 2
N = a / np.sqrt(1 - e2 * sind(lat)**2)
x = (N + h) * (cosd(lat)) * (cosd(lon))
y = (N + h) * (cosd(lat)) * (sind(lon))
# np.ones(np.shape(lon)): Ensure equal shape of x, y, z.
z = (N * (1 - e2) + h) * (sind(lat)) * np.ones(np.shape(lon))
return x, y, z
[docs]def geodetic2geocentric(h, lat, lon, ellipsoid=None, **kwargs):
"""Convert from geodetic to geocentric coordinates.
The geodetic coordinates refer to the reference ellipsoid
specified by input ellipsoid.
See module docstring for a defintion of the geocentric coordinate system.
Parameters:
h: Geodetic height (height above the reference ellipsoid).
lat: Geodetic latitude.
lon: Geodetic longitude.
kwargs: Additional keyword arguments for :func:`cart2geocentric`.
ellipsoid: A tuple with the form (semimajor axis, eccentricity).
Default is 'WGS84' from :class:`ellipsoidmodels`.
Returns:
tuple: Radius, geocentric latiude, geocentric longitude
.. Ported from atmlab. Original author: Bengt Rydberg
"""
if ellipsoid is None:
ellipsoid = ellipsoidmodels()['WGS84']
errtext = 'Invalid excentricity value in ellipsoid model.'
inrange(ellipsoid[1], 0, 1, exclude='upper', text=errtext)
cart = geodetic2cart(h, lat, lon, ellipsoid)
return cart2geocentric(*cart, **kwargs)
[docs]def geocentric2geodetic(r, lat, lon, ellipsoid=None):
"""Convert from geocentric to geodetic coordinates.
The geodetic coordinates refer to the reference ellipsoid
specified by input ellipsoid.
See module docstring for a defintion of the geocentric coordinate system.
Returns:
tuple: Geodetic height, latitude and longitude
Parameters:
r: Radius:
lat: Geocentric latitude.
lon: Geocentric longitude.
ellipsoid: A tuple with the form (semimajor axis, eccentricity).
Default is 'WGS84' from :class:`ellipsoidmodels`.
.. Ported from atmlab. Original author: Bengt Rydberg
"""
if ellipsoid is None:
ellipsoid = ellipsoidmodels()['WGS84']
errtext = 'Invalid excentricity value in ellipsoid model.'
inrange(ellipsoid[1], 0, 1, exclude='upper', text=errtext)
cart = geocentric2cart(r, lat, lon)
return cart2geodetic(*cart, ellipsoid)
[docs]def great_circle_distance(lat1, lon1, lat2, lon2, r=None):
r"""Calculate the distance between two geographical positions
This is the 'As-the-crow-flies' distance between two points, specified by
their latitude and longitude. This uses the so-called haversine formula,
defined by:
.. math::
d = 2 \arcsin \left( \sqrt{
\sin \left(\frac{\varphi_2-\varphi_1}{2} \right)^2
+ \cos(\varphi_1) \cdot \cos(\varphi_1)
\cdot \sin \left(\frac{\lambda_2-\lambda_1}{2} \right)^2} \right)
Args:
lat1: Latitude of position 1.
lon1: Longitude of position 1.
lat2: Latitude of position 2.
lon2: Longitude of position 2.
r (float): The radius (common for both points).
Returns:
If the optional argument *r* is given, the distance in m is returned.
Otherwise the angular distance in degrees is returned.
.. Taken from https://stackoverflow.com/a/29546836/9144990
"""
lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
dlon = lon2 - lon1
dlat = lat2 - lat1
a = np.sin(dlat / 2.0) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(
dlon / 2.0) ** 2
c = 2 * np.arcsin(np.sqrt(a))
if r is None:
return np.rad2deg(c)
else:
return r * c
[docs]def geographic_mean(lat, lon, h=0, ellipsoid=None):
"""Calculate mean position for set of coordinates.
Parameters:
lat (float or ndarray): Latitudes in degrees.
lon (float or ndarray): Longitudes in degrees.
h (float or ndarray): Optiional altitude for
each coordinate (default is).
ellipsoid: A tuple with the form (semimajor axis, eccentricity).
Default is 'WGS84' from :class:`ellipsoidmodels`.
Returns:
tuple: Mean latitudes and longitudes in degrees.
"""
if ellipsoid is None:
ellipsoid = ellipsoidmodels()['WGS84']
x, y, z = geodetic2cart(
h,
lat,
lon,
ellipsoid=ellipsoid)
mh, mlat, mlon = cart2geodetic(
np.mean(x),
np.mean(y),
np.mean(z),
ellipsoid=ellipsoid)
return mlat, mlon
[docs]def cartposlos2geocentric(x, y, z, dx, dy, dz, ppc=None,
lat0=None, lon0=None, za0=None, aa0=None):
"""Convert cartesian POS/LOS to spherical coordinates.
Position is given as (x,y,z), while line-of-sight is given as (dx,dy,dz).
The corresponding quantities in polar coordinates are (r,lat,lon)
and (za,aa), respectively.
See *Contents* for defintion of coordinate systems.
If the optional arguments are given, it is ensured that latitude and
longitude are kept constant for zenith or nadir cases, and the longitude
and azimuth angle for N-S cases. The optional input shall be interpreted
as the [x,y,z] is obtained by moving from [r0,lat0,lon0] in the direction
of [za0,aa0].
This version is different from the atmlab version by normalizing the los-
vector and demanding all or nothing for the optional arguments to work.
Parameters:
x: Coordinate in x dimension.
y: Coordinate in y dimension.
z: Coordinate in z dimension.
dx: LOS component in x dimension.
dy: LOS component in y dimension.
dz: LOS component in z dimension.
ppc: Propagation path constant = r0*sin(za0).
lat0: Original latitude.
lon0: Original longitude.
za0: Orignal zenith angle.
aa0: Orignal azimuth angle.
Returns:
tuple(ndarray): Radius, Latitude, Longitude,
Zenith angle, Azimuth angle
.. Ported from atmlab. Original author: Bengt Rydberg
"""
# Here be dragons!
# Broadcast all input variables to the same shape. Atleast (1)
if(ppc is not None and za0 is not None and lat0 is not None and
aa0 is not None and lon0 is not None):
x, y, z, dx, dy, dz, ppc, lat0, lon0, za0, aa0 = _broadcast(
x, y, z, dx, dy, dz, ppc, lat0, lon0, za0, aa0)
elif ppc is not None:
x, y, z, dx, dy, dz, ppc = _broadcast(x, y, z, dx, dy, dz, ppc)
else:
x, y, z, dx, dy, dz = _broadcast(x, y, z, dx, dy, dz)
r, lat, lon = cart2geocentric(x, y, z, lat0, lon0, za0, aa0)
# Renormalize for length of the variables (not in atmlab)
norm_r = np.sqrt(dx**2 + dy**2 + dz**2)
dx = dx / norm_r
dy = dy / norm_r
dz = dz / norm_r
coslat = np.cos(np.deg2rad(lat))
sinlat = np.sin(np.deg2rad(lat))
coslon = np.cos(np.deg2rad(lon))
sinlon = np.sin(np.deg2rad(lon))
dr = np.clip(coslat * coslon * dx + sinlat * dz + coslat * sinlon * dy,
-1., 1.)
# Get LOS angle
if ppc is None:
za = np.rad2deg(np.arccos(dr))
else:
za = np.rad2deg(np.arcsin(ppc / r))
aa = np.zeros(za.shape)
# Fix zenith and azimuth angle with optional input only when all exists
if(za0 is not None and lat0 is not None and
aa0 is not None and lon0 is not None):
# Determine the type for zenith
noz = np.logical_or(za0 < 1e-06, za0 > 180 - 1e-06)
nan = np.isnan(za)
pre = np.logical_and(~noz, nan)
# Either set or do not
za[noz] = za0[noz]
za[pre] = 90.
# NB: removed check for dr<0 since by putting dr==1 is more sensible
# Determine the type for azimuth
cir1 = abs(aa0) < 1e-06
cir2 = np.logical_or(cir1, abs(aa0 - 180) < 1e-06)
same = np.equal(lon, lon0)
circ = np.logical_and(cir2, same)
left = np.logical_and(cir1, ~same)
right = np.logical_and(~cir1, ~same)
# This should set all cases
aa[circ] = aa0[circ]
aa[left] = 180.
aa[right] = 0.
else:
# Determine the type of calculations to be carried out
noz = np.logical_or(za < 1e-06, za > 180 - 1e-06)
pol = abs(lat) > 90 - 1e-08
pre = np.logical_and(~noz, pol)
non = np.logical_and(~noz, ~pol)
aa[noz] = 0.
aa[pre] = np.rad2deg(np.arctan2(dy[pre], dx[pre]))
dlat = (- sinlat[non] * coslon[non] / r[non] * dx[non] + coslat[non] /
r[non] * dz[non] - sinlat[non] * sinlon[non] / r[non] * dy[non]
)
dlon = (- sinlon[non] / coslat[non] / r[non] * dx[non] + coslon[non] /
coslat[non] / r[non] * dy[non])
aa[non] = (np.rad2deg(np.arccos(r[non] *
dlat / np.sin(np.deg2rad(za[non])))))
fix = np.logical_or(np.isnan(aa), ~np.isreal(aa))
aa[np.logical_and(fix, dlat >= 0)] = 0
aa[np.logical_and(fix, dlat < 0)] = 180
aa[np.logical_and(~fix, dlon < 0)] *= -1
return r, lat, lon, za, aa
[docs]def geocentricposlos2cart(r, lat, lon, za, aa):
"""
Convert from spherical POS/LOS to cartesian POS/LOS
See Contents for a defintion of the geocentric coordinate system.
The local LOS angles are defined following the EAST-NORTH-UP system:
za aa
90 0 points towards north
90 90 points towards east
0 aa points up
FORMAT [x,y,z,dx,dy,dz]=geocentricposlos2cart(r,lat,lon,za,az)
OUT
x Coordinate in x dimension
y Coordinate in y dimension
z Coordinate in z dimension
dx LOS component in x dimension
dy LOS component in y dimension
dz LOS
IN
r Radius
lat Latitude
lon Longitude
za zenith angle
aa azimuth angle
Ported from atmlab. Original author: Bengt Rydberg 2011-10-31
"""
r, lat, lon, za, aa = _broadcast(r, lat, lon, za, aa)
if any(r == 0):
raise Exception("This function is not handling the case of r = 0.")
if any(lat < -90) or any(lat > 90):
raise RuntimeError("The latitude is out of range")
if any(lon < -180) or any(lon > 180):
raise RuntimeError("The longitude is out of range")
if any(za < 0) or any(za > 180):
raise RuntimeError("The zenith angle is out of range")
deg2rad = np.deg2rad(1)
x = np.empty(r.shape)
y = np.empty(r.shape)
z = np.empty(r.shape)
dx = np.empty(r.shape)
dy = np.empty(r.shape)
dz = np.empty(r.shape)
at_pole = abs(lat) > (90 - 1e-8)
if any(at_pole):
s = np.sign(lat[at_pole])
x[at_pole] = 0.
y[at_pole] = 0.
z[at_pole] = s * r[at_pole]
dz[at_pole] = s * np.cos(deg2rad*(za[at_pole]))
dx[at_pole] = np.sin(deg2rad*(za))
dy[at_pole] = dx[at_pole] * np.sin(deg2rad*(aa))
dx[at_pole] *= np.cos(deg2rad*(aa))
not_pole = np.logical_not(at_pole)
if any(not_pole):
latrad = deg2rad * lat[not_pole]
lonrad = deg2rad * lon[not_pole]
zarad = deg2rad * za[not_pole]
aarad = deg2rad * aa[not_pole]
coslat = np.cos(latrad)
sinlat = np.sin(latrad)
coslon = np.cos(lonrad)
sinlon = np.sin(lonrad)
cosza = np.cos(zarad)
sinza = np.sin(zarad)
cosaa = np.cos(aarad)
sinaa = np.sin(aarad)
x[not_pole] = r[not_pole] * coslat
y[not_pole] = x[not_pole] * sinlon
x[not_pole] *= coslon
z[not_pole] = r[not_pole] * sinlat
dr = cosza
dlat = sinza * cosaa
dlon = sinza * sinaa / coslat
dx[not_pole] = (coslat * coslon * dr - sinlat * coslon * dlat -
coslat * sinlon * dlon)
dz[not_pole] = sinlat * dr + coslat * dlat
dy[not_pole] = (coslat * sinlon * dr - sinlat * sinlon * dlat +
coslat * coslon * dlon)
return x, y, z, dx, dy, dz
def get_ellipsoid_semiminor_axis(ellipsoid):
"""Returns the semiminor axis for the ellipsoid
"""
return ellipsoid[0] * np.sqrt(1.0 - ellipsoid[1]**2)
def line_ellipsoid_intersect(x, y, z, dx, dy, dz,
ellipsoid, altitude=0.0):
""" Finds positions of intersection of line with ellipsoid
Solves for d:
$$(X/a)**2 + (Y/a)**2 + (Z/b)**2 - 1 = 0$$
$$X = x + d * dx$$
$$Y = y + d * dy$$
$$Z = z + d * dz$$
where d is a multiplier of sqrt(dx**2 + dy**2 + dz**2), a is the
semimajor axis and b is the semiminor axis. X, Y, Z gives the
point of intersection. If dx, dy, dz is normalized, d is distance
in ellipsoid size units. If no intersect, the code returns np.nan distance
For real d, the halfway point between the points is the tangent point,
though it is below the ellipsoid + altitude surface
Returns:
d: Array of distance parameters, last dimension gives the 2 solutions
Parameters:
x: Geocentric cartesian x-coordinate
y: Geocentric cartesian y-coordinate
z: Geocentric cartesian z-coordinate
dx: Geocentric cartesian dx-view
dy: Geocentric cartesian dy-view
dz: Geocentric cartesian dz-view
elliposid: ellipsoid model [a, e]
altitude: altitude above elipsoid to intersect at [defaults to zero]
Todo:
Refer to line_sphere_intersect if ellipsoid is sphere
"""
# semimajor axis
a = ellipsoid[0] + altitude
# semiminor axis
b = get_ellipsoid_semiminor_axis(ellipsoid) + altitude
# If these are scalars make them arrays
x, y, z, dx, dy, dz = _broadcast(x, y, z, dx, dy, dz)
# A*d**2 + B*d + C = 0, solve for d
A = ((dx**2 + dy**2) / a**2 + dz**2 / b**2).flatten()
B = (2 * ((x * dx + y * dy) / a**2 + z * dz / b**2)).flatten()
C = ((y**2 + x**2) / a**2 + z**2 / b**2 - 1.0).flatten()
d = np.zeros((len(A), 2))
for i in range(len(A)):
roots = np.roots([A[i], B[i], C[i]])
if any(np.isreal(roots)):
d[i, 0] = roots[0]
d[i, 1] = roots[1]
else:
d[i, 0] = np.nan
d[i, 1] = np.nan
sh = []
for i in x.shape:
sh.append(i)
sh.append(2)
return d.reshape(tuple(sh))
def geometric_limb_zenith_angle(ellipsoid, r, lat, lon,
aa=0.0, alt=0.0, za_acc=1e-10):
"""Numerical method to compute the geometric limb zenith angle knowing
other satellite parameters
Assumes the planet is an ellipsoid and calculates the limb by numerically
minimizing the distance between adjacent intersections of the ellipsoid.
The rate of closing in on the right zenith angle should follow 180/30**N,
though numerical errors at very exact accuracies are not guarded against.
Parameters:
ellipsoid: ellipsoid [a, e]
r: geocentric radius of satellite
lat: geocentric latitude of satellite
lon: geocentric longitude of satellite
aa: azimuth angle of satellite
alt: tangent altitude
za_acc: zenith angle accuracy. Warning: Endless loop if too small
Returns:
za: satellite zenith angle for limb view geometry. Shaped as input
"""
if(za_acc == 0):
raise RuntimeError("Zenith accuracy cannot be 0")
r, lat, lon, aa, alt = _broadcast(r, lat, lon, aa, alt)
# Remember shape because we'll flatten these arrays for easy looping
sh = r.shape
r = r.flatten()
lat = lat.flatten()
lon = lon.flatten()
aa = aa.flatten()
alt = alt.flatten()
za = np.empty_like(r)
for i in range(len(r)):
za_min = 0.0
za_max = 180.0
while True: # Until zenith accuracy is reached
# num=31 means 180/30**N is the number of loops to reach accuracy
zenith_angles = np.linspace(za_min, za_max, num=31)
# Get the Cartesian views for the probing zenith accuracies
x, y, z, dx, dy, dz = geocentricposlos2cart(r[i], lat[i], lon[i],
zenith_angles, aa[i])
# Get intersections with the ellipsoid at the tangent altitude
d = line_ellipsoid_intersect(x, y, z, dx, dy, dz,
ellipsoid, alt[i])
# If we have somehow designed a scenario that is impossible or bad
if np.isnan(d).all():
raise RuntimeError("Not possible to get a tangent point. " +
"All tries miss the ellipsoid.")
elif not np.isnan(d).any():
raise RuntimeError("Not possible to get a tangent point. " +
"Inside ellipsoid or atmosphere.")
# Limb zenith from minimizing distance between intersections
D = d[:, 0] - d[:, 1]
# These are the points of interest
these = np.nan_to_num(D) > 0
# The first of the zenith angles hitting the atmosphere is updated
za_max = zenith_angles[these][0]
# The last of the zenith angles missing the atmosphere is updated
za_min = zenith_angles[np.logical_not(these)][-1]
# Are we there yet?
if (za_max-za_min) < za_acc:
break
# Note assignment to flat zenith angle array
za[i] = za_max
# Note that we return with the proper shape
return za.reshape(sh)
def sphere_plane_intersection(pos, r, theta=np.linspace(-180., 180.)):
"""Computes the intersection of a plane with a sphere
Returns cartesian x, y, z at the provided angles. Sphere is centered
at 0, 0, 0
The solution is from
.. math::
\\vec{r}_c = \\vec{r}_p + r_0 \\left( \\vec{v}_1\\cos\\theta +
\\vec{v}_2\\sin\\theta \\right)
rp is the input position of the center of the interesecting circle,
.. math::
r_0 = \\sqrt{r^2 - \\vec{r_p} \\cdot \\vec{r_p}},\\;\\;
\\vec{n} = \\vec{r}_p/||\\vec{r}_p||,\\;\\;
\\vec{n}\\cdot\\vec{v}_1\\approx0,\\;\\;
\\vec{v}_2 = \\vec{n} \\times \\vec{v}_1
where v1 is so defined that it fully overlaps with the azimuth angle at the
equator. At other latitudes, the azimuth and theta are only the same at 0
degrees.
Parameters:
pos (N, 3-array): position at the center of the circle-like shape that
the intersection will form (i.e., the tangent point of a satellite at
the sphere's radius). Cannot be the 0-vector
r (float): radius of sphere in same units as pos. Cannot be 0.0
theta (n-array): determines positions by angle [degrees]
Returns:
(n, N, 3)-array of points intersecting plane and sphere
Examples:
Find possible positions that a satellite at an orbit about
600 km above the surface can see a tangent point 30 km above
the surface of Earth
>>> x, y, z = geodetic2cart(30e3, 30., 30.)
>>> sat = sphere_plane_intersection(np.array([x, y, z]), 6978e3)
"""
assert r > 0.0, "Cannot work with zero radius sphere"
assert len(pos.shape) == 2, "Need (N, 3) as pos-shape"
d = np.linalg.norm(pos, axis=1)
assert all(d > 0.0), "Cannot work with the 0-vector for position"
number_of_angles = len(theta)
number_of_positions = pos.shape[0]
points = np.empty((number_of_angles, number_of_positions, 3))
if number_of_angles == 0:
return points
# If the sphere is under the plane, or tangential to the plane
for ii in range(number_of_positions):
if d[ii] > r:
for i in range(number_of_angles):
points[i, ii, :] = np.array([np.nan, np.nan, np.nan])
continue
elif d[ii] == r:
for i in range(number_of_angles):
points[i, ii, :] = pos[ii, :]
continue
# The (normalized/unit) plane vector
n = pos[ii] / d[ii]
# Southern vector
ones, lat, lon = cart2geocentric(n[0], n[1], n[2])
if(lat == 90.0):
v1 = np.array([cosd(lon), sind(lon), 0.]) * sind(lat)
else:
v1 = np.array([cosd(lon) * sind(lat),
sind(lon) * sind(lat),
cosd(lat-180.0)])
# Eastern vector
v2 = -np.cross(n, v1) # points East
# Radius of the circle
rc = np.sqrt(r**2 - d[ii]**2)
for i in range(number_of_angles):
points[i, ii, :] = pos[ii, :] + rc * (cosd(theta[i]) * v1 +
sind(theta[i]) * v2)
return points
def tunnel_distance(lat1, lon1, lat2, lon2):
"""Calculate the tunnel distance between two points on the Earth's surface
Args:
lat1: Single number or numpy array with latitudes in degrees.
lon1: Single number or numpy array with longitudes in degrees.
lat2: Single number or numpy array with latitudes in degrees.
lon2: Single number or numpy array with longitudes in degrees.
Returns:
Tunnel distance in meters
"""
points1 = np.column_stack(
geocentric2cart(constants.earth_radius, lat1, lon1)
)
points2 = np.column_stack(
geocentric2cart(constants.earth_radius, lat2, lon2)
)
return np.sqrt(np.sum((points2 - points1)**2, axis=1))
def _broadcast(*args):
""" Similar to broadcast_arrays in numpy but with minimum output size (1,)
"""
shape = np.broadcast(*args).shape
if not shape:
shape = (1,)
return [np.broadcast_to(array, shape) for array in args]