Changeset 656


Ignore:
Timestamp:
Dec 2, 2004, 12:47:16 PM (20 years ago)
Author:
ole
Message:

bug fix + Polygon_function

Location:
inundation/ga/storm_surge/pyvolution
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • inundation/ga/storm_surge/pyvolution/test_util.py

    r655 r656  
    332332        assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
    333333       
     334       
     335       
     336    def test_polygon_function(self):
     337        p1 = [[0,0], [10,0], [10,10], [0,10]]
     338        p2 = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
     339       
     340        f = Polygon_function( [(p1, 1.0)] )     
     341        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #Two first inside p1
     342        assert allclose(z, [1,1,0,0])
     343               
     344       
     345        f = Polygon_function( [(p2, 2.0)] )
     346        z = f([5, 5, 27, 35], [5, 9, 8, -5]) #First and last inside p2
     347        assert allclose(z, [2,0,0,2])   
     348               
     349       
     350        #Combined       
     351        f = Polygon_function( [(p1, 1.0), (p2, 2.0)] )
     352        z = f([5, 5, 27, 35], [5, 9, 8, -5])
     353        assert allclose(z, [2,1,0,2])   
     354       
     355       
     356       
    334357    def test_point_on_line(self):
    335358               
     
    355378        polygon = [[0,0], [1,0], [1,1], [0,1]]
    356379 
    357         x = y = 0.5
    358         assert inside_polygon( (x, y), polygon )               
    359        
    360         y = 1.5
    361         assert not inside_polygon( (x, y), polygon )                           
     380        assert inside_polygon( (0.5, 0.5), polygon )                 
     381        assert not inside_polygon( (0.5, 1.5), polygon )   
     382        assert not inside_polygon( (0.5, -0.5), polygon )       
     383        assert not inside_polygon( (-0.5, 0.5), polygon )               
     384        assert not inside_polygon( (1.5, 0.5), polygon )       
     385       
    362386        #Try point on borders   
    363387        assert inside_polygon( (1., 0.5), polygon, closed=True)
     
    381405        assert not inside_polygon( (0.5, -0.5), polygon )               
    382406       
    383        
    384        
    385407        #Now try the vector formulation returning indices
    386408        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
     
    389411        assert allclose( res, [0,1,2] )
    390412       
     413        #Very convoluted polygon
     414        polygon = [[0,0], [10,10], [15,5], [20, 10], [25,0], [30,10], [40,-10]]
     415        assert inside_polygon( (5, 5), polygon )       
     416        assert inside_polygon( (17, 7), polygon )               
     417        assert inside_polygon( (27, 2), polygon )                       
     418        assert inside_polygon( (35, -5), polygon )                     
     419        assert not inside_polygon( (15, 7), polygon )
     420        assert not inside_polygon( (24, 3), polygon )
     421        assert not inside_polygon( (25, -10), polygon )
     422       
     423       
     424       
     425        #Another combination (that failed)
     426        polygon = [[0,0], [10,0], [10,10], [0,10]]
     427        assert inside_polygon( (5, 5), polygon )       
     428        assert inside_polygon( (7, 7), polygon )       
     429        assert not inside_polygon( (-17, 7), polygon )                 
     430        assert not inside_polygon( (7, 17), polygon )                   
     431        assert not inside_polygon( (17, 7), polygon )                   
     432        assert not inside_polygon( (27, 8), polygon )                   
     433        assert not inside_polygon( (35, -5), polygon )                 
     434       
     435       
     436                       
     437       
     438        #Multiple polygons     
     439       
     440        polygon = [[0,0], [1,0], [1,1], [0,1], [0,0],
     441                   [10,10], [11,10], [11,11], [10,11], [10,10]]
     442        assert inside_polygon( (0.5, 0.5), polygon )               
     443        assert inside_polygon( (10.5, 10.5), polygon )                 
     444       
     445        #FIXME: Fails
     446        #assert not inside_polygon( (5.5, 5.5), polygon )                               
     447
     448        #Polygon with a hole   
     449        polygon = [[-1,-1], [2,-1], [2,2], [-1,2], [-1,-1],
     450                   [0,0], [1,0], [1,1], [0,1], [0,0]]
     451                       
     452        assert inside_polygon( (0, -0.5), polygon )
     453       
     454        #FIXME: Fails   
     455        #assert not inside_polygon( (0.5, 0.5), polygon )       
     456       
     457        points = [ [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]                   
    391458       
    392459
  • inundation/ga/storm_surge/pyvolution/util.py

    r655 r656  
    154154    px = polygon[:,0]
    155155    py = polygon[:,1]   
    156                
     156
     157                       
    157158    #Begin algorithm
    158159    indices = []
     
    161162        x = points[k, 0]
    162163        y = points[k, 1]       
    163 
    164            
     164       
    165165        inside = False
    166166        for i in range(N):
     
    174174                    inside = False
    175175                break
    176        
    177             #Check truly inside polygon             
    178             if px[i] < y and py[j] >= y or\
    179                py[j] < y and py[i] >= y:
    180                 if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
    181                     inside = not inside
     176            else:       
     177                #Check if truly inside polygon             
     178                if py[i] < y and py[j] >= y or\
     179                   py[j] < y and py[i] >= y:
     180                    if px[i] + (y-py[i])/(py[j]-py[i])*(px[j]-px[i]) < x:
     181                        inside = not inside
    182182                   
    183183        if inside: indices.append(k)
     184
     185        #print x, y, polygon, inside
     186
    184187   
    185188    if one_point:
     
    189192             
    190193
     194class Polygon_function:
     195    """Create callable object f: x,y -> z, where a,y,z are vectors and
     196    where f will return different values depending on whether x,y belongs
     197    to specified polygons.
     198   
     199    To instantiate:
     200     
     201       Polygon_function(polygons)
     202       
     203    where polygons is a dictionary of the form
     204   
     205      {P0: v0, P1: v1, ...}
     206     
     207      with Pi being lists of vertices defining polygons and vi either
     208      constants or functions of x,y to be applied to points with the polygon.
     209     
     210    The function takes an optional argument, default which is the value
     211    (or function) to used for points not belonging to any polygon.
     212    For example:
     213   
     214       Polygon_function(polygons, default = 0.03)       
     215     
     216    If omitted the default value will be 0.0 
     217   
     218    Note: If two polygons overlap, the one last in the list takes precedence
     219
     220    Note: default can only be a constant.   
     221    """
     222   
     223    def __init__(self, regions, default = 0.0):
     224       
     225        try:
     226            len(regions)
     227        except:
     228            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons           
     229            raise msg
     230
     231
     232        T = regions[0]                                   
     233        try:
     234            a = len(T)
     235        except:   
     236            msg = 'Polygon_function takes a list of pairs (polygon, value). Got %s' %polygons           
     237            raise msg   
     238                   
     239        assert a == 2, 'Must have two component each: %s' %T       
     240
     241        self.regions = regions             
     242        self.default = default
     243         
     244         
     245    def __call__(self, x, y):
     246        from util import inside_polygon
     247        from Numeric import ones, Float, concatenate, array, reshape, choose
     248       
     249        x = array(x).astype(Float)
     250        y = array(y).astype(Float)     
     251        x = reshape(x, (len(x), 1))
     252        y = reshape(y, (len(y), 1))     
     253               
     254        points = concatenate( (x,y), axis=1 )
     255       
     256        z = ones(x.shape[0], Float) * self.default
     257        for polygon, value in self.regions:
     258            indices = inside_polygon(points, polygon)
     259           
     260            for i in indices:
     261                z[i] = value
     262
     263        return z                 
     264             
    191265class File_function:
    192266    """Read time series from file and return a callable object:
Note: See TracChangeset for help on using the changeset viewer.