karney.geodesic =============== .. py:module:: karney.geodesic .. autoapi-nested-parse:: Geodesic module --------------- Functions --------- .. autoapisummary:: karney.geodesic.sphere_distance_rad karney.geodesic.reckon karney.geodesic.distance Module Contents --------------- .. py:function:: sphere_distance_rad(lat1, lon1, lat2, lon2) Returns surface distance between positions A and B as well as the azimuths. :param lat1: latitude(s) and longitude(s) [rad] of position A. :type lat1: real scalars or vectors of length m. :param lon1: latitude(s) and longitude(s) [rad] of position A. :type lon1: real scalars or vectors of length m. :param lat2: latitude(s) and longitude(s) [rad] of position B. :type lat2: real scalars or vectors of length n. :param lon2: latitude(s) and longitude(s) [rad] of position B. :type lon2: real scalars or vectors of length n. :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. .. rubric:: 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. .. rubric:: 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' .. seealso:: :py:obj:`distance` .. py:function:: reckon(lat_a, lon_a, distance, azimuth, a=6378137, f=1.0 / 298.257223563, long_unroll=False, degrees=False) Returns position B computed from position A, distance and azimuth. :param lat_a: latitude(s) and longitude(s) of position A [rad or deg]. :type lat_a: real scalars or vectors of length k. :param lon_b: latitude(s) and longitude(s) of position A [rad or deg]. :type lon_b: real scalars or vectors of length k. :param distance: ellipsoidal distance [m] between position A and B. :type distance: real scalar or vector of length m. :param azimuth: azimuth [rad or deg] of line at position A. :type azimuth: real scalar or vector of length n. :param a: Semi-major half axis of the Earth ellipsoid given in [m]. :type a: real scalar, default WGS-84 ellipsoid. :param f: Flattening [no unit] of the Earth ellipsoid. If f==0 then spherical Earth with radius a is used in stead of WGS-84. :type f: real scalar, default WGS-84 ellipsoid. :param long_unroll: Controls the treatment of longitude. If it is False then the lon2 is reduced to the range [-180, 180) or [-pi, pi). If it is True, then (lon2 - lon1) determines how many times and in what sense the geodesic has encircled the ellipsoid. :type long_unroll: bool :param degrees: True if latitude, longitude and azimuth is given in degress. :type degrees: bool :returns: * **lat_b, lon_b** (*arrays of length max(k,m,n).*) -- latitude(s) and longitude(s) of position B [rad or deg]. * **azimuth_b** (*real scalars or vectors of length max(k,m,n).*) -- azimuth [rad or deg] of line at position B. .. rubric:: Notes This function solves the direct geodesic problem of finding the final point B and azimuth_b given the position of point A and the azimuth and distance from A to B. Keep :math:`|f| <= 1/50` for full double precision accuracy. .. rubric:: Examples **"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 .. seealso:: :py:obj:`distance` .. py:function:: distance(lat_a, lon_a, lat_b, lon_b, a=6378137, f=1.0 / 298.257223563, degrees=False) Returns shortest surface distance between positions A and B on an ellipsoid. :param lat_a: latitude(s) and longitude(s) of position A [deg or rad]. :type lat_a: real scalars or vectors of length m. :param lon_a: latitude(s) and longitude(s) of position A [deg or rad]. :type lon_a: real scalars or vectors of length m. :param lat_b: latitude(s) and longitude(s) of position B [deg or rad]. :type lat_b: real scalars or vectors of length n. :param lon_b: latitude(s) and longitude(s) of position B [deg or rad]. :type lon_b: real scalars or vectors of length n. :param a: Semi-major axis of the Earth ellipsoid given in [m]. :type a: real scalar, default WGS-84 ellipsoid. :param f: Flattening [no unit] of the Earth ellipsoid. If f==0 then spherical Earth with radius a is used in stead of WGS-84. :type f: real scalar, default WGS-84 ellipsoid. :param degrees: True if latitudes, longitudes and azimuths are given in degress :type degrees: bool :returns: * **distance** (*real scalars or vectors of length max(m,n).*) -- Surface distance [m] from A to B on the ellipsoid * **azimuth_a, azimuth_b** (*real scalars or vectors of length max(m,n).*) -- direction [rad or deg] of line at position a and b relative to North, respectively. .. rubric:: Notes Solves the inverse geodesic problem of finding the length and azimuths of the shortest geodesic between points specified by lat_a, lon_a, lat_b, lon_b on an ellipsoid. Keep :math:`|f| <= 1/50` for full double precision accuracy. .. rubric:: 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' .. seealso:: :py:obj:`sphere_distance_rad`, :py:obj:`reckon`