Merge pull request #39 from CrapsJeroen/inverse-haversine-function

Added inver haversine functionality to the module

README.md

@@ -47,6 +47,30 @@ outputs
  <Unit.METERS: 'm'>, <Unit.MILES: 'mi'>, <Unit.NAUTICAL_MILES: 'nmi'>)
 ```
 
+### Inverse Haversine Formula
+Calculates a point from a given vector (distance and direction) and start point.
+Currently explicitly supports both cardinal (north, east, south, west) and intercardinal (northeast, southeast, southwest, northwest) directions.
+But also allows for explicit angles expressed in Radians.
+
+## Example: Finding arbitary point from Paris
+```python
+from haversine import inverse_haversine, Direction
+from math import pi
+paris = (48.8567, 2.3508) # (lat, lon)
+# Finding 32 km west of Paris
+inverse_haversine(paris, 32, Direction.WEST)
+# returns tuple (49.1444, 2.3508)
+# Finding 32 km southwest of Paris
+inverse_haversine(paris, 32, pi * 1.25)
+# returns tuple (48.5377, 1.8705)
+# Finding 50 miles north of Paris
+inverse_haversine(paris, 50, Direction.NORTH, unit=Unit.MILES)
+# returns tuple (49.5803, 2.3508)
+# Finding 10 nautical miles south of Paris
+inverse_haversine(paris, 10, Direction.SOUTH, unit=Unit.NAUTICAL_MILES)
+# returns tuple (48.6901, 2.3508)
+```
+
 ### Performance optimisation for distances between all points in two vectors
 
 You will need to add [numpy](https://pypi.org/project/numpy/) in order to gain performance with vectors.

haversine/__init__.py

@@ -1 +1 @@
-from .haversine import Unit, haversine, haversine_vector
+from .haversine import Unit, haversine, haversine_vector, Direction, inverse_haversine

haversine/haversine.py

@@ -1,5 +1,6 @@
-from math import radians, cos, sin, asin, sqrt
+from math import radians, cos, sin, asin, sqrt, degrees, pi, atan2
 from enum import Enum
+from typing import Union
 
 
 # mean earth radius - https://en.wikipedia.org/wiki/Earth_radius#Mean_radius
@@ -21,18 +22,40 @@ class Unit(Enum):
     INCHES = 'in'
 
 
+class Direction(Enum):
+    """
+    Enumeration of supported directions.
+    The full list can be checked by iterating over the class; e.g.
+    the expression `tuple(Direction)`.
+    Angles expressed in radians.
+    """
+
+    NORTH = 0
+    NORTHEAST = pi * 0.25
+    EAST = pi * 0.5
+    SOUTHEAST = pi * 0.75
+    SOUTH = pi
+    SOUTHWEST = pi * 1.25
+    WEST = pi * 1.5
+    NORTHWEST = pi * 1.75
+
+
 # Unit values taken from http://www.unitconversion.org/unit_converter/length.html
-_CONVERSIONS = {Unit.KILOMETERS:       1.0,
-                Unit.METERS:           1000.0,
-                Unit.MILES:            0.621371192,
-                Unit.NAUTICAL_MILES:   0.539956803,
-                Unit.FEET:             3280.839895013,
-                Unit.INCHES:           39370.078740158}
+_CONVERSIONS = {
+    Unit.KILOMETERS:       1.0,
+    Unit.METERS:           1000.0,
+    Unit.MILES:            0.621371192,
+    Unit.NAUTICAL_MILES:   0.539956803,
+    Unit.FEET:             3280.839895013,
+    Unit.INCHES:           39370.078740158
+}
+
 
 def get_avg_earth_radius(unit):
     unit = Unit(unit)
     return _AVG_EARTH_RADIUS_KM * _CONVERSIONS[unit]
 
+
 def haversine(point1, point2, unit=Unit.KILOMETERS):
     """ Calculate the great-circle distance between two points on the Earth surface.
 
@@ -118,15 +141,30 @@ def haversine_vector(array1, array2, unit=Unit.KILOMETERS, comb=False):
 
     # If in combination mode, turn coordinates of array1 into column vectors for broadcasting
     if comb:
-        lat1 = numpy.expand_dims(lat1,axis=0)
-        lng1 = numpy.expand_dims(lng1,axis=0)
-        lat2 = numpy.expand_dims(lat2,axis=1)
-        lng2 = numpy.expand_dims(lng2,axis=1)
-    
+        lat1 = numpy.expand_dims(lat1, axis=0)
+        lng1 = numpy.expand_dims(lng1, axis=0)
+        lat2 = numpy.expand_dims(lat2, axis=1)
+        lng2 = numpy.expand_dims(lng2, axis=1)
+
     # calculate haversine
     lat = lat2 - lat1
     lng = lng2 - lng1
     d = (numpy.sin(lat * 0.5) ** 2
          + numpy.cos(lat1) * numpy.cos(lat2) * numpy.sin(lng * 0.5) ** 2)
 
     return 2 * get_avg_earth_radius(unit) * numpy.arcsin(numpy.sqrt(d))
+
+
+def inverse_haversine(point, distance, direction: Union[Direction, float], unit=Unit.KILOMETERS):
+
+    lat, lng = point
+    lat, lng = map(radians, (lat, lng))
+    d = distance
+    r = get_avg_earth_radius(unit)
+    brng = direction.value if isinstance(direction, Direction) else direction
+
+    return_lat = asin(sin(lat) * cos(d / r) + cos(lat) * sin(d / r) * cos(brng))
+    return_lng = lng + atan2(sin(brng) * sin(d / r) * cos(lat), cos(d / r) - sin(lat) * sin(return_lat))
+
+    return_lat, return_lng = map(degrees, (return_lat, return_lng))
+    return return_lat, return_lng

tests/test_inverse_haversine.py

@@ -0,0 +1,47 @@
+from haversine import inverse_haversine, haversine, Unit, Direction
+from numpy import isclose
+from math import pi
+import pytest
+
+LYON = (45.7597, 4.8422)
+PARIS = (48.8567, 2.3508)
+LONDON = (51.509865, -0.118092)
+NEW_YORK = (40.7033962, -74.2351462)
+
+
+@pytest.mark.parametrize(
+    "point, dir, dist, result",
+    [
+        (PARIS, Direction.NORTH, 32, (49.144444, 2.3508)),
+        (PARIS, 0, 32, (49.144444, 2.3508)),
+        (LONDON, Direction.WEST, 50, (51.507778, -0.840556)),
+        (LONDON, pi * 1.5, 50, (51.507778, -0.840556)),
+        (NEW_YORK, Direction.SOUTH, 15, (40.568611, -74.235278)),
+        (NEW_YORK, Direction.NORTHWEST, 50, (41.020556, -74.656667)),
+        (NEW_YORK, pi * 1.25, 50, (40.384722, -74.6525)),
+    ],
+)
+def test_inverse_kilometers(point, dir, dist, result):
+    assert isclose(inverse_haversine(point, dist, dir), result, rtol=1e-5).all()
+
+
+@pytest.mark.parametrize(
+    "point, direction, distance, unit",
+    [
+        (PARIS, Direction.NORTH, 10, Unit.KILOMETERS),
+        (LONDON, Direction.WEST, 32, Unit.MILES),
+        (LYON, Direction.NORTHEAST, 45_000, Unit.METERS),
+        (NEW_YORK, Direction.SOUTH, 15, Unit.NAUTICAL_MILES),
+    ],
+)
+def test_back_and_forth(point, direction, distance, unit):
+    new_point = inverse_haversine(point, distance, direction, unit)
+    assert isclose(haversine(new_point, point, unit), distance, rtol=1e-10)
+
+
+def test_inverse_miles():
+    assert isclose(inverse_haversine(PARIS, 50, Direction.NORTH, unit=Unit.MILES), (49.5803579218996, 2.3508), rtol=1e-5).all()
+
+
+def test_nautical_inverse_miles():
+    assert isclose(inverse_haversine(PARIS, 10, Direction.SOUTH, unit=Unit.NAUTICAL_MILES), (48.69014586638915, 2.3508), rtol=1e-5).all()