Add maplabel.py and refernce from pyproj post (+ fix zones link).

diff --git a/posts/geoscript/maplabel.py b/posts/geoscript/maplabel.py
new file mode 100644 (file)
index 0000000..39b5d69
--- /dev/null
+++ b/posts/geoscript/maplabel.py
@@ -0,0 +1,91 @@
+#!/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/>.
+
+"""Label points on a 2D map.
+
+>>> points = [
+...     'New York, NY',
+...     'Miami, FL',
+...     'New Orleans, LA',
+...     'San Antonio, TX',
+...     'Los Angeles, CA',
+...     'San Fransisco, California',
+...     'Seattle, WA',
+...     'Fargo, ND',
+...     'Chicago, IL',
+...     'Bangor, ME',
+...     ]
+>>> rpoints = list(resolve_points(points))
+>>> plot_points(rpoints)
+>>> _pylab.show()
+"""
+
+import geopy as _geopy
+import pylab as _pylab
+import pyproj as _pyproj
+
+
+__version__ = '0.1'
+
+
+def resolve_points(points, geocoder=_geopy.geocoders.Google()):
+    "Resolve points into labels and lat/lng coords."
+    for p in points:
+        label,(lat,lng) = geocoder.geocode(p)
+        yield (label,(lat,lng))
+
+def plot_points(points, proj=None):
+    """Label points 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 labels on the current pylab axis.
+    """
+    if proj is None:
+        mlng = sum([lng for label,(lat,lng) in points])/len(points)
+        zone = int(round((mlng + 180) / 6.))
+        proj = _pyproj.Proj(proj='utm', zone=zone, ellps='WGS84')
+    for label,(lat,lng) in points:
+        x,y = proj(lng, lat)
+        _pylab.annotate(label, xy=(x, y), xycoords='data',
+                        textcoords='offset points', xytext=(6,0),
+                        verticalalignment='center')
+        _pylab.plot([x,x], [y,y], '.')
+
+
+if __name__ == '__main__':
+    from argparse import ArgumentParser
+
+    p = ArgumentParser(
+        description=__doc__.splitlines()[0],)
+    p.add_argument(
+        '--version', action='version', version='%(prog)s '+__version__)
+    p.add_argument('placefiles', metavar='PLACEFILE', type=str, nargs='+',
+                   help='Place name per line file')
+
+    args = p.parse_args()
+
+    _pylab.hold(True)
+
+    for filename in args.coordfiles:
+        points = list(resolve_points(open(filename, 'r').readlines()))
+        plot_points(points)
+
+    _pylab.axis('equal')
+    _pylab.show()

diff --git a/posts/pyproj.mdwn b/posts/pyproj.mdwn
index d669dba1c10bb8ef506f4e21b262fcbb452d8fce..a533fd1f67769428a6f2ed7967f97ccb137cec86 100644 (file)
--- a/posts/pyproj.mdwn
+++ b/posts/pyproj.mdwn
@@ -128,7 +128,8 @@ lat/lng pairs stored in a text file:
     >>> show()
 
 I've written up a simple script using this approach:
-[[geoscript/maproute.py]].
+[[geoscript/maproute.py]].  I've also written up a simple script to
+draw a map with labeled points: [[geoscript/maplabel.py]].
 
 Note that you can easily get lat/lng pairs using [geopy][] (ebuild in
 my [[Gentoo overlay]]):
@@ -156,7 +157,7 @@ conversions, [GeographicLib][] looks promising.
 [awips226]: http://www.nco.ncep.noaa.gov/pmb/docs/on388/tableb.html#GRID218
 [lcc-param]: http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_2sp.html
 [UTM]: http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system
-[zone]: http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#UTM_zone
+[zones]: http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#UTM_zone
 [geopy]: http://code.google.com/p/geopy/
 [GeographicLib]: http://geographiclib.sourceforge.net/html/