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
|
"""
##########################################################################
#
# Generation of boundary representation from arbitrary geophysical
# fields and initialisation for anisotropic, unstructured meshing.
#
# Copyright (C) 2011-2013 Dr Adam S. Candy, adam.candy@imperial.ac.uk
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
##########################################################################
##########################################################################
#
# Generation of boundary representation from arbitrary geophysical
# fields and initialisation for anisotropic, unstructured meshing.
#
# Copyright (C) 2011-2013 Dr Adam S. Candy, adam.candy@imperial.ac.uk
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
##########################################################################
Sets the read and the write file stream according to
the command line arguments given.
The first argument specifies which shape file the user
wants to specify the boundaries of
The second arguments specifies the boundary polygon
The third argument specifies the file path to which the
new shape has to be written
"""
from shapely.ops import cascaded_union
import shapefile
from shapely.geometry import *
import sys
import matplotlib.pyplot as pyplot
"""
Used if there is only one shape in the boundary shapefile.
"""
def getBoundaryPointsList(bounds):
shapes = bounds.shapes()
pointsList = []
for i in range(len(shapes)):
pointsList.append(shapes[i].points)
polygons = []
for j in range(len(pointsList)):
polygons.append(Polygon([pointsList[j][i] for i in range(len(pointsList[j]))]))
return (polygons,pointsList)
"""
When more than one shapefile, works out which objects from the .shp file are overlapping and return the exterior coords of this new shape.
"""
def overlap (bounds, plot = False):
# Read shapefile and work out overlap. Optional plot.
shapes = bounds.shapes()
pointsList = []
for i in range(len(shapes)):
# Check datapoint is valid.
pointsList.append(shapes[i].points)
# Turn the points into polygons.
polygons = []
for j in range(len(pointsList)):
polygons.append(Polygon([pointsList[j][i] for i in range(len(pointsList[j]))]))
# Add overlapping shapes into a list so we know which to join together.
overlapping = []
for n in range(len(polygons) - 1):
if polygons[n].intersects(polygons[n+1]) == True:
# Two if statements to make sure the same polygon isn't being entered more than once.
if polygons[n] not in overlapping: overlapping.append(polygons[n])
if polygons[n + 1] not in overlapping: overlapping.append(polygons[n + 1])
# Create a new shape from the overlapping shapes.
join = cascaded_union(overlapping)
poly = [join]
# Take the coords. of the perimeter of this new shape.
coords = []
for i in range(len(join.exterior.coords)):
coords.append(list(join.exterior.coords[i]))
# Plot results if True. Store x-y coords of the perimeter in two lists to plot.
if plot == True:
x = []; y = []
for i in range(len(coords)):
x.append(coords[i][0]); y.append(coords[i][1])
# Plot results.
pyplot.plot(x, y)
pyplot.xlim(-4, 4)
pyplot.ylim(-4, 4)
pyplot.show()
return join, coords
"""
Output the final results as a .geo file.
"""
def write_geo (coords, filename):
# Write new shape to .geo.
print coords
target = open("%s.geo" % filename, "w") # Creates .geo file to write to.
for i in range(len(coords)):
# Write point.
target.write('Point(%d) = {%.3f, %.3f, 0, %.3f};\n' %(i + 1, coords[i][0], coords[i][1], 1.0))
# Write the lines connecting the sequential points.
if (i + 1 > 1): target.write('Line(%d) = {%d, %d};\n' % (i, i, i + 1))
# Connect first and last points.
target.write('Line(%d) = {%d, %d};\n' % (i + 1, 1, i + 1))
target.close()
return False
assert len(sys.argv)==4, "Incorrect Number of Arguments passed"
readPath = sys.argv[1]
boundaryPath = sys.argv[2]
writePath = sys.argv[3]
#input stream for the given shape
sf = shapefile.Reader(readPath)
#shapes contained in the given file
shapes = sf.shapes();
#boundary = bounds.shapes[0]
boundary = shapefile.Reader(boundaryPath)
print(len(boundary.shapes()))
if (len(boundary.shapes()) > 1): boundaryPolygons, boundaryPointList1 = overlap(boundary); boundaryPointList = [boundaryPointList1]
else: boundaryPolygons, = getBoundaryPointsList(boundary)
print len(shapes[0].points)
"""
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.
"""
for shape in shapes:
x = []; y = []
polygon = Polygon([shape.points[i] for i in range(len(shape.points))])
if (polygon.intersects(boundaryPolygons)):
for i in range(len(polygon.exterior.coords)):
point = Point(polygon.exterior.coords[i])
if point.intersects(boundaryPolygons):
x.append(shape.points[i][0]); y.append(shape.points[i][1])
print shape.points[i]
pyplot.plot(x, y)
"""
Plot boundary.
"""
x=[]
y=[]
for p in boundaryPointList :
for p2 in p :
x.append(p2[0])
y.append(p2[1])
pyplot.plot(x,y)
pyplot.xlim(min(x)-1,max(x)+1)
pyplot.ylim(min(y)-1,max(y)+1)
pyplot.show()
|