add numba and geographic distance

rdp/__init__.py

@@ -13,6 +13,34 @@
 import numpy as np
 import sys
 
+try:
+    from numba import njit
+    USE_NUMBA = True
+
+except ImportError:
+    # https://stackoverflow.com/questions/3888158
+    import functools
+
+    def optional_arg_decorator(fn):
+        @functools.wraps(fn)
+        def wrapped_decorator(*args, **kwargs):
+            # If no arguments were passed...
+            if len(args) == 1 and len(kwargs) == 0 and callable(args[0]):
+                return fn(args[0])
+            else:
+                def real_decorator(decoratee):
+                    return fn(decoratee, *args, **kwargs)
+                return real_decorator
+        return wrapped_decorator
+
+    @optional_arg_decorator
+    def __noop(func, *args, **kwargs):
+        return(func)
+
+    njit = __noop
+    USE_NUMBA = False
+
+
 if sys.version_info[0] >= 3:
     xrange = range
 
@@ -33,8 +61,60 @@ def pldist(point, start, end):
         return np.linalg.norm(point - start)
 
     return np.divide(
-            np.abs(np.linalg.norm(np.cross(end - start, start - point))),
-            np.linalg.norm(end - start))
+        np.abs(np.linalg.norm(np.cross(end - start, start - point))),
+        np.linalg.norm(end - start))
+
+
+@njit
+def dist_ll_m(pt1, pt2):
+    # distance in meters between two points
+    # distance between two coordinates
+    # includes lat correction to longitute
+    # assumes that lon is in index[0]
+    mlat = (pt1[1] + pt2[1])/2
+    dlat = pt2[1] - pt1[1]
+    dlon = (pt2[0] - pt1[0])*np.cos(mlat*np.pi/180)
+    d2 = dlat**2 + dlon**2
+    return 111320 * sqrt(d2)
+
+
+@njit
+def dist_cart(pt1, pt2):
+    # just distance in cartesian coordinates
+    d = pt1 - pt2
+    return sqrt(d[0]**2 + d[1]**2)
+
+
+@njit
+def dist_to_segment(pt, v1, v2, dist=dist_cart):
+    # distance from point pt to line segment v1, v2
+    lls = (v1[0]-v2[0])**2 + (v1[1]-v2[1])**2
+    vm = (v1+v2)/2
+    chk = (pt[0]-vm[0])**2 + (pt[1]-vm[1])**2
+
+    # approximate far points to center
+    if chk > 4*lls:
+        return dist(pt, vm)
+
+    # not needed with approximation
+    # if (lls == 0):
+    #    return dist_m(pt, v1)
+
+    t = np.dot((pt - v1), (v2 - v1))/lls
+    if t <= 0:
+        pp = v1
+    elif t >= 1:
+        pp = v2
+    else:
+        pp = v1 + t * (v2 - v1)
+    # clip it from 0 to 1
+    return dist(pt, pp)
+
+
+@njit
+def dist_to_segment_lonlat_to_m(pt, v1, v2):
+    # numba doesnt' do partial yet
+    return dist_to_segment(pt, v1, v2, dist=dist_ll_m)
 
 
 def rdp_rec(M, epsilon, dist=pldist):
@@ -48,7 +128,8 @@ def rdp_rec(M, epsilon, dist=pldist):
     :param epsilon: epsilon in the rdp algorithm
     :type epsilon: float
     :param dist: distance function
-    :type dist: function with signature ``f(point, start, end)`` -- see :func:`rdp.pldist`
+    :type dist: function with signature ``f(point, start, end)``
+                -- see :func:`rdp.pldist`
     """
     dmax = 0.0
     index = -1
@@ -69,11 +150,12 @@ def rdp_rec(M, epsilon, dist=pldist):
         return np.vstack((M[0], M[-1]))
 
 
+@njit
 def _rdp_iter(M, start_index, last_index, epsilon, dist=pldist):
     stk = []
     stk.append([start_index, last_index])
     global_start_index = start_index
-    indices = np.ones(last_index - start_index + 1, dtype=bool)
+    indices = np.ones(last_index - start_index + 1, dtype=np.byte)
 
     while stk:
         start_index, last_index = stk.pop()
@@ -95,7 +177,7 @@ def _rdp_iter(M, start_index, last_index, epsilon, dist=pldist):
             for i in xrange(start_index + 1, last_index):
                 indices[i - global_start_index] = False
 
-    return indices
+    return indices > 0
 
 
 def rdp_iter(M, epsilon, dist=pldist, return_mask=False):
@@ -109,7 +191,8 @@ def rdp_iter(M, epsilon, dist=pldist, return_mask=False):
     :param epsilon: epsilon in the rdp algorithm
     :type epsilon: float
     :param dist: distance function
-    :type dist: function with signature ``f(point, start, end)`` -- see :func:`rdp.pldist`
+    :type dist: function with signature ``f(point, start, end)``
+                -- see :func:`rdp.pldist`
     :param return_mask: return the mask of points to keep instead
     :type return_mask: bool
     """
@@ -121,7 +204,12 @@ def rdp_iter(M, epsilon, dist=pldist, return_mask=False):
     return M[mask]
 
 
-def rdp(M, epsilon=0, dist=pldist, algo="iter", return_mask=False):
+def rdp_lonlat_to_m(M, dist=dist_to_segment_lonlat_to_m, **kwargs):
+    return rdp(M, dist=dist, **kwargs)
+
+
+def rdp(M, epsilon=0, dist=dist_to_segment if USE_NUMBA else pldist,
+        algo="iter", return_mask=False):
     """
     Simplifies a given array of points using the Ramer-Douglas-Peucker
     algorithm.
@@ -132,7 +220,7 @@ def rdp(M, epsilon=0, dist=pldist, algo="iter", return_mask=False):
     >>> rdp([[1, 1], [2, 2], [3, 3], [4, 4]])
     [[1, 1], [4, 4]]
 
-    This is a convenience wrapper around both :func:`rdp.rdp_iter` 
+    This is a convenience wrapper around both :func:`rdp.rdp_iter`
     and :func:`rdp.rdp_rec` that detects if the input is a numpy array
     in order to adapt the output accordingly. This means that
     when it is called using a Python list as argument, a Python
@@ -158,12 +246,15 @@ def rdp(M, epsilon=0, dist=pldist, algo="iter", return_mask=False):
            [4, 4]])
 
     :param M: a series of points
-    :type M: numpy array with shape ``(n,d)`` where ``n`` is the number of points and ``d`` their dimension
+    :type M: numpy array with shape ``(n,d)`` where ``n`` is the number of
+             points and ``d`` their dimension
     :param epsilon: epsilon in the rdp algorithm
     :type epsilon: float
     :param dist: distance function
-    :type dist: function with signature ``f(point, start, end)`` -- see :func:`rdp.pldist`
-    :param algo: either ``iter`` for an iterative algorithm or ``rec`` for a recursive algorithm
+    :type dist: function with signature ``f(point, start, end)``
+                -- see :func:`rdp.pldist`
+    :param algo: either ``iter`` for an iterative algorithm
+                 or ``rec`` for a recursive algorithm
     :type algo: string
     :param return_mask: return mask instead of simplified array
     :type return_mask: bool
@@ -173,9 +264,10 @@ def rdp(M, epsilon=0, dist=pldist, algo="iter", return_mask=False):
         algo = partial(rdp_iter, return_mask=return_mask)
     elif algo == "rec":
         if return_mask:
-            raise NotImplementedError("return_mask=True not supported with algo=\"rec\"")
+            raise NotImplementedError(
+                "return_mask=True not supported with algo=\"rec\"")
         algo = rdp_rec
-        
+
     if "numpy" in str(type(M)):
         return algo(M, epsilon, dist)