mirror of
https://github.com/TheAlgorithms/Python.git
synced 2024-10-06 05:39:30 +00:00
c909da9b08
* pre-commit: Upgrade psf/black for stable style 2023 Updating https://github.com/psf/black ... updating 22.12.0 -> 23.1.0 for their `2023 stable style`. * https://github.com/psf/black/blob/main/CHANGES.md#2310 > This is the first [psf/black] release of 2023, and following our stability policy, it comes with a number of improvements to our stable style… Also, add https://github.com/tox-dev/pyproject-fmt and https://github.com/abravalheri/validate-pyproject to pre-commit. I only modified `.pre-commit-config.yaml` and all other files were modified by pre-commit.ci and psf/black. * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
86 lines
3.4 KiB
Python
86 lines
3.4 KiB
Python
from math import atan, cos, radians, sin, tan
|
|
|
|
from .haversine_distance import haversine_distance
|
|
|
|
AXIS_A = 6378137.0
|
|
AXIS_B = 6356752.314245
|
|
EQUATORIAL_RADIUS = 6378137
|
|
|
|
|
|
def lamberts_ellipsoidal_distance(
|
|
lat1: float, lon1: float, lat2: float, lon2: float
|
|
) -> float:
|
|
"""
|
|
Calculate the shortest distance along the surface of an ellipsoid between
|
|
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
|
|
|
|
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
|
|
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
|
|
|
|
>>> 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)
|
|
>>> 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'
|
|
>>> f"{lamberts_ellipsoidal_distance(*SAN_FRANCISCO, *NEW_YORK):0,.0f} meters"
|
|
'4,138,992 meters'
|
|
>>> f"{lamberts_ellipsoidal_distance(*SAN_FRANCISCO, *VENICE):0,.0f} meters"
|
|
'9,737,326 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
|
|
sigma = haversine_distance(lat1, lon1, lat2, lon2) / EQUATORIAL_RADIUS
|
|
|
|
# Intermediate P and Q values
|
|
p_value = (b_lat1 + b_lat2) / 2
|
|
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)
|
|
|
|
# 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)
|
|
|
|
return EQUATORIAL_RADIUS * (sigma - ((flattening / 2) * (x_value + y_value)))
|
|
|
|
|
|
if __name__ == "__main__":
|
|
import doctest
|
|
|
|
doctest.testmod()
|