Changeset 5897 for anuga_core/source/anuga/utilities/polygon.py
- Timestamp:
- Nov 6, 2008, 12:17:15 PM (16 years ago)
- File:
-
- 1 edited
-
anuga_core/source/anuga/utilities/polygon.py (modified) (25 diffs)
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/utilities/polygon.py
r5891 r5897 10 10 # #print 'Could not find scipy - using Numeric' 11 11 12 ##from numpy import float, int, zeros, ones, array, concatenate, reshape, dot, allclose, newaxis, ascontiguousarray 13 import numpy 12 from Numeric import Float, Int, zeros, ones, array, concatenate, reshape, dot, allclose 14 13 15 14 … … 77 76 # FIXME (Ole): Write this in C 78 77 79 line0 = ensure_numeric(line0, numpy.float)80 line1 = ensure_numeric(line1, numpy.float)78 line0 = ensure_numeric(line0, Float) 79 line1 = ensure_numeric(line1, Float) 81 80 82 81 x0 = line0[0,0]; y0 = line0[0,1] … … 90 89 u1 = (x2-x0)*(y1-y0) - (y2-y0)*(x1-x0) 91 90 92 if numpy.allclose(denom, 0.0):91 if allclose(denom, 0.0): 93 92 # Lines are parallel - check if they coincide on a shared a segment 94 93 95 if numpy.allclose( [u0, u1], 0.0 ):94 if allclose( [u0, u1], 0.0 ): 96 95 # We now know that the lines if continued coincide 97 96 # The remaining check will establish if the finite lines share a segment … … 121 120 if line0_starts_on_line1 and line0_ends_on_line1: 122 121 # Shared segment is line0 fully included in line1 123 segment = numpy.array([[x0, y0], [x1, y1]])122 segment = array([[x0, y0], [x1, y1]]) 124 123 125 124 if line1_starts_on_line0 and line1_ends_on_line0: 126 125 # Shared segment is line1 fully included in line0 127 segment = numpy.array([[x2, y2], [x3, y3]])126 segment = array([[x2, y2], [x3, y3]]) 128 127 129 128 … … 131 130 if line0_starts_on_line1 and line1_ends_on_line0: 132 131 # Shared segment from line0 start to line 1 end 133 segment = numpy.array([[x0, y0], [x3, y3]])132 segment = array([[x0, y0], [x3, y3]]) 134 133 135 134 if line1_starts_on_line0 and line0_ends_on_line1: 136 135 # Shared segment from line1 start to line 0 end 137 segment = numpy.array([[x2, y2], [x1, y1]])136 segment = array([[x2, y2], [x1, y1]]) 138 137 139 138 … … 141 140 if line0_starts_on_line1 and line1_starts_on_line0: 142 141 # Shared segment from line0 start to line 1 end 143 segment = numpy.array([[x0, y0], [x2, y2]])142 segment = array([[x0, y0], [x2, y2]]) 144 143 145 144 if line0_ends_on_line1 and line1_ends_on_line0: 146 145 # Shared segment from line0 start to line 1 end 147 segment = numpy.array([[x3, y3], [x1, y1]])146 segment = array([[x3, y3], [x1, y1]]) 148 147 149 148 … … 162 161 163 162 # Sanity check - can be removed to speed up if needed 164 assert numpy.allclose(x, x2 + u1*(x3-x2))165 assert numpy.allclose(y, y2 + u1*(y3-y2))163 assert allclose(x, x2 + u1*(x3-x2)) 164 assert allclose(y, y2 + u1*(y3-y2)) 166 165 167 166 # Check if point found lies within given line segments … … 169 168 # We have intersection 170 169 171 return 1, numpy.array([x, y])170 return 1, array([x, y]) 172 171 else: 173 172 # No intersection … … 204 203 205 204 206 line0 = ensure_numeric(line0, numpy.float)207 line1 = ensure_numeric(line1, numpy.float)205 line0 = ensure_numeric(line0, Float) 206 line1 = ensure_numeric(line1, Float) 208 207 209 208 status, value = _intersection(line0[0,0], line0[0,1], … … 257 256 # converted to a numeric array. 258 257 msg = 'Points could not be converted to Numeric array' 259 raise TypeError,msg258 raise msg 260 259 261 260 try: … … 266 265 # If this fails it is going to be because the points can't be 267 266 # converted to a numeric array. 268 msg = 'Polygon %s could not be converted to Numeric array' % 269 raise TypeError,msg267 msg = 'Polygon %s could not be converted to Numeric array' %(str(polygon)) 268 raise msg 270 269 271 270 if len(points.shape) == 1: 272 271 # Only one point was passed in. Convert to array of points 273 points = numpy.reshape(points, (1,2))272 points = reshape(points, (1,2)) 274 273 275 274 indices, count = separate_points_by_polygon(points, polygon, … … 313 312 #if verbose: print 'Checking input to outside_polygon' 314 313 try: 315 points = ensure_numeric(points, numpy.float)314 points = ensure_numeric(points, Float) 316 315 except NameError, e: 317 316 raise NameError, e … … 321 320 322 321 try: 323 polygon = ensure_numeric(polygon, numpy.float)322 polygon = ensure_numeric(polygon, Float) 324 323 except NameError, e: 325 324 raise NameError, e … … 331 330 if len(points.shape) == 1: 332 331 # Only one point was passed in. Convert to array of points 333 points = numpy.reshape(points, (1,2))332 points = reshape(points, (1,2)) 334 333 335 334 indices, count = separate_points_by_polygon(points, polygon, … … 340 339 if count == len(indices): 341 340 # No points are outside 342 return numpy.array([])341 return array([]) 343 342 else: 344 343 return indices[count:][::-1] #return reversed … … 355 354 #if verbose: print 'Checking input to outside_polygon' 356 355 try: 357 points = ensure_numeric(points, numpy.float)356 points = ensure_numeric(points, Float) 358 357 except NameError, e: 359 358 raise NameError, e … … 363 362 364 363 try: 365 polygon = ensure_numeric(polygon, numpy.float)364 polygon = ensure_numeric(polygon, Float) 366 365 except NameError, e: 367 366 raise NameError, e … … 372 371 if len(points.shape) == 1: 373 372 # Only one point was passed in. Convert to array of points 374 points = numpy.reshape(points, (1,2))373 points = reshape(points, (1,2)) 375 374 376 375 … … 440 439 assert isinstance(verbose, bool), 'Keyword argument "verbose" must be boolean' 441 440 442 ## print 'Before: points=%s, flags=%s' % (type(points), str(points.flags)) 441 443 442 try: 444 ## points = numpy.ascontiguousarray(ensure_numeric(points, numpy.float)) 445 points = ensure_numeric(points, numpy.float) 443 points = ensure_numeric(points, Float) 446 444 except NameError, e: 447 445 raise NameError, e 448 446 except: 449 447 msg = 'Points could not be converted to Numeric array' 450 raise TypeError, msg 451 ## print 'After: points=%s, flags=%s' % (type(points), str(points.flags)) 448 raise msg 452 449 453 450 #if verbose: print 'Checking input to separate_points_by_polygon 2' 454 451 try: 455 ## polygon = numpy.ascontiguousarray(ensure_numeric(polygon, numpy.float)) 456 polygon = ensure_numeric(polygon, numpy.float) 452 polygon = ensure_numeric(polygon, Float) 457 453 except NameError, e: 458 454 raise NameError, e 459 455 except: 460 456 msg = 'Polygon could not be converted to Numeric array' 461 raise TypeError,msg457 raise msg 462 458 463 459 msg = 'Polygon array must be a 2d array of vertices' … … 476 472 # Only one point was passed in. 477 473 # Convert to array of points 478 points = numpy.reshape(points, (1,2))474 points = reshape(points, (1,2)) 479 475 480 476 481 477 msg = 'Point array must have two columns (x,y), ' 482 msg += 'I got points.shape[1] == %d' % points.shape[1]478 msg += 'I got points.shape[1] == %d' %points.shape[0] 483 479 assert points.shape[1] == 2, msg 484 480 … … 495 491 496 492 497 indices = numpy.zeros( M, numpy.int )493 indices = zeros( M, Int ) 498 494 499 495 count = _separate_points_by_polygon(points, polygon, indices, … … 622 618 623 619 try: 624 polygon = ensure_numeric(polygon, numpy.float)620 polygon = ensure_numeric(polygon, Float) 625 621 except NameError, e: 626 622 raise NameError, e … … 631 627 x = polygon[:,0] 632 628 y = polygon[:,1] 633 x = numpy.concatenate((x, [polygon[0,0]]), axis = 0)634 y = numpy.concatenate((y, [polygon[0,1]]), axis = 0)629 x = concatenate((x, [polygon[0,0]]), axis = 0) 630 y = concatenate((y, [polygon[0,1]]), axis = 0) 635 631 636 632 return x, y … … 684 680 """ 685 681 686 def __init__(self, regions, default=0.0, geo_reference=None , verbose=False):682 def __init__(self, regions, default=0.0, geo_reference=None): 687 683 688 684 try: … … 727 723 self.regions.append( (P, value) ) 728 724 729 self.verbose = verbose730 731 725 732 726 733 727 734 728 def __call__(self, x, y): 735 x = numpy.array(x).astype(numpy.float) 736 y = numpy.array(y).astype(numpy.float) 737 738 assert len(x.shape) == 1 and len(y.shape) == 1 739 740 N = x.shape[0] 741 assert y.shape[0] == N 742 743 points = numpy.ascontiguousarray(numpy.concatenate( (x[:,numpy.newaxis], y[:,numpy.newaxis]), axis=1 )) 744 729 x = array(x).astype(Float) 730 y = array(y).astype(Float) 731 732 N = len(x) 733 assert len(y) == N 734 735 points = concatenate( (reshape(x, (N, 1)), 736 reshape(y, (N, 1))), axis=1 ) 737 745 738 if callable(self.default): 746 z = self.default(x, 739 z = self.default(x,y) 747 740 else: 748 z = numpy.ones(N, numpy.float) * self.default741 z = ones(N, Float) * self.default 749 742 750 743 for polygon, value in self.regions: 751 indices = inside_polygon(points, polygon , verbose=self.verbose)744 indices = inside_polygon(points, polygon) 752 745 753 746 # FIXME: This needs to be vectorised 754 747 if callable(value): 755 748 for i in indices: 756 xx = numpy.array([x[i]])757 yy = numpy.array([y[i]])749 xx = array([x[i]]) 750 yy = array([y[i]]) 758 751 z[i] = value(xx, yy)[0] 759 752 else:
Note: See TracChangeset
for help on using the changeset viewer.
