Geodesy#

class Geodetic[source]#

A class to manage both transform and tangent point calculations using oblate spheroid geometry. The transform operations allow conversion between geocentric, xyz, coordinates and geodetic, latitude, longitude and height coordinates. We also provide right-handed, geographic-based, unit vectors for any location

The tangent point calculations include finding the location of a tangent point as well as the entrance and exit locations of a ray through altitude shells.

The class is initialized with the default WGS84 oblate spheroid which should be suitable for most work.

The class is fully implemented in pure python and has been written allow multiple calculations to be performed using array-based arithemetic.

Attributes:
semi_major_axis

Semi-major axis of the current ellipsoid in meters

semi_minor_axis

Semi-minor axis of the current ellipsoid in meters

Methods

get_osculating_spheroid(geodeticlatitude, ...)

Calculates the sphere that best fits the ellipsoid surface (height = 0) at the geodetic location.

llh_from_xyz(xyz)

Converts positions specified in xyz geocentric coordinates in meters to the equivalent geodetic (latitude, longitude, height) coordinates.

set_geoid([A, F, geoid])

Defines the shape parameters of the ellipsoid that will be used for future calculations.

xyz_altitude_intercepts(Hnd, observerxyznd, ...)

Calculate the entrance and exit points of a look vector with a given shell at height H.

xyz_from_llh(llh)

Converts positions specified in geodetic coordinates (lat degrees,lon degrees, height meters) to the equivalent geocentric (x,y,z meters ) coordinates.

xyz_lookvector_from_tangent_altitude(h, ...)

Finds the look vector towards a given tangent altitude at a given azimuth at the observer location.

xyz_north_east_down([xyz, llh])

Returns the three right-handed, orthogonal unit vectors, north, east and down at the provided locations.

xyz_observer_sun_look_for_given_sun_and_tangent_point(...)

Generates an observer who is looking along an instrument bore-sight at a given tangent point location.

xyz_tangent_point_location(observerxyz, lookxyz)

Calculates the xyz geocentric coordinates of the tangent point given a look vector away from an observers location towards the Earth.

xyz_tangent_point_location_legacy(...[, ...])

Calculates the xyz geocentric coordinates of the tangent point given a look vector away from an observers location towards the Earth.

xyz_west_south_up([xyz, llh])

Returns the three right-handed, orthogonal unit vectors, west, south and up at the provided locations.

get_osculating_spheroid(geodeticlatitude: float, geodeticlongitude: float)[source]#

Calculates the sphere that best fits the ellipsoid surface (height = 0) at the geodetic location. This algorithm best fits the spheroid in the latitudinal direction (ie North South). A better choice for the future may be to use a look vector. The osculating sphere is passed to radiative transfer models to “best” represent the true earth at the region of interest using a spherical system.

Parameters:
geodeticlatitude: float

The geodetic latitude in degrees of the location where we need to fit the osculating sphere

geodeticlongitude: float

The geodetic longitude in degrees of the location where we need to fit the osculating sphere.

Returns:
Tuple[ radius:float, offset:Tuple[float, float]]:

Returns a two element tuple containing the radius of the best fit osculating spheroid in the first element and the offset of the origin of the osculating sphere from the origin of the true geoid in the second element.

Notes

The radius of curvature is given by:-

\[R = \frac{ \left[1 + \left(\frac{dy}{dx}\right)^2\right]^{\frac{3}{2}}}{\frac{d^2y}{dx^2}}\]

and for an ellipse:

\[R = \frac{1}{ab} \left[ \frac{a^2}{b^2}y^2 + \frac{b^2}{a^2}x^2\right]^{\frac{3}{2}}\]
llh_from_xyz(xyz: ndarray) ndarray[source]#

Converts positions specified in xyz geocentric coordinates in meters to the equivalent geodetic (latitude, longitude, height) coordinates. The conversion is non-trivial and this implementation uses the algorithm developed by Fukushima. In theory, the algorithm is described as an approximation but in practice it is numerically comparable to exact solutions operating, which provide solutions at the nanometer level of precision, but struggle to solve the cubic equations with any higher precision.

The algorithm is based upon the GCONV2H subroutine provided in file xgconv2.txt which is provided by Toshio Fukushima on his ResearchGate site.

The code is implemented to efficiently convert a multi-dimensional array of positions in one pass of the routine.

Parameters:
llh: array[…, 3], optional

An array of position vectors given in geocentric, xyz, meters. The array can be multi-dimensional but the last dimension must be used for the (x,y,z) fields and its size must be equal to 3, [0=X, 1=Y, 2=Z].

Returns:
array[…, 3]:

Returns the geodetic coordinates of the input geocentric locations. The array is the same size and shape as the input array but the last dimensions, which is still size 3, is now used for the geodetic coordinates latitude=0, longitude = 1, height =2. latitude is given in degrees, -90.0 to +90.0, longitude is given in degrees, usually 0 to 360 (but this is not rigorously enforced), and height is given as meters above the geodetic surface (aka sea level)

property semi_major_axis#

Semi-major axis of the current ellipsoid in meters

property semi_minor_axis#

Semi-minor axis of the current ellipsoid in meters

set_geoid(A: float = None, F: float = None, geoid: str = None)[source]#

Defines the shape parameters of the ellipsoid that will be used for future calculations. The shape is defined by the semi-major axis, A, in meters and the Earth flattening parameter , F. The most common usage is to use a standard, pre-defined geoid. This method only needs to be called if the default oblate spheroid used in the constructor, wgs84, is not suitable.

Parameters:
geoid: str, optional

Predefined, standard oblate spheroids:

  • ‘geocentric’: (6371200.0, 0.0)

  • ‘iau1976’: (6378140.0, 0.00335281)

  • ‘grs80’: (6378137.0, 1.0/298.257222101)

  • ‘merit83’: (6378137.0, 1.0/298.257)

  • ‘wgs84’: (6378137.0, 1.0/298.257223563)}

If it is not set the both A and F must be set. A and F are ignored in this option is set.

A: float, optional

The semi-major axis in meters. This argument is ignored if geoid is set. Argument F must also be set if this argument is set.

F: float, optional

The Earth flattening parameter. This argument is ignored if geoid is set. Argument A must also be set if this argument is set

xyz_altitude_intercepts(Hnd: ndarray, observerxyznd: ndarray, lookxyznd) Tuple[ndarray, ndarray, ndarray][source]#

Calculate the entrance and exit points of a look vector with a given shell at height H. The routine allows the user to pass in a multi-dimensional arrays of look vector, observer position and shell height. The calculation is only meaningful in the look vector is towards the Earth and the required shell height is at or above the tangent point, i.e. the shell height must lie between the tangent point and the observer.

The entrance and exit points are approximately symmetric about the tangent point. The entrance point is located towards the observer on the near-side of the tangent point. The exit point is on the far side, furthest from the tangent point.

The algorithm is internally implemented to iterate over the arrays and process the intercepts one case at a time. It currently uses scipy.optimize.zbrentq to find the shell heights which is not conducive to array based operations. Consequently, the method, simply due to its array looping implementation, may be substantially slower than the other methods in this object. This should not be a significant issue for hundreds or thousands of rays but will be an issue for millions of rays.

Parameters:
Hnd: np.ndarray

The height in meters of the target shell intercepts. This is typically given as an array with the same shape as observerxynd and lookxynd except the last dimension is 1 rather than 3. It is also allowed to be a scalar or a simple sequence for situations where that is compatible with the other parameters

observerxyznd: np.ndarray

A multi-dimensional array that provides the geocentric, xyz location of the observer in meters. The array can be any shape but the last dimension must be 3. It is recommended that this array has the same shape as lookxyznd. The observer is typically a spacecraft or high altitude platform.

lookxyznd: np.ndarray

A multi-dimensional array that provides the geocentric, xyz, look vector of each ray as a dimensionless unit vector. The rays are treated as being away from the observer, towards the Earth, which is opposite to the actual direction of light propagation. The array can be any shape but the last dimension must be 3. It is recommended that this array has the same shape as observerxyznd.

Returns:
Tuple[entrypoint: np.ndarray, exitpoint: np.ndarray, ok: np.ndarray]

Three element tuple containing (i) the location of the entry point of the ray into the shell, (ii) the location of the exit point of the ray from the shell and (iii) a boolean status flag array. All variables are the same shape as the input arrays.

xyz_from_llh(llh: ndarray | Tuple[float, float, float] | List[float]) ndarray[source]#

Converts positions specified in geodetic coordinates (lat degrees,lon degrees, height meters) to the equivalent geocentric (x,y,z meters ) coordinates. The conversion is implemented to it can efficiently process a multi-dimensional array of positions in one pass.

Parameters:
llh: array[…, 3], optional

An array of position vectors given in geodetic latitude, longitude, height (degrees, degrees, meters). The array can be multi-dimensional but the last dimension is reserved for the lat, lon and height fields and it size must be equal to 3, [0=lat, 1=long, 2=height].

Returns:
array[…, 3]:

The geocentric, xyz, locations in meters corresponding to the input geodetic locations. The array is the same size as the input array except the last dimension, which is still of size 3, is now reserved for X=0, Y=1, Z=2.

xyz_lookvector_from_tangent_altitude(h: float, observerxyz: ndarray, boresightazimuthxyz: ndarray | float, convergence_meters=0.1)[source]#

Finds the look vector towards a given tangent altitude at a given azimuth at the observer location. The observer is typically a spacecraft position and should be (well) above the requested tangent point (e.g. more than 5 km). The algorithm makes a first approximation using the osculating sphere which may fail if the requested tangent altitude is too close to the observers altitude.

An iterative method is used to find the look vector required for the given tangent altitude. The algorithm uses oblate spheroid Earth and straight line geometry. Solution iterates until calculated tangent altitude is within convergence_meters meters of the requested height.

Parameters:
h: float

The requested tangent height in meters above the surface of the Earth. This must be a scalar

observerxyx: np.ndarray [3,]

The geocentric, xyz, location of the observer in meters.

boresightazimuthxyz: np.ndarray[3,] or float

The horizontal unit vector at the observer’s location that defines the azimuth of the required look vector. The parameter can also be specified as a floating point number that gives the compass bearing in degrees (0= North, 90= East etc.) of the desired unit vector.

convergence_meters: float, optional

Specifies the convergence. The returned solution will be within this height distance in meters from the target tangent point It may prove difficult to get consistency much better than 0.1 meters due to computational round off issues. default is 0.1 meters

Returns:
np.ndarray[3,]

The geocentric, xyz, dimensionless unit vector that looks from the observer toward the target tangent point at the given horizontal azimuth.

xyz_north_east_down(xyz=None, llh=None) Tuple[ndarray, ndarray, ndarray][source]#

Returns the three right-handed, orthogonal unit vectors, north, east and down at the provided locations. It is a small wrapper around xyz_west_south_up and all the parameters are described in that function with the returning variables being changed from (west, south, up) to (north, east, down)

xyz_observer_sun_look_for_given_sun_and_tangent_point(boresight_tanpt_llh: ndarray, omega_ocf: ndarray, lookazi_at_tanpt_deg: float = 75.0, sza: float = 60.0, saa: float = 135.0, observer_height_m: float = 500000.0)[source]#

Generates an observer who is looking along an instrument bore-sight at a given tangent point location. The Sun is chosen so it is at a given solar zenith angle and solar azimuth angle at that tangent point. The observer is located at a given height. The location of the observer and sun are calculated that meet the specified requirements. A set of look vectors are generated using the direction angles \(\Omega\) given in the optical control frame, centered around the boresight.

Parameters:
boresight_tanpt_llh: Tuple[lat, long , height]

THe latitude, longitude and height of the tangent point that the boresight should look at

omega_ocf: np.ndarray [N, 2]

An array of N look directions specified using the \(\Omega_y\) and \(\Omega_x\) angles in the optical control frame. Briefly, \(\Omega_y\) is upward elevation from the boresight and \(\Omega_x\) is azimuth in an anti-clockwise direction from the bore sight, looking away from the instrument towards port (left). All angles are given in degrees. The array is typically a 2-D array with the last dimension set to 2. \(\Omega_y\) corresponds to the […, 0] elements and \(\Omega_x\) corresponds to thge […, 1] elements. The variable can be a list or tuple, or anything else that can be easily coerced into an array. The variable can also be 1 dimensional (eg [N]) in which case it is coerced into a 2 Dimensional array with \(\Omega_x\) equal to zero.

lookazi_at_tanpt_degfloat

The geographic azimuth at the tangent point of the look vector from the observer to the tangent point. It is specified in degrees 0=N, 90=E etc. Default is 75 degrees.

szafloat

Solar zenith angle in degrees at the desired tangent point. Default is 60 degrees.

saafloat

Solar azimuth angle in degrees at the desired tangent point. default is 135 degrees

observer_heightfloat

The height of the observer above the surface in meters. This must be higher than the requested tangent point. Default is 500,000.0

Returns:
Tuple[observerxyz: np.ndarray[3], sunxyz: np.ndarray[3], lookxyz: np.ndarray[…, 3]]

A three element tuple. The first element is the geocentric, xyz, location of the observer (in meters). The second element is a geocentric, xyz, unit vector toward the sun. The third element is the array of look unit vectors away from the instrument corresponding to each of the look directions given in ``omega_ocf`. The look vectors are returned in the same shape as the look directions except the last dimension is now 3 instead of 2 nad contains the geocentric, xyz, unit vector of each look direction.

xyz_tangent_point_location(observerxyz: ndarray, lookxyz: ndarray, checkconvergence=None, check_tangent_behind=False) ndarray[source]#

Calculates the xyz geocentric coordinates of the tangent point given a look vector away from an observers location towards the Earth. The earth is modelled as an oblate spheroid. All calculations assume straight ray paths. This is the second technique we developed to find the tangent point and seems to be the adhoc default. In practice it is neither substantially quicker or more accurate than the original technique. It seems to work well at sub-meter accuracy (typically around 0.1 meter accuracy), but beyond that more careful handling of numeric roundoff may be required.

The code is implemented to process multiple look vectors and observer positions in one pass of the code using multi-dimensional arrays. However both the observer and look arrays must be the same shape as the same shape flattening factor is iteratively applied to both fields; it fails if they have different shapes. The legacy version of this function can handle broadcasting, eg a single observer and multiple lines of sight.

Parameters:
observerxyz: array[…,3]

An array of observer positions given in geocentric xyz coordinates (meters). It can be multidimensional with the last dimension equal to 3, [0=X, 1=Y, 2=Z]. It should have the same shape as variable lookxyz

lookxyz: array[…,3]

An array of look unit vectors given in geocentric xyz coordinates (dimensionless). The look vectors must be away from the observer towards the Earth. It can be multidimensional with the last dimension equal to 3, [0=X, 1=Y, 2=Z]. Its shape must be the same as observerxyz.

checkconvergence: float or None

Diagnostic internal flag used to provide information about the true tangent point assuming its close to this solution. Its reserved for internal usage.

Returns:
array[…,3]

The location of the tangent point in geocentric, xyz coordinates (meters). The final array size is determined from the broadcasting of observerxyz with lookxyz, but the last dimension will be 3. Locations that failed to find a tangent point are set to Nan. This typically occurs if the look vector is not towards the Earths limb.

xyz_tangent_point_location_legacy(observerxyz: ndarray, lookxyz: ndarray, checkconvergence=None)[source]#

Calculates the xyz geocentric coordinates of the tangent point given a look vector away from an observers location towards the Earth. The earth is modelled as an oblate spheroid. All calculations assume straight ray paths. This is the first technique we developed to find the tangent point. It seems to work well at sub-meter accuracy (typically around 0.1 meter accuracy), but beyond that more careful handling of numeric roundoff may be required.

The code is implemented to process multiple look vectors and observer positions in one pass of the code using multi-dimensional arrays and broadcasting.

Parameters:
observerxyz: array[…,3]

An array of observer positions given in geocentric xyz coordinates (meters). It can be multidimensional with the last dimension equal to 3, [0=X, 1=Y, 2=Z]

lookxyz: array[…,3]

An array of look unit vectors given in geocentric xyz coordinates (dimensionless). The look vectors must be away from the observer towards the Earth. It can be multidimensional with the last dimension equal to 3, [0=X, 1=Y, 2=Z]. Its shape must be the same as observerxyz or at least be broadcast compatible.

checkconvergence: float or None

Diagnostic internal flag used to provide information about the true tangent point assuming its close to this solution. Its reserved for internal usage.

Returns:
array[…,3]

The location of the tangent point in geocentric, xyz coordinates (meters). The final array size is determined from the broadcasting of observerxyz with lookxyz, but the last dimension will be 3. Locations that failed to find a tangent point are set to Nan. This typically occurs if the look vector is not towards the Earths limb.

xyz_west_south_up(xyz=None, llh=None) Tuple[ndarray, ndarray, ndarray][source]#

Returns the three right-handed, orthogonal unit vectors, west, south and up at the provided locations. The unit vectors are always given in xyz geocentric coordinates. The locations at which the unit vectors are required can be specified in either geocentric xyz coordinates using the keyword xyz or in geodetic latitude, longitude, height using keyword llh.

The algorithm allows all the locations to be specified using multi-dimensional arrays with the last dimension set to a size of 3 ali_elements. This allows the code to execute all location calculations in one pass of the code.

It is not recommended, but it is possible to specify locations using both keywords in which case the user must ensure they are identically sized arrays and they have the geocentric and geodetic coordinates pertaining to the same exact locations

Parameters:
xyz: array[…, 3], optional

The geocentric, xyz, locations in meters where the three unit vectors will be calculated. It can be multidimensional with the last dimension equal to 3, [0=X, 1=Y, 2=Z]. It is usually left undefined if keyword llh is used.

llh: array[…, 3], optional

An array of look unit vectors given in geodetic lat, lon, height (derees, degrees, meters). It can be multi-dimensional but the last dimension must be equal to 3, [0=lat, 1=long, 2=height]. It is usually left undefined if keyword xyz is used.

Returns:
Tuple[ west:array[…, 3], south:array[…, 3], up:array[…, 3]]

A three element tuple storing the west, south and up unit vectors expressed as dimensional numbers in geocentric, xyz coordinates. The array sizes are the same as the arrays used to define the input locations.