Add maproute.py script.

diff --git a/posts/geoscript/maproute.py b/posts/geoscript/maproute.py
new file mode 100644 (file)
index 0000000..59422be
--- /dev/null
+++ b/posts/geoscript/maproute.py
@@ -0,0 +1,89 @@
+#!/usr/bin/env python
+#
+# Copyright (C) 2011  W. Trevor King
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+"""
+
+>>> from StringIO import StringIO
+>>> stream = StringIO('\\n'.join([
+...     '#lat  lng',
+...     '41.293595 -73.552322',
+...     '44.998213 -73.258667',
+...     '45.021515 -74.752808',
+...     '44.168809 -76.290894',
+...     '43.463286 -76.433716',
+...     '43.231597 -79.004517',
+...     '42.846168 -78.894654',
+...     '42.255357 -79.762573',
+...     '41.994610 -79.773560',
+...     '42.019101 -75.422974',
+...     '41.478130 -74.994507',
+...     '40.996795 -73.901138',
+...     '']))
+>>> latlngs = read_latlngs(stream)
+>>> print('%.2f m' % route_distance(latlngs))
+1679848.96 m
+>>> plot_route(latlngs)
+>>> _pylab.show()
+"""
+
+import numpy as _numpy
+import pylab as _pylab
+import pyproj as _pyproj
+
+def read_latlngs(stream):
+    "Read in data from space-separated lat/lng lines with '#' comments."
+    return _numpy.loadtxt(stream)
+
+def route_distance(latlngs, geod=_pyproj.Geod(ellps='WGS84')):
+    """Calculate the geodesic distance along an array of lat/lng pairs.
+
+    If `geod` is not given, it defaults to the WGS 84 reference
+    ellipsoide.
+    """
+    az12s,az21s,dists = geod.inv(latlngs[:-1,1], latlngs[:-1,0],
+                                 latlngs[1:,1], latlngs[1:,0])
+    return dists.sum()
+
+def plot_route(latlngs, proj=None):
+    """Plot an array of lat/lng pairs using the projection `proj`.
+
+    If `proj` is not given, it defaults to a projection from the WGS
+    84 reference ellipsoide (for lat/lng) to the closest UTM zone (for
+    x/y).  The closes UTM zone is determined by averaging the input
+    longitudes.
+
+    Plots the line on the current pylab axis.
+    """
+    if proj is None:
+        mlng = latlngs[:,1].mean()
+        zone = int(round((mlng + 180) / 6.))
+        proj = _pyproj.Proj(proj='utm', zone=zone, ellps='WGS84')
+    xs,ys = proj(latlngs[:,1], latlngs[:,0])
+    lines = _pylab.plot(xs, ys, '.-')
+
+if __name__ == '__main__':
+    import sys
+
+    _pylab.hold(True)
+
+    for filename in sys.argv[1:]:
+        latlngs = read_latlngs(filename)
+        dist = route_distance(latlngs)
+        print('%s distance: %g m (%g mi)' % (filename, dist, dist/1.6e3))
+        plot_route(latlngs)
+
+    _pylab.axis('equal')
+    _pylab.show()

diff --git a/posts/pyproj.mdwn b/posts/pyproj.mdwn
index 592baaff33f85069f90b95b8e7318f6956f3dd0f..d669dba1c10bb8ef506f4e21b262fcbb452d8fce 100644 (file)
--- a/posts/pyproj.mdwn
+++ b/posts/pyproj.mdwn
@@ -127,6 +127,9 @@ lat/lng pairs stored in a text file:
     >>> lines = plot(xs, ys, 'r.-')
     >>> show()
 
+I've written up a simple script using this approach:
+[[geoscript/maproute.py]].
+
 Note that you can easily get lat/lng pairs using [geopy][] (ebuild in
 my [[Gentoo overlay]]):