from typing import Tuple,List, Union
import numpy as np
#------------------------------------------------------------------------------
# class WGS84
#------------------------------------------------------------------------------
[docs]class WGS84:
"""
Implements a few geodetic conversion functions using the WGS 84 ellipsoid. This is identical to ISKGeodetic WGS84 option in the sasktran libraries but is coded up in pure pytho
"""
semi_major = 6378137.0
flattening = 1.0 / 298.257223563
#------------------------------------------------------------------------------
# geodetic_to_geocentric
#------------------------------------------------------------------------------
@classmethod
[docs] def geodetic_to_geocentric( cls,
latitude : Union[ float, np.ndarray, List[float]],
longitude : Union[ float, np.ndarray, List[float]],
altitude : Union[ float, np.ndarray, List[float]]) ->Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
"""
This function conevrts geodetic latitude, longitude and altitude into geocentric x y and z coordinates using the WGS84 ellipsoid. It also
returns the unit vectors for north, east and down at the location. The implementation will accept scalar or 1-D array based valuesr
:param latitude:
The geodetic latitude in degrees. Either scalar or list/array.
:param longitude:
The geodetic longitude in degrees. Either scalar or list/array, must be same size as latitude
:param altitude:
The altitude above mean sea level in meters. Either scalar or list/array, must be same size as latitude.
:return:
A tuple of 4 arrays ( location, north, east, down)
- location is the geocentric, xyz, vector location in meters. It is an array (3, npts) where 0 is the x component , 1 is y component and 2 is the z components
- north, the unit vector in geocentric coordinates at the location in the direction of north. It is an array of (3, npts) and is same as location except it is a unit vector
- east, the unit vector in geocentric coordinates at the location in the direction of east. It is an array of (3, npts) and is same as location except it is a unit vector
- down, the unit vector in geocentric coordinates at the location in the direction of down. It is an array of (3, npts) and is same as location except it is a unit vector
Author: Daniel Letros
Date Revised: 2017-02-24
"""
npts = 1 if np.isscalar(latitude) else len(latitude)
fm1 = 1.0 - WGS84.flattening # get 1-f
fm12 = fm1*fm1 # get the square of the value
latrad = np.deg2rad(latitude) # convert the latitude to radians
lonrad = np.deg2rad(longitude) # convert the longitude to radians
coslat = np.cos( latrad ) # get the cosine of the latitude as temporary array
sinlat = np.sin( latrad ) # get the sine of the latitude as temporary array
c = 1.0 / np.sqrt( coslat*coslat + (fm12*sinlat*sinlat) ) # calculate c
s = fm12*c # and s
h = (cls.semi_major * c + altitude) * coslat # get the horizontal component
x = h * np.cos(lonrad) # get the x component, array(npts)
y = h * np.sin(lonrad) # get the y component, array(npts)
z = (cls.semi_major * s + altitude) * sinlat # get the z component, array(npts)
# ---- calculate the down, north and east unit vectors
zero = np.zeros( [npts] ) # create an array of zeroes as a temporaray array
horiz = np.squeeze( np.array([x,y,zero])) # Create an array of horizontal vectors ( 3, npts)
hlen = np.linalg.norm(horiz,axis=0, keepdims=True) # get the length of each vector (npts)
bad = np.where( hlen == 0.0) # find out any that are zero (i.e. at the poles)
hlen[bad] = 1.0 # and set the bad ones to a value that does not fail to Nan
horiz = horiz/hlen # Divide each column by the normalized length, ie produce horizontal unit-vector
vertical = np.squeeze( np.array( [zero, zero, zero+1.0] )) # vertical is a vector parallel to spin axis
down = -horiz*coslat - vertical*sinlat # get vertical down unit vector at observer's location, in geocentric corrdinates
north = vertical*coslat - horiz*sinlat # get due north unit vector at observer's location
east = np.cross( down, north, axis=0 ) # get due east unit vector at observer's location
geocentric = np.array( [x,y,z] ) # get the geocentric as a vector
return (geocentric, north, east, down) #