Create a gist now

Instantly share code, notes, and snippets.

Python code that'll use matplotlib to draw polygons extracted from a GIS shapefile.
#!/usr/bin/python
"""
Uses matplotlib to plot WKT polygons in CSV format.
1. install GDAL
2. download a shapefile
3. use ogr2ogr to convert it to CSV with WKT
ogr2ogr -f CSV output.csv input.shp -lco GEOMETRY=AS_WKT
4. run getPolygons against the output csv
"""
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import matplotlib.path as mpath
from matplotlib.collections import PatchCollection
import csv
import numpy as np
import re
def getPolygons(fn, key):
polys = []
with open(fn) as f:
reader = csv.DictReader(f)
polys = []
areas = {}
for row in reader:
print row.keys()
arean = int(row[key]) #
wkt = row['WKT']
# ignore interior rings
if wkt.startswith("POLYGON"):
wktm = re.match("POLYGON \(\(((?:[0-9.,]|\s)+)\)", wkt)
if wktm:
wktt = wktm.group(1)
wktt = [x.split() for x in wktt.split(",")]
points = np.array([map(float, x) for x in wktt])
points = points.astype(int) # ints seem okay so far
areas[arean] = mpatches.Polygon(points)
elif wkt.startswith("MULTIPOLYGON"):
# OHARE dammit
subparts = re.findall("\(\(((?:[0-9.,]|\s)+)\)", wkt)
points, codes = [], []
for sub in subparts:
wktt = [x.split() for x in sub.split(",")]
# first elt of each polygon's code is a MOVETO, rest are LINETO
cs = [mpath.Path.LINETO] * len(wktt)
cs[0] = mpath.Path.MOVETO
codes += cs
points += [map(float, x) for x in wktt]
ppts = np.array(points).astype(int)
p = mpath.Path(ppts, np.array(codes))
areas[arean] = mpatches.PathPatch(p)
polys = areas.items()
polys.sort()
return PatchCollection([p for n,p in polys])
def main():
import matplotlib.cm as cm
import numpy.random as npr
# this setup code is probably necessary for reasonable use
fig = plt.figure(facecolor="white")
ax= plt.axes()
# pick random values for each CA
values = npr.rand(700)
coll = getPolygons("output3.csv", key='TRACTCE10')
# color the polys using the random values
coll.set_cmap(cm.Reds)
coll.set_array(np.array(values))
coll.set_alpha(0.6) # for pretty
ax.add_collection(coll)
plt.axis('equal')
plt.axis('off')
plt.show()
if __name__ == '__main__':
main()
@jens1245
jens1245 commented Dec 5, 2016

Thank you! One small thing:

Use [-+]? in the regex expressions to allow for minus signs

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment