Changeset 7858
- Timestamp:
- Jun 18, 2010, 4:43:10 PM (15 years ago)
- Location:
- trunk/anuga_core/source/anuga
- Files:
-
- 1 added
- 14 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_core/source/anuga/__init__.py
r7851 r7858 42 42 43 43 from anuga.geometry.polygon import read_polygon, plot_polygons, polygon_area 44 from anuga.geometry.polygon import Polygon_function44 from anuga.geometry.polygon_function import Polygon_function 45 45 46 46 from anuga.abstract_2d_finite_volumes.pmesh2domain import \ … … 52 52 from anuga.utilities.file_utils import copy_code_files 53 53 54 from anuga.geometry.polygon import read_polygon , Polygon_function54 from anuga.geometry.polygon import read_polygon 55 55 from anuga.caching import cache 56 56 -
trunk/anuga_core/source/anuga/culvert_flows/test_culvert_class.py
r7778 r7858 7 7 8 8 from anuga.utilities.system_tools import get_pathname_from_package 9 from anuga.geometry.polygon import Polygon_function9 from anuga.geometry.polygon_function import Polygon_function 10 10 11 11 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross -
trunk/anuga_core/source/anuga/culvert_flows/test_new_culvert_class.py
r7778 r7858 9 9 10 10 from anuga.utilities.system_tools import get_pathname_from_package 11 from anuga.geometry.polygon import Polygon_function11 from anuga.geometry.polygon_function import Polygon_function 12 12 13 13 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross -
trunk/anuga_core/source/anuga/file/csv_file.py
r7854 r7858 11 11 """ 12 12 13 13 14 import csv 14 15 import numpy as num 16 import anuga.utilities.log as log 15 17 16 18 … … 168 170 if completed: 169 171 try: 170 file = str(kwargs['file_name'])172 file_name = str(kwargs['file_name']) 171 173 except: 172 raise 'kwargs must have file_name'174 raise Exception('kwargs must have file_name') 173 175 else: 174 176 # write temp file in output directory 175 177 try: 176 file = str(kwargs['output_dir']) + 'detail_temp.csv'178 file_name = str(kwargs['output_dir']) + 'detail_temp.csv' 177 179 except: 178 raise 'kwargs must have output_dir'180 raise Exception('kwargs must have output_dir') 179 181 180 182 # extracts the header info and the new line info … … 199 201 # try to open! 200 202 try: 201 fid = open(file , 'r')203 fid = open(file_name, 'r') 202 204 file_header = fid.readline() 203 205 fid.close() 204 206 if verbose: log.critical('read file header %s' % file_header) 205 except: 206 msg = 'try to create new file: %s' % file 207 if verbose: log.critical(msg) 207 except Exception: 208 msg = 'try to create new file: %s' % file_name 209 if verbose: 210 log.critical(msg) 208 211 #tries to open file, maybe directory is bad 209 212 try: 210 fid = open(file , 'w')213 fid = open(file_name, 'w') 211 214 fid.write(header) 212 215 fid.close() … … 218 221 # if header is same or this is a new file 219 222 if file_header == str(header): 220 fid = open(file , 'a')223 fid = open(file_name, 'a') 221 224 fid.write(line) 222 225 fid.close() … … 225 228 # if header is different and has completed will append info to 226 229 # end of details_temp.cvs file in output directory 227 file = str(kwargs['output_dir']) + 'detail_temp.csv'228 fid = open(file , 'a')230 file_name = str(kwargs['output_dir']) + 'detail_temp.csv' 231 fid = open(file_name, 'a') 229 232 fid.write(header) 230 233 fid.write(line) … … 244 247 245 248 def load_csv_as_building_polygons(file_name, 246 floor_height=3, 247 clipping_polygons=None): 249 floor_height=3): 248 250 """ 249 251 Convert CSV files of the form: … … 285 287 286 288 287 ##288 # @brief Convert CSV file into a dictionary of polygons and associated values.289 # @param filename The path to the file to read, value_name name for the 4th column290 289 def load_csv_as_polygons(file_name, 291 290 value_name='value', … … 338 337 339 338 msg = 'Did not find expected column header: northing' 340 assert 'northing' in X.keys(), northing339 assert 'northing' in X.keys(), msg 341 340 342 341 msg = 'Did not find expected column header: northing' … … 357 356 past_ids = {} 358 357 last_id = None 359 for i, id in enumerate(X['id']):358 for i, poly_id in enumerate(X['id']): 360 359 361 360 # Check for duplicate polygons 362 if id in past_ids:361 if poly_id in past_ids: 363 362 msg = 'Polygon %s was duplicated in line %d' % (id, i) 364 363 raise Exception, msg 365 364 366 if id not in polygons:365 if poly_id not in polygons: 367 366 # Start new polygon 368 polygons[ id] = []367 polygons[poly_id] = [] 369 368 if values is not None: 370 values[ id] = X[value_name][i]369 values[poly_id] = X[value_name][i] 371 370 372 371 # Keep track of previous polygon ids … … 385 384 386 385 if exclude is True: 387 excluded_polygons[ id]=True388 389 polygons[ id].append(point)386 excluded_polygons[poly_id]=True 387 388 polygons[poly_id].append(point) 390 389 391 390 # Check that value is the same across each polygon 392 391 msg = 'Values must be the same across each polygon.' 393 msg += 'I got %s in line %d but it should have been %s' % (X[value_name][i], i, values[id]) 394 assert values[id] == X[value_name][i], msg 395 396 last_id = id 392 msg += 'I got %s in line %d but it should have been %s' % \ 393 (X[value_name][i], i, values[poly_id]) 394 assert values[poly_id] == X[value_name][i], msg 395 396 last_id = poly_id 397 397 398 398 # Weed out polygons that were not wholly inside clipping polygons 399 for id in excluded_polygons:400 del polygons[ id]399 for poly_id in excluded_polygons: 400 del polygons[poly_id] 401 401 402 402 return polygons, values -
trunk/anuga_core/source/anuga/file/sts.py
r7778 r7858 94 94 outfile.variables[q + Write_sts.RANGE][1] = -max_float # Max 95 95 96 # Doing sts_precision instead of Float gives cast errors. 97 outfile.createVariable('time', netcdf_float, ('number_of_timesteps',)) 98 99 for q in Write_sts.sts_quantities: 100 outfile.createVariable(q, sts_precision, ('number_of_timesteps', 101 'number_of_points')) 102 outfile.createVariable(q + Write_sts.RANGE, sts_precision, 103 ('numbers_in_range',)) 104 # Initialise ranges with small and large sentinels. 105 # If this was in pure Python we could have used None sensibly 106 outfile.variables[q + Write_sts.RANGE][0] = max_float # Min 107 outfile.variables[q + Write_sts.RANGE][1] = -max_float # Max 108 109 if isinstance(times, (list, num.ndarray)): 110 outfile.variables['time'][:] = times #Store time relative 111 112 if verbose: 113 log.critical('------------------------------------------------') 114 log.critical('Statistics:') 115 log.critical(' t in [%f, %f], len(t) == %d' 116 % (num.min(times), num.max(times), len(times.flat))) 96 self.write_dynamic_quantities(outfile, Write_sts.sts_quantities, times) 117 97 118 98 ## … … 265 245 266 246 247 def write_dynamic_quantities(self, outfile, quantities, 248 times, precis = netcdf_float32, verbose = False): 249 """ 250 Write out given quantities to file. 251 """ 252 for q in quantities: 253 outfile.createVariable(q, precis, ('number_of_timesteps', 254 'number_of_points')) 255 outfile.createVariable(q + Write_sts.RANGE, precis, 256 ('numbers_in_range',)) 257 258 # Initialise ranges with small and large sentinels. 259 # If this was in pure Python we could have used None sensibly 260 outfile.variables[q+Write_sts.RANGE][0] = max_float # Min 261 outfile.variables[q+Write_sts.RANGE][1] = -max_float # Max 262 263 # Doing sts_precision instead of Float gives cast errors. 264 outfile.createVariable('time', netcdf_float, ('number_of_timesteps',)) 265 266 if isinstance(times, (list, num.ndarray)): 267 outfile.variables['time'][:] = times # Store time relative 268 269 if verbose: 270 log.critical('------------------------------------------------') 271 log.critical('Statistics:') 272 log.critical(' t in [%f, %f], len(t) == %d' 273 % (num.min(times), num.max(times), len(times.flat))) 274 267 275 268 276 -
trunk/anuga_core/source/anuga/file/sww.py
r7841 r7858 9 9 import anuga.utilities.log as log 10 10 from Scientific.IO.NetCDF import NetCDFFile 11 12 from sts import Write_sts 11 13 12 14 from anuga.coordinate_transforms.geo_reference import \ … … 455 457 456 458 457 class Write_sww :459 class Write_sww(Write_sts): 458 460 """ 459 461 A class to write an SWW file. … … 462 464 manually. 463 465 """ 464 RANGE = '_range'465 EXTREMA = ':extrema'466 466 467 467 def __init__(self, static_quantities, dynamic_quantities): … … 567 567 'number_of_vertices')) 568 568 569 # Doing sww_precision instead of Float gives cast errors. 570 outfile.createVariable('time', netcdf_float, 571 ('number_of_timesteps',)) 572 573 569 574 570 for q in self.static_quantities: 575 571 … … 585 581 outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max 586 582 587 #if 'elevation' in self.static_quantities: 588 # # FIXME: Backwards compat - get rid of z once old view has retired 589 # outfile.createVariable('z', sww_precision, 590 # ('number_of_points',)) 591 592 for q in self.dynamic_quantities: 593 outfile.createVariable(q, sww_precision, ('number_of_timesteps', 594 'number_of_points')) 595 outfile.createVariable(q + Write_sww.RANGE, sww_precision, 596 ('numbers_in_range',)) 597 598 # Initialise ranges with small and large sentinels. 599 # If this was in pure Python we could have used None sensibly 600 outfile.variables[q+Write_sww.RANGE][0] = max_float # Min 601 outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max 602 603 if isinstance(times, (list, num.ndarray)): 604 outfile.variables['time'][:] = times # Store time relative 605 606 if verbose: 607 log.critical('------------------------------------------------') 608 log.critical('Statistics:') 609 log.critical(' t in [%f, %f], len(t) == %d' 610 % (num.min(times), num.max(times), len(times.flat))) 583 584 self.write_dynamic_quantities(outfile, self.dynamic_quantities, times) 585 611 586 612 587 -
trunk/anuga_core/source/anuga/file_conversion/csv2sts.py
r7853 r7858 3 3 """ 4 4 Module to convert a .csv file to an .sts file. 5 6 You can use the function from your own scripts, or call this script on the 7 command line. 5 8 6 9 Use this module to convert an arbitrary csv file to an sts file. … … 42 45 data: 43 46 44 stage = 4, 150.66667, 150.83334, 151, 151.16667, -34, -34.16667, -34.33333,45 -34.5, -1, -5, -9, -13 ;47 stage = 4, 150.66667, 150.83334, 151, 151.16667, -34, -34.16667, 48 -34.33333, -34.5, -1, -5, -9, -13 ; 46 49 47 50 time = 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 ; 48 51 } 49 52 53 As of June 2010 this module has a pylint quality rating of 9.30/10. 50 54 """ 51 55 52 import csv53 import os54 56 import sys 55 57 import getopt 56 58 from anuga.utilities import log 57 import numpy as num58 59 from Scientific.IO.NetCDF import NetCDFFile 59 60 from anuga.file.csv_file import load_csv_as_dict … … 108 109 109 110 def usage(): 111 """ Display usage of this module from the comand line. """ 110 112 print 'csv2sts - convert a csv file to an sts file.' 111 113 print 'Usage: csv2sts [-hv] [--help] [--verbose]', … … 161 163 162 164 if __name__ == "__main__": 163 """ Entry point if run from command line """164 165 main(sys.argv[1:]) -
trunk/anuga_core/source/anuga/file_conversion/urs2sww.py
r7800 r7858 24 24 from anuga.file.sww import Write_sww 25 25 26 27 ################################################################################ 28 # CONVERTING UNGRIDDED URS DATA TO AN SWW FILE 29 ################################################################################ 30 31 ## 32 # @brief Convert URS file(s) (wave prop) to an SWW file. 33 # @param basename_in Stem of the input filenames. 34 # @param basename_out Path to the output SWW file. 35 # @param verbose True if this function is to be verbose. 36 # @param mint 37 # @param maxt 38 # @param mean_stage 39 # @param origin Tuple with geo-ref UTM coordinates (zone, easting, northing). 40 # @param hole_points_UTM 41 # @param zscale 42 # @note Also convert latitude and longitude to UTM. All coordinates are 43 # assumed to be given in the GDA94 datum. 44 # @note Input filename stem has suffixes '-z-mux', '-e-mux' and '-n-mux' 45 # added for relative height, x-velocity and y-velocity respectively. 26 ############################################################### 27 46 28 def urs_ungridded2sww(basename_in='o', basename_out=None, verbose=False, 47 29 mint=None, maxt=None, -
trunk/anuga_core/source/anuga/geometry/aabb.py
r7751 r7858 3 3 Contains a class describing a bounding box. It contains methods 4 4 to split itself and return child boxes. 5 6 As of June 2010 this module has a pylint quality rating of 10/10. 5 7 """ 6 8 … … 15 17 """ 16 18 17 def __init__(self, xmin, xmax , ymin, ymax):19 def __init__(self, xmin, xmax=None, ymin=None, ymax=None): 18 20 """ Define axially-algned bounding box. 19 21 xmin is minimum x … … 22 24 ymax is maximum y (absolute coord, ie, not size) 23 25 """ 24 self.xmin = xmin 25 self.xmax = xmax 26 self.ymin = ymin 27 self.ymax = ymax 26 if not xmax: 27 # try treating first arg as a list of points 28 try: 29 xmin[0][0] 30 except: 31 raise Exception('Single parameter to AABB must be point list.') 32 33 self.xmin, self.ymin = self.xmax, self.ymax = xmin[0] 34 self.include(xmin[1:]) 35 else: 36 self.xmin = xmin 37 self.xmax = xmax 38 self.ymin = ymin 39 self.ymax = ymax 28 40 29 41 … … 96 108 and (self.ymin <= point[1] <= self.ymax) 97 109 110 def include(self, point_list): 111 """ Include points in AABB. 112 Bounding box will be expanded to include passed points 113 114 point_list is a list of points. 115 """ 116 for point in point_list: 117 pt_x, pt_y = point 118 if pt_x < self.xmin: 119 self.xmin = pt_x 120 if pt_x > self.xmax: 121 self.xmax = pt_x 122 if pt_y < self.ymin: 123 self.ymin = pt_y 124 if pt_y > self.ymax: 125 self.ymax = pt_y -
trunk/anuga_core/source/anuga/geometry/polygon.py
r7841 r7858 10 10 import anuga.utilities.log as log 11 11 12 from aabb import AABB 12 13 13 14 ## … … 36 37 37 38 res = _point_on_line(point[0], point[1], 38 line[0, 0], line[0,1],39 line[1, 0], line[1,1],39 line[0, 0], line[0, 1], 40 line[1, 0], line[1, 1], 40 41 rtol, atol) 41 42 … … 49 50 50 51 # result functions for possible states 51 def lines_dont_coincide(p0,p1,p2,p3): return (3, None) 52 def lines_0_fully_included_in_1(p0,p1,p2,p3): return (2, 53 num.array([p0,p1])) 54 def lines_1_fully_included_in_0(p0,p1,p2,p3): return (2, 55 num.array([p2,p3])) 56 def lines_overlap_same_direction(p0,p1,p2,p3): return (2, 57 num.array([p0,p3])) 58 def lines_overlap_same_direction2(p0,p1,p2,p3): return (2, 59 num.array([p2,p1])) 60 def lines_overlap_opposite_direction(p0,p1,p2,p3): return (2, 61 num.array([p0,p2])) 62 def lines_overlap_opposite_direction2(p0,p1,p2,p3): return (2, 63 num.array([p3,p1])) 52 def lines_dont_coincide(p0, p1, p2, p3): 53 return (3, None) 54 55 def lines_0_fully_included_in_1(p0, p1, p2, p3): 56 return (2, num.array([p0, p1])) 57 58 def lines_1_fully_included_in_0(p0, p1, p2, p3): 59 return (2, num.array([p2, p3])) 60 61 def lines_overlap_same_direction(p0, p1, p2, p3): 62 return (2, num.array([p0, p3])) 63 64 def lines_overlap_same_direction2(p0, p1, p2, p3): 65 return (2, num.array([p2, p1])) 66 67 def lines_overlap_opposite_direction(p0, p1, p2, p3): 68 return (2, num.array([p0, p2])) 69 70 def lines_overlap_opposite_direction2(p0, p1, p2, p3): 71 return (2, num.array([p3, p1])) 64 72 65 73 # this function called when an impossible state is found … … 146 154 point_on_line([x3, y3], line0, rtol=rtol, atol=atol)) 147 155 148 return collinear_result[state_tuple]([x0, y0], [x1,y1],149 [x2, y2], [x3,y3])156 return collinear_result[state_tuple]([x0, y0], [x1, y1], 157 [x2, y2], [x3, y3]) 150 158 else: 151 159 # Lines are parallel but aren't collinear … … 269 277 return bool(_is_inside_triangle(point, triangle, int(closed), rtol, atol)) 270 278 271 272 273 # FIXME (Ole): The rest of this function has been made274 # obsolete by the C extension.275 276 # Quickly reject points that are clearly outside277 if point[0] < min(triangle[:, 0]):278 return False279 if point[0] > max(triangle[:, 0]):280 return False281 if point[1] < min(triangle[:, 1]):282 return False283 if point[1] > max(triangle[:, 1]):284 return False285 286 287 # Start search288 A = triangle[0, :]289 B = triangle[1, :]290 C = triangle[2, :]291 292 # Now check if point lies wholly inside triangle293 v0 = C-A294 v1 = B-A295 296 a00 = num.inner(v0, v0)297 a10 = a01 = num.inner(v0, v1)298 a11 = num.inner(v1, v1)299 300 denom = a11*a00 - a01*a10301 302 if abs(denom) > 0.0:303 v = point-A304 b0 = num.inner(v0, v)305 b1 = num.inner(v1, v)306 307 alpha = (b0*a11 - b1*a01)/denom308 beta = (b1*a00 - b0*a10)/denom309 310 if (alpha > 0.0) and (beta > 0.0) and (alpha+beta < 1.0):311 return True312 313 314 if closed is True:315 # Check if point lies on one of the edges316 317 for X, Y in [[A, B], [B, C], [C, A]]:318 res = _point_on_line(point[0], point[1],319 X[0], X[1],320 Y[0], Y[1],321 rtol, atol)322 323 if res:324 return True325 326 return False327 328 329 279 def is_complex(polygon, verbose=False): 330 280 """Check if a polygon is complex (self-intersecting). … … 339 289 340 290 def key_xpos(item): 291 """ Return the x coord out of the passed point for sorting key. """ 341 292 return (item[0][0]) 342 293 … … 377 328 print 'Self-intersecting polygon found, type ', type 378 329 print 'point', point, 379 print 'vertices: ', leftmost, ' - ', l_x[cmp] 330 print 'vertices: ', leftmost, ' - ', l_x[cmp] 380 331 return True 381 332 cmp += 1 … … 422 373 try: 423 374 points = ensure_absolute(points) 424 except NameError, e :425 raise NameError, e 375 except NameError, err: 376 raise NameError, err 426 377 except: 427 378 # If this fails it is going to be because the points can't be … … 513 464 if len(points.shape) == 1: 514 465 # Only one point was passed in. Convert to array of points 515 points = num.reshape(points, (1, 2))466 points = num.reshape(points, (1, 2)) 516 467 517 468 indices, count = separate_points_by_polygon(points, polygon, … … 559 510 if len(points.shape) == 1: 560 511 # Only one point was passed in. Convert to array of points 561 points = num.reshape(points, (1, 2))512 points = num.reshape(points, (1, 2)) 562 513 563 514 indices, count = separate_points_by_polygon(points, polygon, … … 632 583 except: 633 584 msg = 'Points could not be converted to numeric array' 634 raise msg585 raise Exception(msg) 635 586 636 587 try: 637 588 polygon = ensure_numeric(polygon, num.float) 638 589 except NameError, e: 639 raise NameError , e590 raise NameError(e) 640 591 except: 641 592 msg = 'Polygon could not be converted to numeric array' 642 raise msg593 raise Exception(msg) 643 594 644 595 msg = 'Polygon array must be a 2d array of vertices' … … 646 597 647 598 msg = 'Polygon array must have two columns' 648 assert polygon.shape[1] ==2, msg599 assert polygon.shape[1] == 2, msg 649 600 650 601 msg = ('Points array must be 1 or 2 dimensional. ' … … 654 605 if len(points.shape) == 1: 655 606 # Only one point was passed in. Convert to array of points. 656 points = num.reshape(points, (1, 2))607 points = num.reshape(points, (1, 2)) 657 608 658 609 msg = ('Point array must have two columns (x,y), ' … … 663 614 msg = ('Points array must be a 2d array. I got %s.' 664 615 % str(points[:30])) 665 assert len(points.shape) ==2, msg616 assert len(points.shape) == 2, msg 666 617 667 618 msg = 'Points array must have two columns' 668 assert points.shape[1] ==2, msg619 assert points.shape[1] == 2, msg 669 620 670 621 N = polygon.shape[0] # Number of vertices in polygon … … 681 632 return indices, count 682 633 683 ## 684 # @brief Determine area of a polygon. 685 # @param input_polygon The polygon to get area of. 686 # @return A scalar value for the polygon area. 634 687 635 def polygon_area(input_polygon): 688 636 """ Determine area of arbitrary polygon. 689 637 690 Reference: http://mathworld.wolfram.com/PolygonArea.html 691 """ 692 638 input_polygon The polygon to get area of. 639 640 return A scalar value for the polygon area. 641 642 Reference: http://mathworld.wolfram.com/PolygonArea.html 643 """ 693 644 # Move polygon to origin (0,0) to avoid rounding errors 694 645 # This makes a copy of the polygon to avoid destroying it 695 646 input_polygon = ensure_numeric(input_polygon) 696 min_x = min(input_polygon[:, 0])697 min_y = min(input_polygon[:, 1])647 min_x = min(input_polygon[:, 0]) 648 min_y = min(input_polygon[:, 1]) 698 649 polygon = input_polygon - [min_x, min_y] 699 650 … … 742 693 Outputs: 743 694 744 - list of min and max of x and y coordinates745 695 - plot of polygons 746 696 """ … … 754 704 ion() 755 705 hold(True) 756 757 minx = 1e10758 maxx = 0.0759 miny = 1e10760 maxy = 0.0761 706 762 707 if label is None: … … 772 717 alpha = max(0.0, min(1.0, alpha)) 773 718 774 n = len(polygons_points)719 num_points = len(polygons_points) 775 720 colour = [] 776 721 if style is None: 777 722 style_type = 'line' 778 723 style = [] 779 for i in range(n ):724 for i in range(num_points): 780 725 style.append(style_type) 781 726 colour.append('b-') 782 727 else: 783 for s in style: 784 if s == 'line': colour.append('b-') 785 if s == 'outside': colour.append('r.') 786 if s == 'point': colour.append('g.') 787 if s != 'line': 788 if s != 'outside': 789 if s != 'point': 790 colour.append(s) 728 for style_name in style: 729 if style_name == 'line': 730 colour.append('b-') 731 if style_name == 'outside': 732 colour.append('r.') 733 if style_name == 'point': 734 colour.append('g.') 735 if style_name not in ['line', 'outside', 'point']: 736 colour.append(style_name) 791 737 792 738 for i, item in enumerate(polygons_points): 793 x, y = poly_xy(item) 794 if min(x) < minx: 795 minx = min(x) 796 if max(x) > maxx: 797 maxx = max(x) 798 if min(y) < miny: 799 miny = min(y) 800 if max(y) > maxy: 801 maxy = max(y) 802 plot(x, y, colour[i]) 739 pt_x, pt_y = _poly_xy(item) 740 plot(pt_x, pt_y, colour[i]) 803 741 if alpha: 804 fill( x,y, colour[i], alpha=alpha)742 fill(pt_x, pt_y, colour[i], alpha=alpha) 805 743 xlabel('x') 806 744 ylabel('y') … … 814 752 close('all') 815 753 816 vec = [minx, maxx, miny, maxy] 817 return vec 818 819 820 def poly_xy(polygon): 754 755 def _poly_xy(polygon): 821 756 """ this is used within plot_polygons so need to duplicate 822 757 the first point so can have closed polygon in plot … … 829 764 try: 830 765 polygon = ensure_numeric(polygon, num.float) 831 except NameError, e :832 raise NameError, e 766 except NameError, err: 767 raise NameError, err 833 768 except: 834 769 msg = ('Polygon %s could not be converted to numeric array' … … 836 771 raise Exception, msg 837 772 838 x = polygon[:, 0] 839 y = polygon[:, 1] 840 x = num.concatenate((x, [polygon[0, 0]]), axis = 0) 841 y = num.concatenate((y, [polygon[0, 1]]), axis = 0) 842 843 return x, y 844 845 846 class Polygon_function: 847 """Create callable object f: x,y -> z, where a,y,z are vectors and 848 where f will return different values depending on whether x,y belongs 849 to specified polygons. 850 851 To instantiate: 852 853 Polygon_function(polygons) 854 855 where polygons is a list of tuples of the form 856 857 [ (P0, v0), (P1, v1), ...] 858 859 with Pi being lists of vertices defining polygons and vi either 860 constants or functions of x,y to be applied to points with the polygon. 861 862 The function takes an optional argument, default which is the value 863 (or function) to used for points not belonging to any polygon. 864 For example: 865 866 Polygon_function(polygons, default = 0.03) 867 868 If omitted the default value will be 0.0 869 870 Note: If two polygons overlap, the one last in the list takes precedence 871 872 Coordinates specified in the call are assumed to be relative to the 873 origin (georeference) e.g. used by domain. 874 By specifying the optional argument georeference, 875 all points are made relative. 876 877 FIXME: This should really work with geo_spatial point sets. 878 """ 879 880 def __init__(self, regions, default=0.0, geo_reference=None): 881 """ 882 # @brief Create instance of a polygon function. 883 # @param regions A list of (x,y) tuples defining a polygon. 884 # @param default Value or function returning value for points outside poly. 885 # @param geo_reference ?? 886 """ 887 888 try: 889 len(regions) 890 except: 891 msg = ('Polygon_function takes a list of pairs (polygon, value).' 892 'Got %s' % str(regions)) 893 raise Exception, msg 894 895 T = regions[0] 896 897 if isinstance(T, basestring): 898 msg = ('You passed in a list of text values into polygon_function ' 899 'instead of a list of pairs (polygon, value): "%s"' 900 % str(T)) 901 raise Exception, msg 902 903 try: 904 a = len(T) 905 except: 906 msg = ('Polygon_function takes a list of pairs (polygon, value). ' 907 'Got %s' % str(T)) 908 raise Exception, msg 909 910 msg = ('Each entry in regions have two components: (polygon, value). ' 911 'I got %s' % str(T)) 912 assert a == 2, msg 913 914 if geo_reference is None: 915 from anuga.coordinate_transforms.geo_reference import Geo_reference 916 geo_reference = Geo_reference() 917 918 self.default = default 919 920 # Make points in polygons relative to geo_reference 921 self.regions = [] 922 for polygon, value in regions: 923 P = geo_reference.change_points_geo_ref(polygon) 924 self.regions.append((P, value)) 925 926 def __call__(self, x, y): 927 """ 928 # @brief Implement the 'callable' property of Polygon_function. 929 # @param x List of x coordinates of points ot interest. 930 # @param y List of y coordinates of points ot interest. 931 """ 932 x = num.array(x, num.float) 933 y = num.array(y, num.float) 934 935 # x and y must be one-dimensional and same length 936 assert len(x.shape) == 1 and len(y.shape) == 1 937 N = x.shape[0] 938 assert y.shape[0] == N 939 940 points = num.ascontiguousarray(num.concatenate((x[:, num.newaxis], 941 y[:, num.newaxis]), 942 axis = 1 )) 943 944 if callable(self.default): 945 z = self.default(x, y) 946 else: 947 z = num.ones(N, num.float) * self.default 948 949 for polygon, value in self.regions: 950 indices = inside_polygon(points, polygon) 951 952 # FIXME: This needs to be vectorised 953 if callable(value): 954 for i in indices: 955 xx = num.array([x[i]]) 956 yy = num.array([y[i]]) 957 z[i] = value(xx, yy)[0] 958 else: 959 for i in indices: 960 z[i] = value 961 962 if len(z) == 0: 963 msg = ('Warning: points provided to Polygon function did not fall ' 964 'within its regions in [%.2f, %.2f], y in [%.2f, %.2f]' 965 % (min(x), max(x), min(y), max(y))) 966 log.critical(msg) 967 968 return z 773 pts_x = num.concatenate((polygon[:, 0], [polygon[0, 0]]), axis = 0) 774 pts_y = num.concatenate((polygon[:, 1], [polygon[0, 1]]), axis = 0) 775 776 return pts_x, pts_y 777 969 778 970 779 ################################################################################ … … 973 782 974 783 def read_polygon(filename, delimiter=','): 975 """Read points assumed to form a polygon. 976 977 # @param filename Path to file containing polygon data. 978 # @param delimiter Delimiter to split polygon data with. 979 # @return A list of point data from the polygon file. 980 981 982 There must be exactly two numbers in each line separated by the delimiter. 983 No header. 784 """ Read points assumed to form a polygon. 785 786 Also checks to make sure polygon is not complex (self-intersecting). 787 788 filename Path to file containing polygon data. 789 delimiter Delimiter to split polygon data with. 790 A list of point data from the polygon file. 791 792 There must be exactly two numbers in each line separated by the delimiter. 793 No header. 984 794 """ 985 795 … … 1002 812 return polygon 1003 813 1004 ## 1005 # @brief Write polygon data to a file. 1006 # @param polygon Polygon points to write to file. 1007 # @param filename Path to file to write. 1008 # @note Delimiter is assumed to be a comma. 814 1009 815 def write_polygon(polygon, filename=None): 1010 816 """Write polygon to csv file. … … 1047 853 1048 854 # Find outer extent of polygon 1049 max_x = min_x = polygon[0][0] 1050 max_y = min_y = polygon[0][1] 1051 for point in polygon[1:]: 1052 x = point[0] 1053 if x > max_x: max_x = x 1054 if x < min_x: min_x = x 1055 y = point[1] 1056 if y > max_y: max_y = y 1057 if y < min_y: min_y = y 1058 855 extents = AABB(polygon) 856 1059 857 while len(points) < number_of_points: 1060 x = uniform(min_x, max_x)1061 y = uniform(min_y, max_y)858 rand_x = uniform(extents.xmin, extents.xmax) 859 rand_y = uniform(extents.ymin, extents.ymax) 1062 860 1063 861 append = False 1064 if is_inside_polygon([ x,y], polygon):862 if is_inside_polygon([rand_x, rand_y], polygon): 1065 863 append = True 1066 864 … … 1068 866 if exclude is not None: 1069 867 for ex_poly in exclude: 1070 if is_inside_polygon([ x,y], ex_poly):868 if is_inside_polygon([rand_x, rand_y], ex_poly): 1071 869 append = False 1072 870 1073 871 if append is True: 1074 points.append([ x,y])872 points.append([rand_x, rand_y]) 1075 873 1076 874 return points … … 1091 889 """ 1092 890 1093 import exceptions1094 1095 class Found(exceptions.Exception):1096 pass1097 1098 891 polygon = ensure_numeric(polygon) 1099 892 1100 point_in = False 1101 while not point_in: 1102 try: 1103 for poly_point in polygon: # [1:]: 1104 for x_mult in range(-1, 2): 1105 for y_mult in range(-1, 2): 1106 x = poly_point[0] 1107 y = poly_point[1] 1108 1109 if x == 0: 1110 x_delta = x_mult * delta 1111 else: 1112 x_delta = x + x_mult*x*delta 1113 1114 if y == 0: 1115 y_delta = y_mult * delta 1116 else: 1117 y_delta = y + y_mult*y*delta 1118 1119 point = [x_delta, y_delta] 1120 1121 if is_inside_polygon(point, polygon, closed=False): 1122 raise Found 1123 except Found: 1124 point_in = True 1125 else: 1126 delta = delta * 0.1 1127 1128 return point 893 while True: 894 for poly_point in polygon: 895 for x_mult in range(-1, 2): 896 for y_mult in range(-1, 2): 897 pt_x, pt_y = poly_point 898 899 if pt_x == 0: 900 x_delta = x_mult * delta 901 else: 902 x_delta = pt_x + x_mult*pt_x*delta 903 904 if pt_y == 0: 905 y_delta = y_mult * delta 906 else: 907 y_delta = pt_y + y_mult*pt_y*delta 908 909 point = [x_delta, y_delta] 910 911 if is_inside_polygon(point, polygon, closed=False): 912 return point 913 delta = delta * 0.1 1129 914 1130 915 … … 1212 997 reduced_polygon = [] 1213 998 for i, point in enumerate(polygon): 1214 x = point[0] 1215 y = point[1] 1216 if x in [min_x, max_x] and y in [min_y, max_y]: 999 if point[0] in [min_x, max_x] and point[1] in [min_y, max_y]: 1217 1000 # Keep 1218 1001 reduced_polygon.append(point) … … 1318 1101 1319 1102 else: 1320 error_msg = 'C implementations could not be accessed by %s.\n ' %__file__1321 error_msg+= 'Make sure compile_all.py has been run as described in '1322 error_msg+= 'the ANUGA installation guide.'1323 raise Exception( error_msg)1103 ERROR_MSG = 'C implementations could not be accessed by %s.\n ' % __file__ 1104 ERROR_MSG += 'Make sure compile_all.py has been run as described in ' 1105 ERROR_MSG += 'the ANUGA installation guide.' 1106 raise Exception(ERROR_MSG) 1324 1107 1325 1108 -
trunk/anuga_core/source/anuga/geometry/quad.py
r7751 r7858 8 8 may then be iterated over with a proper intersection test to detect actual 9 9 geometry intersections. 10 11 As of June 2010 this module has a pylint quality rating of 10/10. 10 12 11 13 """ -
trunk/anuga_core/source/anuga/geometry/test_polygon.py
r7841 r7858 8 8 9 9 from polygon import * 10 from polygon import _poly_xy 11 from polygon_function import Polygon_function 10 12 from anuga.coordinate_transforms.geo_reference import Geo_reference 11 13 from anuga.geospatial_data.geospatial_data import Geospatial_data … … 1570 1572 # Simplest case: Polygon is the unit square 1571 1573 polygon = [[0,0], [1,0], [1,1], [0,1]] 1572 x, y = poly_xy(polygon)1574 x, y = _poly_xy(polygon) 1573 1575 assert len(x) == len(polygon)+1 1574 1576 assert len(y) == len(polygon)+1 … … 1584 1586 # Arbitrary polygon 1585 1587 polygon = [[1,5], [1,1], [100,10], [1,10], [3,6]] 1586 x, y = poly_xy(polygon)1588 x, y = _poly_xy(polygon) 1587 1589 assert len(x) == len(polygon)+1 1588 1590 assert len(y) == len(polygon)+1 … … 1605 1607 polygon1 = [[0,0], [1,0], [1,1], [0,1]] 1606 1608 polygon2 = [[1,1], [2,1], [3,2], [2,2]] 1607 v = plot_polygons([polygon1, polygon2], figname='test1') 1608 assert len(v) == 4 1609 assert v[0] == 0 1610 assert v[1] == 3 1611 assert v[2] == 0 1612 assert v[3] == 2 1609 plot_polygons([polygon1, polygon2], figname='test1') 1613 1610 1614 1611 # Another case 1615 1612 polygon3 = [[1,5], [10,1], [100,10], [50,10], [3,6]] 1616 1613 v = plot_polygons([polygon2,polygon3], figname='test2') 1617 assert len(v) == 4 1618 assert v[0] == 1 1619 assert v[1] == 100 1620 assert v[2] == 1 1621 assert v[3] == 10 1622 1623 os.remove('test1.png') 1624 os.remove('test2.png') 1614 1615 for file in ['test1.png', 'test2.png']: 1616 assert os.access(file, os.R_OK) 1617 os.remove(file) 1618 1625 1619 1626 1620 def test_inside_polygon_geospatial(self): -
trunk/anuga_core/source/anuga/pmesh/view_tsh.py
r7711 r7858 10 10 from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance 11 11 from anuga.pyvolution.util import file_function 12 from anuga.geometry.polygon import Polygon_function, read_polygon 12 from anuga.geometry.polygon import read_polygon 13 from anuga.geometry.polygon_function import Polygon_function 13 14 from realtime_visualisation_new import * 14 15 -
trunk/anuga_core/source/anuga/shallow_water/sww_interrogate.py
r7778 r7858 4 4 5 5 import os 6 import anuga.utilities.log as log 6 7 import numpy as num 7 8
Note: See TracChangeset
for help on using the changeset viewer.