Changeset 6304 for branches/numpy/anuga/utilities/polygon.py
- Timestamp:
- Feb 10, 2009, 11:11:04 AM (15 years ago)
- Location:
- branches/numpy
- Files:
-
- 1 edited
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
branches/numpy/anuga/utilities/polygon.py
r6223 r6304 3 3 """ 4 4 5 import Numericas num5 import numpy as num 6 6 7 7 from math import sqrt 8 8 from anuga.utilities.numerical_tools import ensure_numeric 9 9 from anuga.geospatial_data.geospatial_data import ensure_absolute, Geospatial_data 10 from anuga.config import netcdf_float 10 11 11 12 … … 15 16 Input: 16 17 point is given by [x, y] 17 18 the equivalent 2x2 Numeric array with each row corresponding to a point.18 line is given by [x0, y0], [x1, y1]] or 19 the equivalent 2x2 numeric array with each row corresponding to a point. 19 20 20 21 Output: … … 100 101 # FIXME (Ole): Write this in C 101 102 102 line0 = ensure_numeric(line0, num. Float)103 line1 = ensure_numeric(line1, num. Float)103 line0 = ensure_numeric(line0, num.float) 104 line1 = ensure_numeric(line1, num.float) 104 105 105 106 x0 = line0[0,0]; y0 = line0[0,1] … … 176 177 177 178 178 line0 = ensure_numeric(line0, num. Float)179 line1 = ensure_numeric(line1, num. Float)179 line0 = ensure_numeric(line0, num.float) 180 line1 = ensure_numeric(line1, num.float) 180 181 181 182 status, value = _intersection(line0[0,0], line0[0,1], … … 203 204 else: 204 205 msg = 'is_inside_polygon must be invoked with one point only' 205 raise msg206 raise Exception, msg 206 207 207 208 … … 228 229 # If this fails it is going to be because the points can't be 229 230 # converted to a numeric array. 230 msg = 'Points could not be converted to Numeric array'231 raisemsg231 msg = 'Points could not be converted to numeric array' 232 raise Exception, msg 232 233 233 234 try: … … 238 239 # If this fails it is going to be because the points can't be 239 240 # converted to a numeric array. 240 msg = 'Polygon %s could not be converted to Numeric array' %(str(polygon))241 raisemsg241 msg = 'Polygon %s could not be converted to numeric array' %(str(polygon)) 242 raise Exception, msg 242 243 243 244 if len(points.shape) == 1: 244 245 # Only one point was passed in. Convert to array of points 245 246 points = num.reshape(points, (1,2)) 246 247 247 248 indices, count = separate_points_by_polygon(points, polygon, … … 270 271 else: 271 272 msg = 'is_outside_polygon must be invoked with one point only' 272 raise msg273 raise Exception, msg 273 274 274 275 … … 285 286 #if verbose: print 'Checking input to outside_polygon' 286 287 try: 287 points = ensure_numeric(points, num. Float)288 points = ensure_numeric(points, num.float) 288 289 except NameError, e: 289 290 raise NameError, e 290 291 except: 291 msg = 'Points could not be converted to Numeric array'292 raisemsg292 msg = 'Points could not be converted to numeric array' 293 raise Exception, msg 293 294 294 295 try: 295 polygon = ensure_numeric(polygon, num. Float)296 polygon = ensure_numeric(polygon, num.float) 296 297 except NameError, e: 297 298 raise NameError, e 298 299 except: 299 msg = 'Polygon could not be converted to Numeric array'300 raisemsg300 msg = 'Polygon could not be converted to numeric array' 301 raise Exception, msg 301 302 302 303 303 304 if len(points.shape) == 1: 304 305 # Only one point was passed in. Convert to array of points 305 306 points = num.reshape(points, (1,2)) 306 307 307 308 indices, count = separate_points_by_polygon(points, polygon, … … 327 328 #if verbose: print 'Checking input to outside_polygon' 328 329 try: 329 points = ensure_numeric(points, num. Float)330 points = ensure_numeric(points, num.float) 330 331 except NameError, e: 331 332 raise NameError, e 332 333 except: 333 msg = 'Points could not be converted to Numeric array'334 raisemsg334 msg = 'Points could not be converted to numeric array' 335 raise Exception, msg 335 336 336 337 try: 337 polygon = ensure_numeric(polygon, num. Float)338 polygon = ensure_numeric(polygon, num.float) 338 339 except NameError, e: 339 340 raise NameError, e 340 341 except: 341 msg = 'Polygon could not be converted to Numeric array'342 raisemsg342 msg = 'Polygon could not be converted to numeric array' 343 raise Exception, msg 343 344 344 345 if len(points.shape) == 1: 345 346 # Only one point was passed in. Convert to array of points 346 347 points = num.reshape(points, (1,2)) 347 348 348 349 … … 413 414 414 415 416 # points = ensure_numeric(points, num.float) 415 417 try: 416 points = ensure_numeric(points, num. Float)418 points = ensure_numeric(points, num.float) 417 419 except NameError, e: 418 420 raise NameError, e 419 421 except: 420 msg = 'Points could not be converted to Numeric array'421 raisemsg422 msg = 'Points could not be converted to numeric array' 423 raise Exception, msg 422 424 423 425 #if verbose: print 'Checking input to separate_points_by_polygon 2' 424 426 try: 425 polygon = ensure_numeric(polygon, num. Float)427 polygon = ensure_numeric(polygon, num.float) 426 428 except NameError, e: 427 429 raise NameError, e 428 430 except: 429 msg = 'Polygon could not be converted to Numeric array'430 raisemsg431 msg = 'Polygon could not be converted to numeric array' 432 raise Exception, msg 431 433 432 434 msg = 'Polygon array must be a 2d array of vertices' … … 449 451 450 452 msg = 'Point array must have two columns (x,y), ' 451 msg += 'I got points.shape[1] == %d' % points.shape[0]453 msg += 'I got points.shape[1] == %d' % points.shape[1] 452 454 assert points.shape[1] == 2, msg 453 455 … … 464 466 465 467 466 indices = num.zeros( M, num. Int )468 indices = num.zeros( M, num.int ) 467 469 468 470 count = _separate_points_by_polygon(points, polygon, indices, … … 599 601 600 602 try: 601 polygon = ensure_numeric(polygon, num. Float)603 polygon = ensure_numeric(polygon, num.float) 602 604 except NameError, e: 603 605 raise NameError, e 604 606 except: 605 msg = 'Polygon %s could not be converted to Numeric array' %(str(polygon))606 raise msg607 msg = 'Polygon %s could not be converted to numeric array' %(str(polygon)) 608 raise Exception, msg 607 609 608 610 x = polygon[:,0] … … 666 668 geo_reference=None): 667 669 668 669 670 670 try: 671 len(regions) 672 except: 671 673 msg = 'Polygon_function takes a list of pairs (polygon, value).' 672 674 msg += 'Got %s' %regions 673 raise msg675 raise Exception, msg 674 676 675 677 … … 682 684 raise Exception, msg 683 685 684 686 try: 685 687 a = len(T) 686 688 except: 687 689 msg = 'Polygon_function takes a list of pairs (polygon, value).' 688 690 msg += 'Got %s' %str(T) 689 raise msg691 raise Exception, msg 690 692 691 693 msg = 'Each entry in regions have two components: (polygon, value).' 692 694 msg +='I got %s' %str(T) 693 695 assert a == 2, msg 694 696 695 697 … … 711 713 712 714 def __call__(self, x, y): 713 x = num.array(x, num.Float) 714 y = num.array(y, num.Float) 715 716 N = len(x) 717 assert len(y) == N 718 719 points = num.concatenate((num.reshape(x, (N, 1)), 720 num.reshape(y, (N, 1))), axis=1) 721 722 if callable(self.default): 723 z = self.default(x,y) 724 else: 725 z = num.ones(N, num.Float) * self.default 726 727 for polygon, value in self.regions: 728 indices = inside_polygon(points, polygon) 729 730 # FIXME: This needs to be vectorised 731 if callable(value): 732 for i in indices: 733 xx = num.array([x[i]]) 734 yy = num.array([y[i]]) 715 x = num.array(x, num.float) 716 y = num.array(y, num.float) 717 718 # x and y must be one-dimensional and same length 719 assert len(x.shape) == 1 and len(y.shape) == 1 720 N = x.shape[0] 721 assert y.shape[0] == N 722 723 points = num.ascontiguousarray(num.concatenate((x[:,num.newaxis], 724 y[:,num.newaxis]), 725 axis=1 )) 726 727 if callable(self.default): 728 z = self.default(x,y) 729 else: 730 z = num.ones(N, num.float) * self.default 731 732 for polygon, value in self.regions: 733 indices = inside_polygon(points, polygon) 734 735 # FIXME: This needs to be vectorised 736 if callable(value): 737 for i in indices: 738 xx = num.array([x[i]]) 739 yy = num.array([y[i]]) 735 740 z[i] = value(xx, yy)[0] 736 737 738 741 else: 742 for i in indices: 743 z[i] = value 739 744 740 745 if len(z) == 0: … … 998 1003 interpolation_points: Interpolate polyline data to these positions. 999 1004 List of coordinate pairs [x, y] of 1000 data points or an nx2 Numeric array or a Geospatial_data object1005 data points or an nx2 numeric array or a Geospatial_data object 1001 1006 rtol, atol: Used to determine whether a point is on the polyline or not. See point_on_line. 1002 1007 … … 1008 1013 interpolation_points = interpolation_points.get_data_points(absolute=True) 1009 1014 1010 interpolated_values = num.zeros(len(interpolation_points), num. Float)1011 1012 data = ensure_numeric(data, num. Float)1013 polyline_nodes = ensure_numeric(polyline_nodes, num. Float)1014 interpolation_points = ensure_numeric(interpolation_points, num. Float)1015 gauge_neighbour_id = ensure_numeric(gauge_neighbour_id, num. Int)1015 interpolated_values = num.zeros(len(interpolation_points), num.float) 1016 1017 data = ensure_numeric(data, num.float) 1018 polyline_nodes = ensure_numeric(polyline_nodes, num.float) 1019 interpolation_points = ensure_numeric(interpolation_points, num.float) 1020 gauge_neighbour_id = ensure_numeric(gauge_neighbour_id, num.int) 1016 1021 1017 1022 n = polyline_nodes.shape[0] # Number of nodes in polyline
Note: See TracChangeset
for help on using the changeset viewer.