Changeset 7858


Ignore:
Timestamp:
Jun 18, 2010, 4:43:10 PM (14 years ago)
Author:
hudson
Message:

Refactorings to increase code quality, fixed missing log import from sww_interrogate.

Location:
trunk/anuga_core/source/anuga
Files:
1 added
14 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/source/anuga/__init__.py

    r7851 r7858  
    4242
    4343from anuga.geometry.polygon import read_polygon, plot_polygons, polygon_area
    44 from anuga.geometry.polygon import Polygon_function
     44from anuga.geometry.polygon_function import Polygon_function
    4545
    4646from anuga.abstract_2d_finite_volumes.pmesh2domain import \
     
    5252from anuga.utilities.file_utils import copy_code_files
    5353
    54 from anuga.geometry.polygon import read_polygon, Polygon_function
     54from anuga.geometry.polygon import read_polygon
    5555from anuga.caching import cache
    5656
  • trunk/anuga_core/source/anuga/culvert_flows/test_culvert_class.py

    r7778 r7858  
    77
    88from anuga.utilities.system_tools import get_pathname_from_package
    9 from anuga.geometry.polygon import Polygon_function
     9from anuga.geometry.polygon_function import Polygon_function
    1010       
    1111from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
  • trunk/anuga_core/source/anuga/culvert_flows/test_new_culvert_class.py

    r7778 r7858  
    99
    1010from anuga.utilities.system_tools import get_pathname_from_package
    11 from anuga.geometry.polygon import Polygon_function
     11from anuga.geometry.polygon_function import Polygon_function
    1212       
    1313from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
  • trunk/anuga_core/source/anuga/file/csv_file.py

    r7854 r7858  
    1111"""
    1212
     13
    1314import csv
    1415import numpy as num
     16import anuga.utilities.log as log
    1517
    1618
     
    168170    if completed:
    169171        try:
    170             file = str(kwargs['file_name'])
     172            file_name = str(kwargs['file_name'])
    171173        except:
    172             raise 'kwargs must have file_name'
     174            raise Exception('kwargs must have file_name')
    173175    else:
    174176        # write temp file in output directory
    175177        try:
    176             file = str(kwargs['output_dir']) + 'detail_temp.csv'
     178            file_name = str(kwargs['output_dir']) + 'detail_temp.csv'
    177179        except:
    178             raise 'kwargs must have output_dir'
     180            raise Exception('kwargs must have output_dir')
    179181
    180182    # extracts the header info and the new line info
     
    199201    # try to open!
    200202    try:
    201         fid = open(file, 'r')
     203        fid = open(file_name, 'r')
    202204        file_header = fid.readline()
    203205        fid.close()
    204206        if verbose: log.critical('read file header %s' % file_header)
    205     except:
    206         msg = 'try to create new file: %s' % file
    207         if verbose: log.critical(msg)
     207    except Exception:
     208        msg = 'try to create new file: %s' % file_name
     209        if verbose:
     210            log.critical(msg)
    208211        #tries to open file, maybe directory is bad
    209212        try:
    210             fid = open(file, 'w')
     213            fid = open(file_name, 'w')
    211214            fid.write(header)
    212215            fid.close()
     
    218221    # if header is same or this is a new file
    219222    if file_header == str(header):
    220         fid = open(file, 'a')
     223        fid = open(file_name, 'a')
    221224        fid.write(line)
    222225        fid.close()
     
    225228        # if header is different and has completed will append info to
    226229        # end of details_temp.cvs file in output directory
    227         file = str(kwargs['output_dir']) + 'detail_temp.csv'
    228         fid = open(file, 'a')
     230        file_name = str(kwargs['output_dir']) + 'detail_temp.csv'
     231        fid = open(file_name, 'a')
    229232        fid.write(header)
    230233        fid.write(line)
     
    244247
    245248def load_csv_as_building_polygons(file_name,
    246                           floor_height=3,
    247                           clipping_polygons=None):
     249                          floor_height=3):
    248250    """
    249251    Convert CSV files of the form:
     
    285287           
    286288
    287 ##
    288 # @brief Convert CSV file into a dictionary of polygons and associated values.
    289 # @param filename The path to the file to read, value_name name for the 4th column
    290289def load_csv_as_polygons(file_name,
    291290                 value_name='value',
     
    338337   
    339338    msg = 'Did not find expected column header: northing'   
    340     assert 'northing' in X.keys(), northing
     339    assert 'northing' in X.keys(), msg
    341340   
    342341    msg = 'Did not find expected column header: northing'       
     
    357356    past_ids = {}
    358357    last_id = None
    359     for i, id in enumerate(X['id']):
     358    for i, poly_id in enumerate(X['id']):
    360359
    361360        # Check for duplicate polygons
    362         if id in past_ids:
     361        if poly_id in past_ids:
    363362            msg = 'Polygon %s was duplicated in line %d' % (id, i)
    364363            raise Exception, msg
    365364       
    366         if id not in polygons:
     365        if poly_id not in polygons:
    367366            # Start new polygon
    368             polygons[id] = []
     367            polygons[poly_id] = []
    369368            if values is not None:
    370                 values[id] = X[value_name][i]
     369                values[poly_id] = X[value_name][i]
    371370
    372371            # Keep track of previous polygon ids
     
    385384               
    386385            if exclude is True:
    387                 excluded_polygons[id]=True
    388 
    389         polygons[id].append(point)   
     386                excluded_polygons[poly_id]=True
     387
     388        polygons[poly_id].append(point)   
    390389           
    391390        # Check that value is the same across each polygon
    392391        msg = 'Values must be the same across each polygon.'
    393         msg += 'I got %s in line %d but it should have been %s' % (X[value_name][i], i, values[id])
    394         assert values[id] == X[value_name][i], msg
    395 
    396         last_id = id
     392        msg += 'I got %s in line %d but it should have been %s' % \
     393                            (X[value_name][i], i, values[poly_id])
     394        assert values[poly_id] == X[value_name][i], msg
     395
     396        last_id = poly_id
    397397
    398398    # Weed out polygons that were not wholly inside clipping polygons
    399     for id in excluded_polygons:
    400         del polygons[id]
     399    for poly_id in excluded_polygons:
     400        del polygons[poly_id]
    401401       
    402402    return polygons, values
  • trunk/anuga_core/source/anuga/file/sts.py

    r7778 r7858  
    9494        outfile.variables[q + Write_sts.RANGE][1] = -max_float # Max
    9595
    96         # Doing sts_precision instead of Float gives cast errors.
    97         outfile.createVariable('time', netcdf_float, ('number_of_timesteps',))
    98 
    99         for q in Write_sts.sts_quantities:
    100             outfile.createVariable(q, sts_precision, ('number_of_timesteps',
    101                                                       'number_of_points'))
    102             outfile.createVariable(q + Write_sts.RANGE, sts_precision,
    103                                    ('numbers_in_range',))
    104             # Initialise ranges with small and large sentinels.
    105             # If this was in pure Python we could have used None sensibly
    106             outfile.variables[q + Write_sts.RANGE][0] = max_float  # Min
    107             outfile.variables[q + Write_sts.RANGE][1] = -max_float # Max
    108 
    109         if isinstance(times, (list, num.ndarray)):
    110             outfile.variables['time'][:] = times    #Store time relative
    111 
    112         if verbose:
    113             log.critical('------------------------------------------------')
    114             log.critical('Statistics:')
    115             log.critical('    t in [%f, %f], len(t) == %d'
    116                          % (num.min(times), num.max(times), len(times.flat)))
     96        self.write_dynamic_quantities(outfile, Write_sts.sts_quantities, times)
    11797
    11898    ##
     
    265245
    266246
     247    def write_dynamic_quantities(self, outfile, quantities,
     248                    times, precis = netcdf_float32, verbose = False):   
     249        """
     250            Write out given quantities to file.
     251        """
     252        for q in quantities:
     253            outfile.createVariable(q, precis, ('number_of_timesteps',
     254                                                      'number_of_points'))
     255            outfile.createVariable(q + Write_sts.RANGE, precis,
     256                                   ('numbers_in_range',))
     257
     258            # Initialise ranges with small and large sentinels.
     259            # If this was in pure Python we could have used None sensibly
     260            outfile.variables[q+Write_sts.RANGE][0] = max_float  # Min
     261            outfile.variables[q+Write_sts.RANGE][1] = -max_float # Max
     262
     263        # Doing sts_precision instead of Float gives cast errors.
     264        outfile.createVariable('time', netcdf_float, ('number_of_timesteps',))
     265
     266        if isinstance(times, (list, num.ndarray)):
     267            outfile.variables['time'][:] = times    # Store time relative
     268
     269        if verbose:
     270            log.critical('------------------------------------------------')
     271            log.critical('Statistics:')
     272            log.critical('    t in [%f, %f], len(t) == %d'
     273                         % (num.min(times), num.max(times), len(times.flat)))
     274
    267275
    268276
  • trunk/anuga_core/source/anuga/file/sww.py

    r7841 r7858  
    99import anuga.utilities.log as log
    1010from Scientific.IO.NetCDF import NetCDFFile
     11
     12from sts import Write_sts
    1113
    1214from anuga.coordinate_transforms.geo_reference import \
     
    455457
    456458
    457 class Write_sww:
     459class Write_sww(Write_sts):
    458460    """
    459461        A class to write an SWW file.
     
    462464        manually.
    463465    """
    464     RANGE = '_range'
    465     EXTREMA = ':extrema'
    466466
    467467    def __init__(self, static_quantities, dynamic_quantities):
     
    567567                                                       'number_of_vertices'))
    568568
    569         # Doing sww_precision instead of Float gives cast errors.
    570         outfile.createVariable('time', netcdf_float,
    571                                ('number_of_timesteps',))
    572 
    573                                
     569
    574570        for q in self.static_quantities:
    575571           
     
    585581            outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
    586582
    587         #if 'elevation' in self.static_quantities:   
    588         #    # FIXME: Backwards compat - get rid of z once old view has retired
    589         #    outfile.createVariable('z', sww_precision,
    590         #                           ('number_of_points',))
    591                                
    592         for q in self.dynamic_quantities:
    593             outfile.createVariable(q, sww_precision, ('number_of_timesteps',
    594                                                       'number_of_points'))
    595             outfile.createVariable(q + Write_sww.RANGE, sww_precision,
    596                                    ('numbers_in_range',))
    597 
    598             # Initialise ranges with small and large sentinels.
    599             # If this was in pure Python we could have used None sensibly
    600             outfile.variables[q+Write_sww.RANGE][0] = max_float  # Min
    601             outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
    602 
    603         if isinstance(times, (list, num.ndarray)):
    604             outfile.variables['time'][:] = times    # Store time relative
    605 
    606         if verbose:
    607             log.critical('------------------------------------------------')
    608             log.critical('Statistics:')
    609             log.critical('    t in [%f, %f], len(t) == %d'
    610                          % (num.min(times), num.max(times), len(times.flat)))
     583       
     584        self.write_dynamic_quantities(outfile, self.dynamic_quantities, times)
     585
    611586
    612587
  • trunk/anuga_core/source/anuga/file_conversion/csv2sts.py

    r7853 r7858  
    33"""
    44    Module to convert a .csv file to an .sts file.
     5   
     6    You can use the function from your own scripts, or call this script on the
     7    command line.
    58   
    69    Use this module to convert an arbitrary csv file to an sts file.
     
    4245        data:
    4346
    44          stage = 4, 150.66667, 150.83334, 151, 151.16667, -34, -34.16667, -34.33333,
    45             -34.5, -1, -5, -9, -13 ;
     47         stage = 4, 150.66667, 150.83334, 151, 151.16667, -34, -34.16667,
     48                -34.33333, -34.5, -1, -5, -9, -13 ;
    4649
    4750         time = 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 ;
    4851        }
    4952   
     53    As of June 2010 this module has a pylint quality rating of 9.30/10.
    5054"""
    5155
    52 import csv
    53 import os
    5456import sys
    5557import getopt
    5658from anuga.utilities import log
    57 import numpy as num
    5859from Scientific.IO.NetCDF import NetCDFFile
    5960from anuga.file.csv_file import load_csv_as_dict
     
    108109
    109110def usage():
     111    """ Display usage of this module from the comand line. """
    110112    print 'csv2sts - convert a csv file to an sts file.'
    111113    print 'Usage: csv2sts [-hv] [--help] [--verbose]',
     
    161163   
    162164if __name__ == "__main__":
    163     """ Entry point if run from command line """
    164165    main(sys.argv[1:])       
  • trunk/anuga_core/source/anuga/file_conversion/urs2sww.py

    r7800 r7858  
    2424from anuga.file.sww import Write_sww 
    2525
    26 
    27 ################################################################################
    28 # CONVERTING UNGRIDDED URS DATA TO AN SWW FILE
    29 ################################################################################
    30 
    31 ##
    32 # @brief Convert URS file(s) (wave prop) to an SWW file.
    33 # @param basename_in Stem of the input filenames.
    34 # @param basename_out Path to the output SWW file.
    35 # @param verbose True if this function is to be verbose.
    36 # @param mint
    37 # @param maxt
    38 # @param mean_stage
    39 # @param origin Tuple with geo-ref UTM coordinates (zone, easting, northing).
    40 # @param hole_points_UTM
    41 # @param zscale
    42 # @note Also convert latitude and longitude to UTM. All coordinates are
    43 #       assumed to be given in the GDA94 datum.
    44 # @note Input filename stem has suffixes '-z-mux', '-e-mux' and '-n-mux'
    45 #       added for relative height, x-velocity and y-velocity respectively.
     26###############################################################
     27
    4628def urs_ungridded2sww(basename_in='o', basename_out=None, verbose=False,
    4729                      mint=None, maxt=None,
  • trunk/anuga_core/source/anuga/geometry/aabb.py

    r7751 r7858  
    33    Contains a class describing a bounding box. It contains methods
    44    to split itself and return child boxes.
     5   
     6    As of June 2010 this module has a pylint quality rating of 10/10.
    57"""
    68
     
    1517    """
    1618   
    17     def __init__(self, xmin, xmax, ymin, ymax):
     19    def __init__(self, xmin, xmax=None, ymin=None, ymax=None):
    1820        """ Define axially-algned bounding box.
    1921            xmin is minimum x
     
    2224            ymax is maximum y (absolute coord, ie, not size)
    2325        """
    24         self.xmin = xmin   
    25         self.xmax = xmax
    26         self.ymin = ymin   
    27         self.ymax = ymax
     26        if not xmax:
     27            # try treating first arg as a list of points
     28            try:
     29                xmin[0][0]
     30            except:
     31                raise Exception('Single parameter to AABB must be point list.')
     32               
     33            self.xmin, self.ymin = self.xmax, self.ymax = xmin[0]
     34            self.include(xmin[1:])
     35        else:
     36            self.xmin = xmin   
     37            self.xmax = xmax
     38            self.ymin = ymin   
     39            self.ymax = ymax
    2840
    2941
     
    96108                and (self.ymin <= point[1] <= self.ymax)
    97109       
     110    def include(self, point_list):
     111        """ Include points in AABB.
     112            Bounding box will be expanded to include passed points
     113       
     114            point_list is a list of points.
     115        """
     116        for point in point_list:
     117            pt_x, pt_y = point
     118            if pt_x < self.xmin:
     119                self.xmin = pt_x
     120            if pt_x > self.xmax:
     121                self.xmax = pt_x
     122            if pt_y < self.ymin:
     123                self.ymin = pt_y               
     124            if pt_y > self.ymax:
     125                self.ymax = pt_y               
  • trunk/anuga_core/source/anuga/geometry/polygon.py

    r7841 r7858  
    1010import anuga.utilities.log as log
    1111
     12from aabb import AABB
    1213
    1314##
     
    3637
    3738    res = _point_on_line(point[0], point[1],
    38                          line[0,0], line[0,1],
    39                          line[1,0], line[1,1],
     39                         line[0, 0], line[0, 1],
     40                         line[1, 0], line[1, 1],
    4041                         rtol, atol)
    4142
     
    4950
    5051# result functions for possible states
    51 def lines_dont_coincide(p0,p1,p2,p3):               return (3, None)
    52 def lines_0_fully_included_in_1(p0,p1,p2,p3):       return (2,
    53                                                             num.array([p0,p1]))
    54 def lines_1_fully_included_in_0(p0,p1,p2,p3):       return (2,
    55                                                             num.array([p2,p3]))
    56 def lines_overlap_same_direction(p0,p1,p2,p3):      return (2,
    57                                                             num.array([p0,p3]))
    58 def lines_overlap_same_direction2(p0,p1,p2,p3):     return (2,
    59                                                             num.array([p2,p1]))
    60 def lines_overlap_opposite_direction(p0,p1,p2,p3):  return (2,
    61                                                             num.array([p0,p2]))
    62 def lines_overlap_opposite_direction2(p0,p1,p2,p3): return (2,
    63                                                             num.array([p3,p1]))
     52def lines_dont_coincide(p0, p1, p2, p3):
     53    return (3, None)
     54   
     55def lines_0_fully_included_in_1(p0, p1, p2, p3):
     56    return (2, num.array([p0, p1]))
     57   
     58def lines_1_fully_included_in_0(p0, p1, p2, p3):
     59    return (2, num.array([p2, p3]))
     60   
     61def lines_overlap_same_direction(p0, p1, p2, p3):
     62    return (2, num.array([p0, p3]))
     63   
     64def lines_overlap_same_direction2(p0, p1, p2, p3):
     65    return (2, num.array([p2, p1]))
     66   
     67def lines_overlap_opposite_direction(p0, p1, p2, p3):
     68    return (2, num.array([p0, p2]))
     69   
     70def lines_overlap_opposite_direction2(p0, p1, p2, p3):
     71    return (2, num.array([p3, p1]))
    6472
    6573# this function called when an impossible state is found
     
    146154                           point_on_line([x3, y3], line0, rtol=rtol, atol=atol))
    147155
    148             return collinear_result[state_tuple]([x0,y0], [x1,y1],
    149                                                  [x2,y2], [x3,y3])
     156            return collinear_result[state_tuple]([x0, y0], [x1, y1],
     157                                                 [x2, y2], [x3, y3])
    150158        else:
    151159            # Lines are parallel but aren't collinear
     
    269277    return bool(_is_inside_triangle(point, triangle, int(closed), rtol, atol))
    270278   
    271    
    272 
    273     # FIXME (Ole): The rest of this function has been made
    274     # obsolete by the C extension.
    275    
    276     # Quickly reject points that are clearly outside
    277     if point[0] < min(triangle[:, 0]):
    278         return False
    279     if point[0] > max(triangle[:, 0]):
    280         return False   
    281     if point[1] < min(triangle[:, 1]):
    282         return False
    283     if point[1] > max(triangle[:, 1]):
    284         return False       
    285 
    286 
    287     # Start search   
    288     A = triangle[0, :]
    289     B = triangle[1, :]
    290     C = triangle[2, :]
    291    
    292     # Now check if point lies wholly inside triangle
    293     v0 = C-A
    294     v1 = B-A       
    295        
    296     a00 = num.inner(v0, v0)
    297     a10 = a01 = num.inner(v0, v1)
    298     a11 = num.inner(v1, v1)
    299    
    300     denom = a11*a00 - a01*a10
    301    
    302     if abs(denom) > 0.0:
    303         v = point-A
    304         b0 = num.inner(v0, v)       
    305         b1 = num.inner(v1, v)           
    306    
    307         alpha = (b0*a11 - b1*a01)/denom
    308         beta = (b1*a00 - b0*a10)/denom       
    309 
    310         if (alpha > 0.0) and (beta > 0.0) and (alpha+beta < 1.0):
    311             return True
    312 
    313 
    314     if closed is True:
    315         # Check if point lies on one of the edges
    316        
    317         for X, Y in [[A, B], [B, C], [C, A]]:
    318             res = _point_on_line(point[0], point[1],
    319                                  X[0], X[1],
    320                                  Y[0], Y[1],
    321                                  rtol, atol)
    322            
    323             if res:
    324                 return True
    325                
    326     return False
    327 
    328    
    329279def is_complex(polygon, verbose=False):
    330280    """Check if a polygon is complex (self-intersecting).
     
    339289           
    340290    def key_xpos(item):
     291        """ Return the x coord out of the passed point for sorting key. """
    341292        return (item[0][0])
    342293   
     
    377328                        print 'Self-intersecting polygon found, type ', type
    378329                        print 'point', point,
    379                         print 'vertices: ', leftmost, ' - ', l_x[cmp]               
     330                        print 'vertices: ', leftmost, ' - ', l_x[cmp] 
    380331                    return True           
    381332            cmp += 1
     
    422373    try:
    423374        points = ensure_absolute(points)
    424     except NameError, e:
    425         raise NameError, e
     375    except NameError, err:
     376        raise NameError, err
    426377    except:
    427378        # If this fails it is going to be because the points can't be
     
    513464    if len(points.shape) == 1:
    514465        # Only one point was passed in. Convert to array of points
    515         points = num.reshape(points, (1,2))
     466        points = num.reshape(points, (1, 2))
    516467
    517468    indices, count = separate_points_by_polygon(points, polygon,
     
    559510    if len(points.shape) == 1:
    560511        # Only one point was passed in. Convert to array of points
    561         points = num.reshape(points, (1,2))
     512        points = num.reshape(points, (1, 2))
    562513
    563514    indices, count = separate_points_by_polygon(points, polygon,
     
    632583        except:
    633584            msg = 'Points could not be converted to numeric array'
    634             raise msg
     585            raise Exception(msg)
    635586
    636587        try:
    637588            polygon = ensure_numeric(polygon, num.float)
    638589        except NameError, e:
    639             raise NameError, e
     590            raise NameError(e)
    640591        except:
    641592            msg = 'Polygon could not be converted to numeric array'
    642             raise msg
     593            raise Exception(msg)
    643594
    644595        msg = 'Polygon array must be a 2d array of vertices'
     
    646597
    647598        msg = 'Polygon array must have two columns'
    648         assert polygon.shape[1]==2, msg
     599        assert polygon.shape[1] == 2, msg
    649600
    650601        msg = ('Points array must be 1 or 2 dimensional. '
     
    654605        if len(points.shape) == 1:
    655606            # Only one point was passed in.  Convert to array of points.
    656             points = num.reshape(points, (1,2))
     607            points = num.reshape(points, (1, 2))
    657608   
    658609            msg = ('Point array must have two columns (x,y), '
     
    663614            msg = ('Points array must be a 2d array. I got %s.'
    664615                   % str(points[:30]))
    665             assert len(points.shape)==2, msg
     616            assert len(points.shape) == 2, msg
    666617
    667618            msg = 'Points array must have two columns'
    668             assert points.shape[1]==2, msg
     619            assert points.shape[1] == 2, msg
    669620
    670621    N = polygon.shape[0] # Number of vertices in polygon
     
    681632    return indices, count
    682633
    683 ##
    684 # @brief Determine area of a polygon.
    685 # @param input_polygon The polygon to get area of.
    686 # @return A scalar value for the polygon area.
     634
    687635def polygon_area(input_polygon):
    688636    """ Determine area of arbitrary polygon.
    689637
    690     Reference:     http://mathworld.wolfram.com/PolygonArea.html
    691     """
    692 
     638        input_polygon The polygon to get area of.
     639       
     640        return A scalar value for the polygon area.
     641
     642        Reference:     http://mathworld.wolfram.com/PolygonArea.html
     643    """
    693644    # Move polygon to origin (0,0) to avoid rounding errors
    694645    # This makes a copy of the polygon to avoid destroying it
    695646    input_polygon = ensure_numeric(input_polygon)
    696     min_x = min(input_polygon[:,0])
    697     min_y = min(input_polygon[:,1])
     647    min_x = min(input_polygon[:, 0])
     648    min_y = min(input_polygon[:, 1])
    698649    polygon = input_polygon - [min_x, min_y]
    699650
     
    742693    Outputs:
    743694
    744     - list of min and max of x and y coordinates
    745695    - plot of polygons
    746696    """
     
    754704    ion()
    755705    hold(True)
    756 
    757     minx = 1e10
    758     maxx = 0.0
    759     miny = 1e10
    760     maxy = 0.0
    761706
    762707    if label is None:
     
    772717            alpha = max(0.0, min(1.0, alpha))
    773718
    774     n = len(polygons_points)
     719    num_points = len(polygons_points)
    775720    colour = []
    776721    if style is None:
    777722        style_type = 'line'
    778723        style = []
    779         for i in range(n):
     724        for i in range(num_points):
    780725            style.append(style_type)
    781726            colour.append('b-')
    782727    else:
    783         for s in style:
    784             if s == 'line': colour.append('b-')
    785             if s == 'outside': colour.append('r.')
    786             if s == 'point': colour.append('g.')
    787             if s != 'line':
    788                 if s != 'outside':
    789                     if s != 'point':
    790                         colour.append(s)
     728        for style_name in style:
     729            if style_name == 'line':
     730                colour.append('b-')
     731            if style_name == 'outside':
     732                colour.append('r.')
     733            if style_name == 'point':
     734                colour.append('g.')
     735            if style_name not in ['line', 'outside', 'point']:
     736                colour.append(style_name)
    791737
    792738    for i, item in enumerate(polygons_points):
    793         x, y = poly_xy(item)
    794         if min(x) < minx:
    795             minx = min(x)
    796         if max(x) > maxx:
    797             maxx = max(x)
    798         if min(y) < miny:
    799             miny = min(y)
    800         if max(y) > maxy:
    801             maxy = max(y)
    802         plot(x, y, colour[i])
     739        pt_x, pt_y = _poly_xy(item)
     740        plot(pt_x, pt_y, colour[i])
    803741        if alpha:
    804             fill(x, y, colour[i], alpha=alpha)
     742            fill(pt_x, pt_y, colour[i], alpha=alpha)
    805743        xlabel('x')
    806744        ylabel('y')
     
    814752    close('all')
    815753
    816     vec = [minx, maxx, miny, maxy]
    817     return vec
    818 
    819 
    820 def poly_xy(polygon):
     754
     755def _poly_xy(polygon):
    821756    """ this is used within plot_polygons so need to duplicate
    822757        the first point so can have closed polygon in plot
     
    829764    try:
    830765        polygon = ensure_numeric(polygon, num.float)
    831     except NameError, e:
    832         raise NameError, e
     766    except NameError, err:
     767        raise NameError, err
    833768    except:
    834769        msg = ('Polygon %s could not be converted to numeric array'
     
    836771        raise Exception, msg
    837772
    838     x = polygon[:, 0]
    839     y = polygon[:, 1]
    840     x = num.concatenate((x, [polygon[0, 0]]), axis = 0)
    841     y = num.concatenate((y, [polygon[0, 1]]), axis = 0)
    842 
    843     return x, y
    844 
    845 
    846 class Polygon_function:
    847     """Create callable object f: x,y -> z, where a,y,z are vectors and
    848     where f will return different values depending on whether x,y belongs
    849     to specified polygons.
    850 
    851     To instantiate:
    852 
    853        Polygon_function(polygons)
    854 
    855     where polygons is a list of tuples of the form
    856 
    857       [ (P0, v0), (P1, v1), ...]
    858 
    859       with Pi being lists of vertices defining polygons and vi either
    860       constants or functions of x,y to be applied to points with the polygon.
    861 
    862     The function takes an optional argument, default which is the value
    863     (or function) to used for points not belonging to any polygon.
    864     For example:
    865 
    866        Polygon_function(polygons, default = 0.03)
    867 
    868     If omitted the default value will be 0.0
    869 
    870     Note: If two polygons overlap, the one last in the list takes precedence
    871 
    872     Coordinates specified in the call are assumed to be relative to the
    873     origin (georeference) e.g. used by domain.
    874     By specifying the optional argument georeference,
    875     all points are made relative.
    876 
    877     FIXME: This should really work with geo_spatial point sets.
    878     """
    879 
    880     def __init__(self, regions, default=0.0, geo_reference=None):
    881         """
    882         # @brief Create instance of a polygon function.
    883         # @param regions A list of (x,y) tuples defining a polygon.
    884         # @param default Value or function returning value for points outside poly.
    885         # @param geo_reference ??
    886         """
    887 
    888         try:
    889             len(regions)
    890         except:
    891             msg = ('Polygon_function takes a list of pairs (polygon, value).'
    892                    'Got %s' % str(regions))
    893             raise Exception, msg
    894 
    895         T = regions[0]
    896 
    897         if isinstance(T, basestring):
    898             msg = ('You passed in a list of text values into polygon_function '
    899                    'instead of a list of pairs (polygon, value): "%s"'
    900                    % str(T))
    901             raise Exception, msg
    902 
    903         try:
    904             a = len(T)
    905         except:
    906             msg = ('Polygon_function takes a list of pairs (polygon, value). '
    907                    'Got %s' % str(T))
    908             raise Exception, msg
    909 
    910         msg = ('Each entry in regions have two components: (polygon, value). '
    911                'I got %s' % str(T))
    912         assert a == 2, msg
    913 
    914         if geo_reference is None:
    915             from anuga.coordinate_transforms.geo_reference import Geo_reference
    916             geo_reference = Geo_reference()
    917 
    918         self.default = default
    919 
    920         # Make points in polygons relative to geo_reference
    921         self.regions = []
    922         for polygon, value in regions:
    923             P = geo_reference.change_points_geo_ref(polygon)
    924             self.regions.append((P, value))
    925 
    926     def __call__(self, x, y):
    927         """
    928         # @brief Implement the 'callable' property of Polygon_function.
    929         # @param x List of x coordinates of points ot interest.
    930         # @param y List of y coordinates of points ot interest.
    931         """
    932         x = num.array(x, num.float)
    933         y = num.array(y, num.float)
    934 
    935         # x and y must be one-dimensional and same length
    936         assert len(x.shape) == 1 and len(y.shape) == 1
    937         N = x.shape[0]
    938         assert y.shape[0] == N
    939 
    940         points = num.ascontiguousarray(num.concatenate((x[:, num.newaxis],
    941                                                         y[:, num.newaxis]),
    942                                                        axis = 1 ))
    943 
    944         if callable(self.default):
    945             z = self.default(x, y)
    946         else:
    947             z = num.ones(N, num.float) * self.default
    948 
    949         for polygon, value in self.regions:
    950             indices = inside_polygon(points, polygon)
    951 
    952             # FIXME: This needs to be vectorised
    953             if callable(value):
    954                 for i in indices:
    955                     xx = num.array([x[i]])
    956                     yy = num.array([y[i]])
    957                     z[i] = value(xx, yy)[0]
    958             else:
    959                 for i in indices:
    960                     z[i] = value
    961 
    962         if len(z) == 0:
    963             msg = ('Warning: points provided to Polygon function did not fall '
    964                    'within its regions in [%.2f, %.2f], y in [%.2f, %.2f]'
    965                    % (min(x), max(x), min(y), max(y)))
    966             log.critical(msg)
    967 
    968         return z
     773    pts_x = num.concatenate((polygon[:, 0], [polygon[0, 0]]), axis = 0)
     774    pts_y = num.concatenate((polygon[:, 1], [polygon[0, 1]]), axis = 0)
     775
     776    return pts_x, pts_y
     777
    969778
    970779################################################################################
     
    973782
    974783def read_polygon(filename, delimiter=','):
    975     """Read points assumed to form a polygon.
    976 
    977 # @param filename Path to file containing polygon data.
    978 # @param delimiter Delimiter to split polygon data with.
    979 # @return A list of point data from the polygon file.
    980 
    981 
    982     There must be exactly two numbers in each line separated by the delimiter.
    983     No header.
     784    """ Read points assumed to form a polygon.
     785
     786        Also checks to make sure polygon is not complex (self-intersecting).
     787
     788        filename Path to file containing polygon data.
     789        delimiter Delimiter to split polygon data with.
     790        A list of point data from the polygon file.
     791
     792        There must be exactly two numbers in each line separated by the delimiter.
     793        No header.
    984794    """
    985795
     
    1002812    return polygon
    1003813
    1004 ##
    1005 # @brief Write polygon data to a file.
    1006 # @param polygon Polygon points to write to file.
    1007 # @param filename Path to file to write.
    1008 # @note Delimiter is assumed to be a comma.
     814
    1009815def write_polygon(polygon, filename=None):
    1010816    """Write polygon to csv file.
     
    1047853
    1048854    # Find outer extent of polygon
    1049     max_x = min_x = polygon[0][0]
    1050     max_y = min_y = polygon[0][1]
    1051     for point in polygon[1:]:
    1052         x = point[0]
    1053         if x > max_x: max_x = x
    1054         if x < min_x: min_x = x
    1055         y = point[1]
    1056         if y > max_y: max_y = y
    1057         if y < min_y: min_y = y
    1058 
     855    extents = AABB(polygon)
     856   
    1059857    while len(points) < number_of_points:
    1060         x = uniform(min_x, max_x)
    1061         y = uniform(min_y, max_y)
     858        rand_x = uniform(extents.xmin, extents.xmax)
     859        rand_y = uniform(extents.ymin, extents.ymax)
    1062860
    1063861        append = False
    1064         if is_inside_polygon([x,y], polygon):
     862        if is_inside_polygon([rand_x, rand_y], polygon):
    1065863            append = True
    1066864
     
    1068866            if exclude is not None:
    1069867                for ex_poly in exclude:
    1070                     if is_inside_polygon([x, y], ex_poly):
     868                    if is_inside_polygon([rand_x, rand_y], ex_poly):
    1071869                        append = False
    1072870
    1073871        if append is True:
    1074             points.append([x, y])
     872            points.append([rand_x, rand_y])
    1075873
    1076874    return points
     
    1091889    """
    1092890
    1093     import exceptions
    1094 
    1095     class Found(exceptions.Exception):
    1096         pass
    1097 
    1098891    polygon = ensure_numeric(polygon)
    1099892   
    1100     point_in = False
    1101     while not point_in:
    1102         try:
    1103             for poly_point in polygon:     # [1:]:
    1104                 for x_mult in range(-1, 2):
    1105                     for y_mult in range(-1, 2):
    1106                         x = poly_point[0]
    1107                         y = poly_point[1]
    1108 
    1109                         if x == 0:
    1110                             x_delta = x_mult * delta
    1111                         else:
    1112                             x_delta = x + x_mult*x*delta
    1113 
    1114                         if y == 0:
    1115                             y_delta = y_mult * delta
    1116                         else:
    1117                             y_delta = y + y_mult*y*delta
    1118 
    1119                         point = [x_delta, y_delta]
    1120 
    1121                         if is_inside_polygon(point, polygon, closed=False):
    1122                             raise Found
    1123         except Found:
    1124             point_in = True
    1125         else:
    1126             delta = delta * 0.1
    1127 
    1128     return point
     893    while True:
     894        for poly_point in polygon:
     895            for x_mult in range(-1, 2):
     896                for y_mult in range(-1, 2):
     897                    pt_x, pt_y = poly_point
     898
     899                    if pt_x == 0:
     900                        x_delta = x_mult * delta
     901                    else:
     902                        x_delta = pt_x + x_mult*pt_x*delta
     903
     904                    if pt_y == 0:
     905                        y_delta = y_mult * delta
     906                    else:
     907                        y_delta = pt_y + y_mult*pt_y*delta
     908
     909                    point = [x_delta, y_delta]
     910
     911                    if is_inside_polygon(point, polygon, closed=False):
     912                        return point
     913        delta = delta * 0.1
    1129914
    1130915
     
    1212997    reduced_polygon = []
    1213998    for i, point in enumerate(polygon):
    1214         x = point[0]
    1215         y = point[1]
    1216         if x in [min_x, max_x] and y in [min_y, max_y]:
     999        if point[0] in [min_x, max_x] and point[1] in [min_y, max_y]:
    12171000            # Keep
    12181001            reduced_polygon.append(point)
     
    13181101
    13191102else:
    1320     error_msg = 'C implementations could not be accessed by %s.\n ' %__file__
    1321     error_msg += 'Make sure compile_all.py has been run as described in '
    1322     error_msg += 'the ANUGA installation guide.'
    1323     raise Exception(error_msg)
     1103    ERROR_MSG = 'C implementations could not be accessed by %s.\n ' % __file__
     1104    ERROR_MSG += 'Make sure compile_all.py has been run as described in '
     1105    ERROR_MSG += 'the ANUGA installation guide.'
     1106    raise Exception(ERROR_MSG)
    13241107
    13251108
  • trunk/anuga_core/source/anuga/geometry/quad.py

    r7751 r7858  
    88may then be iterated over with a proper intersection test to detect actual
    99geometry intersections.
     10
     11As of June 2010 this module has a pylint quality rating of 10/10.
    1012
    1113"""
  • trunk/anuga_core/source/anuga/geometry/test_polygon.py

    r7841 r7858  
    88
    99from polygon import *
     10from polygon import _poly_xy
     11from polygon_function import Polygon_function
    1012from anuga.coordinate_transforms.geo_reference import Geo_reference
    1113from anuga.geospatial_data.geospatial_data import Geospatial_data
     
    15701572        # Simplest case: Polygon is the unit square
    15711573        polygon = [[0,0], [1,0], [1,1], [0,1]]
    1572         x, y = poly_xy(polygon)
     1574        x, y = _poly_xy(polygon)
    15731575        assert len(x) == len(polygon)+1
    15741576        assert len(y) == len(polygon)+1
     
    15841586    # Arbitrary polygon
    15851587        polygon = [[1,5], [1,1], [100,10], [1,10], [3,6]]
    1586         x, y = poly_xy(polygon)
     1588        x, y = _poly_xy(polygon)
    15871589        assert len(x) == len(polygon)+1
    15881590        assert len(y) == len(polygon)+1
     
    16051607        polygon1 = [[0,0], [1,0], [1,1], [0,1]]
    16061608        polygon2 = [[1,1], [2,1], [3,2], [2,2]]
    1607         v = plot_polygons([polygon1, polygon2], figname='test1')
    1608         assert len(v) == 4
    1609         assert v[0] == 0
    1610         assert v[1] == 3
    1611         assert v[2] == 0
    1612         assert v[3] == 2
     1609        plot_polygons([polygon1, polygon2], figname='test1')
    16131610
    16141611        # Another case
    16151612        polygon3 = [[1,5], [10,1], [100,10], [50,10], [3,6]]
    16161613        v = plot_polygons([polygon2,polygon3], figname='test2')
    1617         assert len(v) == 4
    1618         assert v[0] == 1
    1619         assert v[1] == 100
    1620         assert v[2] == 1
    1621         assert v[3] == 10
    1622 
    1623         os.remove('test1.png')
    1624         os.remove('test2.png')
     1614
     1615        for file in ['test1.png', 'test2.png']:
     1616            assert os.access(file, os.R_OK)
     1617            os.remove(file)
     1618
    16251619
    16261620    def test_inside_polygon_geospatial(self):
  • trunk/anuga_core/source/anuga/pmesh/view_tsh.py

    r7711 r7858  
    1010from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
    1111from anuga.pyvolution.util import file_function
    12 from anuga.geometry.polygon import Polygon_function, read_polygon
     12from anuga.geometry.polygon import read_polygon
     13from anuga.geometry.polygon_function import Polygon_function
    1314from realtime_visualisation_new import *
    1415
  • trunk/anuga_core/source/anuga/shallow_water/sww_interrogate.py

    r7778 r7858  
    44
    55import os
     6import anuga.utilities.log as log
    67import numpy as num
    78
Note: See TracChangeset for help on using the changeset viewer.