Changeset 2655
- Timestamp:
- Apr 4, 2006, 3:49:42 PM (18 years ago)
- Location:
- inundation
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/fit_interpolate/interpolate.py
r2577 r2655 23 23 24 24 from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \ 25 ArrayType, allclose, take, NewAxis 25 ArrayType, allclose, take, NewAxis, arange 26 26 27 27 from pyvolution.mesh import Mesh … … 31 31 from pyvolution.quad import build_quadtree 32 32 from utilities.numerical_tools import ensure_numeric, mean, INF 33 from utilities.polygon import in side_polygon33 from utilities.polygon import in_and_outside_polygon 34 34 from geospatial_data.geospatial_data import Geospatial_data 35 35 … … 128 128 point_coordinates = ensure_numeric(point_coordinates, Float) 129 129 130 #Remove points falling outside mesh boundary131 # do this bit later - that sorta means this becomes an object132 # get a list of what indices are outside the boundary133 # maybe fill these rows with n/a?134 130 if verbose: print 'Getting indices inside mesh boundary' 131 self.inside_poly_indices, self.outside_poly_indices = \ 132 in_and_outside_polygon(point_coordinates, 133 self.mesh.get_boundary_polygon(), 134 closed = True, verbose = verbose) 135 135 136 136 #Build n x m interpolation matrix … … 143 143 A = Sparse(n,m) 144 144 145 #Compute matrix elements 146 for i in range(n):145 #Compute matrix elements for points inside the mesh 146 for i in self.inside_poly_indices: 147 147 #For each data_coordinate point 148 148 if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n) … … 164 164 A[i,j] = sigmas[j] 165 165 else: 166 # Hack so this message only appears when the test is 167 #run in fit_interpoate directory. 168 #Error message reminds me to do something 169 #with points outside the polygon 170 if os.getcwd()[-15:] == "fit_interpolate": 171 print 'Could not find triangle for point', x 166 msg = 'Could not find triangle for point', x 167 raise Exception(msg) 172 168 return A 173 169 … … 232 228 self._A_can_be_reused = False 233 229 start=0 234 z = self.interpolate_block(f, point_coordinates[0:0]) 230 # creating a dummy array to concatenate to. 231 z = zeros((0,2)) # 2 for points in 2D space 235 232 for end in range(start_blocking_len 236 233 ,len(point_coordinates) … … 251 248 doesn't occur. 252 249 250 Return the point data, z. 251 253 252 See interpolate for doc info. 254 253 """ … … 259 258 260 259 def _get_point_data_z(self, f, verbose=False): 261 return self._A * f 262 260 """ 261 Return the point data, z. 262 263 Precondition, 264 The _A matrix has been created 265 """ 266 z = self._A * f 267 for i in self.outside_poly_indices: # Taking into account points outside the 268 # mesh. 269 z[i,:] = INF 270 return z 263 271 264 272 -
inundation/fit_interpolate/test_interpolate.py
r2651 r2655 20 20 from coordinate_transforms.geo_reference import Geo_reference 21 21 from shallow_water import Domain, Transmissive_boundary 22 from utilities.numerical_tools import mean 22 from utilities.numerical_tools import mean, INF 23 23 from data_manager import get_dataobject 24 24 from geospatial_data.geospatial_data import Geospatial_data … … 237 237 assert allclose(sum(results, axis=1), 1.0) 238 238 239 #FIXME - have to change this test to check default info 240 def NO_test_arbitrary_datapoints_some_outside(self): 239 def test_arbitrary_datapoints_some_outside(self): 241 240 """Try arbitrary datapoints one outside the triangle. 242 241 That one should be ignored … … 253 252 data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44], [5.0, 7.0]] 254 253 255 256 interp = Interpolate(points, vertices, data, precrop = True) 257 258 results = interp._build_interpolation_matrix_A(data).todense() 259 assert allclose(sum(results, axis=1), 1.0) 260 261 interp = Interpolate(points, vertices, data, precrop = False) 254 interp = Interpolate(points, vertices, data) 262 255 results = interp._build_interpolation_matrix_A(data).todense() 263 256 assert allclose(sum(results, axis=1), [1,1,1,0]) … … 301 294 assert allclose(A, answer) 302 295 303 296 #FIXME - tests are hopefully failing due to points outsdie the 297 # mesh 304 298 def test_interpolate_attributes_to_points(self): 305 299 v0 = [0.0, 0.0] … … 442 436 443 437 z = interp.interpolate(f, point_coords) 444 answer = [ [0., 0., 0., 0.]]# (-1,-1)438 answer = array([ [INF, INF, INF, INF]]) # (-1,-1) 445 439 446 440 #print "***********" … … 452 446 # of the mesh? Not currently. 453 447 454 assert allclose(z, answer) 455 448 for i in range(4): 449 self.failUnless( z[0,i] == answer[0,i], 'Fail!') 450 456 451 def test_interpolate_attributes_to_pointsIV(self): 457 452 a = [-1.0, 0.0] … … 594 589 595 590 def test_interpolation_interface_time_only(self): 596 """Test spatio-temporal interpolation 597 Test that spatio temporal function performs the correct 598 interpolations in both time and space 599 """ 591 592 # Test spatio-temporal interpolation 593 # Test that spatio temporal function performs the correct 594 # interpolations in both time and space 595 600 596 601 597 … … 667 663 668 664 def test_interpolation_interface_spatial_only(self): 669 """Test spatio-temporal interpolation with constant time 670 """ 671 665 # Test spatio-temporal interpolation with constant time 666 672 667 #Three timesteps 673 668 time = [1.0, 5.0, 6.0] … … 728 723 729 724 def test_interpolation_interface(self): 730 """Test spatio-temporal interpolation 731 Test that spatio temporal function performs the correct 732 interpolations in both time and space 733 """ 734 725 # Test spatio-temporal interpolation 726 # Test that spatio temporal function performs the correct 727 # interpolations in both time and space 728 735 729 736 730 #Three timesteps … … 793 787 794 788 795 def qtest_interpolation_interface(self):789 def test_points_outside_the_polygon(self): 796 790 a = [-1.0, 0.0] 797 791 b = [3.0, 4.0] … … 815 809 [-0.9, -1.5], 816 810 [0.5, -1.9], 817 [999999, 9999999]] 811 [999999, 9999999]] # point Outside poly 818 812 geo_data = Geospatial_data(data_points = point_coords) 819 813 … … 827 821 2*linear_function(point_coords) ] 828 822 answer = transpose(answer) 823 answer[11,:] = [INF, INF] 829 824 #print "z",z 830 #print "answer",answer 831 assert allclose(z, answer) 832 825 #print "answer _ fixed",answer 826 assert allclose(z[0:10], answer[0:10]) 827 self.failUnless( z[11,1] == answer[11,1], 'Fail!') 828 self.failUnless( z[11,0] == answer[11,0], 'Fail!') 829 830 831 def test_points_outside_the_polygon(self): 832 a = [-1.0, 0.0] 833 b = [3.0, 4.0] 834 c = [4.0, 1.0] 835 d = [-3.0, 2.0] #3 836 e = [-1.0, -2.0] 837 f = [1.0, -2.0] #5 838 839 vertices = [a, b, c, d,e,f] 840 triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc 841 842 point_coords = [[-2.0, 2.0], 843 [-1.0, 1.0], 844 [9999.0, 9999.0], # point Outside poly 845 [-9999.0, 1.0], # point Outside poly 846 [2.0, 1.0], 847 [0.0, 0.0], 848 [1.0, 0.0], 849 [0.0, -1.0], 850 [-0.2, -0.5], 851 [-0.9, -1.5], 852 [0.5, -1.9], 853 [999999, 9999999]] # point Outside poly 854 geo_data = Geospatial_data(data_points = point_coords) 855 856 interp = Interpolate(vertices, triangles) 857 f = array([linear_function(vertices),2*linear_function(vertices) ]) 858 f = transpose(f) 859 #print "f",f 860 z = interp.interpolate(f, geo_data) 861 #z = interp.interpolate(f, point_coords) 862 answer = [linear_function(point_coords), 863 2*linear_function(point_coords) ] 864 answer = transpose(answer) 865 answer[2,:] = [INF, INF] 866 answer[3,:] = [INF, INF] 867 answer[11,:] = [INF, INF] 868 #print "z",z 869 #print "answer _ fixed",answer 870 assert allclose(z[0:1], answer[0:1]) 871 assert allclose(z[4:10], answer[4:10]) 872 for i in [2,3,11]: 873 self.failUnless( z[i,1] == answer[11,1], 'Fail!') 874 self.failUnless( z[i,0] == answer[11,0], 'Fail!') 833 875 834 876 #------------------------------------------------------------- 835 877 if __name__ == "__main__": 878 879 #suite = unittest.makeSuite(Test_Interpolate,'test') 836 880 suite = unittest.makeSuite(Test_Interpolate,'test') 881 #suite = unittest.makeSuite(Test_Interpolate,'test_points_outside_the_polygon') 837 882 runner = unittest.TextTestRunner(verbosity=1) 838 883 runner.run(suite) -
inundation/utilities/polygon.py
r2633 r2655 127 127 verbose=verbose) 128 128 129 129 130 if one_point: 130 131 return count != 1 … … 135 136 else: 136 137 return indices[count:][::-1] #return reversed 138 139 140 def in_and_outside_polygon(points, polygon, closed = True, verbose = False): 141 """Determine points inside and outside a polygon 142 143 See separate_points_by_polygon for documentation 144 145 Returns an array of points inside and an array of points outside the polygon 146 """ 147 148 if verbose: print 'Checking input to outside_polygon' 149 try: 150 points = ensure_numeric(points, Float) 151 except NameError, e: 152 raise NameError, e 153 except: 154 msg = 'Points could not be converted to Numeric array' 155 raise msg 156 157 try: 158 polygon = ensure_numeric(polygon, Float) 159 except NameError, e: 160 raise NameError, e 161 except: 162 msg = 'Polygon could not be converted to Numeric array' 163 raise msg 164 165 166 167 if len(points.shape) == 1: 168 one_point = True 169 points = reshape(points, (1,2)) 170 else: 171 one_point = False 172 173 indices, count = separate_points_by_polygon(points, polygon, 174 closed=closed, 175 verbose=verbose) 176 # Returns an array of points inside and an array of points outside 177 # the polygon 178 179 if count == len(indices): 180 # No points are outside 181 return indices[:count],[] 182 else: 183 return indices[:count], indices[count:][::-1] #return reversed 137 184 138 185 -
inundation/utilities/test_polygon.py
r2442 r2655 477 477 assert inside_polygon(point, polygon) 478 478 479 def test_in_and_outside_polygon_main(self): 480 481 482 #Simplest case: Polygon is the unit square 483 polygon = [[0,0], [1,0], [1,1], [0,1]] 484 485 inside, outside = in_and_outside_polygon( (0.5, 0.5), polygon ) 486 assert inside[0] == 0 487 assert len(outside) == 0 488 489 inside, outside = in_and_outside_polygon( (1., 0.5), polygon, closed=True) 490 assert inside[0] == 0 491 assert len(outside) == 0 492 493 inside, outside = in_and_outside_polygon( (1., 0.5), polygon, closed=False) 494 assert len(inside) == 0 495 assert outside[0] == 0 496 497 points = [(1., 0.25),(1., 0.75) ] 498 inside, outside = in_and_outside_polygon( points, polygon, closed=True) 499 assert (inside, [0,1]) 500 assert len(outside) == 0 501 502 inside, outside = in_and_outside_polygon( points, polygon, closed=False) 503 assert len(inside) == 0 504 assert (outside, [0,1]) 505 506 507 points = [(100., 0.25),(0.5, 0.5) ] 508 inside, outside = in_and_outside_polygon( points, polygon) 509 assert (inside, [1]) 510 assert outside[0] == 0 511 512 points = [(100., 0.25),(0.5, 0.5), (39,20), (0.6,0.7),(56,43),(67,90) ] 513 inside, outside = in_and_outside_polygon( points, polygon) 514 assert (inside, [1,3]) 515 assert (outside, [0,2,4,5]) 516 517 def zzztest_inside_polygon_main(self): 518 print "inside",inside 519 print "outside",outside 520 521 assert not inside_polygon( (0.5, 1.5), polygon ) 522 assert not inside_polygon( (0.5, -0.5), polygon ) 523 assert not inside_polygon( (-0.5, 0.5), polygon ) 524 assert not inside_polygon( (1.5, 0.5), polygon ) 525 526 #Try point on borders 527 assert inside_polygon( (1., 0.5), polygon, closed=True) 528 assert inside_polygon( (0.5, 1), polygon, closed=True) 529 assert inside_polygon( (0., 0.5), polygon, closed=True) 530 assert inside_polygon( (0.5, 0.), polygon, closed=True) 531 532 assert not inside_polygon( (0.5, 1), polygon, closed=False) 533 assert not inside_polygon( (0., 0.5), polygon, closed=False) 534 assert not inside_polygon( (0.5, 0.), polygon, closed=False) 535 assert not inside_polygon( (1., 0.5), polygon, closed=False) 536 537 538 539 #From real example (that failed) 540 polygon = [[20,20], [40,20], [40,40], [20,40]] 541 points = [ [40, 50] ] 542 res = inside_polygon(points, polygon) 543 assert len(res) == 0 544 545 polygon = [[20,20], [40,20], [40,40], [20,40]] 546 points = [ [25, 25], [30, 20], [40, 50], [90, 20], [40, 90] ] 547 res = inside_polygon(points, polygon) 548 assert len(res) == 2 549 assert allclose(res, [0,1]) 550 479 551 #------------------------------------------------------------- 480 552 if __name__ == "__main__":
Note: See TracChangeset
for help on using the changeset viewer.