karney.geodesic

Functions

`sphere_distance_rad`(lat1, lon1, lat2, lon2)	Returns surface distance between positions A and B as well as the azimuths.
`reckon`(lat_a, lon_a, distance, azimuth[, a, f, ...])	Returns position B computed from position A, distance and azimuth.
`distance`(lat_a, lon_a, lat_b, lon_b[, a, f, degrees])	Returns shortest surface distance between positions A and B on an ellipsoid.

Module Contents

karney.geodesic.sphere_distance_rad(lat1, lon1, lat2, lon2)

Returns surface distance between positions A and B as well as the azimuths.

Parameters:

lat1 (real scalars or vectors of length m.) – latitude(s) and longitude(s) [rad] of position A.
lon1 (real scalars or vectors of length m.) – latitude(s) and longitude(s) [rad] of position A.
lat2 (real scalars or vectors of length n.) – latitude(s) and longitude(s) [rad] of position B.
lon2 (real scalars or vectors of length n.) – latitude(s) and longitude(s) [rad] of position B.

Returns:

distance (real scalars or vectors of length max(m,n).) – Surface distance [rad] from A to B on the sphere
azimuth_a, azimuth_b (real scalars or vectors of length max(m,n).) – direction [rad] of line at position A and B relative to North, respectively.

Notes

Solves the inverse geodesic problem of finding the length and azimuths of the shortest geodesic between points specified by lat1, lon1, lat2, lon2 on a sphere.

Examples

“Surface distance”

Find the surface distance sAB (i.e. great circle distance) between two positions A and B. The heights of A and B are ignored, i.e. if they don’t have zero height, we seek the distance between the points that are at the surface of the Earth, directly above/below A and B. Use Earth radius 6371e3 m. Compare the results with exact calculations for the WGS-84 ellipsoid.

In geodesy this is known as “The second geodetic problem” or “The inverse geodetic problem” for a sphere/ellipsoid.

Solution for a sphere:

>>> import numpy as np
>>> from karney.geodesic import rad, sphere_distance_rad, distance

>>> latlon_a = (88, 0)
>>> latlon_b = (89, -170)
>>> latlons = latlon_a + latlon_b
>>> r_Earth = 6371e3  # m, mean Earth radius
>>> s_AB = sphere_distance_rad(*rad(latlons))[0]*r_Earth
>>> s_AB = distance(*latlons, a=r_Earth, f=0, degrees=True)[0] # or alternatively

>>> "Ex5: Great circle = {:5.2f} km".format(s_AB / 1000)
'Ex5: Great circle = 332.46 km'

Exact solution for the WGS84 ellipsoid:

>>> s_12, azi1, azi2 = distance(*latlons, degrees=True)
>>> "Ellipsoidal distance = {:5.2f} km".format(s_12 / 1000)
'Ellipsoidal distance = 333.95 km'

“A and azimuth/distance to B”

We have an initial position A, direction of travel given as an azimuth (bearing) relative to north (clockwise), and finally the distance to travel along a great circle given as sAB. Use Earth radius 6371e3 m to find the destination point B.

In geodesy this is known as “The first geodetic problem” or “The direct geodetic problem” for a sphere.

>>> import numpy as np
>>> from karney.geodesic import reckon
>>> lat, lon = 80, -90
>>> msg = "Ex8, Destination: lat, lon = {:4.4f} deg, {:4.4f} deg"

Greatcircle solution:

>>> datum = dict(a=6371e3, f=0)
>>> lat2, lon2, azi_b = reckon(lat, lon, distance=1000, azimuth=200, degrees=True, **datum)

>>> msg.format(lat2, lon2)
'Ex8, Destination: lat, lon = 79.9915 deg, -90.0177 deg'

>>> np.allclose(azi_b, -160.0174292682187)
True

Exact solution:

>>> lat_b, lon_b, azi_b = reckon(lat, lon, distance=1000, azimuth=200, degrees=True)
>>> msg.format(lat_b, lon_b)
'Ex8, Destination: lat, lon = 79.9916 deg, -90.0176 deg'

>>> np.allclose(azi_b, -160.01742926820506)
True

“Surface distance”

Find the surface distance sAB (i.e. great circle distance) between two positions A and B. The heights of A and B are ignored, i.e. if they don’t have zero height, we seek the distance between the points that are at the surface of the Earth, directly above/below A and B. Use Earth radius 6371e3 m. Compare the results with exact calculations for the WGS-84 ellipsoid.

In geodesy this is known as “The second geodetic problem” or “The inverse geodetic problem” for a sphere/ellipsoid.

Solution for a sphere:

>>> import numpy as np
>>> from karney.geodesic import rad, sphere_distance_rad, distance

>>> latlon_a = (88, 0)
>>> latlon_b = (89, -170)
>>> latlons = latlon_a + latlon_b
>>> r_Earth = 6371e3  # m, mean Earth radius
>>> s_AB = sphere_distance_rad(*rad(latlons))[0]*r_Earth
>>> s_AB = distance(*latlons, a=r_Earth, f=0, degrees=True)[0] # or alternatively

>>> "Ex5: Great circle = {:5.2f} km".format(s_AB / 1000)
'Ex5: Great circle = 332.46 km'

Exact solution for the WGS84 ellipsoid:

>>> s_12, azi1, azi2 = distance(*latlons, degrees=True)
>>> "Ellipsoidal distance = {:5.2f} km".format(s_12 / 1000)
'Ellipsoidal distance = 333.95 km'