17
by Adam Candy
Import of Varun's branch lp:~varun-verma11/+junk/gshhstoshp into the folder 'shape'. |
1 |
"""
|
2 |
Sets the read and the write file stream according to
|
|
3 |
the command line arguments given.
|
|
4 |
||
5 |
The first argument specifies which shape file the user
|
|
6 |
wants to specify the boundaries of
|
|
7 |
||
8 |
The second arguments specifies the boundary polygon
|
|
9 |
||
10 |
The third argument specifies the file path to which the
|
|
11 |
new shape has to be written
|
|
12 |
"""
|
|
13 |
||
14 |
||
15 |
from shapely.ops import cascaded_union |
|
16 |
import shapefile |
|
17 |
from shapely.geometry import * |
|
18 |
import sys |
|
19 |
import matplotlib.pyplot as pyplot |
|
20 |
||
21 |
"""
|
|
22 |
Used if there is only one shape in the boundary shapefile.
|
|
23 |
"""
|
|
24 |
def getBoundaryPointsList(bounds): |
|
25 |
shapes = bounds.shapes() |
|
26 |
pointsList = [] |
|
27 |
for i in range(len(shapes)): |
|
28 |
pointsList.append(shapes[i].points) |
|
29 |
polygons = [] |
|
30 |
for j in range(len(pointsList)): |
|
31 |
polygons.append(Polygon([pointsList[j][i] for i in range(len(pointsList[j]))])) |
|
32 |
return (polygons,pointsList) |
|
33 |
||
34 |
"""
|
|
35 |
When more than one shapefile, works out which objects from the .shp file are overlapping and return the exterior coords of this new shape.
|
|
36 |
"""
|
|
37 |
def overlap (bounds, plot = False): |
|
38 |
# Read shapefile and work out overlap. Optional plot.
|
|
39 |
shapes = bounds.shapes() |
|
40 |
pointsList = [] |
|
41 |
for i in range(len(shapes)): |
|
42 |
# Check datapoint is valid.
|
|
43 |
pointsList.append(shapes[i].points) |
|
44 |
||
45 |
# Turn the points into polygons.
|
|
46 |
polygons = [] |
|
47 |
for j in range(len(pointsList)): |
|
48 |
polygons.append(Polygon([pointsList[j][i] for i in range(len(pointsList[j]))])) |
|
49 |
||
50 |
# Add overlapping shapes into a list so we know which to join together.
|
|
51 |
overlapping = [] |
|
52 |
for n in range(len(polygons) - 1): |
|
53 |
if polygons[n].intersects(polygons[n+1]) == True: |
|
54 |
# Two if statements to make sure the same polygon isn't being entered more than once.
|
|
55 |
if polygons[n] not in overlapping: overlapping.append(polygons[n]) |
|
56 |
if polygons[n + 1] not in overlapping: overlapping.append(polygons[n + 1]) |
|
57 |
||
58 |
# Create a new shape from the overlapping shapes.
|
|
59 |
join = cascaded_union(overlapping) |
|
60 |
poly = [join] |
|
61 |
||
62 |
# Take the coords. of the perimeter of this new shape.
|
|
63 |
coords = [] |
|
64 |
for i in range(len(join.exterior.coords)): |
|
65 |
coords.append(list(join.exterior.coords[i])) |
|
66 |
||
67 |
||
68 |
# Plot results if True. Store x-y coords of the perimeter in two lists to plot.
|
|
69 |
if plot == True: |
|
70 |
x = []; y = [] |
|
71 |
for i in range(len(coords)): |
|
72 |
x.append(coords[i][0]); y.append(coords[i][1]) |
|
73 |
||
74 |
# Plot results.
|
|
75 |
pyplot.plot(x, y) |
|
76 |
pyplot.xlim(-4, 4) |
|
77 |
pyplot.ylim(-4, 4) |
|
78 |
pyplot.show() |
|
79 |
return join, coords |
|
80 |
||
81 |
||
82 |
||
83 |
"""
|
|
84 |
Output the final results as a .geo file.
|
|
85 |
"""
|
|
86 |
def write_geo (coords, filename): |
|
87 |
# Write new shape to .geo.
|
|
88 |
print coords |
|
89 |
||
90 |
target = open("%s.geo" % filename, "w") # Creates .geo file to write to. |
|
91 |
for i in range(len(coords)): |
|
92 |
# Write point.
|
|
93 |
target.write('Point(%d) = {%.3f, %.3f, 0, %.3f};\n' %(i + 1, coords[i][0], coords[i][1], 1.0)) |
|
94 |
# Write the lines connecting the sequential points.
|
|
95 |
if (i + 1 > 1): target.write('Line(%d) = {%d, %d};\n' % (i, i, i + 1)) |
|
96 |
# Connect first and last points.
|
|
97 |
target.write('Line(%d) = {%d, %d};\n' % (i + 1, 1, i + 1)) |
|
98 |
target.close() |
|
99 |
return False |
|
100 |
||
101 |
||
102 |
assert len(sys.argv)==4, "Incorrect Number of Arguments passed" |
|
103 |
||
104 |
readPath = sys.argv[1] |
|
105 |
boundaryPath = sys.argv[2] |
|
106 |
writePath = sys.argv[3] |
|
107 |
||
108 |
#input stream for the given shape
|
|
109 |
sf = shapefile.Reader(readPath) |
|
110 |
||
111 |
#shapes contained in the given file
|
|
112 |
shapes = sf.shapes(); |
|
113 |
||
114 |
#boundary = bounds.shapes[0]
|
|
115 |
boundary = shapefile.Reader(boundaryPath) |
|
116 |
||
117 |
print(len(boundary.shapes())) |
|
118 |
if (len(boundary.shapes()) > 1): boundaryPolygons, boundaryPointList1 = overlap(boundary); boundaryPointList = [boundaryPointList1] |
|
119 |
else: boundaryPolygons, boundaryPointList = getBoundaryPointsList(boundary) |
|
120 |
print len(shapes[0].points) |
|
121 |
||
122 |
"""
|
|
123 |
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 the exterior cooords are taken and turned into Shapely Points and checked to see which components of the shape lie within the boundary then plots result.
|
|
124 |
"""
|
|
125 |
for shape in shapes: |
|
126 |
x = []; y = [] |
|
127 |
polygon = Polygon([shape.points[i] for i in range(len(shape.points))]) |
|
128 |
if (polygon.intersects(boundaryPolygons)): |
|
129 |
for i in range(len(polygon.exterior.coords)): |
|
130 |
point = Point(polygon.exterior.coords[i]) |
|
131 |
if point.intersects(boundaryPolygons): |
|
132 |
x.append(shape.points[i][0]); y.append(shape.points[i][1]) |
|
133 |
print shape.points[i] |
|
134 |
pyplot.plot(x, y) |
|
135 |
"""
|
|
136 |
Plot boundary.
|
|
137 |
"""
|
|
138 |
x=[] |
|
139 |
y=[] |
|
140 |
for p in boundaryPointList : |
|
141 |
for p2 in p : |
|
142 |
x.append(p2[0]) |
|
143 |
y.append(p2[1]) |
|
144 |
pyplot.plot(x,y) |
|
145 |
pyplot.xlim(min(x)-1,max(x)+1) |
|
146 |
pyplot.ylim(min(y)-1,max(y)+1) |
|
147 |
pyplot.show() |