Changeset 2655


Ignore:
Timestamp:
Apr 4, 2006, 3:49:42 PM (18 years ago)
Author:
duncan
Message:

Adding ability to remove points outside of the mesh when interpolating

Location:
inundation
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • inundation/fit_interpolate/interpolate.py

    r2577 r2655  
    2323
    2424from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \
    25      ArrayType, allclose, take, NewAxis
     25     ArrayType, allclose, take, NewAxis, arange
    2626
    2727from pyvolution.mesh import Mesh
     
    3131from pyvolution.quad import build_quadtree
    3232from utilities.numerical_tools import ensure_numeric, mean, INF
    33 from utilities.polygon import inside_polygon
     33from utilities.polygon import in_and_outside_polygon
    3434from geospatial_data.geospatial_data import Geospatial_data
    3535
     
    128128        point_coordinates = ensure_numeric(point_coordinates, Float)
    129129       
    130         #Remove points falling outside mesh boundary
    131         # do this bit later - that sorta means this becomes an object
    132         # get a list of what indices are outside the boundary
    133         # 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)
    135135       
    136136        #Build n x m interpolation matrix
     
    143143        A = Sparse(n,m)
    144144
    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:
    147147            #For each data_coordinate point
    148148            if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n)
     
    164164                    A[i,j] = sigmas[j]
    165165            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)
    172168        return A
    173169
     
    232228                self._A_can_be_reused = False
    233229                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
    235232                for end in range(start_blocking_len
    236233                                 ,len(point_coordinates)
     
    251248        doesn't occur.
    252249
     250        Return the point data, z.
     251       
    253252        See interpolate for doc info.
    254253        """
     
    259258
    260259    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
    263271
    264272
  • inundation/fit_interpolate/test_interpolate.py

    r2651 r2655  
    2020from coordinate_transforms.geo_reference import Geo_reference
    2121from shallow_water import Domain, Transmissive_boundary
    22 from utilities.numerical_tools import mean
     22from utilities.numerical_tools import mean, INF
    2323from data_manager import get_dataobject
    2424from geospatial_data.geospatial_data import Geospatial_data
     
    237237        assert allclose(sum(results, axis=1), 1.0)
    238238
    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):
    241240        """Try arbitrary datapoints one outside the triangle.
    242241        That one should be ignored
     
    253252        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44], [5.0, 7.0]]
    254253
    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)
    262255        results = interp._build_interpolation_matrix_A(data).todense()
    263256        assert allclose(sum(results, axis=1), [1,1,1,0])
     
    301294        assert allclose(A, answer)
    302295
    303 
     296        #FIXME -  tests are hopefully failing due to points outsdie the
     297        # mesh
    304298    def test_interpolate_attributes_to_points(self):
    305299        v0 = [0.0, 0.0]
     
    442436
    443437        z = interp.interpolate(f, point_coords)
    444         answer = [ [0., 0., 0., 0.]] # (-1,-1)
     438        answer = array([ [INF, INF, INF, INF]]) # (-1,-1)
    445439
    446440        #print "***********"
     
    452446        # of the mesh? Not currently. 
    453447
    454         assert allclose(z, answer)
    455 
     448        for i in range(4):
     449            self.failUnless( z[0,i] == answer[0,i], 'Fail!')
     450       
    456451    def test_interpolate_attributes_to_pointsIV(self):
    457452        a = [-1.0, 0.0]
     
    594589
    595590    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       
    600596
    601597
     
    667663
    668664    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       
    672667        #Three timesteps
    673668        time = [1.0, 5.0, 6.0]
     
    728723
    729724    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   
    735729
    736730        #Three timesteps
     
    793787
    794788
    795     def qtest_interpolation_interface(self):
     789    def test_points_outside_the_polygon(self):
    796790        a = [-1.0, 0.0]
    797791        b = [3.0, 4.0]
     
    815809                        [-0.9, -1.5],
    816810                        [0.5, -1.9],
    817                         [999999, 9999999]]
     811                        [999999, 9999999]] # point Outside poly
    818812        geo_data = Geospatial_data(data_points = point_coords)
    819813
     
    827821                  2*linear_function(point_coords) ]
    828822        answer = transpose(answer)
     823        answer[11,:] = [INF, INF]
    829824        #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!')
    833875       
    834876#-------------------------------------------------------------
    835877if __name__ == "__main__":
     878
     879    #suite = unittest.makeSuite(Test_Interpolate,'test')
    836880    suite = unittest.makeSuite(Test_Interpolate,'test')
     881    #suite = unittest.makeSuite(Test_Interpolate,'test_points_outside_the_polygon')
    837882    runner = unittest.TextTestRunner(verbosity=1)
    838883    runner.run(suite)
  • inundation/utilities/polygon.py

    r2633 r2655  
    127127                                                verbose=verbose)
    128128
     129   
    129130    if one_point:
    130131        return count != 1
     
    135136        else:
    136137            return indices[count:][::-1]  #return reversed
     138       
     139
     140def 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
    137184
    138185
  • inundation/utilities/test_polygon.py

    r2442 r2655  
    477477        assert inside_polygon(point, polygon)
    478478
     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       
    479551#-------------------------------------------------------------
    480552if __name__ == "__main__":
Note: See TracChangeset for help on using the changeset viewer.