-
Notifications
You must be signed in to change notification settings - Fork 3
/
simplify.py
342 lines (269 loc) · 12.6 KB
/
simplify.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
#! /usr/bin/env python
# encoding: utf-8
__author__ = 'asimmons'
## This script is to created to simplify
## shapefile geometry using the
## Visvalingam algorithm found here
## http://www2.dcs.hull.ac.uk/CISRG/publications/DPs/DP10/DP10.html
## Threshold is the area of the largest allowed triangle
import fiona
from shapely.geometry import shape, mapping, Polygon, MultiPolygon, LineString, MultiLineString
import heapq
import shapely
from shapely.geometry.polygon import LinearRing
import sys
class TriangleCalculator:
def __init__(self, point, index):
# Need to add better validation
# Save instance variables
self.point = point
self.ringIndex = index
self.prevTriangle = None
self.nextTriangle = None
# enables the instantiation of 'TriangleCalculator' to be compared
# by the calcArea().
def __cmp__(self, other):
return cmp(self.calcArea(), other.calcArea())
## calculate the effective area of a triangle given
## its vertices -- using the cross product
def calcArea(self):
# Add validation
if not self.prevTriangle or not self.nextTriangle:
print "ERROR:"
p1 = self.point
p2 = self.prevTriangle.point
p3 = self.nextTriangle.point
area = abs(p1[0] * (p2[1] - p3[1]) + p2[0] * (p3[1] - p1[1]) + p3[0] * (p1[1] - p2[1])) / 2.0
#print "area = " + str(area) + ", point = " + str(self.point)
return area
class GeomSimplify:
def simplify_line(self, line, threshold):
# unlike rings: we need to keep beginning and end points static throughout the simplification process
# Build list of Triangles from the line points
triangleArray = []
## each triangle contains an index and a point (x,y)
# handle line 'interior' (i.e. the vertices
# between start and end) first -- explicitly
# defined using the below slice notation
# i.e. [1:-1]
for index, point in enumerate(line.coords[1:-1]):
triangleArray.append(TriangleCalculator(point, index))
# then create start/end points separate from the triangleArray (meaning
# we cannot have the start/end points included in the heap sort)
startIndex = 0
endIndex = len(line.coords)-1
startTriangle = TriangleCalculator(line.coords[startIndex], startIndex)
endTriangle = TriangleCalculator(line.coords[endIndex], endIndex)
# Hook up triangles with next and prev references (doubly-linked list)
# NOTE: linked list are composed of nodes, which have at
# least one link to another node (and this is a doubly-linked list..pointing at
# both our prevTriangle & our nextTriangle)
# NOTE: in this code block the 'triangle' is our 'triangle node'
for index, triangle in enumerate(triangleArray):
# set prevIndex to be the adjacent point to index
prevIndex = index - 1
nextIndex = index + 1
if prevIndex >= 0:
triangle.prevTriangle = triangleArray[prevIndex]
else:
triangle.prevTriangle = startTriangle
if nextIndex < len(triangleArray):
triangle.nextTriangle = triangleArray[nextIndex]
else:
triangle.nextTriangle = endTriangle
# Build a min-heap from the TriangleCalculator list
# print "heapify"
heapq.heapify(triangleArray)
# Simplify steps...
# Note: in contrast
# to our function 'simplify_ring'
# we can allow our array to go down to 0 and STILL have a valid line
# because we will still have the start and end points
while len(triangleArray) > 0:
# if the smallest triangle is greater than the threshold, we can stop
# i.e. loop to point where the heap head is >= threshold
if triangleArray[0].calcArea() >= threshold:
#print "break"
break
else:
# print statement for debugging - prints area's and coords of deleted/simplified pts
#print "simplify...triangle area's and their corresponding points that were less then the threshold"
#print "area = " + str(triangleArray[0].calcArea()) + ", point = " + str(triangleArray[0].point)
prev = triangleArray[0].prevTriangle
next = triangleArray[0].nextTriangle
prev.nextTriangle = next
next.prevTriangle = prev
# This has to be done after updating the linked list
# in order for the areas to be correct when the
# heap re-sorts
# print "popping (i.e. re-measuring area & comparing)"
heapq.heappop(triangleArray)
#print "area = " + str(triangle.calcArea()) + ", point = " + str(triangle.point)
#print "done popping (i.e. area that is less than threshold, and will have point removed)"
# Create an list of indices from the triangleRing heap
indexList = []
for triangle in triangleArray:
# add 1 b/c the triangle array's first index is actually the second point
indexList.append(triangle.ringIndex + 1)
# Append start and end points back into the array
indexList.append(startTriangle.ringIndex)
indexList.append(endTriangle.ringIndex)
# Sort the index list
indexList.sort()
# Create a new simplified ring
simpleLine = []
for index in indexList:
simpleLine.append(line.coords[index])
# Convert list into LineString
simpleLine = LineString(simpleLine)
# print statements for debugging to check if points are being reduced...
#print "Starting size (incl. beginning/end point): " + str(len(line.coords))
#print "Ending size (incl. beginning/end point): " + str(len(simpleLine.coords))
#print "Starting Coord: " + str(line.coords[startIndex])
#print "End Coord: " + str(line.coords[endIndex])
#print list(simpleLine.coords)
return simpleLine
def simplify_ring(self, ring, threshold):
# Build list of TriangleCalculators
triangleRing = []
## each triangle contains an index and a point (x,y)
## because rings have a point on top of a point
## we are skipping the last point by using slice notation[:-1]
## *i.e. 'a[:-1]' # everything except the last item*
for index, point in enumerate(ring.coords[:-1]):
triangleRing.append(TriangleCalculator(point, index))
# Hook up triangles with next and prev references (doubly-linked list)
for index, triangle in enumerate(triangleRing):
# set prevIndex to be the adjacent point to index
# these steps are necessary for dealing with
# closed rings
prevIndex = index - 1
if prevIndex < 0:
# if prevIndex is less than 0, then it means index = 0, and
# the prevIndex is set to last value in the index
# (i.e. adjacent to index[0])
prevIndex = len(triangleRing) - 1
# set nextIndex adjacent to index
nextIndex = index + 1
if nextIndex == len(triangleRing):
# if nextIndex is equivalent to the length of the array
# set nextIndex to 0
nextIndex = 0
triangle.prevTriangle = triangleRing[prevIndex]
triangle.nextTriangle = triangleRing[nextIndex]
# Build a min-heap from the TriangleCalculator list
heapq.heapify(triangleRing)
# Simplify
while len(triangleRing) > 2:
# if the smallest triangle is greater than the threshold, we can stop
# i.e. loop to point where the heap head is >= threshold
if triangleRing[0].calcArea() >= threshold:
break
else:
prev = triangleRing[0].prevTriangle
next = triangleRing[0].nextTriangle
prev.nextTriangle = next
next.prevTriangle = prev
# This has to be done after updating the linked list
# in order for the areas to be correct when the
# heap re-sorts
heapq.heappop(triangleRing)
# Handle case where we've removed too many points for the ring to be a polygon
if len(triangleRing) < 3:
return None
# Create an list of indices from the triangleRing heap
indexList = []
for triangle in triangleRing:
indexList.append(triangle.ringIndex)
# Sort the index list
indexList.sort()
# Create a new simplified ring
simpleRing = []
for index in indexList:
simpleRing.append(ring.coords[index])
# Convert list into LinearRing
simpleRing = LinearRing(simpleRing)
# print statements for debugging to check if points are being reduced...
#print "Starting size: " + str(len(ring.coords))
#print "Ending size: " + str(len(simpleRing.coords))
return simpleRing
def simplify_multipolygon(self, mpoly, threshold):
# break multipolygon into polys
polyList = mpoly.geoms
simplePolyList = []
# call simplify_polygon() on each
for poly in polyList:
simplePoly = self.simplify_polygon(poly, threshold)
#if not none append to list
if simplePoly:
simplePolyList.append(simplePoly)
# check that polygon count > 0, otherwise return None
if not simplePolyList:
return None
# put back into multipolygon
return MultiPolygon(simplePolyList)
def simplify_polygon(self, poly, threshold):
# Get exterior ring
simpleExtRing = self.simplify_ring(poly.exterior, threshold)
# If the exterior ring was removed by simplification, return None
if simpleExtRing is None:
return None
simpleIntRings = []
for ring in poly.interiors:
simpleRing = self.simplify_ring(ring, threshold)
if simpleRing is not None:
simpleIntRings.append(simpleRing)
return shapely.geometry.Polygon(simpleExtRing, simpleIntRings)
def simplify_multiline(self, mline, threshold):
# break MultiLineString into lines
lineList = mline.geoms
simpleLineList = []
# call simplify_line on each
for line in lineList:
simpleLine = self.simplify_line(line, threshold)
#if not none append to list
if simpleLine:
simpleLineList.append(simpleLine)
# check that line count > 0, otherwise return None
if not simpleLineList:
return None
# put back into multilinestring
return MultiLineString(simpleLineList)
def process_file(self, inFile, outFile, threshold):
with fiona.open(inFile, 'r') as input:
meta = input.meta
# The outFile has the same crs, schema as inFile
with fiona.open(outFile, 'w', **meta) as output:
# Read shapely geometries from file
# Loop through all shapely objects
for myGeom in input:
myShape = shape(myGeom['geometry'])
if isinstance(myShape, Polygon):
myShape = self.simplify_polygon(myShape, threshold)
elif isinstance(myShape, MultiPolygon):
myShape = self.simplify_multipolygon(myShape, threshold)
elif isinstance(myShape, LineString):
myShape = self.simplify_line(myShape, threshold)
elif isinstance(myShape, MultiLineString):
myShape = self.simplify_multiline(myShape, threshold)
else:
raise ValueError('Unhandled geometry type: ' + repr(myShape.type))
# write to outfile
if myShape is not None:
output.write({'geometry':mapping(myShape), 'properties': myGeom['properties']})
def main():
print "number of arguments (incl. py file name): " + str(len(sys.argv))
if len(sys.argv) != 4:
print "Wrong amount of arguments!"
usage()
exit()
inputFile = sys.argv[1]
outputFile = sys.argv[2]
threshold = sys.argv[3]
geomSimplifyObject = GeomSimplify()
geomSimplifyObject.process_file(inputFile, outputFile, float(threshold))
print "Finished simplifying file!"
def usage():
print "python simplify.py <input file> <output file> <threshold>"
if __name__ == "__main__":
main()