import fiona
import shapely.geometry as geometry
input_shapefile = 'data/concave_demo_points.shp'
shapefile = fiona.open(input_shapefile)
points = [geometry.shape(point['geometry']) for point in shapefile]
print("We have {0} points!".format(len(points)))
We have 143 points!
We just loaded a "shapefile." Shapefiles are one of the most commonly used formats for saving geometries. Developed by ESRI, and open.
We can use pylab now to plot those points.
%matplotlib inline
import pylab as pl
# The following proves that I don't know numpy in the slightest
x = [p.coords.xy[0] for p in points]
y = [p.coords.xy[1] for p in points]
pl.figure(figsize=(10,10))
_ = pl.plot(x,y,'o', color='#f16824')
But now we have points. I used QGIS to draw these points in the rough shape of a capital "H". I highly recommend that you brew install qgis-22
(or equivalent) if you want to play with geometric features.
Now that you have a collection of points, you can interrogate the collection. Many shapely operations result in a different kind of geometry than the one you're working with. Since our geometry is a collection of points, I can instantiate a MultiPoint, and then ask that MultiPoint for its envelope, which is a Polygon.
point_collection = geometry.MultiPoint(list(points))
point_collection.envelope
<shapely.geometry.polygon.Polygon at 0x7f3734463b00>
We should take a look at that envelope. matplotlib can help us out, but polygons aren't functions, so we need to use PolygonPatch.
from descartes import PolygonPatch
def plot_polygon(polygon):
fig = pl.figure(figsize=(10,10))
ax = fig.add_subplot(111)
margin = .3
x_min, y_min, x_max, y_max = polygon.bounds
ax.set_xlim([x_min-margin, x_max+margin])
ax.set_ylim([y_min-margin, y_max+margin])
patch = PolygonPatch(polygon, fc='#999999', ec='#000000', fill=True, zorder=-1)
ax.add_patch(patch)
return fig
_ = plot_polygon(point_collection.envelope)
_ = pl.plot(x,y,'o', color='#f16824')
Shapely enables tons of geometric set theoretic operations in python. It might be interesting to compute the convex hull, for instance.
convex_hull_polygon = point_collection.convex_hull
_ = plot_polygon(convex_hull_polygon)
_ = pl.plot(x,y,'o', color='#f16824')
Convex Hull does a neat thing for us. It draws a polygon around the original points, but without any margin (contrast with envelope.)
It doesn't look like an H though. What if we want to derive a polygon that's more representative of the density of these points?
From Wikipedia:
In computational geometry, an alpha shape, or α-shape, is a family of piecewise linear simple curves in the Euclidean plane associated with the shape of a finite set of points. They were first defined by Edelsbrunner, Kirkpatrick & Seidel (1983). The alpha-shape associated with a set of points is a generalization of the concept of the convex hull, i.e. every convex hull is an alpha-shape but not every alpha shape is a convex hull.
But shapely doesn't have an implementation of concave hull.
Sean Gillies has some alpha shapes code. The gist is, compute the Delaunay triangles for the points, and then only add the triangles that meet a certain heuristic.
from shapely.ops import cascaded_union, polygonize
from scipy.spatial import Delaunay
import numpy as np
import math
def alpha_shape(points, alpha):
"""
Compute the alpha shape (concave hull) of a set of points.
@param points: Iterable container of points.
@param alpha: alpha value to influence the gooeyness of the border. Smaller
numbers don't fall inward as much as larger numbers. Too large,
and you lose everything!
"""
if len(points) < 4:
# When you have a triangle, there is no sense in computing an alpha
# shape.
return geometry.MultiPoint(list(points)).convex_hull
def add_edge(edges, edge_points, coords, i, j):
"""Add a line between the i-th and j-th points, if not in the list already"""
if (i, j) in edges or (j, i) in edges:
# already added
return
edges.add( (i, j) )
edge_points.append(coords[ [i, j] ])
coords = np.array([point.coords[0] for point in points])
tri = Delaunay(coords)
edges = set()
edge_points = []
# loop over triangles:
# ia, ib, ic = indices of corner points of the triangle
for ia, ib, ic in tri.vertices:
pa = coords[ia]
pb = coords[ib]
pc = coords[ic]
# Lengths of sides of triangle
a = math.sqrt((pa[0]-pb[0])**2 + (pa[1]-pb[1])**2)
b = math.sqrt((pb[0]-pc[0])**2 + (pb[1]-pc[1])**2)
c = math.sqrt((pc[0]-pa[0])**2 + (pc[1]-pa[1])**2)
# Semiperimeter of triangle
s = (a + b + c)/2.0
# Area of triangle by Heron's formula
area = math.sqrt(s*(s-a)*(s-b)*(s-c))
circum_r = a*b*c/(4.0*area)
# Here's the radius filter.
#print circum_r
if circum_r < 1.0/alpha:
add_edge(edges, edge_points, coords, ia, ib)
add_edge(edges, edge_points, coords, ib, ic)
add_edge(edges, edge_points, coords, ic, ia)
m = geometry.MultiLineString(edge_points)
triangles = list(polygonize(m))
return cascaded_union(triangles), edge_points
concave_hull, edge_points = alpha_shape(points, alpha=1.87)
_ = plot_polygon(concave_hull)
_ = pl.plot(x,y,'o', color='#f16824')
It turns out that the alpha value and the scale of the points matters a lot when it comes to how well the Delaunay triangulation method will work. You can usually play with the alpha value to find a suitable response. I'll do both: scale up the "H" and try some different alpha values.
To get more points, I opened up QGIS, drew an "H" like polygon, generated regular points within the bounds of the "H", and then spatially joined them to remove any points outside the "H".
input_shapefile = 'data/demo_poly_scaled_points_join.shp'
new_shapefile = fiona.open(input_shapefile)
new_points = [geometry.shape(point['geometry']) for point in new_shapefile]
print("We have {0} points!".format(len(new_points)))
x = [p.coords.xy[0] for p in new_points]
y = [p.coords.xy[1] for p in new_points]
pl.figure(figsize=(10,10))
_ = pl.plot(x,y,'o', color='#f16824')
We have 497 points!
from matplotlib.collections import LineCollection
for i in range(9):
alpha = (i+1)*.1
concave_hull, edge_points = alpha_shape(new_points, alpha=alpha)
#print concave_hull
lines = LineCollection(edge_points)
pl.figure(figsize=(10,10))
pl.title('Alpha={0} Delaunay triangulation'.format(alpha))
pl.gca().add_collection(lines)
delaunay_points = np.array([point.coords[0] for point in new_points])
pl.plot(delaunay_points[:,0], delaunay_points[:,1], 'o', color='#f16824')
_ = plot_polygon(concave_hull)
_ = pl.plot(x,y,'o', color='#f16824')
So in this case, alpha of about 0.4 looks pretty good. We can use shapely's buffer operation to clean up that polygon a bit and smooth out any of the jagged edges.
alpha = .4
concave_hull, edge_points = alpha_shape(new_points, alpha=alpha)
plot_polygon(concave_hull)
_ = pl.plot(x,y,'o', color='#f16824')
plot_polygon(concave_hull.buffer(1))
_ = pl.plot(x,y,'o', color='#f16824')