from IPython.core.display import Image Image(filename='files/conway.png') from IPython.core.display import Image Image(filename='files/GA.png') Image(filename='files/luna.jpg') Image(filename='files/luna_real.JPG') import os import numpy as np from matplotlib.patches import Polygon import matplotlib.pyplot as plt import pysal import pylab as pl import shapely import shapely.geometry from shapely.geometry import shape from shapely.geometry import Point, LineString, Polygon, MultiPolygon from shapely.wkt import dumps, loads import random #grab a random color for coloring a figure def get_random_color(): r = lambda: random.randint(0,255) return('#%02X%02X%02X' % (r(),r(),r())) #initialize the figure and draw the shape def plot(shapelyGeometries, color_dict={"fill":"#AADDCC", "line":"#666666","hole_fill":"#ffffff", "hole_line":"#999999" }): 'Plot shapelyGeometries' figure = pl.figure(num=None, figsize=(10, 10), dpi=180) axes = pl.axes() axes.set_aspect('equal', 'datalim') axes.xaxis.set_visible(False) axes.yaxis.set_visible(False) draw(shapelyGeometries, color_dict) #Check the type and break up multipolygons def draw(gs, color_dict): 'Draw shapelyGeometries' # Handle single and lists of geometries try: gs = iter(gs) except TypeError: gs = [gs] #Route polygons and multipolygons to the right place for g in gs: gType = g.geom_type if gType.startswith('Multi') or gType == 'GeometryCollection': draw(g.geoms, color_dict) else: draw_(g, color_dict) #Break the shape into its interior and exterior rings def draw_(g, color_dict): 'Draw a shapelyGeometry; thanks to Sean Gilles' gType = g.geom_type if gType == 'Point': pl.plot(g.x, g.y, 'k,') elif gType == 'LineString': x, y = g.xy pl.plot(x, y, 'b-', color=color_dict["line"]) elif gType == 'Polygon': #can draw parts as multiple colors if not color_dict: color_dict={"fill":get_random_color(), "line":"#666666", "hole_fill":"#FFFFFF", "hole_line":"#999999" } x, y = g.exterior.xy pl.fill(x, y, color=color_dict["fill"], aa=True) pl.plot(x, y, color=color_dict["line"], aa=True, lw=1.0) for hole in g.interiors: x, y = hole.xy pl.fill(x, y, color=color_dict["hole_fill"], aa=True) pl.plot(x, y, color=color_dict["hole_line"], aa=True, lw=1.0) pt = Point(0, 0) print pt plot(pt) circle = pt.buffer(.5) plot([circle, pt]) print circle small_circle = pt.buffer(.25) plot([circle, small_circle, pt], color_dict=None) donut = circle.difference(small_circle) plot(donut) lineString1 = LineString(((-3, -5), (3, 5))) print lineString1 plot(lineString1) lineString2 = LineString(((-3, -5), (3, 5), (-3, 5), (3, -5))) print lineString2 plot(lineString2) pt1 = Point(-3, 6.5) pt2 = Point(3, 6.5) left = pt1.buffer(.5) right = pt2.buffer(.5) poly = shapely.geometry.Polygon( ((-3, 5), (3, 5), (0, 0), (-3, 5))) plot([lineString2, left, right, poly], {"fill":"#003300", "line":"#000000","hole_fill":"#000000", "hole_line":"#000000" }) Image(filename='files/emoticon_bunny.jpeg') polygon1 = shapely.geometry.Polygon( ((1, 1), (5, 1), (5, 5), (1, 5), (1, 1)), # Shell [ ((2, 2), (2, 3), (3, 3), (3, 2), (2, 2)), # Hole ((3, 3), (3, 4), (4, 4), (4, 3), (3, 3)), # Hole ] ) print polygon1 plot(polygon1) exterior = polygon1.exterior print "Exterior Ring: " + str(exterior) print "\nInterior Rings: " for hole in polygon1.interiors: print "\t" + str(hole) circle = Point(5,2.5).buffer(2) plot([circle, polygon1], color_dict=None) new_polygon = polygon1.union(circle) plot(new_polygon) print new_polygon.type circle2 = Point(7,7).buffer(2) plot([circle2, polygon1], color_dict=None) new_multipolygon = polygon1.union(circle2) plot(new_multipolygon) print new_multipolygon.type print "CIRCLE AREA: %f" % (circle2.area) print "POLYGON AREA: %f" % (polygon1.area) print "MULTIPOLYGON AREA: %f" % (new_multipolygon.area) exteriors = new_multipolygon.geoms for exterior in exteriors: print "Exterior Ring:" print "\t" + str(exterior) if not exterior.interiors: continue print "\n\tInterior Rings: " for hole in exterior.interiors: print "\t\t" + str(hole) import fiona def get_from_fiona(filename, id_attr): from shapely.geometry import shape data = {} c = fiona.open(filename, 'r') for shape_data in c: name = shape_data['properties'][id_attr] parsed_shape = shape(shape_data['geometry']) data[name] = parsed_shape return data OrderedDict([(u'REGION', u'3'), (u'DIVISION', u'5'), (u'STATEFP', u'12'), (u'STATENS', u'00294478'), (u'GEOID', u'12'), (u'STUSPS', u'FL'), (u'NAME', u'Florida'), (u'LSAD', u'00'), (u'MTFCC', u'G4000'), (u'FUNCSTAT', u'A'), (u'ALAND', 138897453172.0), (u'AWATER', 31413676956.0), (u'INTPTLAT', u'+28.4574302'), (u'INTPTLON', u'-082.4091478')]) OrderedDict([(u'REGION', u'2'), (u'DIVISION', u'3'), (u'STATEFP', u'17'), (u'STATENS', u'01779784'), (u'GEOID', u'17'), (u'STUSPS', u'IL'), (u'NAME', u'Illinois'), (u'LSAD', u'00'), (u'MTFCC', u'G4000'), (u'FUNCSTAT', u'A'), (u'ALAND', 143793994610.0), (u'AWATER', 6201680290.0), (u'INTPTLAT', u'+40.1028754'), (u'INTPTLON', u'-089.1526108')]) from shapely.wkt import dumps, loads state_dict = get_from_fiona("files/tl_2013_us_state/tl_2013_us_state.shp", "NAME") #print(shapely.wkt.dumps(state_dict["Ohio"])) ohio = state_dict["Ohio"] print ohio.convex_hull plot(ohio) plot([state_dict["Pennsylvania"], state_dict["Ohio"], state_dict["West Virginia"], state_dict["Maryland"], state_dict["Kentucky"], state_dict["Indiana"], state_dict["Michigan"]], color_dict=None) print Image(filename='files/MI.png') #This function will help total up a count of all exterior ring points in a polygon or multipolygon def points(s): if s.type == 'Polygon': return len(s.exterior.coords) if s.type == 'Point': return s else: x = 0 for g in s.geoms: if g.type != 'Point' and g.type != 'LineString': x += len(g.exterior.coords) return x import time hawaii = (state_dict["Hawaii"]) now = time.time() hawaii_simple = hawaii.simplify(.1) simple_time = time.time() - now now = time.time() hawaii_simple_not_preserved = hawaii.simplify(.1, preserve_topology=False) simple_not_preserved_time = time.time() - now print "Points with preserving topology " + str((points(hawaii_simple))) print "Time cost: " + str(simple_time) print print "Points without preserving topology " + str((points(hawaii_simple_not_preserved))) print "Time cost: " + str(simple_not_preserved_time) plot(hawaii_simple) plot(hawaii_simple_preserved) texas = state_dict["Texas"] print "TOTAL POINTS: %d" % (points(texas)) plot(state_dict["Texas"]) texas_simple = texas.simplify(.001) plot(texas_simple) print "TOTAL POINTS: %d" % (points(texas_simple)) texas_simplest = texas.simplify(.1) plot(texas_simplest) print "TOTAL POINTS: %d" % (points(texas_simplest)) not_texas = texas.simplify(1.1) plot(not_texas) print "TOTAL POINTS: %d" % (points(not_texas)) plot(texas.convex_hull) plot(texas.envelope) print(texas.bounds) print(texas.centroid) center_circle1 = texas.centroid.buffer(.2) center_circle2 = texas.centroid.buffer(.5) center_circle3 = texas.centroid.buffer(1) plot([texas, center_circle3, center_circle2, center_circle1, texas.centroid], color_dict=None) from shapely.ops import cascaded_union, unary_union buffered_points_list = [] x_ys = range(5); for x in x_ys: buffered_points_list.append(Point(x,x).buffer(1)) union_shape = cascaded_union(buffered_points_list) plot(buffered_points_list, color_dict=None) plot(union_shape) intersection_shape = buffered_points_list[0].intersection(buffered_points_list[1]) plot([buffered_points_list[0], buffered_points_list[1], intersection_shape], color_dict=None) difference_shape1 = buffered_points_list[0].difference(buffered_points_list[1]) difference_shape2 = buffered_points_list[1].difference(buffered_points_list[0]) plot([difference_shape1,difference_shape2], color_dict=None) print "Do the full circle shapes intersect?" print(buffered_points_list[0].intersects(buffered_points_list[1])) print "Do the difference shapes intersect?" print(difference_shape1.intersects(difference_shape2)) print "Are the difference shapes overlapping?" print(difference_shape1.overlaps(difference_shape2)) difference_shape3 = buffered_points_list[1].symmetric_difference(buffered_points_list[0]) plot([difference_shape3]) sad_polygon1 = Polygon(((-3, -5),(-3,0),(-5,0),(-3,0),(-3, 5),(3, 5), (3, -5), (-3, -5))) plot(sad_polygon1) print explain_validity(sad_polygon1) print "SAD POLYGON AREA: %f" % (sad_polygon1.area) fixed_polygon1 = sad_polygon1.buffer(.03) plot(fixed_polygon1) print explain_validity(fixed_polygon1) print "HAPPY POLYGON AREA: %f" % (fixed_polygon1.area) from shapely.validation import explain_validity sad_polygon2 = Polygon(((-3, -5), (3, 5), (-3, 5), (3, -5), (-3, -5))) plot(sad_polygon2) print explain_validity(sad_polygon2) print "The bowtie's area: " + str(sad_polygon2.area) fixed_polygon2 = sad_polygon2.buffer(0) plot(fixed_polygon2) print explain_validity(fixed_polygon2) print fixed_polygon2 print "The area of our fixed shape:" + str(fixed_polygon2.area) from shapely.validation import explain_validity print explain_validity(sad_polygon2) pt = Point(0,0) band_aid = pt.buffer(.2) plot([sad_polygon2, band_aid], color_dict=None) try: trythis = sad_polygon2.union(band_aid) except(shapely.geos.TopologicalError): print "ERROR" --------------------------------------------------------------------------- TopologicalError Traceback (most recent call last) in () 4 band_aid = pt.buffer(.2) 5 plot([sad_polygon2, band_aid]) ----> 6 trythis = sad_polygon2.union(band_aid) 7 8 try: /Users/alisonalvarez/anaconda/lib/python2.7/site-packages/shapely/geometry/base.pyc in union(self, other) 489 def union(self, other): 490 """Returns the union of the geometries (Shapely geometry)""" --> 491 return geom_factory(self.impl['union'](self, other)) 492 493 # Unary predicates /Users/alisonalvarez/anaconda/lib/python2.7/site-packages/shapely/topology.pyc in __call__(self, this, other, *args) 45 if not this.is_valid: 46 raise TopologicalError( ---> 47 "The operation '%s' produced a null geometry. Likely cause is invalidity of the geometry %s" % (self.fn.__name__, repr(this))) 48 elif not other.is_valid: 49 raise TopologicalError( TopologicalError: The operation 'GEOSUnion_r' produced a null geometry. Likely cause is invalidity of the geometry sad_polygon3 = shapely.geometry.Polygon( ((1, 1), (5, 1), (5, 5), (1, 5), (1, 1)), # Shell [ ((2, 2), (2, 3), (3, 3), (3, 2), (2, 2)), # Hole ((6, 6), (6, 7), (7, 7), (7, 6), (6, 6)), # Hole ] ) plot(sad_polygon3) print explain_validity(sad_polygon3) print "SAD POLYGON AREA: %f" % (sad_polygon3.area) fixed_polygon3 = sad_polygon3.buffer(0) plot(fixed_polygon3) print explain_validity(fixed_polygon3) print "FIXED POLYGON AREA: %f" % (fixed_polygon3.area) print fixed_polygon3 sad_polygon4 = shapely.geometry.Polygon( ((1, 1), (5, 1), (5, 5), (1, 5), (1, 1)), # Shell [ ((2, 2), (2, 4), (4, 4), (4, 2), (2, 2)), # Hole ((3, 3), (3, 3.5), (3.5, 3.5), (3.5, 3), (3, 3)), # Hole ] ) plot(sad_polygon4) print explain_validity(sad_polygon4) print "SAD POLYGON AREA: %f" % (sad_polygon4.area) fixed_polygon4 = sad_polygon4.buffer(0) plot(fixed_polygon4) print explain_validity(fixed_polygon4) print "FIXED POLYGON AREA: %f" % (fixed_polygon4.area) print fixed_polygon4 sad_polygon5 = shapely.geometry.Polygon( ((1, 1), (5, 1), (5, 5), (1, 5), (1, 1)), # Shell [ ((2, 2), (2, 3.5), (3, 3.5), (3, 2), (2, 2)), # Hole ((3, 3), (3, 4), (4, 4), (4, 3), (3, 3)), # Hole ] ) print explain_validity(sad_polygon5) print sad_polygon5.is_valid plot(sad_polygon5) print "SAD POLYGON AREA: %f" % (sad_polygon5.area) fixed_polygon5 = sad_polygon5.buffer(0) plot(fixed_polygon5) print explain_validity(fixed_polygon5) print "FIXED POLYGON AREA: %f" % (fixed_polygon5.area) print fixed_polygon5 sad_polygon6 = shapely.geometry.Polygon( ((1, 1), (5, 1), (5, 5), (1, 5), (1, 1)), # Shell [ ((2, 2), (2, 2), (3, 3), (3, 3), (2, 2)), # Hole ] ) plot(sad_polygon6) print explain_validity(sad_polygon6) def remove_bad_shapes(shape): if shape.type == 'Multipolygon': shape_list = shape.geoms elif shape.type == 'Polygon': shape_list = [shape] else: return None polygon_list = [] for g in shape_list: print g print explain_validity(g) hole_list = [] for hole in g.interiors: hole_shape = Polygon(hole) if "Too few points" not in explain_validity(hole_shape): hole_list.append(hole) else: print "Removed hole: " print "\t" + str(hole) if "Too few points" not in explain_validity(Polygon(g.exterior.coords)): print "\t" + explain_validity(Polygon(g.exterior.coords)) polygon_list.append(Polygon(g.exterior.coords, hole_list)) else: print "Removed: " print "\t" + str(g.exterior.coords) return MultiPolygon(polygon_list) fixed_polygon6 = remove_bad_shapes(sad_polygon6) plot(fixed_polygon6) print "\nThe validity of the fixed shape:" print explain_validity(fixed_polygon6) Image(filename='files/zip_union.png') Image(filename='files/separate_zips.png') Image(filename='files/not_cleaned.png') Image(filename='files/cleaned_up.png') radius_earth = 3963.207 # miles def degreesToRadians(d): if (type(d) == type('')): print "ERROR: d is of type string....value = ", d try: d = float(d) except: return 0.0 return (d / 360.0) * 2.0 * math.pi def greatCircleDistance(A, B): (lat1, lon1) = map(degreesToRadians, A) (lat2, lon2) = map(degreesToRadians, B) global radius_earth x1 = math.cos(lat1) * math.cos(lon1) y1 = math.cos(lat1) * math.sin(lon1) z1 = math.sin(lat1) x2 = math.cos(lat2) * math.cos(lon2) y2 = math.cos(lat2) * math.sin(lon2) z2 = math.sin(lat2) C = math.sqrt( ((x1-x2) * (x1-x2)) + ((y1-y2) * (y1-y2)) + ((z1-z2) * (z1-z2)) ) alpha = math.asin(C/2.0) return (2.0 * alpha * radius_earth) #use the bounding box of a shape to set the degree of curve correction def get_sq_mi(shape): (minx, miny, maxx, maxy) = shape.bounds top_left = [maxy, maxx] top_right = [miny, maxx] bottom_left = [maxy, minx] bottom_right = [miny, minx] width = greatCircleDistance(top_left, top_right) height = greatCircleDistance(bottom_right, top_right) return (width * height) #this code will traverse the components of a list of polygons or multipolygons and remove any pieces not meeting the area threshold def simplify_shapes(to_denoise_list, polygon_threshold=0.15, hole_threshold=0.1): denoised_list = [] for shape in to_denoise_list: #handle polgyon vs. multipolygon case if shape.type == 'Polygon': before_shapes = 1 shape_list = [shape] if shape.type == 'MultiPolygon': before_shapes = len(shape.geoms) shape_list = shape.geoms polygon_list = [] for g in shape_list: polygon_sq_mi = get_sq_mi(g) if polygon_sq_mi > polygon_threshold: if len(g.interiors) > 0: hole_list = [] for hole in g.interiors: hole_shape = shapely_Polygon(hole) hole_sq_mi = get_sq_mi(hole_shape) if hole_sq_mi > hole_threshold: hole_list.append(hole) else: g = Polygon(g.exterior.coords, hole_list) polygon_list.append(g) else: #this polygon is below the threshold and should be left out continue denoised_list.append(MultiPolygon(polygon_list)) print "Total polygons before: " + str(before_shapes) print "Total polygons after: " + str(len(polygon_list)) return denoised_list hawaii_list = simplify_shapes([hawaii], polygon_threshold=250, hole_threshold=250) plot(hawaii_list)