diff --git a/geodesy/haversine_distance.py b/geodesy/haversine_distance.py index 93e625770f9d..824a4a255a86 100644 --- a/geodesy/haversine_distance.py +++ b/geodesy/haversine_distance.py @@ -1,54 +1,44 @@ -from math import asin, atan, cos, radians, sin, sqrt, tan - -AXIS_A = 6378137.0 -AXIS_B = 6356752.314245 -RADIUS = 6378137 +from math import asin, cos, radians, sin, sqrt def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> float: """ - Calculate great circle distance between two points in a sphere, - given longitudes and latitudes https://en.wikipedia.org/wiki/Haversine_formula - - We know that the globe is "sort of" spherical, so a path between two points - isn't exactly a straight line. We need to account for the Earth's curvature - when calculating distance from point A to B. This effect is negligible for - small distances but adds up as distance increases. The Haversine method treats - the earth as a sphere which allows us to "project" the two points A and B - onto the surface of that sphere and approximate the spherical distance between - them. Since the Earth is not a perfect sphere, other methods which model the - Earth's ellipsoidal nature are more accurate but a quick and modifiable - computation like Haversine can be handy for shorter range distances. + Calculate the great-circle distance between two points + on the Earth specified by latitude and longitude using + the Haversine formula. Args: - lat1, lon1: latitude and longitude of coordinate 1 - lat2, lon2: latitude and longitude of coordinate 2 + lat1, lon1: Latitude and longitude of point 1 in decimal degrees. + lat2, lon2: Latitude and longitude of point 2 in decimal degrees. + + Returns: - geographical distance between two points in metres + Distance between the two points in meters. + >>> from collections import namedtuple >>> point_2d = namedtuple("point_2d", "lat lon") >>> SAN_FRANCISCO = point_2d(37.774856, -122.424227) >>> YOSEMITE = point_2d(37.864742, -119.537521) >>> f"{haversine_distance(*SAN_FRANCISCO, *YOSEMITE):0,.0f} meters" - '254,352 meters' + '254,033 meters' """ - # CONSTANTS per WGS84 https://en.wikipedia.org/wiki/World_Geodetic_System - # Distance in metres(m) - # Equation parameters - # Equation https://en.wikipedia.org/wiki/Haversine_formula#Formulation - flattening = (AXIS_A - AXIS_B) / AXIS_A - phi_1 = atan((1 - flattening) * tan(radians(lat1))) - phi_2 = atan((1 - flattening) * tan(radians(lat2))) - lambda_1 = radians(lon1) - lambda_2 = radians(lon2) - # Equation - sin_sq_phi = sin((phi_2 - phi_1) / 2) - sin_sq_lambda = sin((lambda_2 - lambda_1) / 2) - # Square both values - sin_sq_phi *= sin_sq_phi - sin_sq_lambda *= sin_sq_lambda - h_value = sqrt(sin_sq_phi + (cos(phi_1) * cos(phi_2) * sin_sq_lambda)) - return 2 * RADIUS * asin(h_value) + + radius = 6378137 # earth radius (meters) + + lat1_rad = radians(lat1) + lat2_rad = radians(lat2) + delta_lat = radians(lat2 - lat1) + delta_lon = radians(lon2 - lon1) + + # Haversine formula + a = ( + sin(delta_lat / 2) ** 2 + + cos(lat1_rad) * cos(lat2_rad) * sin(delta_lon / 2) ** 2 + ) + c = 2 * asin(sqrt(a)) + + # Great-Circle Distance + return radius * c if __name__ == "__main__": diff --git a/geodesy/lamberts_ellipsoidal_distance.py b/geodesy/lamberts_ellipsoidal_distance.py index 4805674e51ab..79ce62ed7a09 100644 --- a/geodesy/lamberts_ellipsoidal_distance.py +++ b/geodesy/lamberts_ellipsoidal_distance.py @@ -15,22 +15,21 @@ def lamberts_ellipsoidal_distance( two points on the surface of earth given longitudes and latitudes https://en.wikipedia.org/wiki/Geographical_distance#Lambert's_formula_for_long_lines - NOTE: This algorithm uses geodesy/haversine_distance.py to compute central angle, - sigma - + NOTE: This algorithm uses geodesy/haversine_distance.py to compute the + central angle, sigma Representing the earth as an ellipsoid allows us to approximate distances between points on the surface much better than a sphere. Ellipsoidal formulas treat the Earth as an oblate ellipsoid which means accounting for the flattening that happens - at the North and South poles. Lambert's formulae provide accuracy on the order of - 10 meteres over thousands of kilometeres. Other methods can provide - millimeter-level accuracy but this is a simpler method to calculate long range + at the North and South poles. Lambert's formulas provide accuracy on the order of + 10 meters over thousands of kilometers. Other methods can provide + millimeter-level accuracy, but this is a simpler method to calculate long-range distances without increasing computational intensity. Args: lat1, lon1: latitude and longitude of coordinate 1 lat2, lon2: latitude and longitude of coordinate 2 Returns: - geographical distance between two points in metres + geographical distance between two points in meters >>> from collections import namedtuple >>> point_2d = namedtuple("point_2d", "lat lon") @@ -39,25 +38,20 @@ def lamberts_ellipsoidal_distance( >>> NEW_YORK = point_2d(40.713019, -74.012647) >>> VENICE = point_2d(45.443012, 12.313071) >>> f"{lamberts_ellipsoidal_distance(*SAN_FRANCISCO, *YOSEMITE):0,.0f} meters" - '254,351 meters' + '254,032 meters' >>> f"{lamberts_ellipsoidal_distance(*SAN_FRANCISCO, *NEW_YORK):0,.0f} meters" - '4,138,992 meters' + '4,133,295 meters' >>> f"{lamberts_ellipsoidal_distance(*SAN_FRANCISCO, *VENICE):0,.0f} meters" - '9,737,326 meters' + '9,719,525 meters' """ # CONSTANTS per WGS84 https://en.wikipedia.org/wiki/World_Geodetic_System - # Distance in metres(m) - # Equation Parameters - # https://en.wikipedia.org/wiki/Geographical_distance#Lambert's_formula_for_long_lines flattening = (AXIS_A - AXIS_B) / AXIS_A # Parametric latitudes - # https://en.wikipedia.org/wiki/Latitude#Parametric_(or_reduced)_latitude b_lat1 = atan((1 - flattening) * tan(radians(lat1))) b_lat2 = atan((1 - flattening) * tan(radians(lat2))) - # Compute central angle between two points - # using haversine theta. sigma = haversine_distance / equatorial radius + # Compute the central angle between two points using the haversine function sigma = haversine_distance(lat1, lon1, lat2, lon2) / EQUATORIAL_RADIUS # Intermediate P and Q values @@ -65,13 +59,11 @@ def lamberts_ellipsoidal_distance( q_value = (b_lat2 - b_lat1) / 2 # Intermediate X value - # X = (sigma - sin(sigma)) * sin^2Pcos^2Q / cos^2(sigma/2) x_numerator = (sin(p_value) ** 2) * (cos(q_value) ** 2) - x_demonimator = cos(sigma / 2) ** 2 - x_value = (sigma - sin(sigma)) * (x_numerator / x_demonimator) + x_denominator = cos(sigma / 2) ** 2 + x_value = (sigma - sin(sigma)) * (x_numerator / x_denominator) # Intermediate Y value - # Y = (sigma + sin(sigma)) * cos^2Psin^2Q / sin^2(sigma/2) y_numerator = (cos(p_value) ** 2) * (sin(q_value) ** 2) y_denominator = sin(sigma / 2) ** 2 y_value = (sigma + sin(sigma)) * (y_numerator / y_denominator) diff --git a/geodesy/temp_code_runner_file.py b/geodesy/temp_code_runner_file.py new file mode 100644 index 000000000000..289e85db5e43 --- /dev/null +++ b/geodesy/temp_code_runner_file.py @@ -0,0 +1,3 @@ +# if __name__ == "__main__": +# import doctest +# doctest.testmod()