1 | |
---|
2 | |
---|
3 | |
---|
4 | |
---|
5 | ##################################### |
---|
6 | #POLYGON STUFF |
---|
7 | # |
---|
8 | #FIXME: All these should be put into new module polygon.py |
---|
9 | |
---|
10 | |
---|
11 | |
---|
12 | def point_on_line(x, y, x0, y0, x1, y1): |
---|
13 | """Determine whether a point is on a line segment |
---|
14 | |
---|
15 | Input: x, y, x0, x0, x1, y1: where |
---|
16 | point is given by x, y |
---|
17 | line is given by (x0, y0) and (x1, y1) |
---|
18 | |
---|
19 | """ |
---|
20 | |
---|
21 | from Numeric import array, dot, allclose |
---|
22 | from math import sqrt |
---|
23 | |
---|
24 | a = array([x - x0, y - y0]) |
---|
25 | a_normal = array([a[1], -a[0]]) |
---|
26 | |
---|
27 | b = array([x1 - x0, y1 - y0]) |
---|
28 | |
---|
29 | if dot(a_normal, b) == 0: |
---|
30 | #Point is somewhere on the infinite extension of the line |
---|
31 | |
---|
32 | len_a = sqrt(sum(a**2)) |
---|
33 | len_b = sqrt(sum(b**2)) |
---|
34 | if dot(a, b) >= 0 and len_a <= len_b: |
---|
35 | return True |
---|
36 | else: |
---|
37 | return False |
---|
38 | else: |
---|
39 | return False |
---|
40 | |
---|
41 | |
---|
42 | |
---|
43 | #These have been redefined to use separate_by_polygon which will |
---|
44 | #put all inside indices in the first part of the indices array and |
---|
45 | #outside indices in the last |
---|
46 | |
---|
47 | def inside_polygon(points, polygon, closed = True, verbose = False): |
---|
48 | """See separate_points_by_polygon for documentation |
---|
49 | """ |
---|
50 | |
---|
51 | from Numeric import array, Float, reshape |
---|
52 | |
---|
53 | if verbose: print 'Checking input to inside_polygon' |
---|
54 | try: |
---|
55 | points = ensure_numeric(points, Float) |
---|
56 | except: |
---|
57 | msg = 'Points could not be converted to Numeric array' |
---|
58 | raise msg |
---|
59 | |
---|
60 | try: |
---|
61 | polygon = ensure_numeric(polygon, Float) |
---|
62 | except: |
---|
63 | msg = 'Polygon could not be converted to Numeric array' |
---|
64 | raise msg |
---|
65 | |
---|
66 | |
---|
67 | |
---|
68 | if len(points.shape) == 1: |
---|
69 | one_point = True |
---|
70 | points = reshape(points, (1,2)) |
---|
71 | else: |
---|
72 | one_point = False |
---|
73 | |
---|
74 | indices, count = separate_points_by_polygon(points, polygon, |
---|
75 | closed, verbose) |
---|
76 | |
---|
77 | if one_point: |
---|
78 | return count == 1 |
---|
79 | else: |
---|
80 | return indices[:count] |
---|
81 | |
---|
82 | def outside_polygon(points, polygon, closed = True, verbose = False): |
---|
83 | """See separate_points_by_polygon for documentation |
---|
84 | """ |
---|
85 | |
---|
86 | from Numeric import array, Float, reshape |
---|
87 | |
---|
88 | if verbose: print 'Checking input to outside_polygon' |
---|
89 | try: |
---|
90 | points = ensure_numeric(points, Float) |
---|
91 | except: |
---|
92 | msg = 'Points could not be converted to Numeric array' |
---|
93 | raise msg |
---|
94 | |
---|
95 | try: |
---|
96 | polygon = ensure_numeric(polygon, Float) |
---|
97 | except: |
---|
98 | msg = 'Polygon could not be converted to Numeric array' |
---|
99 | raise msg |
---|
100 | |
---|
101 | |
---|
102 | |
---|
103 | if len(points.shape) == 1: |
---|
104 | one_point = True |
---|
105 | points = reshape(points, (1,2)) |
---|
106 | else: |
---|
107 | one_point = False |
---|
108 | |
---|
109 | indices, count = separate_points_by_polygon(points, polygon, |
---|
110 | closed, verbose) |
---|
111 | |
---|
112 | if one_point: |
---|
113 | return count != 1 |
---|
114 | else: |
---|
115 | return indices[count:][::-1] #return reversed |
---|
116 | |
---|
117 | |
---|
118 | def separate_points_by_polygon(points, polygon, |
---|
119 | closed = True, verbose = False): |
---|
120 | """Determine whether points are inside or outside a polygon |
---|
121 | |
---|
122 | Input: |
---|
123 | points - Tuple of (x, y) coordinates, or list of tuples |
---|
124 | polygon - list of vertices of polygon |
---|
125 | closed - (optional) determine whether points on boundary should be |
---|
126 | regarded as belonging to the polygon (closed = True) |
---|
127 | or not (closed = False) |
---|
128 | |
---|
129 | Outputs: |
---|
130 | indices: array of same length as points with indices of points falling |
---|
131 | inside the polygon listed from the beginning and indices of points |
---|
132 | falling outside listed from the end. |
---|
133 | |
---|
134 | count: count of points falling inside the polygon |
---|
135 | |
---|
136 | The indices of points inside are obtained as indices[:count] |
---|
137 | The indices of points outside are obtained as indices[count:] |
---|
138 | |
---|
139 | |
---|
140 | Examples: |
---|
141 | U = [[0,0], [1,0], [1,1], [0,1]] #Unit square |
---|
142 | |
---|
143 | separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U) |
---|
144 | will return the indices [0, 2, 1] and count == 2 as only the first |
---|
145 | and the last point are inside the unit square |
---|
146 | |
---|
147 | Remarks: |
---|
148 | The vertices may be listed clockwise or counterclockwise and |
---|
149 | the first point may optionally be repeated. |
---|
150 | Polygons do not need to be convex. |
---|
151 | Polygons can have holes in them and points inside a hole is |
---|
152 | regarded as being outside the polygon. |
---|
153 | |
---|
154 | Algorithm is based on work by Darel Finley, |
---|
155 | http://www.alienryderflex.com/polygon/ |
---|
156 | |
---|
157 | """ |
---|
158 | |
---|
159 | from Numeric import array, Float, reshape, Int, zeros |
---|
160 | |
---|
161 | |
---|
162 | #Input checks |
---|
163 | try: |
---|
164 | points = ensure_numeric(points, Float) |
---|
165 | except: |
---|
166 | msg = 'Points could not be converted to Numeric array' |
---|
167 | raise msg |
---|
168 | |
---|
169 | try: |
---|
170 | polygon = ensure_numeric(polygon, Float) |
---|
171 | except: |
---|
172 | msg = 'Polygon could not be converted to Numeric array' |
---|
173 | raise msg |
---|
174 | |
---|
175 | assert len(polygon.shape) == 2,\ |
---|
176 | 'Polygon array must be a 2d array of vertices' |
---|
177 | |
---|
178 | assert polygon.shape[1] == 2,\ |
---|
179 | 'Polygon array must have two columns' |
---|
180 | |
---|
181 | assert len(points.shape) == 2,\ |
---|
182 | 'Points array must be a 2d array' |
---|
183 | |
---|
184 | assert points.shape[1] == 2,\ |
---|
185 | 'Points array must have two columns' |
---|
186 | |
---|
187 | N = polygon.shape[0] #Number of vertices in polygon |
---|
188 | M = points.shape[0] #Number of points |
---|
189 | |
---|
190 | px = polygon[:,0] |
---|
191 | py = polygon[:,1] |
---|
192 | |
---|
193 | #Used for an optimisation when points are far away from polygon |
---|
194 | minpx = min(px); maxpx = max(px) |
---|
195 | minpy = min(py); maxpy = max(py) |
---|
196 | |
---|
197 | |
---|
198 | #Begin main loop |
---|
199 | indices = zeros(M, Int) |
---|
200 | |
---|
201 | inside_index = 0 #Keep track of points inside |
---|
202 | outside_index = M-1 #Keep track of points outside (starting from end) |
---|
203 | |
---|
204 | for k in range(M): |
---|
205 | |
---|
206 | if verbose: |
---|
207 | if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M) |
---|
208 | |
---|
209 | #for each point |
---|
210 | x = points[k, 0] |
---|
211 | y = points[k, 1] |
---|
212 | |
---|
213 | inside = False |
---|
214 | |
---|
215 | if not x > maxpx or x < minpx or y > maxpy or y < minpy: |
---|
216 | #Check polygon |
---|
217 | for i in range(N): |
---|
218 | j = (i+1)%N |
---|
219 | |
---|
220 | #Check for case where point is contained in line segment |
---|
221 | if point_on_line(x, y, px[i], py[i], px[j], py[j]): |
---|
222 | if closed: |
---|
223 | inside = True |
---|
224 | else: |
---|
225 | inside = False |
---|
226 | break |
---|
227 | else: |
---|
228 | #Check if truly inside polygon |
---|
229 | if py[i] < y and py[j] >= y or\ |
---|
230 | py[j] < y and py[i] >= y: |
---|
231 | if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x: |
---|
232 | inside = not inside |
---|
233 | |
---|
234 | if inside: |
---|
235 | indices[inside_index] = k |
---|
236 | inside_index += 1 |
---|
237 | else: |
---|
238 | indices[outside_index] = k |
---|
239 | outside_index -= 1 |
---|
240 | |
---|
241 | return indices, inside_index |
---|
242 | |
---|
243 | |
---|
244 | def separate_points_by_polygon_c(points, polygon, |
---|
245 | closed = True, verbose = False): |
---|
246 | """Determine whether points are inside or outside a polygon |
---|
247 | |
---|
248 | C-wrapper |
---|
249 | """ |
---|
250 | |
---|
251 | from Numeric import array, Float, reshape, zeros, Int |
---|
252 | |
---|
253 | |
---|
254 | if verbose: print 'Checking input to separate_points_by_polygon' |
---|
255 | #Input checks |
---|
256 | try: |
---|
257 | points = ensure_numeric(points, Float) |
---|
258 | except: |
---|
259 | msg = 'Points could not be converted to Numeric array' |
---|
260 | raise msg |
---|
261 | |
---|
262 | #if verbose: print 'Checking input to separate_points_by_polygon 2' |
---|
263 | try: |
---|
264 | polygon = ensure_numeric(polygon, Float) |
---|
265 | except: |
---|
266 | msg = 'Polygon could not be converted to Numeric array' |
---|
267 | raise msg |
---|
268 | |
---|
269 | if verbose: print 'check' |
---|
270 | |
---|
271 | assert len(polygon.shape) == 2,\ |
---|
272 | 'Polygon array must be a 2d array of vertices' |
---|
273 | |
---|
274 | assert polygon.shape[1] == 2,\ |
---|
275 | 'Polygon array must have two columns' |
---|
276 | |
---|
277 | assert len(points.shape) == 2,\ |
---|
278 | 'Points array must be a 2d array' |
---|
279 | |
---|
280 | assert points.shape[1] == 2,\ |
---|
281 | 'Points array must have two columns' |
---|
282 | |
---|
283 | N = polygon.shape[0] #Number of vertices in polygon |
---|
284 | M = points.shape[0] #Number of points |
---|
285 | |
---|
286 | from util_ext import separate_points_by_polygon |
---|
287 | |
---|
288 | if verbose: print 'Allocating array for indices' |
---|
289 | |
---|
290 | indices = zeros( M, Int ) |
---|
291 | |
---|
292 | if verbose: print 'Calling C-version of inside poly' |
---|
293 | count = separate_points_by_polygon(points, polygon, indices, |
---|
294 | int(closed), int(verbose)) |
---|
295 | |
---|
296 | return indices, count |
---|
297 | |
---|
298 | |
---|
299 | |
---|
300 | class Polygon_function: |
---|
301 | """Create callable object f: x,y -> z, where a,y,z are vectors and |
---|
302 | where f will return different values depending on whether x,y belongs |
---|
303 | to specified polygons. |
---|
304 | |
---|
305 | To instantiate: |
---|
306 | |
---|
307 | Polygon_function(polygons) |
---|
308 | |
---|
309 | where polygons is a list of tuples of the form |
---|
310 | |
---|
311 | [ (P0, v0), (P1, v1), ...] |
---|
312 | |
---|
313 | with Pi being lists of vertices defining polygons and vi either |
---|
314 | constants or functions of x,y to be applied to points with the polygon. |
---|
315 | |
---|
316 | The function takes an optional argument, default which is the value |
---|
317 | (or function) to used for points not belonging to any polygon. |
---|
318 | For example: |
---|
319 | |
---|
320 | Polygon_function(polygons, default = 0.03) |
---|
321 | |
---|
322 | If omitted the default value will be 0.0 |
---|
323 | |
---|
324 | Note: If two polygons overlap, the one last in the list takes precedence |
---|
325 | |
---|
326 | FIXME? : Currently, coordinates specified here are assumed to be relative to the origin (georeference) used by domain. This function is more general than domain so perhaps it'd be an idea to allow the optional argument georeference??? |
---|
327 | |
---|
328 | """ |
---|
329 | |
---|
330 | def __init__(self, regions, default = 0.0): |
---|
331 | |
---|
332 | try: |
---|
333 | len(regions) |
---|
334 | except: |
---|
335 | msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons |
---|
336 | raise msg |
---|
337 | |
---|
338 | |
---|
339 | T = regions[0] |
---|
340 | try: |
---|
341 | a = len(T) |
---|
342 | except: |
---|
343 | msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons |
---|
344 | raise msg |
---|
345 | |
---|
346 | assert a == 2, 'Must have two component each: %s' %T |
---|
347 | |
---|
348 | self.regions = regions |
---|
349 | self.default = default |
---|
350 | |
---|
351 | |
---|
352 | def __call__(self, x, y): |
---|
353 | from util import inside_polygon |
---|
354 | from Numeric import ones, Float, concatenate, array, reshape, choose |
---|
355 | |
---|
356 | x = array(x).astype(Float) |
---|
357 | y = array(y).astype(Float) |
---|
358 | |
---|
359 | N = len(x) |
---|
360 | assert len(y) == N |
---|
361 | |
---|
362 | points = concatenate( (reshape(x, (N, 1)), |
---|
363 | reshape(y, (N, 1))), axis=1 ) |
---|
364 | |
---|
365 | if callable(self.default): |
---|
366 | z = self.default(x,y) |
---|
367 | else: |
---|
368 | z = ones(N, Float) * self.default |
---|
369 | |
---|
370 | for polygon, value in self.regions: |
---|
371 | indices = inside_polygon(points, polygon) |
---|
372 | |
---|
373 | #FIXME: This needs to be vectorised |
---|
374 | if callable(value): |
---|
375 | for i in indices: |
---|
376 | xx = array([x[i]]) |
---|
377 | yy = array([y[i]]) |
---|
378 | z[i] = value(xx, yy)[0] |
---|
379 | else: |
---|
380 | for i in indices: |
---|
381 | z[i] = value |
---|
382 | |
---|
383 | return z |
---|
384 | |
---|
385 | def read_polygon(filename): |
---|
386 | """Read points assumed to form a polygon |
---|
387 | There must be exactly two numbers in each line |
---|
388 | """ |
---|
389 | |
---|
390 | #Get polygon |
---|
391 | fid = open(filename) |
---|
392 | lines = fid.readlines() |
---|
393 | fid.close() |
---|
394 | polygon = [] |
---|
395 | for line in lines: |
---|
396 | fields = line.split(',') |
---|
397 | polygon.append( [float(fields[0]), float(fields[1])] ) |
---|
398 | |
---|
399 | return polygon |
---|
400 | |
---|
401 | def populate_polygon(polygon, number_of_points, seed = None, exclude = None): |
---|
402 | """Populate given polygon with uniformly distributed points. |
---|
403 | |
---|
404 | Input: |
---|
405 | polygon - list of vertices of polygon |
---|
406 | number_of_points - (optional) number of points |
---|
407 | seed - seed for random number generator (default=None) |
---|
408 | exclude - list of polygons (inside main polygon) from where points should be excluded |
---|
409 | |
---|
410 | Output: |
---|
411 | points - list of points inside polygon |
---|
412 | |
---|
413 | Examples: |
---|
414 | populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 ) |
---|
415 | will return five randomly selected points inside the unit square |
---|
416 | """ |
---|
417 | |
---|
418 | from random import uniform, seed |
---|
419 | |
---|
420 | seed(seed) |
---|
421 | |
---|
422 | points = [] |
---|
423 | |
---|
424 | #Find outer extent of polygon |
---|
425 | max_x = min_x = polygon[0][0] |
---|
426 | max_y = min_y = polygon[0][1] |
---|
427 | for point in polygon[1:]: |
---|
428 | x = point[0] |
---|
429 | if x > max_x: max_x = x |
---|
430 | if x < min_x: min_x = x |
---|
431 | y = point[1] |
---|
432 | if y > max_y: max_y = y |
---|
433 | if y < min_y: min_y = y |
---|
434 | |
---|
435 | |
---|
436 | while len(points) < number_of_points: |
---|
437 | x = uniform(min_x, max_x) |
---|
438 | y = uniform(min_y, max_y) |
---|
439 | |
---|
440 | append = False |
---|
441 | if inside_polygon( [x,y], polygon ): |
---|
442 | |
---|
443 | append = True |
---|
444 | |
---|
445 | #Check exclusions |
---|
446 | if exclude is not None: |
---|
447 | for ex_poly in exclude: |
---|
448 | if inside_polygon( [x,y], ex_poly ): |
---|
449 | append = False |
---|
450 | |
---|
451 | |
---|
452 | if append is True: |
---|
453 | points.append([x,y]) |
---|
454 | |
---|
455 | return points |
---|
456 | |
---|
457 | |
---|
458 | |
---|
459 | |
---|
460 | |
---|
461 | ############################################## |
---|
462 | #Initialise module |
---|
463 | |
---|
464 | import compile |
---|
465 | if compile.can_use_C_extension('polygon_ext.c'): |
---|
466 | from polygon_ext import point_on_line |
---|
467 | separate_points_by_polygon = separate_points_by_polygon_c |
---|
468 | |
---|
469 | |
---|
470 | if __name__ == "__main__": |
---|
471 | pass |
---|