1 | """This module contains various auxiliary function used by pyvolution. |
---|
2 | |
---|
3 | It is also a clearing house for functions that may later earn a module |
---|
4 | of their own. |
---|
5 | """ |
---|
6 | |
---|
7 | #FIXME: Move polygon stuff out |
---|
8 | |
---|
9 | def angle(v): |
---|
10 | """Compute angle between e1 (the unit vector in the x-direction) |
---|
11 | and the specified vector |
---|
12 | """ |
---|
13 | |
---|
14 | from math import acos, pi, sqrt |
---|
15 | from Numeric import sum, array |
---|
16 | |
---|
17 | l = sqrt( sum (array(v)**2)) |
---|
18 | v1 = v[0]/l |
---|
19 | v2 = v[1]/l |
---|
20 | |
---|
21 | |
---|
22 | theta = acos(v1) |
---|
23 | |
---|
24 | #try: |
---|
25 | # theta = acos(v1) |
---|
26 | #except ValueError, e: |
---|
27 | # print 'WARNING (util.py): Angle acos(%s) failed: %s' %(str(v1), e) |
---|
28 | # |
---|
29 | # #FIXME, hack to avoid acos(1.0) Value error |
---|
30 | # # why is it happening? |
---|
31 | # # is it catching something we should avoid? |
---|
32 | # #FIXME (Ole): When did this happen? We need a unit test to |
---|
33 | # #reveal this condition or otherwise remove the hack |
---|
34 | # |
---|
35 | # s = 1e-6 |
---|
36 | # if (v1+s > 1.0) and (v1-s < 1.0) : |
---|
37 | # theta = 0.0 |
---|
38 | # elif (v1+s > -1.0) and (v1-s < -1.0): |
---|
39 | # theta = 3.1415926535897931 |
---|
40 | # print 'WARNING (util.py): angle v1 is %f, setting acos to %f '\ |
---|
41 | # %(v1, theta) |
---|
42 | |
---|
43 | if v2 < 0: |
---|
44 | #Quadrant 3 or 4 |
---|
45 | theta = 2*pi-theta |
---|
46 | |
---|
47 | return theta |
---|
48 | |
---|
49 | |
---|
50 | def anglediff(v0, v1): |
---|
51 | """Compute difference between angle of vector x0, y0 and x1, y1. |
---|
52 | This is used for determining the ordering of vertices, |
---|
53 | e.g. for checking if they are counter clockwise. |
---|
54 | |
---|
55 | Always return a positive value |
---|
56 | """ |
---|
57 | |
---|
58 | from math import pi |
---|
59 | |
---|
60 | a0 = angle(v0) |
---|
61 | a1 = angle(v1) |
---|
62 | |
---|
63 | #Ensure that difference will be positive |
---|
64 | if a0 < a1: |
---|
65 | a0 += 2*pi |
---|
66 | |
---|
67 | return a0-a1 |
---|
68 | |
---|
69 | |
---|
70 | def mean(x): |
---|
71 | from Numeric import sum |
---|
72 | return sum(x)/len(x) |
---|
73 | |
---|
74 | |
---|
75 | def point_on_line(x, y, x0, y0, x1, y1): |
---|
76 | """Determine whether a point is on a line segment |
---|
77 | |
---|
78 | Input: x, y, x0, x0, x1, y1: where |
---|
79 | point is given by x, y |
---|
80 | line is given by (x0, y0) and (x1, y1) |
---|
81 | |
---|
82 | """ |
---|
83 | |
---|
84 | from Numeric import array, dot, allclose |
---|
85 | from math import sqrt |
---|
86 | |
---|
87 | a = array([x - x0, y - y0]) |
---|
88 | a_normal = array([a[1], -a[0]]) |
---|
89 | |
---|
90 | b = array([x1 - x0, y1 - y0]) |
---|
91 | |
---|
92 | if dot(a_normal, b) == 0: |
---|
93 | #Point is somewhere on the infinite extension of the line |
---|
94 | |
---|
95 | len_a = sqrt(sum(a**2)) |
---|
96 | len_b = sqrt(sum(b**2)) |
---|
97 | if dot(a, b) >= 0 and len_a <= len_b: |
---|
98 | return True |
---|
99 | else: |
---|
100 | return False |
---|
101 | else: |
---|
102 | return False |
---|
103 | |
---|
104 | def ensure_numeric(A, typecode = None): |
---|
105 | """Ensure that sequence is a Numeric array. |
---|
106 | Inputs: |
---|
107 | A: Sequance. If A is already a Numeric array it will be returned |
---|
108 | unaltered |
---|
109 | If not, an attempt is made to convert it to a Numeric |
---|
110 | array |
---|
111 | typecode: Numeric type. If specified, use this in the conversion. |
---|
112 | If not, let Numeric decide |
---|
113 | |
---|
114 | This function is necessary as array(A) can cause memory overflow. |
---|
115 | """ |
---|
116 | |
---|
117 | from Numeric import ArrayType, array |
---|
118 | |
---|
119 | if typecode is None: |
---|
120 | if type(A) == ArrayType: |
---|
121 | return A |
---|
122 | else: |
---|
123 | return array(A) |
---|
124 | else: |
---|
125 | if type(A) == ArrayType: |
---|
126 | if A.typecode == typecode: |
---|
127 | return array(A) #FIXME: Shouldn't this be just return A? |
---|
128 | else: |
---|
129 | return A.astype(typecode) |
---|
130 | else: |
---|
131 | return array(A).astype(typecode) |
---|
132 | |
---|
133 | |
---|
134 | |
---|
135 | def file_function(filename, |
---|
136 | domain = None, |
---|
137 | quantities = None, |
---|
138 | interpolation_points = None, verbose = False): |
---|
139 | """Read time history of spatial data from NetCDF file and return |
---|
140 | a callable object. |
---|
141 | |
---|
142 | If the file has extension 'sww' then it is assumed to be spatio-temporal |
---|
143 | or temporal and the callable object will have the form f(t,x,y) or f(t) |
---|
144 | depending on whether the file contains spatial data |
---|
145 | |
---|
146 | If the file has extension 'tms' then it is assumed to be temporal only |
---|
147 | and the callable object will have the form f(t) |
---|
148 | |
---|
149 | Either form will return interpolated values based on the input file |
---|
150 | using the underlying interpolation_function. |
---|
151 | |
---|
152 | If domain is specified, model time (domain.starttime) |
---|
153 | will be checked and possibly modified. |
---|
154 | |
---|
155 | All times are assumed to be in UTC |
---|
156 | |
---|
157 | All spatial information is assumed to be in UTM coordinates. |
---|
158 | |
---|
159 | See Interpolation function for further documentation |
---|
160 | """ |
---|
161 | |
---|
162 | |
---|
163 | #FIXME (OLE): Should check origin of domain against that of file |
---|
164 | #In fact, this is where origin should be converted to that of domain |
---|
165 | #Also, check that file covers domain fully. |
---|
166 | #If we use the suggested Point_set class for interpolation points |
---|
167 | #here it would be easier |
---|
168 | |
---|
169 | #Take into account: |
---|
170 | #- domain's georef |
---|
171 | #- sww file's georef |
---|
172 | #- interpolation points georef |
---|
173 | |
---|
174 | |
---|
175 | assert type(filename) == type(''),\ |
---|
176 | 'First argument to File_function must be a string' |
---|
177 | |
---|
178 | try: |
---|
179 | fid = open(filename) |
---|
180 | except Exception, e: |
---|
181 | msg = 'File "%s" could not be opened: Error="%s"'\ |
---|
182 | %(filename, e) |
---|
183 | raise msg |
---|
184 | |
---|
185 | line = fid.readline() |
---|
186 | fid.close() |
---|
187 | |
---|
188 | if quantities is None: |
---|
189 | if domain is not None: |
---|
190 | quantities = domain.conserved_quantities |
---|
191 | |
---|
192 | |
---|
193 | |
---|
194 | if line[:3] == 'CDF': |
---|
195 | return get_netcdf_file_function(filename, domain, quantities, |
---|
196 | interpolation_points, verbose = verbose) |
---|
197 | else: |
---|
198 | raise 'Must be a NetCDF File' |
---|
199 | |
---|
200 | |
---|
201 | |
---|
202 | def get_netcdf_file_function(filename, domain=None, quantity_names=None, |
---|
203 | interpolation_points=None, verbose = False): |
---|
204 | """Read time history of spatial data from NetCDF sww file and |
---|
205 | return a callable object f(t,x,y) |
---|
206 | which will return interpolated values based on the input file. |
---|
207 | |
---|
208 | If domain is specified, model time (domain.starttime) |
---|
209 | will be checked and possibly modified |
---|
210 | |
---|
211 | All times are assumed to be in UTC |
---|
212 | |
---|
213 | See Interpolation function for further documetation |
---|
214 | |
---|
215 | |
---|
216 | """ |
---|
217 | |
---|
218 | #verbose = True |
---|
219 | |
---|
220 | #FIXME: Check that model origin is the same as file's origin |
---|
221 | #(both in UTM coordinates) |
---|
222 | #If not - modify those from file to match domain |
---|
223 | #Take this code from e.g. dem2pts in data_manager.py |
---|
224 | #FIXME: Use geo_reference to read and write xllcorner... |
---|
225 | |
---|
226 | |
---|
227 | import time, calendar, types |
---|
228 | from config import time_format |
---|
229 | from Scientific.IO.NetCDF import NetCDFFile |
---|
230 | from Numeric import array, zeros, Float, alltrue, concatenate, reshape |
---|
231 | |
---|
232 | #Open NetCDF file |
---|
233 | if verbose: print 'Reading', filename |
---|
234 | fid = NetCDFFile(filename, 'r') |
---|
235 | |
---|
236 | if type(quantity_names) == types.StringType: |
---|
237 | quantity_names = [quantity_names] |
---|
238 | |
---|
239 | if quantity_names is None or len(quantity_names) < 1: |
---|
240 | #If no quantities are specified get quantities from file |
---|
241 | #x, y, time are assumed as the independent variables so |
---|
242 | #they are excluded as potentiol quantities |
---|
243 | quantity_names = [] |
---|
244 | for name in fid.variables.keys(): |
---|
245 | if name not in ['x', 'y', 'time']: |
---|
246 | quantity_names.append(name) |
---|
247 | |
---|
248 | if len(quantity_names) < 1: |
---|
249 | msg = 'ERROR: At least one independent value must be specified' |
---|
250 | raise msg |
---|
251 | |
---|
252 | |
---|
253 | #Now assert that requested quantitites (and the independent ones) |
---|
254 | #are present in file |
---|
255 | missing = [] |
---|
256 | for quantity in ['time'] + quantity_names: |
---|
257 | if not fid.variables.has_key(quantity): |
---|
258 | missing.append(quantity) |
---|
259 | |
---|
260 | if len(missing) > 0: |
---|
261 | msg = 'Quantities %s could not be found in file %s'\ |
---|
262 | %(str(missing), filename) |
---|
263 | fid.close() |
---|
264 | raise msg |
---|
265 | |
---|
266 | #Decide whether this data has a spatial dimension |
---|
267 | spatial = True |
---|
268 | for quantity in ['x', 'y']: |
---|
269 | if not fid.variables.has_key(quantity): |
---|
270 | spatial = False |
---|
271 | |
---|
272 | if filename[-3:] == 'tms' and spatial is True: |
---|
273 | msg = 'Files of type tms must not contain spatial information' |
---|
274 | raise msg |
---|
275 | |
---|
276 | if filename[-3:] == 'sww' and spatial is False: |
---|
277 | msg = 'Files of type sww must contain spatial information' |
---|
278 | raise msg |
---|
279 | |
---|
280 | #Get first timestep |
---|
281 | try: |
---|
282 | starttime = fid.starttime[0] |
---|
283 | except ValueError: |
---|
284 | msg = 'Could not read starttime from file %s' %filename |
---|
285 | raise msg |
---|
286 | |
---|
287 | #Get variables |
---|
288 | if verbose: print 'Get variables' |
---|
289 | time = fid.variables['time'][:] |
---|
290 | |
---|
291 | if spatial: |
---|
292 | #Get origin |
---|
293 | xllcorner = fid.xllcorner[0] |
---|
294 | yllcorner = fid.yllcorner[0] |
---|
295 | zone = fid.zone[0] |
---|
296 | |
---|
297 | x = fid.variables['x'][:] |
---|
298 | y = fid.variables['y'][:] |
---|
299 | triangles = fid.variables['volumes'][:] |
---|
300 | |
---|
301 | x = reshape(x, (len(x),1)) |
---|
302 | y = reshape(y, (len(y),1)) |
---|
303 | vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array |
---|
304 | |
---|
305 | |
---|
306 | if domain is not None: |
---|
307 | #Update domain.startime if it is *earlier* than starttime |
---|
308 | if starttime > domain.starttime: |
---|
309 | msg = 'WARNING: Start time as specified in domain (%f)'\ |
---|
310 | %domain.starttime |
---|
311 | msg += ' is earlier than the starttime of file %s (%f).'\ |
---|
312 | %(filename, starttime) |
---|
313 | msg += ' Modifying domain starttime accordingly.' |
---|
314 | |
---|
315 | if verbose: print msg |
---|
316 | domain.starttime = starttime #Modifying model time |
---|
317 | if verbose: print 'Domain starttime is now set to %f'\ |
---|
318 | %domain.starttime |
---|
319 | |
---|
320 | |
---|
321 | #If domain.startime is *later* than starttime, |
---|
322 | #move time back - relative to domain's time |
---|
323 | if domain.starttime > starttime: |
---|
324 | time = time - domain.starttime + starttime |
---|
325 | |
---|
326 | #FIXME Use method in geo to reconcile |
---|
327 | #if spatial: |
---|
328 | #assert domain.geo_reference.xllcorner == xllcorner |
---|
329 | #assert domain.geo_reference.yllcorner == yllcorner |
---|
330 | #assert domain.geo_reference.zone == zone |
---|
331 | |
---|
332 | if verbose: |
---|
333 | print 'File_function data obtained from: %s' %filename |
---|
334 | print ' References:' |
---|
335 | #print ' Datum ....' #FIXME |
---|
336 | if spatial: |
---|
337 | print ' Lower left corner: [%f, %f]'\ |
---|
338 | %(xllcorner, yllcorner) |
---|
339 | print ' Start time: %f' %starttime |
---|
340 | |
---|
341 | |
---|
342 | #Produce values for desired data points at |
---|
343 | #each timestep |
---|
344 | |
---|
345 | quantities = {} |
---|
346 | for i, name in enumerate(quantity_names): |
---|
347 | quantities[name] = fid.variables[name][:] |
---|
348 | fid.close() |
---|
349 | |
---|
350 | |
---|
351 | from least_squares import Interpolation_function |
---|
352 | |
---|
353 | if not spatial: |
---|
354 | vertex_coordinates = triangles = interpolation_points = None |
---|
355 | |
---|
356 | return Interpolation_function(time, |
---|
357 | quantities, |
---|
358 | quantity_names, |
---|
359 | vertex_coordinates, |
---|
360 | triangles, |
---|
361 | interpolation_points, |
---|
362 | verbose = verbose) |
---|
363 | |
---|
364 | |
---|
365 | |
---|
366 | |
---|
367 | |
---|
368 | def read_xya(filename, format = 'netcdf'): |
---|
369 | """Read simple xya file, possibly with a header in the first line, just |
---|
370 | x y [attributes] |
---|
371 | separated by whitespace |
---|
372 | |
---|
373 | Format can be either ASCII or NetCDF |
---|
374 | |
---|
375 | Return list of points, list of attributes and |
---|
376 | attribute names if present otherwise None |
---|
377 | """ |
---|
378 | #FIXME: Probably obsoleted by similar function in load_ASCII |
---|
379 | #FIXME: Name read_data_points (or see 'load' in pylab) |
---|
380 | |
---|
381 | from Scientific.IO.NetCDF import NetCDFFile |
---|
382 | |
---|
383 | if format.lower() == 'netcdf': |
---|
384 | #Get NetCDF |
---|
385 | |
---|
386 | fid = NetCDFFile(filename, 'r') |
---|
387 | |
---|
388 | # Get the variables |
---|
389 | points = fid.variables['points'][:] |
---|
390 | keys = fid.variables.keys() |
---|
391 | attributes = {} |
---|
392 | for key in keys: |
---|
393 | if key != 'points': |
---|
394 | attributes[key] = fid.variables[key][:] |
---|
395 | |
---|
396 | fid.close() |
---|
397 | else: |
---|
398 | #Read as ASCII file assuming that it is separated by whitespace |
---|
399 | fid = open(filename) |
---|
400 | lines = fid.readlines() |
---|
401 | fid.close() |
---|
402 | |
---|
403 | #Check if there is a header line |
---|
404 | fields = lines[0].strip().split() |
---|
405 | try: |
---|
406 | float(fields[0]) |
---|
407 | except: |
---|
408 | #This must be a header line |
---|
409 | attribute_names = fields |
---|
410 | lines = lines[1:] |
---|
411 | else: |
---|
412 | attribute_names = ['elevation'] #HACK |
---|
413 | |
---|
414 | attributes = {} |
---|
415 | for key in attribute_names: |
---|
416 | attributes[key] = [] |
---|
417 | |
---|
418 | points = [] |
---|
419 | for line in lines: |
---|
420 | fields = line.strip().split() |
---|
421 | points.append( (float(fields[0]), float(fields[1])) ) |
---|
422 | for i, key in enumerate(attribute_names): |
---|
423 | attributes[key].append( float(fields[2+i]) ) |
---|
424 | |
---|
425 | return points, attributes |
---|
426 | |
---|
427 | |
---|
428 | ##################################### |
---|
429 | #POLYGON STUFF |
---|
430 | # |
---|
431 | #FIXME: All these should be put into new module polygon.py |
---|
432 | |
---|
433 | |
---|
434 | |
---|
435 | #These have been redefined to use separate_by_polygon which will |
---|
436 | #put all inside indices in the first part of the indices array and |
---|
437 | #outside indices in the last |
---|
438 | |
---|
439 | def inside_polygon(points, polygon, closed = True, verbose = False): |
---|
440 | """See separate_points_by_polygon for documentation |
---|
441 | """ |
---|
442 | |
---|
443 | from Numeric import array, Float, reshape |
---|
444 | |
---|
445 | if verbose: print 'Checking input to inside_polygon' |
---|
446 | try: |
---|
447 | points = ensure_numeric(points, Float) |
---|
448 | except: |
---|
449 | msg = 'Points could not be converted to Numeric array' |
---|
450 | raise msg |
---|
451 | |
---|
452 | try: |
---|
453 | polygon = ensure_numeric(polygon, Float) |
---|
454 | except: |
---|
455 | msg = 'Polygon could not be converted to Numeric array' |
---|
456 | raise msg |
---|
457 | |
---|
458 | |
---|
459 | |
---|
460 | if len(points.shape) == 1: |
---|
461 | one_point = True |
---|
462 | points = reshape(points, (1,2)) |
---|
463 | else: |
---|
464 | one_point = False |
---|
465 | |
---|
466 | indices, count = separate_points_by_polygon(points, polygon, |
---|
467 | closed, verbose) |
---|
468 | |
---|
469 | if one_point: |
---|
470 | return count == 1 |
---|
471 | else: |
---|
472 | return indices[:count] |
---|
473 | |
---|
474 | def outside_polygon(points, polygon, closed = True, verbose = False): |
---|
475 | """See separate_points_by_polygon for documentation |
---|
476 | """ |
---|
477 | |
---|
478 | from Numeric import array, Float, reshape |
---|
479 | |
---|
480 | if verbose: print 'Checking input to outside_polygon' |
---|
481 | try: |
---|
482 | points = ensure_numeric(points, Float) |
---|
483 | except: |
---|
484 | msg = 'Points could not be converted to Numeric array' |
---|
485 | raise msg |
---|
486 | |
---|
487 | try: |
---|
488 | polygon = ensure_numeric(polygon, Float) |
---|
489 | except: |
---|
490 | msg = 'Polygon could not be converted to Numeric array' |
---|
491 | raise msg |
---|
492 | |
---|
493 | |
---|
494 | |
---|
495 | if len(points.shape) == 1: |
---|
496 | one_point = True |
---|
497 | points = reshape(points, (1,2)) |
---|
498 | else: |
---|
499 | one_point = False |
---|
500 | |
---|
501 | indices, count = separate_points_by_polygon(points, polygon, |
---|
502 | closed, verbose) |
---|
503 | |
---|
504 | if one_point: |
---|
505 | return count != 1 |
---|
506 | else: |
---|
507 | return indices[count:][::-1] #return reversed |
---|
508 | |
---|
509 | |
---|
510 | def separate_points_by_polygon(points, polygon, |
---|
511 | closed = True, verbose = False): |
---|
512 | """Determine whether points are inside or outside a polygon |
---|
513 | |
---|
514 | Input: |
---|
515 | points - Tuple of (x, y) coordinates, or list of tuples |
---|
516 | polygon - list of vertices of polygon |
---|
517 | closed - (optional) determine whether points on boundary should be |
---|
518 | regarded as belonging to the polygon (closed = True) |
---|
519 | or not (closed = False) |
---|
520 | |
---|
521 | Outputs: |
---|
522 | indices: array of same length as points with indices of points falling |
---|
523 | inside the polygon listed from the beginning and indices of points |
---|
524 | falling outside listed from the end. |
---|
525 | |
---|
526 | count: count of points falling inside the polygon |
---|
527 | |
---|
528 | The indices of points inside are obtained as indices[:count] |
---|
529 | The indices of points outside are obtained as indices[count:] |
---|
530 | |
---|
531 | |
---|
532 | Examples: |
---|
533 | U = [[0,0], [1,0], [1,1], [0,1]] #Unit square |
---|
534 | |
---|
535 | separate_points_by_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]], U) |
---|
536 | will return the indices [0, 2, 1] and count == 2 as only the first |
---|
537 | and the last point are inside the unit square |
---|
538 | |
---|
539 | Remarks: |
---|
540 | The vertices may be listed clockwise or counterclockwise and |
---|
541 | the first point may optionally be repeated. |
---|
542 | Polygons do not need to be convex. |
---|
543 | Polygons can have holes in them and points inside a hole is |
---|
544 | regarded as being outside the polygon. |
---|
545 | |
---|
546 | Algorithm is based on work by Darel Finley, |
---|
547 | http://www.alienryderflex.com/polygon/ |
---|
548 | |
---|
549 | """ |
---|
550 | |
---|
551 | from Numeric import array, Float, reshape, Int, zeros |
---|
552 | |
---|
553 | |
---|
554 | #Input checks |
---|
555 | try: |
---|
556 | points = ensure_numeric(points, Float) |
---|
557 | except: |
---|
558 | msg = 'Points could not be converted to Numeric array' |
---|
559 | raise msg |
---|
560 | |
---|
561 | try: |
---|
562 | polygon = ensure_numeric(polygon, Float) |
---|
563 | except: |
---|
564 | msg = 'Polygon could not be converted to Numeric array' |
---|
565 | raise msg |
---|
566 | |
---|
567 | assert len(polygon.shape) == 2,\ |
---|
568 | 'Polygon array must be a 2d array of vertices' |
---|
569 | |
---|
570 | assert polygon.shape[1] == 2,\ |
---|
571 | 'Polygon array must have two columns' |
---|
572 | |
---|
573 | assert len(points.shape) == 2,\ |
---|
574 | 'Points array must be a 2d array' |
---|
575 | |
---|
576 | assert points.shape[1] == 2,\ |
---|
577 | 'Points array must have two columns' |
---|
578 | |
---|
579 | N = polygon.shape[0] #Number of vertices in polygon |
---|
580 | M = points.shape[0] #Number of points |
---|
581 | |
---|
582 | px = polygon[:,0] |
---|
583 | py = polygon[:,1] |
---|
584 | |
---|
585 | #Used for an optimisation when points are far away from polygon |
---|
586 | minpx = min(px); maxpx = max(px) |
---|
587 | minpy = min(py); maxpy = max(py) |
---|
588 | |
---|
589 | |
---|
590 | #Begin main loop |
---|
591 | indices = zeros(M, Int) |
---|
592 | |
---|
593 | inside_index = 0 #Keep track of points inside |
---|
594 | outside_index = M-1 #Keep track of points outside (starting from end) |
---|
595 | |
---|
596 | for k in range(M): |
---|
597 | |
---|
598 | if verbose: |
---|
599 | if k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M) |
---|
600 | |
---|
601 | #for each point |
---|
602 | x = points[k, 0] |
---|
603 | y = points[k, 1] |
---|
604 | |
---|
605 | inside = False |
---|
606 | |
---|
607 | if not x > maxpx or x < minpx or y > maxpy or y < minpy: |
---|
608 | #Check polygon |
---|
609 | for i in range(N): |
---|
610 | j = (i+1)%N |
---|
611 | |
---|
612 | #Check for case where point is contained in line segment |
---|
613 | if point_on_line(x, y, px[i], py[i], px[j], py[j]): |
---|
614 | if closed: |
---|
615 | inside = True |
---|
616 | else: |
---|
617 | inside = False |
---|
618 | break |
---|
619 | else: |
---|
620 | #Check if truly inside polygon |
---|
621 | if py[i] < y and py[j] >= y or\ |
---|
622 | py[j] < y and py[i] >= y: |
---|
623 | if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x: |
---|
624 | inside = not inside |
---|
625 | |
---|
626 | if inside: |
---|
627 | indices[inside_index] = k |
---|
628 | inside_index += 1 |
---|
629 | else: |
---|
630 | indices[outside_index] = k |
---|
631 | outside_index -= 1 |
---|
632 | |
---|
633 | return indices, inside_index |
---|
634 | |
---|
635 | |
---|
636 | def separate_points_by_polygon_c(points, polygon, |
---|
637 | closed = True, verbose = False): |
---|
638 | """Determine whether points are inside or outside a polygon |
---|
639 | |
---|
640 | C-wrapper |
---|
641 | """ |
---|
642 | |
---|
643 | from Numeric import array, Float, reshape, zeros, Int |
---|
644 | |
---|
645 | |
---|
646 | if verbose: print 'Checking input to separate_points_by_polygon' |
---|
647 | #Input checks |
---|
648 | try: |
---|
649 | points = ensure_numeric(points, Float) |
---|
650 | except: |
---|
651 | msg = 'Points could not be converted to Numeric array' |
---|
652 | raise msg |
---|
653 | |
---|
654 | #if verbose: print 'Checking input to separate_points_by_polygon 2' |
---|
655 | try: |
---|
656 | polygon = ensure_numeric(polygon, Float) |
---|
657 | except: |
---|
658 | msg = 'Polygon could not be converted to Numeric array' |
---|
659 | raise msg |
---|
660 | |
---|
661 | if verbose: print 'check' |
---|
662 | |
---|
663 | assert len(polygon.shape) == 2,\ |
---|
664 | 'Polygon array must be a 2d array of vertices' |
---|
665 | |
---|
666 | assert polygon.shape[1] == 2,\ |
---|
667 | 'Polygon array must have two columns' |
---|
668 | |
---|
669 | assert len(points.shape) == 2,\ |
---|
670 | 'Points array must be a 2d array' |
---|
671 | |
---|
672 | assert points.shape[1] == 2,\ |
---|
673 | 'Points array must have two columns' |
---|
674 | |
---|
675 | N = polygon.shape[0] #Number of vertices in polygon |
---|
676 | M = points.shape[0] #Number of points |
---|
677 | |
---|
678 | from util_ext import separate_points_by_polygon |
---|
679 | |
---|
680 | if verbose: print 'Allocating array for indices' |
---|
681 | |
---|
682 | indices = zeros( M, Int ) |
---|
683 | |
---|
684 | if verbose: print 'Calling C-version of inside poly' |
---|
685 | count = separate_points_by_polygon(points, polygon, indices, |
---|
686 | int(closed), int(verbose)) |
---|
687 | |
---|
688 | return indices, count |
---|
689 | |
---|
690 | |
---|
691 | |
---|
692 | class Polygon_function: |
---|
693 | """Create callable object f: x,y -> z, where a,y,z are vectors and |
---|
694 | where f will return different values depending on whether x,y belongs |
---|
695 | to specified polygons. |
---|
696 | |
---|
697 | To instantiate: |
---|
698 | |
---|
699 | Polygon_function(polygons) |
---|
700 | |
---|
701 | where polygons is a list of tuples of the form |
---|
702 | |
---|
703 | [ (P0, v0), (P1, v1), ...] |
---|
704 | |
---|
705 | with Pi being lists of vertices defining polygons and vi either |
---|
706 | constants or functions of x,y to be applied to points with the polygon. |
---|
707 | |
---|
708 | The function takes an optional argument, default which is the value |
---|
709 | (or function) to used for points not belonging to any polygon. |
---|
710 | For example: |
---|
711 | |
---|
712 | Polygon_function(polygons, default = 0.03) |
---|
713 | |
---|
714 | If omitted the default value will be 0.0 |
---|
715 | |
---|
716 | Note: If two polygons overlap, the one last in the list takes precedence |
---|
717 | |
---|
718 | 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??? |
---|
719 | |
---|
720 | """ |
---|
721 | |
---|
722 | def __init__(self, regions, default = 0.0): |
---|
723 | |
---|
724 | try: |
---|
725 | len(regions) |
---|
726 | except: |
---|
727 | msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons |
---|
728 | raise msg |
---|
729 | |
---|
730 | |
---|
731 | T = regions[0] |
---|
732 | try: |
---|
733 | a = len(T) |
---|
734 | except: |
---|
735 | msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons |
---|
736 | raise msg |
---|
737 | |
---|
738 | assert a == 2, 'Must have two component each: %s' %T |
---|
739 | |
---|
740 | self.regions = regions |
---|
741 | self.default = default |
---|
742 | |
---|
743 | |
---|
744 | def __call__(self, x, y): |
---|
745 | from util import inside_polygon |
---|
746 | from Numeric import ones, Float, concatenate, array, reshape, choose |
---|
747 | |
---|
748 | x = array(x).astype(Float) |
---|
749 | y = array(y).astype(Float) |
---|
750 | |
---|
751 | N = len(x) |
---|
752 | assert len(y) == N |
---|
753 | |
---|
754 | points = concatenate( (reshape(x, (N, 1)), |
---|
755 | reshape(y, (N, 1))), axis=1 ) |
---|
756 | |
---|
757 | if callable(self.default): |
---|
758 | z = self.default(x,y) |
---|
759 | else: |
---|
760 | z = ones(N, Float) * self.default |
---|
761 | |
---|
762 | for polygon, value in self.regions: |
---|
763 | indices = inside_polygon(points, polygon) |
---|
764 | |
---|
765 | #FIXME: This needs to be vectorised |
---|
766 | if callable(value): |
---|
767 | for i in indices: |
---|
768 | xx = array([x[i]]) |
---|
769 | yy = array([y[i]]) |
---|
770 | z[i] = value(xx, yy)[0] |
---|
771 | else: |
---|
772 | for i in indices: |
---|
773 | z[i] = value |
---|
774 | |
---|
775 | return z |
---|
776 | |
---|
777 | def read_polygon(filename): |
---|
778 | """Read points assumed to form a polygon |
---|
779 | There must be exactly two numbers in each line |
---|
780 | """ |
---|
781 | |
---|
782 | #Get polygon |
---|
783 | fid = open(filename) |
---|
784 | lines = fid.readlines() |
---|
785 | fid.close() |
---|
786 | polygon = [] |
---|
787 | for line in lines: |
---|
788 | fields = line.split(',') |
---|
789 | polygon.append( [float(fields[0]), float(fields[1])] ) |
---|
790 | |
---|
791 | return polygon |
---|
792 | |
---|
793 | def populate_polygon(polygon, number_of_points, seed = None, exclude = None): |
---|
794 | """Populate given polygon with uniformly distributed points. |
---|
795 | |
---|
796 | Input: |
---|
797 | polygon - list of vertices of polygon |
---|
798 | number_of_points - (optional) number of points |
---|
799 | seed - seed for random number generator (default=None) |
---|
800 | exclude - list of polygons (inside main polygon) from where points should be excluded |
---|
801 | |
---|
802 | Output: |
---|
803 | points - list of points inside polygon |
---|
804 | |
---|
805 | Examples: |
---|
806 | populate_polygon( [[0,0], [1,0], [1,1], [0,1]], 5 ) |
---|
807 | will return five randomly selected points inside the unit square |
---|
808 | """ |
---|
809 | |
---|
810 | from random import uniform, seed |
---|
811 | |
---|
812 | seed(seed) |
---|
813 | |
---|
814 | points = [] |
---|
815 | |
---|
816 | #Find outer extent of polygon |
---|
817 | max_x = min_x = polygon[0][0] |
---|
818 | max_y = min_y = polygon[0][1] |
---|
819 | for point in polygon[1:]: |
---|
820 | x = point[0] |
---|
821 | if x > max_x: max_x = x |
---|
822 | if x < min_x: min_x = x |
---|
823 | y = point[1] |
---|
824 | if y > max_y: max_y = y |
---|
825 | if y < min_y: min_y = y |
---|
826 | |
---|
827 | |
---|
828 | while len(points) < number_of_points: |
---|
829 | x = uniform(min_x, max_x) |
---|
830 | y = uniform(min_y, max_y) |
---|
831 | |
---|
832 | append = False |
---|
833 | if inside_polygon( [x,y], polygon ): |
---|
834 | |
---|
835 | append = True |
---|
836 | |
---|
837 | #Check exclusions |
---|
838 | if exclude is not None: |
---|
839 | for ex_poly in exclude: |
---|
840 | if inside_polygon( [x,y], ex_poly ): |
---|
841 | append = False |
---|
842 | |
---|
843 | |
---|
844 | if append is True: |
---|
845 | points.append([x,y]) |
---|
846 | |
---|
847 | return points |
---|
848 | |
---|
849 | |
---|
850 | |
---|
851 | def tsh2sww(filename, verbose=True): |
---|
852 | """ |
---|
853 | to check if a tsh/msh file 'looks' good. |
---|
854 | """ |
---|
855 | |
---|
856 | #FIXME Move to data_manager |
---|
857 | from shallow_water import Domain |
---|
858 | from pmesh2domain import pmesh_to_domain_instance |
---|
859 | import time, os |
---|
860 | from data_manager import get_dataobject |
---|
861 | from os import sep, path |
---|
862 | #from util import mean |
---|
863 | |
---|
864 | if verbose == True:print 'Creating domain from', filename |
---|
865 | domain = pmesh_to_domain_instance(filename, Domain) |
---|
866 | if verbose == True:print "Number of triangles = ", len(domain) |
---|
867 | |
---|
868 | domain.smooth = True |
---|
869 | domain.format = 'sww' #Native netcdf visualisation format |
---|
870 | file_path, filename = path.split(filename) |
---|
871 | filename, ext = path.splitext(filename) |
---|
872 | domain.filename = filename |
---|
873 | domain.reduction = mean |
---|
874 | if verbose == True:print "file_path",file_path |
---|
875 | if file_path == "":file_path = "." |
---|
876 | domain.set_datadir(file_path) |
---|
877 | |
---|
878 | if verbose == True: |
---|
879 | print "Output written to " + domain.get_datadir() + sep + \ |
---|
880 | domain.filename + "." + domain.format |
---|
881 | sww = get_dataobject(domain) |
---|
882 | sww.store_connectivity() |
---|
883 | sww.store_timestep('stage') |
---|
884 | |
---|
885 | #################################################################### |
---|
886 | #Python versions of function that are also implemented in util_gateway.c |
---|
887 | # |
---|
888 | |
---|
889 | def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2): |
---|
890 | """ |
---|
891 | """ |
---|
892 | |
---|
893 | det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0) |
---|
894 | a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0) |
---|
895 | a /= det |
---|
896 | |
---|
897 | b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0) |
---|
898 | b /= det |
---|
899 | |
---|
900 | return a, b |
---|
901 | |
---|
902 | |
---|
903 | def gradient2_python(x0, y0, x1, y1, q0, q1): |
---|
904 | """Compute radient based on two points and enforce zero gradient |
---|
905 | in the direction orthogonal to (x1-x0), (y1-y0) |
---|
906 | """ |
---|
907 | |
---|
908 | #Old code |
---|
909 | #det = x0*y1 - x1*y0 |
---|
910 | #if det != 0.0: |
---|
911 | # a = (y1*q0 - y0*q1)/det |
---|
912 | # b = (x0*q1 - x1*q0)/det |
---|
913 | |
---|
914 | #Correct code (ON) |
---|
915 | det = (x1-x0)**2 + (y1-y0)**2 |
---|
916 | if det != 0.0: |
---|
917 | a = (x1-x0)*(q1-q0)/det |
---|
918 | b = (y1-y0)*(q1-q0)/det |
---|
919 | |
---|
920 | return a, b |
---|
921 | |
---|
922 | |
---|
923 | ############################################## |
---|
924 | #Initialise module |
---|
925 | |
---|
926 | import compile |
---|
927 | if compile.can_use_C_extension('util_ext.c'): |
---|
928 | from util_ext import gradient, gradient2, point_on_line |
---|
929 | separate_points_by_polygon = separate_points_by_polygon_c |
---|
930 | else: |
---|
931 | gradient = gradient_python |
---|
932 | gradient2 = gradient2_python |
---|
933 | |
---|
934 | |
---|
935 | if __name__ == "__main__": |
---|
936 | pass |
---|
937 | |
---|
938 | |
---|
939 | |
---|
940 | |
---|