Skip to content

Created geodesy section with one algorithm #1757

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 7 commits into from
Feb 16, 2020
Merged

Created geodesy section with one algorithm #1757

Changes from 2 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
36c1792
implemented haversine
eightysixth Feb 16, 2020
395a806
updated docstring
eightysixth Feb 16, 2020
e440ed2
added type hints
eightysixth Feb 16, 2020
3dfd133
added type hints
eightysixth Feb 16, 2020
e037a9f
Merge branch 'master' of https://github.com/eightysixth/Python
eightysixth Feb 16, 2020
73e8401
improved docstring and math usage
eightysixth Feb 16, 2020
b5d09f3
f"{haversine_distance(*SAN_FRANCISCO, *YOSEMITE):0,.0f} meters"
cclauss Feb 16, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
48 changes: 48 additions & 0 deletions geodesy/haversine_distance.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,48 @@
from math import asin, atan, cos, sin, sqrt, tan, pow, radians


def haversine_distance(lat1, lon1, lat2, lon2):
"""
Calculate great circle distance between two points in a sphere,
given longitudes and latitudes.
(https://en.wikipedia.org/wiki/Haversine_formula)

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

>>> int(haversine_distance(37.774856, -122.424227, 37.864742, -119.537521)) # From SF to Yosemite
254352
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is the unit of measure in great circle calculations? Google maps says 164 miles or 264 kilometers but that is road distance, not as-the-crow-flies distance. Helping readers connect what they are learning with what they already know will improve the value of this submission.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm planning on adding two more algorithms which improve approximation so haversine is a good starter method. Perhaps, a readme file can follow with comparisons and detailed explanations. For now:

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.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is the unit of measure? Is it meters, kilometres, miles, lightyears?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Line 15 - metres.. Clarified that in the latest commit as well.

unit of measure meant is irrelevant for explanation.


"""

# CONSTANTS per WGS84 https://en.wikipedia.org/wiki/World_Geodetic_System
AXIS_A = 6378137.0
AXIS_B = 6356752.314245
RADIUS = 6378137

# 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 = pow(sin((phi_2 - phi_1) / 2), 2)
sin_sq_lambda = pow(sin((lambda_2 - lambda_1) / 2), 2)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Using math.pow() is very slow compared to pythons ** notation for powers.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What about sqrt vs ** .5?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

According to this stackoverflow thread you should use math.sqrt.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And as a side note, we can simplify the above lines even more: a * a is a lot faster than a ** 2.

So you can change the above lines to (sin(phi_2 - phi_1) / 2) * (sin(phi_2 - phi_1) / 2) etc

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Did that and made it a little more readable

h_value = sqrt(sin_sq_phi + (cos(phi_1) * cos(phi_2) * sin_sq_lambda))

distance = 2 * RADIUS * asin(h_value)

return distance


if __name__ == "__main__":
import doctest

doctest.testmod()