Ignore:
Timestamp:
Feb 25, 2009, 9:37:22 AM (15 years ago)
Author:
rwilson
Message:

numpy changes.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/numpy/anuga/utilities/polygon.py

    r6360 r6410  
    11#!/usr/bin/env python
    2 """Polygon manipulations
    3 """
     2
     3"""Polygon manipulations"""
    44
    55import numpy as num
     
    1111
    1212
     13##
     14# @brief Determine whether a point is on a line segment.
     15# @param point (x, y) of point in question (tuple, list or array).
     16# @param line ((x1,y1), (x2,y2)) for line (tuple, list or array).
     17# @param rtol Relative error for 'close'.
     18# @param atol Absolute error for 'close'.
     19# @return True or False.
    1320def point_on_line(point, line, rtol=1.0e-5, atol=1.0e-8):
    1421    """Determine whether a point is on a line segment
    1522
    16     Input: 
     23    Input:
    1724        point is given by [x, y]
    1825        line is given by [x0, y0], [x1, y1]] or
     
    2128    Output:
    2229
    23     Note: Line can be degenerate and function still works to discern coinciding points from non-coinciding.
     30    Note: Line can be degenerate and function still works to discern coinciding
     31          points from non-coinciding.
    2432    """
    2533
     
    3139                         line[1,0], line[1,1],
    3240                         rtol, atol)
    33    
     41
    3442    return bool(res)
    3543
     
    4250# result functions for possible states
    4351def lines_dont_coincide(p0,p1,p2,p3):               return (3, None)
    44 def lines_0_fully_included_in_1(p0,p1,p2,p3):       return (2, num.array([p0,p1]))
    45 def lines_1_fully_included_in_0(p0,p1,p2,p3):       return (2, num.array([p2,p3]))
    46 def lines_overlap_same_direction(p0,p1,p2,p3):      return (2, num.array([p0,p3]))
    47 def lines_overlap_same_direction2(p0,p1,p2,p3):     return (2, num.array([p2,p1]))
    48 def lines_overlap_opposite_direction(p0,p1,p2,p3):  return (2, num.array([p0,p2]))
    49 def lines_overlap_opposite_direction2(p0,p1,p2,p3): return (2, num.array([p3,p1]))
     52def lines_0_fully_included_in_1(p0,p1,p2,p3):       return (2,
     53                                                            num.array([p0,p1]))
     54def lines_1_fully_included_in_0(p0,p1,p2,p3):       return (2,
     55                                                            num.array([p2,p3]))
     56def lines_overlap_same_direction(p0,p1,p2,p3):      return (2,
     57                                                            num.array([p0,p3]))
     58def lines_overlap_same_direction2(p0,p1,p2,p3):     return (2,
     59                                                            num.array([p2,p1]))
     60def lines_overlap_opposite_direction(p0,p1,p2,p3):  return (2,
     61                                                            num.array([p0,p2]))
     62def lines_overlap_opposite_direction2(p0,p1,p2,p3): return (2,
     63                                                            num.array([p3,p1]))
    5064
    5165# this function called when an impossible state is found
    52 def lines_error(p1, p2, p3, p4): raise RuntimeError, ("INTERNAL ERROR: p1=%s, p2=%s, p3=%s, p4=%s" %
    53                                                       (str(p1), str(p2), str(p3), str(p4)))
     66def lines_error(p1, p2, p3, p4):
     67    raise RuntimeError, ('INTERNAL ERROR: p1=%s, p2=%s, p3=%s, p4=%s'
     68                         % (str(p1), str(p2), str(p3), str(p4)))
    5469
    5570#                     0s1    0e1    1s0    1e0   # line 0 starts on 1, 0 ends 1, 1 starts 0, 1 ends 0
     
    7287                   }
    7388
     89##
     90# @brief Finds intersection point of two line segments.
     91# @param line0 First line ((x1,y1), (x2,y2)).
     92# @param line1 Second line ((x1,y1), (x2,y2)).
     93# @param rtol Relative error for 'close'.
     94# @param atol Absolute error for 'close'.
     95# @return (status, value) where:
     96#             status = 0 - no intersection, value set to None
     97#                      1 - intersection found, value=(x,y)
     98#                      2 - lines collienar, overlap, value=overlap segment
     99#                      3 - lines collinear, no overlap, value is None
     100#                      4 - lines parallel, value is None
    74101def intersection(line0, line1, rtol=1.0e-5, atol=1.0e-8):
    75     """Returns intersecting point between two line segments or None
    76     (if parallel or no intersection is found).
    77 
    78     However, if parallel lines coincide partly (i.e. shara a common segment,
     102    """Returns intersecting point between two line segments.
     103
     104    However, if parallel lines coincide partly (i.e. share a common segment),
    79105    the line segment where lines coincide is returned
    80    
    81106
    82107    Inputs:
     
    85110                      corresponding to a point.
    86111
    87 
    88112    Output:
    89         status, value
    90 
    91         where status and value is interpreted as follows
    92        
     113        status, value - where status and value is interpreted as follows:
    93114        status == 0: no intersection, value set to None.
    94115        status == 1: intersection point found and returned in value as [x,y].
    95         status == 2: Collinear overlapping lines found. Value takes the form [[x0,y0], [x1,y1]].
     116        status == 2: Collinear overlapping lines found.
     117                     Value takes the form [[x0,y0], [x1,y1]].
    96118        status == 3: Collinear non-overlapping lines. Value set to None.
    97         status == 4: Lines are parallel with a fixed distance apart. Value set to None.
    98    
     119        status == 4: Lines are parallel. Value set to None.
    99120    """
    100121
     
    102123
    103124    line0 = ensure_numeric(line0, num.float)
    104     line1 = ensure_numeric(line1, num.float)   
     125    line1 = ensure_numeric(line1, num.float)
    105126
    106127    x0 = line0[0,0]; y0 = line0[0,1]
     
    113134    u0 = (x3-x2)*(y0-y2) - (y3-y2)*(x0-x2)
    114135    u1 = (x2-x0)*(y1-y0) - (y2-y0)*(x1-x0)
    115        
     136
    116137    if num.allclose(denom, 0.0, rtol=rtol, atol=atol):
    117138        # Lines are parallel - check if they are collinear
     
    123144                           point_on_line([x3, y3], line0, rtol=rtol, atol=atol))
    124145
    125             return collinear_result[state_tuple]([x0,y0],[x1,y1],[x2,y2],[x3,y3])
     146            return collinear_result[state_tuple]([x0,y0], [x1,y1],
     147                                                 [x2,y2], [x3,y3])
    126148        else:
    127149            # Lines are parallel but aren't collinear
    128             return 4, None #FIXME (Ole): Add distance here instead of None 
     150            return 4, None #FIXME (Ole): Add distance here instead of None
    129151    else:
    130152        # Lines are not parallel, check if they intersect
    131153        u0 = u0/denom
    132         u1 = u1/denom       
     154        u1 = u1/denom
    133155
    134156        x = x0 + u0*(x1-x0)
     
    137159        # Sanity check - can be removed to speed up if needed
    138160        assert num.allclose(x, x2 + u1*(x3-x2), rtol=rtol, atol=atol)
    139         assert num.allclose(y, y2 + u1*(y3-y2), rtol=rtol, atol=atol)       
     161        assert num.allclose(y, y2 + u1*(y3-y2), rtol=rtol, atol=atol)
    140162
    141163        # Check if point found lies within given line segments
    142         if 0.0 <= u0 <= 1.0 and 0.0 <= u1 <= 1.0: 
     164        if 0.0 <= u0 <= 1.0 and 0.0 <= u1 <= 1.0:
    143165            # We have intersection
    144166            return 1, num.array([x, y])
     
    147169            return 0, None
    148170
    149 
     171##
     172# @brief Finds intersection point of two line segments.
     173# @param line0 First line ((x1,y1), (x2,y2)).
     174# @param line1 Second line ((x1,y1), (x2,y2)).
     175# @return (status, value) where:
     176#             status = 0 - no intersection, value set to None
     177#                      1 - intersection found, value=(x,y)
     178#                      2 - lines collienar, overlap, value=overlap segment
     179#                      3 - lines collinear, no overlap, value is None
     180#                      4 - lines parallel, value is None
     181# @note Wrapper for C function.
    150182def NEW_C_intersection(line0, line1):
    151     #FIXME(Ole): To write in C
    152     """Returns intersecting point between two line segments or None
    153     (if parallel or no intersection is found).
    154 
    155     However, if parallel lines coincide partly (i.e. shara a common segment,
     183    """Returns intersecting point between two line segments.
     184
     185    However, if parallel lines coincide partly (i.e. share a common segment),
    156186    the line segment where lines coincide is returned
    157    
    158187
    159188    Inputs:
    160189        line0, line1: Each defined by two end points as in: [[x0, y0], [x1, y1]]
    161                       A line can also be a 2x2 numeric array with each row
     190                      A line can also be a 2x2 numpy array with each row
    162191                      corresponding to a point.
    163192
    164 
    165193    Output:
    166         status, value
    167 
    168         where status is interpreted as follows
    169        
    170         status == 0: no intersection with value set to None
    171         status == 1: One intersection point found and returned in value as [x,y]
    172         status == 2: Coinciding line segment found. Value taks the form [[x0,y0], [x1,y1]]
    173         status == 3: Lines would coincide but only if extended. Value set to None
    174         status == 4: Lines are parallel with a fixed distance apart. Value set to None.
    175    
    176     """
    177 
     194        status, value - where status and value is interpreted as follows:
     195        status == 0: no intersection, value set to None.
     196        status == 1: intersection point found and returned in value as [x,y].
     197        status == 2: Collinear overlapping lines found.
     198                     Value takes the form [[x0,y0], [x1,y1]].
     199        status == 3: Collinear non-overlapping lines. Value set to None.
     200        status == 4: Lines are parallel. Value set to None.
     201    """
    178202
    179203    line0 = ensure_numeric(line0, num.float)
    180     line1 = ensure_numeric(line1, num.float)   
     204    line1 = ensure_numeric(line1, num.float)
    181205
    182206    status, value = _intersection(line0[0,0], line0[0,1],
     
    187211    return status, value
    188212
    189 
    190 
    191 
     213##
     214# @brief Determine if one point is inside a polygon.
     215# @param point The point of interest.
     216# @param polygon The polygon to test inclusion in.
     217# @param closed True if points on boundary are considered 'inside' polygon.
     218# @param verbose True if this function is to be verbose.
     219# @return True if point is inside the polygon.
     220# @note Uses inside_polygon() to do the work.
     221# @note Raises Exception if more than one point supplied.
    192222def is_inside_polygon(point, polygon, closed=True, verbose=False):
    193223    """Determine if one point is inside a polygon
     
    195225    See inside_polygon for more details
    196226    """
    197 
    198 #    print 'polygon.py: 198, is_inside_polygon: point=%s' % str(point)
    199 #    print 'polygon.py: 198, is_inside_polygon: polygon=%s' % str(polygon)
    200227
    201228    indices = inside_polygon(point, polygon, closed, verbose)
     
    208235        msg = 'is_inside_polygon must be invoked with one point only'
    209236        raise Exception, msg
    210    
    211 
     237
     238##
     239# @brief Determine which of a set of points are inside a polygon.
     240# @param points A set of points (tuple, list or array).
     241# @param polygon A set of points defining a polygon (tuple, list or array).
     242# @param closed True if points on boundary are considered 'inside' polygon.
     243# @param verbose True if this function is to be verbose.
     244# @return A list of indices of points inside the polygon.
    212245def inside_polygon(points, polygon, closed=True, verbose=False):
    213246    """Determine points inside a polygon
    214247
    215248       Functions inside_polygon and outside_polygon have been defined in
    216        terms af separate_by_polygon which will put all inside indices in
     249       terms of separate_by_polygon which will put all inside indices in
    217250       the first part of the indices array and outside indices in the last
    218251
     
    222255       a list or a numeric array
    223256    """
    224 
    225     #if verbose: print 'Checking input to inside_polygon'
    226257
    227258    try:
     
    232263        # If this fails it is going to be because the points can't be
    233264        # converted to a numeric array.
    234         msg = 'Points could not be converted to numeric array' 
     265        msg = 'Points could not be converted to numeric array'
    235266        raise Exception, msg
    236267
     
    242273        # If this fails it is going to be because the points can't be
    243274        # converted to a numeric array.
    244         msg = 'Polygon %s could not be converted to numeric array' %(str(polygon))
     275        msg = ('Polygon %s could not be converted to numeric array'
     276               % (str(polygon)))
    245277        raise Exception, msg
    246278
     
    253285                                                verbose=verbose)
    254286
    255 #    print 'polygon.py: 255, inside_polygon: indices=%s' % str(indices)
    256 #    print 'polygon.py: 255, inside_polygon: count=%s' % str(count)
    257287    # Return indices of points inside polygon
    258288    return indices[:count]
    259289
    260 
    261 
     290##
     291# @brief Determine if one point is outside a polygon.
     292# @param point The point of interest.
     293# @param polygon The polygon to test inclusion in.
     294# @param closed True if points on boundary are considered 'inside' polygon.
     295# @param verbose True if this function is to be verbose.
     296# @return True if point is outside the polygon.
     297# @note Uses inside_polygon() to do the work.
    262298def is_outside_polygon(point, polygon, closed=True, verbose=False,
    263299                       points_geo_ref=None, polygon_geo_ref=None):
     
    268304
    269305    indices = outside_polygon(point, polygon, closed, verbose)
    270                               #points_geo_ref, polygon_geo_ref)
    271306
    272307    if indices.shape[0] == 1:
     
    277312        msg = 'is_outside_polygon must be invoked with one point only'
    278313        raise Exception, msg
    279    
    280 
     314
     315##
     316# @brief Determine which of a set of points are outside a polygon.
     317# @param points A set of points (tuple, list or array).
     318# @param polygon A set of points defining a polygon (tuple, list or array).
     319# @param closed True if points on boundary are considered 'inside' polygon.
     320# @param verbose True if this function is to be verbose.
     321# @return A list of indices of points outside the polygon.
    281322def outside_polygon(points, polygon, closed = True, verbose = False):
    282323    """Determine points outside a polygon
    283324
    284325       Functions inside_polygon and outside_polygon have been defined in
    285        terms af separate_by_polygon which will put all inside indices in
     326       terms of separate_by_polygon which will put all inside indices in
    286327       the first part of the indices array and outside indices in the last
    287328
     
    289330    """
    290331
    291     #if verbose: print 'Checking input to outside_polygon'
    292332    try:
    293333        points = ensure_numeric(points, num.float)
     
    306346        raise Exception, msg
    307347
    308 
    309348    if len(points.shape) == 1:
    310349        # Only one point was passed in. Convert to array of points
     
    321360    else:
    322361        return indices[count:][::-1]  #return reversed
    323        
    324 
    325 def in_and_outside_polygon(points, polygon, closed = True, verbose = False):
     362
     363##
     364# @brief Separate a list of points into two sets inside+outside a polygon.
     365# @param points A set of points (tuple, list or array).
     366# @param polygon A set of points defining a polygon (tuple, list or array).
     367# @param closed True if points on boundary are considered 'inside' polygon.
     368# @param verbose True if this function is to be verbose.
     369# @return A tuple (in, out) of point indices for poinst inside amd outside.
     370def in_and_outside_polygon(points, polygon, closed=True, verbose=False):
    326371    """Determine points inside and outside a polygon
    327372
    328373       See separate_points_by_polygon for documentation
    329374
    330        Returns an array of points inside and an array of points outside the polygon
    331     """
    332 
    333     #if verbose: print 'Checking input to outside_polygon'
     375       Returns an array of points inside and array of points outside the polygon
     376    """
     377
    334378    try:
    335379        points = ensure_numeric(points, num.float)
     
    352396        points = num.reshape(points, (1,2))
    353397
    354 
    355398    indices, count = separate_points_by_polygon(points, polygon,
    356399                                                closed=closed,
    357400                                                verbose=verbose)
    358    
     401
    359402    # Returns indices of points inside and indices of points outside
    360403    # the polygon
    361 
    362404    if count == len(indices):
    363405        # No points are outside
     
    366408        return  indices[:count], indices[count:][::-1]  #return reversed
    367409
    368 
    369 def separate_points_by_polygon(points, polygon,
    370                                closed = True, verbose = False):
     410##
     411# @brief Sort a list of points into contiguous points inside+outside a polygon.
     412# @param points A set of points (tuple, list or array).
     413# @param polygon A set of points defining a polygon (tuple, list or array).
     414# @param closed True if points on boundary are considered 'inside' polygon.
     415# @param verbose True if this function is to be verbose.
     416# @return (indices, count) where indices are point indices and count is the
     417#          delimiter index between point inside (on left) and others.
     418def separate_points_by_polygon(points, polygon, closed=True, verbose=False):
    371419    """Determine whether points are inside or outside a polygon
    372420
     
    388436       The indices of points outside are obtained as indices[count:]
    389437
    390 
    391438    Examples:
    392439       U = [[0,0], [1,0], [1,1], [0,1]] #Unit square
     
    409456    """
    410457
    411 
    412     #if verbose: print 'Checking input to separate_points_by_polygon'
    413 
    414 
    415458    #Input checks
    416 
    417     assert isinstance(closed, bool), 'Keyword argument "closed" must be boolean'
    418     assert isinstance(verbose, bool), 'Keyword argument "verbose" must be boolean'
    419 
    420 
    421 #    points = ensure_numeric(points, num.float)
     459    assert isinstance(closed, bool), 'Keyword arg "closed" must be boolean'
     460    assert isinstance(verbose, bool), 'Keyword arg "verbose" must be boolean'
     461
    422462    try:
    423463        points = ensure_numeric(points, num.float)
     
    428468        raise Exception, msg
    429469
    430     #if verbose: print 'Checking input to separate_points_by_polygon 2'
    431470    try:
    432471        polygon = ensure_numeric(polygon, num.float)
     
    440479    assert len(polygon.shape) == 2, msg
    441480
    442     msg = 'Polygon array must have two columns' 
     481    msg = 'Polygon array must have two columns'
    443482    assert polygon.shape[1] == 2, msg
    444483
    445 
    446     msg = 'Points array must be 1 or 2 dimensional.'
    447     msg += ' I got %d dimensions' %len(points.shape)
     484    msg = ('Points array must be 1 or 2 dimensional. I got %d dimensions'
     485           % len(points.shape))
    448486    assert 0 < len(points.shape) < 3, msg
    449 
    450487
    451488    if len(points.shape) == 1:
     
    454491        points = num.reshape(points, (1,2))
    455492
    456    
    457     msg = 'Point array must have two columns (x,y), '
    458     msg += 'I got points.shape[1] == %d' % points.shape[1]
     493    msg = ('Point array must have two columns (x,y). I got points.shape[1]=%d'
     494           % points.shape[1])
    459495    assert points.shape[1] == 2, msg
    460496
    461        
    462     msg = 'Points array must be a 2d array. I got %s' %str(points[:30])
     497    msg = 'Points array must be a 2d array. I got %s' % str(points[:30])
    463498    assert len(points.shape) == 2, msg
    464499
     
    466501    assert points.shape[1] == 2, msg
    467502
    468 
    469     N = polygon.shape[0] #Number of vertices in polygon
    470     M = points.shape[0]  #Number of points
    471 
     503    N = polygon.shape[0]    # Number of vertices in polygon
     504    M = points.shape[0]     # Number of points
    472505
    473506    indices = num.zeros( M, num.int )
     
    476509                                        int(closed), int(verbose))
    477510
    478     if verbose: print 'Found %d points (out of %d) inside polygon'\
    479        %(count, M)
     511    if verbose: print 'Found %d points (out of %d) inside polygon' % (count, M)
     512
    480513    return indices, count
    481514
    482 
     515##
     516# @brief Determine area of a polygon.
     517# @param input_polygon The polygon to get area of.
     518# @return A scalar value for the polygon area.
    483519def polygon_area(input_polygon):
    484     """ Determin area of arbitrary polygon
    485     Reference
    486     http://mathworld.wolfram.com/PolygonArea.html
    487     """
    488    
     520    """ Determine area of arbitrary polygon.
     521
     522    Reference:     http://mathworld.wolfram.com/PolygonArea.html
     523    """
     524
    489525    # Move polygon to origin (0,0) to avoid rounding errors
    490526    # This makes a copy of the polygon to avoid destroying it
    491527    input_polygon = ensure_numeric(input_polygon)
    492528    min_x = min(input_polygon[:,0])
    493     min_y = min(input_polygon[:,1])   
     529    min_y = min(input_polygon[:,1])
    494530    polygon = input_polygon - [min_x, min_y]
    495531
    496     # Compute area   
     532    # Compute area
    497533    n = len(polygon)
    498534    poly_area = 0.0
     
    509545        yi = pti[1]
    510546        poly_area += xi*yi1 - xi1*yi
    511        
     547
    512548    return abs(poly_area/2)
    513549
    514 def plot_polygons(polygons_points, style=None,
    515                   figname=None, label=None, verbose=False):
    516    
     550##
     551# @brief Plot a set of polygons.
     552# @param polygons_points List of polygons to plot.
     553# @param style List of styles for each polygon.
     554# @param figname Name to save figure to.
     555# @param label Title for the plot.
     556# @param verbose True if this function is to be verbose.
     557# @return A list of min/max x and y values [minx, maxx, miny, maxy].
     558# @note A style value is 'line' for polygons, 'outside' for points outside.
     559def plot_polygons(polygons_points,
     560                  style=None,
     561                  figname=None,
     562                  label=None,
     563                  verbose=False):
    517564    """ Take list of polygons and plot.
    518565
     
    524571                     - for a polygon, use 'line'
    525572                     - for points falling outside a polygon, use 'outside'
    526                        
     573
    527574    figname          - name to save figure to
    528575
     
    533580    - list of min and max of x and y coordinates
    534581    - plot of polygons
    535     """
    536 
    537     from pylab import ion, hold, plot, axis, figure, legend, savefig, xlabel, ylabel, title, close, title
    538 
    539     assert type(polygons_points) == list,\
    540                'input must be a list of polygons and/or points'
    541                
     582    """
     583
     584    from pylab import ion, hold, plot, axis, figure, legend, savefig, xlabel, \
     585                      ylabel, title, close, title
     586
     587    assert type(polygons_points) == list, \
     588                'input must be a list of polygons and/or points'
     589
    542590    ion()
    543591    hold(True)
     
    548596    maxy = 0.0
    549597
    550     if label is None: label = ''
     598    if label is None:
     599        label = ''
    551600
    552601    n = len(polygons_points)
    553602    colour = []
    554603    if style is None:
    555         style_type = 'line' 
     604        style_type = 'line'
    556605        style = []
    557606        for i in range(n):
     
    560609    else:
    561610        for s in style:
    562             if s == 'line': colour.append('b-')           
     611            if s == 'line': colour.append('b-')
    563612            if s == 'outside': colour.append('r.')
    564613            if s <> 'line':
    565614                if s <> 'outside':
    566615                    colour.append('g.')
    567            
     616
    568617    for i, item in enumerate(polygons_points):
    569618        x, y = poly_xy(item)
     
    594643    close('all')
    595644
    596     vec = [minx,maxx,miny,maxy]
    597 
     645    vec = [minx, maxx, miny, maxy]
    598646    return vec
    599647
     648##
     649# @brief
     650# @param polygon A set of points defining a polygon.
     651# @param verbose True if this function is to be verbose.
     652# @return A tuple (x, y) of X and Y coordinates of the polygon.
     653# @note We duplicate the first point so can have closed polygon in plot.
    600654def poly_xy(polygon, verbose=False):
    601655    """ this is used within plot_polygons so need to duplicate
    602656        the first point so can have closed polygon in plot
    603657    """
    604 
    605     #if verbose: print 'Checking input to poly_xy'
    606658
    607659    try:
     
    610662        raise NameError, e
    611663    except:
    612         msg = 'Polygon %s could not be converted to numeric array' %(str(polygon))
     664        msg = ('Polygon %s could not be converted to numeric array'
     665               % (str(polygon)))
    613666        raise Exception, msg
    614667
     
    617670    x = num.concatenate((x, [polygon[0,0]]), axis = 0)
    618671    y = num.concatenate((y, [polygon[0,1]]), axis = 0)
    619    
     672
    620673    return x, y
    621    
    622 #    x = []
    623 #    y = []
    624 #    n = len(poly)
    625 #    firstpt = poly[0]
    626 #    for i in range(n):
    627 #        thispt = poly[i]
    628 #        x.append(thispt[0])
    629 #        y.append(thispt[1])
    630 
    631 #    x.append(firstpt[0])
    632 #    y.append(firstpt[1])
    633    
    634 #    return x, y
    635 
     674
     675
     676##
     677# @brief Define a class that defines a callable object for a polygon.
     678# @note Object created is function: f: x,y -> z
     679#       where x, y and z are vectors and z depends on whether x,y belongs
     680#       to specified polygons.
    636681class Polygon_function:
    637682    """Create callable object f: x,y -> z, where a,y,z are vectors and
     
    668713    """
    669714
    670     def __init__(self,
    671                  regions,
    672                  default=0.0,
    673                  geo_reference=None):
    674 
     715    ##
     716    # @brief Create instance of a polygon function.
     717    # @param regions A list of (x,y) tuples defining a polygon.
     718    # @param default Value or function returning value for points outside poly.
     719    # @param geo_reference ??
     720    def __init__(self, regions, default=0.0, geo_reference=None):
    675721        try:
    676722            len(regions)
    677723        except:
    678             msg = 'Polygon_function takes a list of pairs (polygon, value).'
    679             msg += 'Got %s' %regions
     724            msg = ('Polygon_function takes a list of pairs (polygon, value).'
     725                   'Got %s' % str(regions))
    680726            raise Exception, msg
    681727
    682 
    683728        T = regions[0]
    684729
    685730        if isinstance(T, basestring):
    686             msg = 'You passed in a list of text values into polygon_function'
    687             msg += ' instead of a list of pairs (polygon, value): "%s"' %T
    688            
     731            msg = ('You passed in a list of text values into polygon_function '
     732                   'instead of a list of pairs (polygon, value): "%s"'
     733                   % str(T))
    689734            raise Exception, msg
    690        
     735
    691736        try:
    692737            a = len(T)
    693738        except:
    694             msg = 'Polygon_function takes a list of pairs (polygon, value).'
    695             msg += 'Got %s' %str(T)
     739            msg = ('Polygon_function takes a list of pairs (polygon, value). '
     740                   'Got %s' % str(T))
    696741            raise Exception, msg
    697742
    698         msg = 'Each entry in regions have two components: (polygon, value).'
    699         msg +='I got %s' %str(T)
     743        msg = ('Each entry in regions have two components: (polygon, value). '
     744               'I got %s' % str(T))
    700745        assert a == 2, msg
    701 
    702746
    703747        if geo_reference is None:
    704748            from anuga.coordinate_transforms.geo_reference import Geo_reference
    705749            geo_reference = Geo_reference()
    706 
    707750
    708751        self.default = default
     
    714757            self.regions.append((P, value))
    715758
    716 
    717 
    718 
     759    ##
     760    # @brief Implement the 'callable' property of Polygon_function.
     761    # @param x List of x coordinates of points ot interest.
     762    # @param y List of y coordinates of points ot interest.
    719763    def __call__(self, x, y):
    720764        x = num.array(x, num.float)
     
    731775
    732776        if callable(self.default):
    733             z = self.default(x,y)
     777            z = self.default(x, y)
    734778        else:
    735779            z = num.ones(N, num.float) * self.default
     
    749793
    750794        if len(z) == 0:
    751             msg = 'Warning: points provided to Polygon function did not fall within'
    752             msg += 'its regions'
    753             msg += 'x in [%.2f, %.2f], y in [%.2f, %.2f]' % (min(x), max(x),
    754                                                              min(y), max(y))
     795            msg = ('Warning: points provided to Polygon function did not fall '
     796                   'within its regions in [%.2f, %.2f], y in [%.2f, %.2f]'
     797                   % (min(x), max(x), min(y), max(y)))
    755798            print msg
    756799
    757            
    758800        return z
    759801
    760 
    761        
    762 # Functions to read and write polygon information       
    763 def read_polygon(filename, split=','):
     802################################################################################
     803# Functions to read and write polygon information
     804################################################################################
     805
     806##
     807# @brief Read polygon data from a file.
     808# @param filename Path to file containing polygon data.
     809# @param delimiter Delimiter to split polygon data with.
     810# @return A list of point data from the polygon file.
     811def read_polygon(filename, delimiter=','):
    764812    """Read points assumed to form a polygon.
    765        There must be exactly two numbers in each line separated by a comma.
    766        No header.
     813
     814    There must be exactly two numbers in each line separated by the delimiter.
     815    No header.
    767816    """
    768817
     
    772821    polygon = []
    773822    for line in lines:
    774         fields = line.split(split)
    775         polygon.append( [float(fields[0]), float(fields[1])] )
     823        fields = line.split(delimiter)
     824        polygon.append([float(fields[0]), float(fields[1])])
    776825
    777826    return polygon
    778827
    779 
     828##
     829# @brief Write polygon data to a file.
     830# @param polygon Polygon points to write to file.
     831# @param filename Path to file to write.
     832# @note Delimiter is assumed to be a comma.
    780833def write_polygon(polygon, filename=None):
    781834    """Write polygon to csv file.
    782        There will be exactly two numbers, easting and northing,
    783        in each line separated by a comma.
    784        
    785        No header.   
     835
     836    There will be exactly two numbers, easting and northing, in each line
     837    separated by a comma.
     838
     839    No header.
    786840    """
    787841
    788842    fid = open(filename, 'w')
    789843    for point in polygon:
    790         fid.write('%f, %f\n' %point)
     844        fid.write('%f, %f\n' % point)
    791845    fid.close()
    792    
    793 
     846
     847##
     848# @brief Unimplemented.
    794849def read_tagged_polygons(filename):
    795850    """
    796851    """
    797852    pass
    798    
     853
     854##
     855# @brief Populate given polygon with uniformly distributed points.
     856# @param polygon Polygon to uniformly fill.
     857# @param number_of_points Number of points required in polygon.
     858# @param seed Seed for random number generator.
     859# @param exclude List of polygons inside main where points should be excluded.
     860# @return List of random points inside input polygon.
     861# @note Delimiter is assumed to be a comma.
    799862def populate_polygon(polygon, number_of_points, seed=None, exclude=None):
    800863    """Populate given polygon with uniformly distributed points.
     
    804867       number_of_points - (optional) number of points
    805868       seed - seed for random number generator (default=None)
    806        exclude - list of polygons (inside main polygon) from where points should be excluded
     869       exclude - list of polygons (inside main polygon) from where points
     870                 should be excluded
    807871
    808872    Output:
     
    831895        if y < min_y: min_y = y
    832896
    833 
    834897    while len(points) < number_of_points:
    835898        x = uniform(min_x, max_x)
     
    838901        append = False
    839902        if is_inside_polygon([x,y], polygon):
    840 
    841903            append = True
    842904
     
    847909                        append = False
    848910
    849 
    850911        if append is True:
    851912            points.append([x,y])
     
    853914    return points
    854915
    855 
     916##
     917# @brief Get a point inside a polygon that is close to an edge.
     918# @param polygon List  of vertices of polygon.
     919# @param delta Maximum distance from an edge is delta * sqrt(2).
     920# @return A point that is inside polgon and close to the polygon edge.
    856921def point_in_polygon(polygon, delta=1e-8):
    857922    """Return a point inside a given polygon which will be close to the
     
    865930       points - a point inside polygon
    866931
    867        searches in all diagonals and up and down (not left and right)
    868     """
     932       searches in all diagonals and up and down (not left and right).
     933    """
     934
    869935    import exceptions
     936
    870937    class Found(exceptions.Exception): pass
    871938
     
    873940    while not point_in:
    874941        try:
    875             for poly_point in polygon: #[1:]:
    876                 for x_mult in range (-1,2):
    877                     for y_mult in range (-1,2):
     942            for poly_point in polygon:     # [1:]:
     943                for x_mult in range(-1, 2):
     944                    for y_mult in range(-1, 2):
    878945                        x = poly_point[0]
    879946                        y = poly_point[1]
     947
    880948                        if x == 0:
    881                             x_delta = x_mult*delta
     949                            x_delta = x_mult * delta
    882950                        else:
    883                             x_delta = x+x_mult*x*delta
     951                            x_delta = x + x_mult*x*delta
    884952
    885953                        if y == 0:
    886                             y_delta = y_mult*delta
     954                            y_delta = y_mult * delta
    887955                        else:
    888                             y_delta = y+y_mult*y*delta
     956                            y_delta = y + y_mult*y*delta
    889957
    890958                        point = [x_delta, y_delta]
    891                         #print "point",point
     959
    892960                        if is_inside_polygon(point, polygon, closed=False):
    893961                            raise Found
     
    895963            point_in = True
    896964        else:
    897             delta = delta*0.1
     965            delta = delta * 0.1
     966
    898967    return point
    899968
    900 
     969##
     970# @brief Calculate approximate number of triangles inside a bounding polygon.
     971# @param interior_regions
     972# @param bounding_poly
     973# @param remainder_res
     974# @return The number of triangles.
    901975def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
    902976    """Calculate the approximate number of triangles inside the
    903977    bounding polygon and the other interior regions
    904978
    905     Polygon areas are converted to square Kms 
     979    Polygon areas are converted to square Kms
    906980
    907981    FIXME: Add tests for this function
    908982    """
    909    
     983
    910984    from anuga.utilities.polygon import polygon_area
    911985
    912 
    913986    # TO DO check if any of the regions fall inside one another
    914987
    915     print '----------------------------------------------------------------------------'
    916     print 'Polygon   Max triangle area (m^2)   Total area (km^2)   Estimated #triangles'
    917     print '----------------------------------------------------------------------------'   
    918        
     988    print '--------------------------------------------------------------------'
     989    print ('Polygon   Max triangle area (m^2)   Total area (km^2)   '
     990           'Estimated #triangles')
     991    print '--------------------------------------------------------------------'
     992
    919993    no_triangles = 0.0
    920994    area = polygon_area(bounding_poly)
    921    
     995
    922996    for poly, resolution in interior_regions:
    923997        this_area = polygon_area(poly)
     
    925999        no_triangles += this_triangles
    9261000        area -= this_area
    927        
     1001
    9281002        print 'Interior ',
    929         print ('%.0f' %resolution).ljust(25),
    930         print ('%.2f' %(this_area/1000000)).ljust(19),
    931         print '%d' %(this_triangles)
    932        
     1003        print ('%.0f' % resolution).ljust(25),
     1004        print ('%.2f' % (this_area/1000000)).ljust(19),
     1005        print '%d' % (this_triangles)
     1006
    9331007    bound_triangles = area/remainder_res
    9341008    no_triangles += bound_triangles
    9351009
    9361010    print 'Bounding ',
    937     print ('%.0f' %remainder_res).ljust(25),
    938     print ('%.2f' %(area/1000000)).ljust(19),
    939     print '%d' %(bound_triangles)   
     1011    print ('%.0f' % remainder_res).ljust(25),
     1012    print ('%.2f' % (area/1000000)).ljust(19),
     1013    print '%d' % (bound_triangles)
    9401014
    9411015    total_number_of_triangles = no_triangles/0.7
    9421016
    9431017    print 'Estimated total number of triangles: %d' %total_number_of_triangles
    944     print 'Note: This is generally about 20% less than the final amount'   
     1018    print 'Note: This is generally about 20% less than the final amount'
    9451019
    9461020    return int(total_number_of_triangles)
    9471021
    948 
     1022##
     1023# @brief Reduce number of points in polygon by the specified factor.
     1024# @param polygon The polygon to reduce.
     1025# @param factor The factor to reduce polygon points by (default 10).
     1026# @return The reduced polygon points list.
     1027# @note The extrema of both axes are preserved.
    9491028def decimate_polygon(polygon, factor=10):
    9501029    """Reduce number of points in polygon by the specified
     
    9631042    max_y = max(num_polygon[:,1])
    9641043    min_x = min(num_polygon[:,0])
    965     min_y = min(num_polygon[:,1])       
     1044    min_y = min(num_polygon[:,1])
    9661045
    9671046    # Keep only some points making sure extrema are kept
    968     reduced_polygon = []   
     1047    reduced_polygon = []
    9691048    for i, point in enumerate(polygon):
    9701049        x = point[0]
    971         y = point[1]       
     1050        y = point[1]
    9721051        if x in [min_x, max_x] and y in [min_y, max_y]:
    9731052            # Keep
     
    9751054        else:
    9761055            if len(reduced_polygon)*factor < i:
    977                 reduced_polygon.append(point)               
     1056                reduced_polygon.append(point)
    9781057
    9791058    return reduced_polygon
    980    
    981    
    982        
    9831059
    9841060##
     
    9971073                         verbose=False):
    9981074    """Interpolate linearly between values data on polyline nodes
    999     of a polyline to list of interpolation points. 
     1075    of a polyline to list of interpolation points.
    10001076
    10011077    data is the data on the polyline nodes.
    1002 
    10031078
    10041079    Inputs:
    10051080      data: Vector or array of data at the polyline nodes.
    1006       polyline_nodes: Location of nodes where data is available.     
     1081      polyline_nodes: Location of nodes where data is available.
    10071082      gauge_neighbour_id: ?
    10081083      interpolation_points: Interpolate polyline data to these positions.
    10091084          List of coordinate pairs [x, y] of
    10101085          data points or an nx2 numeric array or a Geospatial_data object
    1011       rtol, atol: Used to determine whether a point is on the polyline or not. See point_on_line.
     1086      rtol, atol: Used to determine whether a point is on the polyline or not.
     1087                  See point_on_line.
    10121088
    10131089    Output:
    10141090      Interpolated values at interpolation points
    10151091    """
    1016    
     1092
    10171093    if isinstance(interpolation_points, Geospatial_data):
    1018         interpolation_points = interpolation_points.get_data_points(absolute=True)
     1094        interpolation_points = interpolation_points.\
     1095                                    get_data_points(absolute=True)
    10191096
    10201097    interpolated_values = num.zeros(len(interpolation_points), num.float)
     
    10251102    gauge_neighbour_id = ensure_numeric(gauge_neighbour_id, num.int)
    10261103
    1027     n = polyline_nodes.shape[0] # Number of nodes in polyline       
     1104    n = polyline_nodes.shape[0]    # Number of nodes in polyline
     1105
    10281106    # Input sanity check
    10291107    msg = 'interpolation_points are not given (interpolate.py)'
    10301108    assert interpolation_points is not None, msg
     1109
    10311110    msg = 'function value must be specified at every interpolation node'
    1032     assert data.shape[0]==polyline_nodes.shape[0], msg
     1111    assert data.shape[0] == polyline_nodes.shape[0], msg
     1112
    10331113    msg = 'Must define function value at one or more nodes'
    1034     assert data.shape[0]>0, msg
     1114    assert data.shape[0] > 0, msg
    10351115
    10361116    if n == 1:
     
    10411121                              polyline_nodes,
    10421122                              gauge_neighbour_id,
    1043                               interpolation_points,                               
     1123                              interpolation_points,
    10441124                              interpolated_values,
    10451125                              rtol,
    10461126                              atol)
    1047        
     1127
    10481128    return interpolated_values
    10491129
    1050        
     1130##
     1131# @brief
     1132# @param data
     1133# @param polyline_nodes
     1134# @param gauge_neighbour_id
     1135# @param interpolation_points
     1136# @param interpolated_values
     1137# @param rtol
     1138# @param atol
     1139# @return
     1140# @note OBSOLETED BY C-EXTENSION
    10511141def _interpolate_polyline(data,
    1052                           polyline_nodes, 
    1053                           gauge_neighbour_id, 
    1054                           interpolation_points, 
     1142                          polyline_nodes,
     1143                          gauge_neighbour_id,
     1144                          interpolation_points,
    10551145                          interpolated_values,
    10561146                          rtol=1.0e-6,
    10571147                          atol=1.0e-8):
    10581148    """Auxiliary function used by interpolate_polyline
    1059    
     1149
    10601150    NOTE: OBSOLETED BY C-EXTENSION
    10611151    """
    1062    
    1063     number_of_nodes = len(polyline_nodes)               
     1152
     1153    number_of_nodes = len(polyline_nodes)
    10641154    number_of_points = len(interpolation_points)
    1065    
    1066     for j in range(number_of_nodes):               
     1155
     1156    for j in range(number_of_nodes):
    10671157        neighbour_id = gauge_neighbour_id[j]
    1068        
    1069         # FIXME(Ole): I am convinced that gauge_neighbour_id can be discarded, but need to check with John J.
     1158
     1159        # FIXME(Ole): I am convinced that gauge_neighbour_id can be discarded,
     1160        #             but need to check with John J.
    10701161        # Keep it for now (17 Jan 2009)
    1071         # When gone, we can simply interpolate between neighbouring nodes, i.e. neighbour_id = j+1.
    1072         # and the test below becomes something like: if j < number_of_nodes... 
    1073        
     1162        # When gone, we can simply interpolate between neighbouring nodes,
     1163        # i.e. neighbour_id = j+1.
     1164        # and the test below becomes something like: if j < number_of_nodes...
     1165
    10741166        if neighbour_id >= 0:
    10751167            x0, y0 = polyline_nodes[j,:]
    10761168            x1, y1 = polyline_nodes[neighbour_id,:]
    1077            
     1169
    10781170            segment_len = sqrt((x1-x0)**2 + (y1-y0)**2)
    1079             segment_delta = data[neighbour_id] - data[j]           
     1171            segment_delta = data[neighbour_id] - data[j]
    10801172            slope = segment_delta/segment_len
    1081            
    1082                
    1083             for i in range(number_of_points):               
    1084                
     1173
     1174            for i in range(number_of_points):
    10851175                x, y = interpolation_points[i,:]
    1086                 if point_on_line([x, y],
    1087                                  [[x0, y0], [x1, y1]],
    1088                                  rtol=rtol,
    1089                                  atol=atol):
    1090                                  
    1091 
     1176                if point_on_line([x, y], [[x0, y0], [x1, y1]],
     1177                                 rtol=rtol, atol=atol):
    10921178                    dist = sqrt((x-x0)**2 + (y-y0)**2)
    10931179                    interpolated_values[i] = slope*dist + data[j]
    1094      
    1095 
    1096    
    1097 
    1098 ##############################################
    1099 #Initialise module
     1180
     1181
     1182################################################################################
     1183# Initialise module
     1184################################################################################
    11001185
    11011186from anuga.utilities import compile
     
    11041189    from polygon_ext import _point_on_line
    11051190    from polygon_ext import _separate_points_by_polygon
    1106     from polygon_ext import _interpolate_polyline   
     1191    from polygon_ext import _interpolate_polyline
    11071192    #from polygon_ext import _intersection
    11081193
Note: See TracChangeset for help on using the changeset viewer.