2
Sets the read and the write file stream according to
3
the command line arguments given.
5
The first argument specifies which shape file the user
6
wants to specify the boundaries of
8
The second arguments specifies the boundary polygon
10
The third argument specifies the file path to which the
11
new shape has to be written
15
from shapely.ops import cascaded_union
17
from shapely.geometry import *
19
import matplotlib.pyplot as pyplot
22
Write the results to .shp.
24
def writeShapeFile(points, filepath) :
25
#begin the instance pf writer class
26
w = shapefile.Writer()
27
#ensure shape and records the balance
34
w.point(pList[0][0],pList[0][1])
35
w.field("%d_FLD"%i,"C","40")
39
w.field("%d_FLD"%i,"C","40")
43
w.field("%d_FLD"%i,"C","40")
46
print("Number of shapes Written %d" %i)
50
Used if there is only one shape in the boundary shapefile.
52
def getBoundaryPointsList(bounds):
53
shapes = bounds.shapes()
54
pointsList = (shapes[0].points)
55
polygon = Polygon(pointsList)
56
return (polygon,pointsList)
59
When more than one shapefile, works out which objects from the .shp file are overlapping and return the exterior coords of this new shape.
61
def overlap (bounds, plot = False):
62
# Read shapefile and work out overlap. Optional plot.
63
shapes = bounds.shapes()
65
for i in range(len(shapes)):
66
# Check datapoint is valid.
67
pointsList.append(shapes[i].points)
69
# Turn the points into polygons.
71
for j in range(len(pointsList)):
72
polygons.append(Polygon([pointsList[j][i] for i in range(len(pointsList[j]))]))
74
# Add overlapping shapes into a list so we know which to join together.
76
for n in range(len(polygons) - 1):
77
if polygons[n].intersects(polygons[n+1]) == True:
78
# Two if statements to make sure the same polygon isn't being entered more than once.
79
if polygons[n] not in overlapping: overlapping.append(polygons[n])
80
if polygons[n + 1] not in overlapping: overlapping.append(polygons[n + 1])
82
# Create a new shape from the overlapping shapes.
83
join = cascaded_union(overlapping)
86
# Take the coords. of the perimeter of this new shape.
88
for i in range(len(join.exterior.coords)):
89
coords.append(list(join.exterior.coords[i]))
92
# Plot results if True. Store x-y coords of the perimeter in two lists to plot.
95
for i in range(len(coords)):
96
x.append(coords[i][0]); y.append(coords[i][1])
108
Output the final results as a .geo file.
110
def write_geo (coords, filename):
111
# Write new shape to .geo.
114
target = open("%s.geo" % filename, "w") # Creates .geo file to write to.
115
for i in range(len(coords)):
117
target.write('Point(%d) = {%.3f, %.3f, 0, %.3f};\n' %(i + 1, coords[i][0], coords[i][1], 1.0))
118
# Write the lines connecting the sequential points.
119
if (i + 1 > 1): target.write('Line(%d) = {%d, %d};\n' % (i, i, i + 1))
120
# Connect first and last points.
121
target.write('Line(%d) = {%d, %d};\n' % (i + 1, 1, i + 1))
126
assert len(sys.argv)==5, "Incorrect Number of Arguments passed"
128
readPath = sys.argv[1]
129
boundaryPath = sys.argv[2]
130
writePath = sys.argv[3]
131
areaThreshold = float(sys.argv[4])
133
#input stream for the given shape
134
sf = shapefile.Reader(readPath)
136
#shapes contained in the given file
137
shapes = sf.shapes();
139
#boundary = bounds.shapes[0]
140
boundary = shapefile.Reader(boundaryPath)
142
if (len(boundary.shapes()) > 1): boundaryPolygons, boundaryPointList1 = overlap(boundary); boundaryPointList = [boundaryPointList1]
143
else: boundaryPolygons, boundaryPointList = getBoundaryPointsList(boundary)
146
Takes shape from shapefile and converts to a Shapely Polygon. Checks if this polygon lies within the boundary using a.intersect(b). If it does it will perform a.intersection(b) operation returning a Polygon/MultiPolygon which lies within the boundary and then plots result.
150
x = []; y = []; shp = []
151
polygon = Polygon([shape.points[i] for i in range(len(shape.points))])
152
if (polygon.intersects(boundaryPolygons)):
153
intersection = boundaryPolygons.intersection(polygon)
154
if intersection.area >= areaThreshold:
156
if intersection.geom_type == 'Polygon':
157
for i in range(len(list(intersection.exterior.coords))):
158
x.append(intersection.exterior.coords[i][0]); y.append(intersection.exterior.coords[i][1])
160
shp.append([intersection.exterior.coords[i][0], intersection.exterior.coords[i][1]])
162
if intersection.geom_type == 'MultiPolygon':
163
for j in range(len(intersection)):
164
for i in range(len(list(intersection[j].exterior.coords))):
165
x.append(intersection[j].exterior.coords[i][0]); y.append(intersection[j].exterior.coords[i][1])
167
shp.append([intersection[j].exterior.coords[i][0], intersection[j].exterior.coords[i][1]])
168
shapeList.append(shp)
169
writeShapeFile(shapeList, writePath)
175
for i in range(len(boundaryPointList)):
176
x.append(boundaryPointList[i][0]); y.append(boundaryPointList[i][1])
178
pyplot.xlim(min(x)-1,max(x)+1)
179
pyplot.ylim(min(y)-1,max(y)+1)