Add support for Rhumb and RhumbLine

MANIFEST.in

@@ -1,2 +1,3 @@
 include LICENSE.txt
 include README.md
+include cython_extern_src/*.cpp

README.md

@@ -5,6 +5,7 @@
 Pythonic [`geographiclib` package](https://pypi.python.org/pypi/geographiclib) reimplements geodesic math in Python. Using Cython bindings for C++ version of `geographiclib` allows for **60-90x** speedup for geodesic functions.
 
 **Warning:** only `Inverse()`, `Direct()`, and `InverseLine()` for `Geodesic.WGS84` are implemented. For `GeodesicLine`, only `Position()` is implemented.
+For `Rhumb.WGS84`, `Direct()`, `Inverse()`, `Line()`, and `InverseLine()` are implemented. For `RhumbLine`, `Position()` is implemented, as well as an extra method `Positions()` as an efficient workaround to Cython limitations.
 
 ```bash
 # for Ubuntu
@@ -59,10 +60,13 @@ pip install cython
 
 There are two Cython files: `cgeographiclib.pxd` describing the C++ API of the `libGeographic` library, and `geographiclib_cython.pyx` describing classes that will be visible from Python user code. 
 The `.pyx` imports `.pxd` to learn about C++ classes and functions available to be called.
+An additional file `cython_extern_src/rhumb_line_positions.cpp` contains C++ code that is directly included in the `.pyx` file.
 
 We wrap C++ classes `Geodesic` and `GeodesicLine` into Cython classes `PyGeodesic` and `PyGeodesicLine`. 
 Additionally, we define a pure Python class `Geodesic` with a single field `WGS84` to mimic the behavior of the official `geographiclib` package.
 
+Similarly, we wrap C++ classes `Rhumb` and `RhumbLine` into Cython classes `PyRhumb` and `PyRhumbLine`. (Note: Due to limitations in Cython, the `PyRhumbLine` class is not able to directly wrap the C++ `RhumbLine`, so a workaround is used.) We also define a pure Python class `Rhumb` with a single field `WGS84` to mimic the behavior of the official `geographiclib` package.
+
 There is also `setup.py` file. 
 This file describes how to build the extension module, using `distutils`.
 In there, we specify the library to link with as `libraries=['Geographic']`. The `Geographic` stands for `libGeographic.so` that we previously installed.

cgeographiclib.pxd

@@ -14,4 +14,15 @@ cdef extern from "GeographicLib/GeodesicLine.hpp" namespace "GeographicLib":
         double Distance()
         double Position(double s12, double &lat2, double &lon2)
 
+cdef extern from "GeographicLib/Rhumb.hpp" namespace "GeographicLib":
+    cdef cppclass Rhumb:
+        Rhumb(double a, double f)  except +
+        double Direct(double lat1, double lon1, double azi12, double s12, double &lat2, double &lon2)
+        double Inverse(double lat1, double lon1, double lat2, double lon2, double &s12, double &azi12)
+        RhumbLine Line(double lat1, double lon1, double azi12)
+        @staticmethod
+        Rhumb& WGS84()
 
+cdef extern from "GeographicLib/Rhumb.hpp" namespace "GeographicLib":
+    cdef cppclass RhumbLine:
+        double Position(double s12, double &lat2, double &lon2)

cython_extern_src/rhumb_line_positions.cpp

@@ -0,0 +1,9 @@
+#include <GeographicLib/Rhumb.hpp>
+
+static void rhumb_line_positions(const GeographicLib::Rhumb *rhumb, double lat1, double lon1, double azi12, double s12[], size_t num_positions, double lat2[], double lon2[] ) {
+    GeographicLib::RhumbLine line = rhumb->Line(lat1, lon1, azi12);
+
+    for (size_t i = 0; i < num_positions; i++) {
+        line.Position(s12[i], lat2[i], lon2[i]);
+    }
+}

geographiclib_cython.pyx

@@ -2,6 +2,9 @@
 
 cimport cgeographiclib
 from libcpp cimport bool as cbool
+from cpython cimport array
+import array
+from collections.abc import Iterable
 
 class Geodesic:
     WGS84 = PyGeodesic()
@@ -64,3 +67,114 @@ cdef class PyGeodesicLine:
         return {
             'lat2': lat2, 'lon2': lon2, 's12': s12
         }
+
+class Rhumb:
+    WGS84 = PyRhumb()
+
+cdef class PyRhumb:
+
+    cdef const cgeographiclib.Rhumb* rhumb
+    cdef cbool needs_delete
+
+    def __cinit__(self, double radius=0.0):
+        if radius == 0.0:
+            self.rhumb = &(cgeographiclib.Rhumb.WGS84())
+            self.needs_delete = False
+        else:
+            self.rhumb = new cgeographiclib.Rhumb(radius, 0.0)
+            self.needs_delete = True
+
+    def __dealloc__(self):
+        if self.needs_delete:
+            del self.rhumb
+
+    cpdef Direct(self, double lat1, double lon1, double azi12, double s12):
+        cdef double lat2 = 0.0, lon2 = 0.0
+        self.rhumb.Direct(lat1, lon1, azi12, s12, lat2, lon2)
+        return {
+             'lat1': lat1, 'lon1': lon1,
+             'lat2': lat2, 'lon2': lon2,
+             's12': s12, 'azi12': azi12,
+        }
+
+    cpdef Inverse(self, double lat1, double lon1, double lat2, double lon2):
+        cdef double s12 = 0.0, azi12 = 0.0
+        self.rhumb.Inverse(lat1, lon1, lat2, lon2, s12, azi12)
+        return {
+            'lat1': lat1, 'lon1': lon1,
+            'lat2': lat2, 'lon2': lon2,
+            's12': s12, 'azi12': azi12
+        }
+
+    cpdef Line(self, double lat1, double lon1, double azi12):
+        return PyRhumbLine(self, lat1, lon1, azi12)
+
+    cpdef InverseLine(self, double lat1, double lon1, double lat2, double lon2):
+        inverse = self.Inverse(lat1, lon1, lat2, lon2)
+        cdef PyRhumbLine line = self.Line(lat1, lon1, inverse['azi12'])
+        return line
+
+    cpdef _line_position(self, double lat1, double lon1, double azi12, double s12):
+        cdef double lat2 = 0.0, lon2 = 0.0
+        self.rhumb.Line(lat1, lon1, azi12).Position(s12, lat2, lon2)
+        return {
+            'lat2': lat2, 'lon2': lon2, 's12': s12
+        }
+
+    cpdef _line_positions(self, double lat1, double lon1, double azi12, array.array[double] s12):
+
+        cdef int num_positions = len(s12)
+        cdef array.array lat2 = array.clone(s12, num_positions, zero=False)
+        cdef array.array lon2 = array.clone(s12, num_positions, zero=False)
+
+        cdef double *s12_ptr = s12.data.as_doubles
+        cdef double *lat2_ptr = lat2.data.as_doubles
+        cdef double *lon2_ptr = lon2.data.as_doubles
+
+        rhumb_line_positions(self.rhumb, lat1, lon1, azi12, s12_ptr , num_positions, lat2_ptr, lon2_ptr)
+
+        output = [
+            {
+                'lat2': lat2[i],
+                'lon2': lon2[i],
+                's12': s12[i]
+            }
+            for i in range(num_positions)
+        ]
+
+        return output
+
+cdef class PyRhumbLine:
+    """
+    Because Cython cannot stack allocate classes lacking a nullary constructor,
+    and the RhumbLine class is not copy-assignable, this wrapper must construct
+    a new RhumbLine instance with each call to a position method.
+    """
+
+    cdef readonly PyRhumb rhumb
+    cdef readonly double lat1
+    cdef readonly double lon1
+    cdef readonly double azi12
+
+
+    def __cinit__(self, PyRhumb rhumb, double lat1, double lon1, double azi12):
+        self.rhumb = rhumb
+        self.lat1 = lat1
+        self.lon1 = lon1
+        self.azi12 = azi12
+
+    cpdef Position(self, double s12):
+        return self.rhumb._line_position(self.lat1, self.lon1, self.azi12, s12)
+
+    cpdef Positions(self, s12: Iterable[float]):
+        """
+        Because we can't hold on to a RhumbLine instance in Cython code, this method is provided
+        to calculate multiple positions from a given line, only instantiating it once.
+        """
+        s12_array = array.array('d', s12)
+        return self.rhumb._line_positions(self.lat1, self.lon1, self.azi12, s12_array)
+
+cdef extern from "cython_extern_src/rhumb_line_positions.cpp":
+    # Source C++ definitions from an external file.
+    # Note that this file must be covered by a MANIFEST.in entry to be included in the source distribution by setuptools!
+    void rhumb_line_positions(const cgeographiclib.Rhumb *rhumb, double lat1, double lon1, double azi12, double s12[], size_t num_positions, double lat2[], double lon2[])

setup.py

@@ -11,8 +11,10 @@
 extensions = [Extension(
     "geographiclib_cython",
     ["geographiclib_cython"+file_ext],
-    libraries=["Geographic"]
-)]
+    libraries=["Geographic"],
+    extra_compile_args=["-std=c++11"]
+)
+]
 
 if USE_CYTHON:
     from Cython.Build import cythonize
@@ -21,7 +23,7 @@
 setup(
     name="geographiclib-cython-bindings",
     description="""Cython extension module for C++ geographiclib functions""",
-    version="0.2.0",
+    version="0.3.0",
     author="Sergey Serebryakov",
     author_email="[email protected]",
     url="https://github.com/megaserg/geographiclib-cython-bindings",

tests.py

@@ -1,7 +1,9 @@
 import unittest
 import pytest
+import timeit
 from geographiclib.geodesic import Geodesic
 from geographiclib_cython import Geodesic as CythonGeodesic
+from geographiclib_cython import Rhumb, PyRhumbLine
 from geopy.distance import great_circle
 
 # Run with: python -m pytest tests.py
@@ -39,3 +41,122 @@ def test_sphere_distance(self):
         expected = great_circle((10, 20), (30, 40))
 
         assert actual['s12'] == pytest.approx(expected.meters, 1e-10)
+
+
+class TestRhumb(unittest.TestCase):
+    # Test case from the manual for RhumbSolve
+    LatLonJFK = 40 + 38/60 + 23/3600, -1*(73 + 46/60 + 44/3600)
+    LatLonSIN = 1 + 21/60 + 33/3600, 103 + 59/60 + 22/3600
+    azi_JFK_SIN = 103.5828330034
+    s_JFK_SIN = 18523563.04238
+
+    def test_direct(self):
+        lat1, lon1 = self.LatLonJFK
+        azi12 = self.azi_JFK_SIN
+        s12 = self.s_JFK_SIN
+
+        expected = {
+            'lat2': self.LatLonSIN[0],
+            'lon2': self.LatLonSIN[1]
+        }
+
+        actual = Rhumb.WGS84.Direct(lat1, lon1, azi12, s12)
+
+        assert actual['lat2'] == pytest.approx(expected['lat2'], 1e-10)
+        assert actual['lon2'] == pytest.approx(expected['lon2'], 1e-10)
+
+    def test_inverse(self):
+        lat1, lon1 = self.LatLonJFK
+        lat2, lon2 = self.LatLonSIN
+
+        expected = {
+            'azi12': self.azi_JFK_SIN,
+            's12': self.s_JFK_SIN
+        }
+
+        actual = Rhumb.WGS84.Inverse(lat1, lon1, lat2, lon2)
+
+        assert actual['s12'] == pytest.approx(expected['s12'], 1e-10)
+        assert actual['azi12'] == pytest.approx(expected['azi12'], 1e-10)
+
+    def test_line_position(self):
+        lat1, lon1 = self.LatLonJFK
+        azi12 = self.azi_JFK_SIN
+        s12 = self.s_JFK_SIN
+
+        line = Rhumb.WGS84.Line(lat1, lon1, azi12)
+
+        expected = {
+            'lat2': self.LatLonSIN[0],
+            'lon2': self.LatLonSIN[1]
+        }
+
+        actual = line.Position(s12)
+
+        assert actual['lat2'] == pytest.approx(expected['lat2'], 1e-10)
+        assert actual['lon2'] == pytest.approx(expected['lon2'], 1e-10)
+
+    def test_line_positions(self):
+        lat1, lon1 = self.LatLonJFK
+        azi12 = self.azi_JFK_SIN
+        s12 = self.s_JFK_SIN
+
+        line = Rhumb.WGS84.Line(lat1, lon1, azi12)
+
+        distances = [0, 1*s12/4, 2*s12/4, 3*s12/4, s12]
+        num_distances = len(distances)
+
+        actual = line.Positions(distances)
+
+        assert len(actual) == num_distances
+
+        assert actual[0]['lat2'] == pytest.approx(self.LatLonJFK[0], 1e-10)
+        assert actual[0]['lon2'] == pytest.approx(self.LatLonJFK[1], 1e-10)
+
+        assert actual[num_distances-1]['lat2'] == pytest.approx(self.LatLonSIN[0], 1e-10)
+        assert actual[num_distances-1]['lon2'] == pytest.approx(self.LatLonSIN[1], 1e-10)
+
+    def test_line_positions_type(self):
+        lat1, lon1 = self.LatLonJFK
+        azi12 = self.azi_JFK_SIN
+
+        line = Rhumb.WGS84.Line(lat1, lon1, azi12)
+
+        with pytest.raises(TypeError):
+            line.Positions(2)
+
+    def test_line_positions_methods(self):
+        lat1, lon1 = self.LatLonJFK
+        azi12 = self.azi_JFK_SIN
+        s12 = self.s_JFK_SIN
+
+        line = Rhumb.WGS84.Line(lat1, lon1, azi12)
+        num_distances = 1000
+        distances = [s12*(x/num_distances) for x in range(num_distances)]
+
+        def compute_single():
+            return [line.Position(distance) for distance in distances]
+
+        def compute_multiple():
+            return line.Positions(distances)
+
+        # Compare equality of results
+        output_single = compute_single()
+        output_multiple = compute_multiple()
+        assert output_single == output_multiple
+
+        # Compare timing
+        num_repetitions = 1000
+        time_single = timeit.timeit(compute_single, number=num_repetitions)
+        time_multiple = timeit.timeit(compute_multiple, number=num_repetitions)
+
+        assert time_multiple < 0.95*time_single  # make sure our optimization is actually helping!
+
+    def test_inverse_line(self):
+        lat1, lon1 = self.LatLonJFK
+        lat2, lon2 = self.LatLonSIN
+
+        inverse_line = Rhumb.WGS84.InverseLine(lat1, lon1, lat2, lon2)
+
+        assert inverse_line.azi12 == pytest.approx(self.azi_JFK_SIN, 1e-10)
+