Optimize and fix errors in plot interpolation.

Use NumPy functions a bit more and avoid using np.append too much.
The previous code did not ensure that interpolated points were less than
the tolerance apart.

1 files changed, 20 insertions, 17 deletions

diff --git a/docs/plot/plot.py b/docs/plot/plot.py
index 9d388375..52bd7fd2 100644
--- a/docs/plot/plot.py
+++ b/docs/plot/plot.py
@@ -107,27 +107,29 @@ def interp_coords(coords, tol):
     x = np.array(x)
     y = np.array(y)
 
-    l = len(x)
-    dsts = np.sqrt((x[0:-1]-x[1:l])**2 + (y[0:-1]-y[1:l])**2)
-    I = dsts > tol
+    dsts = np.hypot(np.diff(x), np.diff(y))
+    # Points to the *beginning* of the segment.
+    I = np.where(dsts > tol)[0]
 
     offset = 0
-    xy = np.ndarray(shape=(2, 0))
-    for i, b in enumerate(I):
-        xy = np.append(xy, (x[i], y[i]))
-        if not b:
-            continue
+    xy = []
+    for i in I:
+        # Add points that are already below tolerance.
+        xy.append((x[offset:i], y[offset:i]))
+
+        # Interpolate between points above tolerance.
+        n = np.ceil(dsts[i] / tol)
+        x1 = np.linspace(x[i], x[i + 1], n + 1)
+        y1 = np.linspace(y[i], y[i + 1], n + 1)
+        xy.append((x1[:-1], y1[:-1]))
+
+        offset = i + 1
 
-        n = round(dsts[i] / tol/2)
-        x1 = np.linspace(x[i], x[i+1], n+2)
-        y1 = np.linspace(y[i], y[i+1], n+2)
-        xy = np.append(xy, zip(x1[1:], y1[1:]))
+    xy.append((x[offset:], y[offset:]))
 
-    # the last coordinate is not added in the loop above
-    #xy = np.append(xy, (x[i+1], y[i+1])).reshape((-1,2)).tolist()
-    xy = xy.reshape((-1, 2)).tolist()
+    xy = np.concatenate(xy, axis=1)
 
-    return xy
+    return xy.T
 
 
 def project(coordinates, proj_string, in_radians=False):
@@ -149,6 +151,7 @@ def project(coordinates, proj_string, in_radians=False):
     # proj expects binary input to be in radians
     if not in_radians:
         coordinates = np.deg2rad(coordinates)
+    coordinates = coordinates.astype(np.double, copy=False)
 
     # set up cmd call. -b for binary in/out
     args = [PROJ, '-b']
@@ -156,7 +159,7 @@ def project(coordinates, proj_string, in_radians=False):
 
     proc = subprocess.Popen(args, stdin=subprocess.PIPE, stdout=subprocess.PIPE,
                             env={'PROJ_LIB': os.path.abspath(PROJ_LIB)})
-    stdout, _ = proc.communicate(coordinates.tobytes())
+    stdout, _ = proc.communicate(coordinates.tobytes(order='C'))
 
     out = np.frombuffer(stdout, dtype=np.double)
     return np.reshape(out, (-1, 2))