Index: /trunk/anuga_core/source/anuga/__init__.py
===================================================================
--- /trunk/anuga_core/source/anuga/__init__.py	(revision 7857)
+++ /trunk/anuga_core/source/anuga/__init__.py	(revision 7858)
@@ -42,5 +42,5 @@
 
 from anuga.geometry.polygon import read_polygon, plot_polygons, polygon_area
-from anuga.geometry.polygon import Polygon_function
+from anuga.geometry.polygon_function import Polygon_function
 
 from anuga.abstract_2d_finite_volumes.pmesh2domain import \
@@ -52,5 +52,5 @@
 from anuga.utilities.file_utils import copy_code_files
 
-from anuga.geometry.polygon import read_polygon, Polygon_function
+from anuga.geometry.polygon import read_polygon
 from anuga.caching import cache
 
Index: /trunk/anuga_core/source/anuga/culvert_flows/test_culvert_class.py
===================================================================
--- /trunk/anuga_core/source/anuga/culvert_flows/test_culvert_class.py	(revision 7857)
+++ /trunk/anuga_core/source/anuga/culvert_flows/test_culvert_class.py	(revision 7858)
@@ -7,5 +7,5 @@
 
 from anuga.utilities.system_tools import get_pathname_from_package
-from anuga.geometry.polygon import Polygon_function
+from anuga.geometry.polygon_function import Polygon_function
         
 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
Index: /trunk/anuga_core/source/anuga/culvert_flows/test_new_culvert_class.py
===================================================================
--- /trunk/anuga_core/source/anuga/culvert_flows/test_new_culvert_class.py	(revision 7857)
+++ /trunk/anuga_core/source/anuga/culvert_flows/test_new_culvert_class.py	(revision 7858)
@@ -9,5 +9,5 @@
 
 from anuga.utilities.system_tools import get_pathname_from_package
-from anuga.geometry.polygon import Polygon_function
+from anuga.geometry.polygon_function import Polygon_function
         
 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
Index: /trunk/anuga_core/source/anuga/file/csv_file.py
===================================================================
--- /trunk/anuga_core/source/anuga/file/csv_file.py	(revision 7857)
+++ /trunk/anuga_core/source/anuga/file/csv_file.py	(revision 7858)
@@ -11,6 +11,8 @@
 """
 
+
 import csv
 import numpy as num
+import anuga.utilities.log as log
 
 
@@ -168,13 +170,13 @@
     if completed:
         try:
-            file = str(kwargs['file_name'])
+            file_name = str(kwargs['file_name'])
         except:
-            raise 'kwargs must have file_name'
+            raise Exception('kwargs must have file_name')
     else:
         # write temp file in output directory
         try:
-            file = str(kwargs['output_dir']) + 'detail_temp.csv'
+            file_name = str(kwargs['output_dir']) + 'detail_temp.csv'
         except:
-            raise 'kwargs must have output_dir'
+            raise Exception('kwargs must have output_dir')
 
     # extracts the header info and the new line info
@@ -199,14 +201,15 @@
     # try to open!
     try:
-        fid = open(file, 'r')
+        fid = open(file_name, 'r')
         file_header = fid.readline()
         fid.close()
         if verbose: log.critical('read file header %s' % file_header)
-    except:
-        msg = 'try to create new file: %s' % file
-        if verbose: log.critical(msg)
+    except Exception:
+        msg = 'try to create new file: %s' % file_name
+        if verbose:
+            log.critical(msg)
         #tries to open file, maybe directory is bad
         try:
-            fid = open(file, 'w')
+            fid = open(file_name, 'w')
             fid.write(header)
             fid.close()
@@ -218,5 +221,5 @@
     # if header is same or this is a new file
     if file_header == str(header):
-        fid = open(file, 'a')
+        fid = open(file_name, 'a')
         fid.write(line)
         fid.close()
@@ -225,6 +228,6 @@
         # if header is different and has completed will append info to
         # end of details_temp.cvs file in output directory
-        file = str(kwargs['output_dir']) + 'detail_temp.csv'
-        fid = open(file, 'a')
+        file_name = str(kwargs['output_dir']) + 'detail_temp.csv'
+        fid = open(file_name, 'a')
         fid.write(header)
         fid.write(line)
@@ -244,6 +247,5 @@
 
 def load_csv_as_building_polygons(file_name,
-                          floor_height=3,
-                          clipping_polygons=None):
+                          floor_height=3):
     """
     Convert CSV files of the form:
@@ -285,7 +287,4 @@
             
 
-##
-# @brief Convert CSV file into a dictionary of polygons and associated values.
-# @param filename The path to the file to read, value_name name for the 4th column
 def load_csv_as_polygons(file_name,
                  value_name='value',
@@ -338,5 +337,5 @@
     
     msg = 'Did not find expected column header: northing'    
-    assert 'northing' in X.keys(), northing
+    assert 'northing' in X.keys(), msg
     
     msg = 'Did not find expected column header: northing'        
@@ -357,16 +356,16 @@
     past_ids = {}
     last_id = None
-    for i, id in enumerate(X['id']):
+    for i, poly_id in enumerate(X['id']):
 
         # Check for duplicate polygons
-        if id in past_ids:
+        if poly_id in past_ids:
             msg = 'Polygon %s was duplicated in line %d' % (id, i)
             raise Exception, msg
         
-        if id not in polygons:
+        if poly_id not in polygons:
             # Start new polygon
-            polygons[id] = []
+            polygons[poly_id] = []
             if values is not None:
-                values[id] = X[value_name][i]
+                values[poly_id] = X[value_name][i]
 
             # Keep track of previous polygon ids
@@ -385,18 +384,19 @@
                 
             if exclude is True:
-                excluded_polygons[id]=True
-
-        polygons[id].append(point)    
+                excluded_polygons[poly_id]=True
+
+        polygons[poly_id].append(point)    
             
         # Check that value is the same across each polygon
         msg = 'Values must be the same across each polygon.'
-        msg += 'I got %s in line %d but it should have been %s' % (X[value_name][i], i, values[id])
-        assert values[id] == X[value_name][i], msg
-
-        last_id = id
+        msg += 'I got %s in line %d but it should have been %s' % \
+                            (X[value_name][i], i, values[poly_id])
+        assert values[poly_id] == X[value_name][i], msg
+
+        last_id = poly_id
 
     # Weed out polygons that were not wholly inside clipping polygons
-    for id in excluded_polygons:
-        del polygons[id]
+    for poly_id in excluded_polygons:
+        del polygons[poly_id]
         
     return polygons, values
Index: /trunk/anuga_core/source/anuga/file/sts.py
===================================================================
--- /trunk/anuga_core/source/anuga/file/sts.py	(revision 7857)
+++ /trunk/anuga_core/source/anuga/file/sts.py	(revision 7858)
@@ -94,25 +94,5 @@
         outfile.variables[q + Write_sts.RANGE][1] = -max_float # Max
 
-        # Doing sts_precision instead of Float gives cast errors.
-        outfile.createVariable('time', netcdf_float, ('number_of_timesteps',))
-
-        for q in Write_sts.sts_quantities:
-            outfile.createVariable(q, sts_precision, ('number_of_timesteps',
-                                                      'number_of_points'))
-            outfile.createVariable(q + Write_sts.RANGE, sts_precision,
-                                   ('numbers_in_range',))
-            # Initialise ranges with small and large sentinels.
-            # If this was in pure Python we could have used None sensibly
-            outfile.variables[q + Write_sts.RANGE][0] = max_float  # Min
-            outfile.variables[q + Write_sts.RANGE][1] = -max_float # Max
-
-        if isinstance(times, (list, num.ndarray)):
-            outfile.variables['time'][:] = times    #Store time relative
-
-        if verbose:
-            log.critical('------------------------------------------------')
-            log.critical('Statistics:')
-            log.critical('    t in [%f, %f], len(t) == %d'
-                         % (num.min(times), num.max(times), len(times.flat)))
+        self.write_dynamic_quantities(outfile, Write_sts.sts_quantities, times)
 
     ##
@@ -265,4 +245,32 @@
 
 
+    def write_dynamic_quantities(self, outfile, quantities,
+                    times, precis = netcdf_float32, verbose = False):   
+        """
+            Write out given quantities to file.
+        """
+        for q in quantities:
+            outfile.createVariable(q, precis, ('number_of_timesteps',
+                                                      'number_of_points'))
+            outfile.createVariable(q + Write_sts.RANGE, precis,
+                                   ('numbers_in_range',))
+
+            # Initialise ranges with small and large sentinels.
+            # If this was in pure Python we could have used None sensibly
+            outfile.variables[q+Write_sts.RANGE][0] = max_float  # Min
+            outfile.variables[q+Write_sts.RANGE][1] = -max_float # Max
+
+        # Doing sts_precision instead of Float gives cast errors.
+        outfile.createVariable('time', netcdf_float, ('number_of_timesteps',))
+
+        if isinstance(times, (list, num.ndarray)):
+            outfile.variables['time'][:] = times    # Store time relative
+
+        if verbose:
+            log.critical('------------------------------------------------')
+            log.critical('Statistics:')
+            log.critical('    t in [%f, %f], len(t) == %d'
+                         % (num.min(times), num.max(times), len(times.flat)))
+
 
 
Index: /trunk/anuga_core/source/anuga/file/sww.py
===================================================================
--- /trunk/anuga_core/source/anuga/file/sww.py	(revision 7857)
+++ /trunk/anuga_core/source/anuga/file/sww.py	(revision 7858)
@@ -9,4 +9,6 @@
 import anuga.utilities.log as log
 from Scientific.IO.NetCDF import NetCDFFile
+
+from sts import Write_sts
 
 from anuga.coordinate_transforms.geo_reference import \
@@ -455,5 +457,5 @@
 
 
-class Write_sww:
+class Write_sww(Write_sts):
     """
         A class to write an SWW file.
@@ -462,6 +464,4 @@
         manually.
     """
-    RANGE = '_range'
-    EXTREMA = ':extrema'
 
     def __init__(self, static_quantities, dynamic_quantities):
@@ -567,9 +567,5 @@
                                                        'number_of_vertices'))
 
-        # Doing sww_precision instead of Float gives cast errors.
-        outfile.createVariable('time', netcdf_float,
-                               ('number_of_timesteps',))
-
-                               
+
         for q in self.static_quantities:
             
@@ -585,28 +581,7 @@
             outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
 
-        #if 'elevation' in self.static_quantities:    
-        #    # FIXME: Backwards compat - get rid of z once old view has retired
-        #    outfile.createVariable('z', sww_precision,
-        #                           ('number_of_points',))
-                               
-        for q in self.dynamic_quantities:
-            outfile.createVariable(q, sww_precision, ('number_of_timesteps',
-                                                      'number_of_points'))
-            outfile.createVariable(q + Write_sww.RANGE, sww_precision,
-                                   ('numbers_in_range',))
-
-            # Initialise ranges with small and large sentinels.
-            # If this was in pure Python we could have used None sensibly
-            outfile.variables[q+Write_sww.RANGE][0] = max_float  # Min
-            outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
-
-        if isinstance(times, (list, num.ndarray)):
-            outfile.variables['time'][:] = times    # Store time relative
-
-        if verbose:
-            log.critical('------------------------------------------------')
-            log.critical('Statistics:')
-            log.critical('    t in [%f, %f], len(t) == %d'
-                         % (num.min(times), num.max(times), len(times.flat)))
+        
+        self.write_dynamic_quantities(outfile, self.dynamic_quantities, times)
+
 
 
Index: /trunk/anuga_core/source/anuga/file_conversion/csv2sts.py
===================================================================
--- /trunk/anuga_core/source/anuga/file_conversion/csv2sts.py	(revision 7857)
+++ /trunk/anuga_core/source/anuga/file_conversion/csv2sts.py	(revision 7858)
@@ -3,4 +3,7 @@
 """
     Module to convert a .csv file to an .sts file.
+    
+    You can use the function from your own scripts, or call this script on the 
+    command line.
     
     Use this module to convert an arbitrary csv file to an sts file.
@@ -42,18 +45,16 @@
         data:
 
-         stage = 4, 150.66667, 150.83334, 151, 151.16667, -34, -34.16667, -34.33333, 
-            -34.5, -1, -5, -9, -13 ;
+         stage = 4, 150.66667, 150.83334, 151, 151.16667, -34, -34.16667,
+                -34.33333, -34.5, -1, -5, -9, -13 ;
 
          time = 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 ;
         }
     
+    As of June 2010 this module has a pylint quality rating of 9.30/10.
 """
 
-import csv
-import os
 import sys
 import getopt
 from anuga.utilities import log
-import numpy as num
 from Scientific.IO.NetCDF import NetCDFFile
 from anuga.file.csv_file import load_csv_as_dict
@@ -108,4 +109,5 @@
 
 def usage():
+    """ Display usage of this module from the comand line. """
     print 'csv2sts - convert a csv file to an sts file.'
     print 'Usage: csv2sts [-hv] [--help] [--verbose]',
@@ -161,4 +163,3 @@
     
 if __name__ == "__main__":
-    """ Entry point if run from command line """
     main(sys.argv[1:])       
Index: /trunk/anuga_core/source/anuga/file_conversion/urs2sww.py
===================================================================
--- /trunk/anuga_core/source/anuga/file_conversion/urs2sww.py	(revision 7857)
+++ /trunk/anuga_core/source/anuga/file_conversion/urs2sww.py	(revision 7858)
@@ -24,24 +24,6 @@
 from anuga.file.sww import Write_sww  
 
-
-################################################################################
-# CONVERTING UNGRIDDED URS DATA TO AN SWW FILE
-################################################################################
-
-##
-# @brief Convert URS file(s) (wave prop) to an SWW file.
-# @param basename_in Stem of the input filenames.
-# @param basename_out Path to the output SWW file.
-# @param verbose True if this function is to be verbose.
-# @param mint
-# @param maxt
-# @param mean_stage
-# @param origin Tuple with geo-ref UTM coordinates (zone, easting, northing).
-# @param hole_points_UTM
-# @param zscale
-# @note Also convert latitude and longitude to UTM. All coordinates are
-#       assumed to be given in the GDA94 datum.
-# @note Input filename stem has suffixes '-z-mux', '-e-mux' and '-n-mux'
-#       added for relative height, x-velocity and y-velocity respectively.
+###############################################################
+
 def urs_ungridded2sww(basename_in='o', basename_out=None, verbose=False,
                       mint=None, maxt=None,
Index: /trunk/anuga_core/source/anuga/geometry/aabb.py
===================================================================
--- /trunk/anuga_core/source/anuga/geometry/aabb.py	(revision 7857)
+++ /trunk/anuga_core/source/anuga/geometry/aabb.py	(revision 7858)
@@ -3,4 +3,6 @@
     Contains a class describing a bounding box. It contains methods
     to split itself and return child boxes.
+    
+    As of June 2010 this module has a pylint quality rating of 10/10.
 """
 
@@ -15,5 +17,5 @@
     """
     
-    def __init__(self, xmin, xmax, ymin, ymax):
+    def __init__(self, xmin, xmax=None, ymin=None, ymax=None):
         """ Define axially-algned bounding box.
             xmin is minimum x
@@ -22,8 +24,18 @@
             ymax is maximum y (absolute coord, ie, not size)
         """
-        self.xmin = xmin    
-        self.xmax = xmax
-        self.ymin = ymin    
-        self.ymax = ymax
+        if not xmax:
+            # try treating first arg as a list of points
+            try:
+                xmin[0][0]
+            except:
+                raise Exception('Single parameter to AABB must be point list.')
+                
+            self.xmin, self.ymin = self.xmax, self.ymax = xmin[0]
+            self.include(xmin[1:])
+        else:
+            self.xmin = xmin    
+            self.xmax = xmax
+            self.ymin = ymin    
+            self.ymax = ymax
 
 
@@ -96,2 +108,18 @@
                 and (self.ymin <= point[1] <= self.ymax)
         
+    def include(self, point_list):
+        """ Include points in AABB.
+            Bounding box will be expanded to include passed points
+        
+            point_list is a list of points.
+        """
+        for point in point_list:
+            pt_x, pt_y = point
+            if pt_x < self.xmin:
+                self.xmin = pt_x
+            if pt_x > self.xmax:
+                self.xmax = pt_x
+            if pt_y < self.ymin:
+                self.ymin = pt_y                
+            if pt_y > self.ymax:
+                self.ymax = pt_y                
Index: /trunk/anuga_core/source/anuga/geometry/polygon.py
===================================================================
--- /trunk/anuga_core/source/anuga/geometry/polygon.py	(revision 7857)
+++ /trunk/anuga_core/source/anuga/geometry/polygon.py	(revision 7858)
@@ -10,4 +10,5 @@
 import anuga.utilities.log as log
 
+from aabb import AABB
 
 ##
@@ -36,6 +37,6 @@
 
     res = _point_on_line(point[0], point[1],
-                         line[0,0], line[0,1],
-                         line[1,0], line[1,1],
+                         line[0, 0], line[0, 1],
+                         line[1, 0], line[1, 1],
                          rtol, atol)
 
@@ -49,17 +50,24 @@
 
 # result functions for possible states
-def lines_dont_coincide(p0,p1,p2,p3):               return (3, None)
-def lines_0_fully_included_in_1(p0,p1,p2,p3):       return (2,
-                                                            num.array([p0,p1]))
-def lines_1_fully_included_in_0(p0,p1,p2,p3):       return (2,
-                                                            num.array([p2,p3]))
-def lines_overlap_same_direction(p0,p1,p2,p3):      return (2,
-                                                            num.array([p0,p3]))
-def lines_overlap_same_direction2(p0,p1,p2,p3):     return (2,
-                                                            num.array([p2,p1]))
-def lines_overlap_opposite_direction(p0,p1,p2,p3):  return (2,
-                                                            num.array([p0,p2]))
-def lines_overlap_opposite_direction2(p0,p1,p2,p3): return (2,
-                                                            num.array([p3,p1]))
+def lines_dont_coincide(p0, p1, p2, p3):
+    return (3, None)
+    
+def lines_0_fully_included_in_1(p0, p1, p2, p3):
+    return (2, num.array([p0, p1]))
+    
+def lines_1_fully_included_in_0(p0, p1, p2, p3):
+    return (2, num.array([p2, p3]))
+    
+def lines_overlap_same_direction(p0, p1, p2, p3):
+    return (2, num.array([p0, p3]))
+    
+def lines_overlap_same_direction2(p0, p1, p2, p3):
+    return (2, num.array([p2, p1]))
+    
+def lines_overlap_opposite_direction(p0, p1, p2, p3):
+    return (2, num.array([p0, p2]))
+    
+def lines_overlap_opposite_direction2(p0, p1, p2, p3):
+    return (2, num.array([p3, p1]))
 
 # this function called when an impossible state is found
@@ -146,6 +154,6 @@
                            point_on_line([x3, y3], line0, rtol=rtol, atol=atol))
 
-            return collinear_result[state_tuple]([x0,y0], [x1,y1],
-                                                 [x2,y2], [x3,y3])
+            return collinear_result[state_tuple]([x0, y0], [x1, y1],
+                                                 [x2, y2], [x3, y3])
         else:
             # Lines are parallel but aren't collinear
@@ -269,62 +277,4 @@
     return bool(_is_inside_triangle(point, triangle, int(closed), rtol, atol))
     
-    
-
-    # FIXME (Ole): The rest of this function has been made 
-    # obsolete by the C extension.
-    
-    # Quickly reject points that are clearly outside
-    if point[0] < min(triangle[:, 0]):
-        return False 
-    if point[0] > max(triangle[:, 0]): 
-        return False    
-    if point[1] < min(triangle[:, 1]):
-        return False
-    if point[1] > max(triangle[:, 1]):
-        return False        
-
-
-    # Start search    
-    A = triangle[0, :]
-    B = triangle[1, :]
-    C = triangle[2, :]
-    
-    # Now check if point lies wholly inside triangle
-    v0 = C-A 
-    v1 = B-A        
-        
-    a00 = num.inner(v0, v0)
-    a10 = a01 = num.inner(v0, v1)
-    a11 = num.inner(v1, v1)
-    
-    denom = a11*a00 - a01*a10
-    
-    if abs(denom) > 0.0:
-        v = point-A
-        b0 = num.inner(v0, v)        
-        b1 = num.inner(v1, v)            
-    
-        alpha = (b0*a11 - b1*a01)/denom
-        beta = (b1*a00 - b0*a10)/denom        
-
-        if (alpha > 0.0) and (beta > 0.0) and (alpha+beta < 1.0):
-            return True 
-
-
-    if closed is True:
-        # Check if point lies on one of the edges
-        
-        for X, Y in [[A, B], [B, C], [C, A]]:
-            res = _point_on_line(point[0], point[1],
-                                 X[0], X[1],
-                                 Y[0], Y[1],
-                                 rtol, atol)
-            
-            if res:
-                return True
-                
-    return False
-
-    
 def is_complex(polygon, verbose=False):
     """Check if a polygon is complex (self-intersecting).
@@ -339,4 +289,5 @@
             
     def key_xpos(item):
+        """ Return the x coord out of the passed point for sorting key. """
         return (item[0][0])
     
@@ -377,5 +328,5 @@
                         print 'Self-intersecting polygon found, type ', type
                         print 'point', point,
-                        print 'vertices: ', leftmost, ' - ', l_x[cmp]                
+                        print 'vertices: ', leftmost, ' - ', l_x[cmp]  
                     return True            
             cmp += 1
@@ -422,6 +373,6 @@
     try:
         points = ensure_absolute(points)
-    except NameError, e:
-        raise NameError, e
+    except NameError, err:
+        raise NameError, err
     except:
         # If this fails it is going to be because the points can't be
@@ -513,5 +464,5 @@
     if len(points.shape) == 1:
         # Only one point was passed in. Convert to array of points
-        points = num.reshape(points, (1,2))
+        points = num.reshape(points, (1, 2))
 
     indices, count = separate_points_by_polygon(points, polygon,
@@ -559,5 +510,5 @@
     if len(points.shape) == 1:
         # Only one point was passed in. Convert to array of points
-        points = num.reshape(points, (1,2))
+        points = num.reshape(points, (1, 2))
 
     indices, count = separate_points_by_polygon(points, polygon,
@@ -632,13 +583,13 @@
         except:
             msg = 'Points could not be converted to numeric array'
-            raise msg
+            raise Exception(msg)
 
         try:
             polygon = ensure_numeric(polygon, num.float)
         except NameError, e:
-            raise NameError, e
+            raise NameError(e)
         except:
             msg = 'Polygon could not be converted to numeric array'
-            raise msg
+            raise Exception(msg)
 
         msg = 'Polygon array must be a 2d array of vertices'
@@ -646,5 +597,5 @@
 
         msg = 'Polygon array must have two columns' 
-        assert polygon.shape[1]==2, msg
+        assert polygon.shape[1] == 2, msg
 
         msg = ('Points array must be 1 or 2 dimensional. '
@@ -654,5 +605,5 @@
         if len(points.shape) == 1:
             # Only one point was passed in.  Convert to array of points.
-            points = num.reshape(points, (1,2))
+            points = num.reshape(points, (1, 2))
     
             msg = ('Point array must have two columns (x,y), '
@@ -663,8 +614,8 @@
             msg = ('Points array must be a 2d array. I got %s.'
                    % str(points[:30]))
-            assert len(points.shape)==2, msg
+            assert len(points.shape) == 2, msg
 
             msg = 'Points array must have two columns'
-            assert points.shape[1]==2, msg
+            assert points.shape[1] == 2, msg
 
     N = polygon.shape[0] # Number of vertices in polygon
@@ -681,19 +632,19 @@
     return indices, count
 
-##
-# @brief Determine area of a polygon.
-# @param input_polygon The polygon to get area of.
-# @return A scalar value for the polygon area.
+
 def polygon_area(input_polygon):
     """ Determine area of arbitrary polygon.
 
-    Reference:     http://mathworld.wolfram.com/PolygonArea.html
-    """
-
+        input_polygon The polygon to get area of.
+        
+        return A scalar value for the polygon area.
+
+        Reference:     http://mathworld.wolfram.com/PolygonArea.html
+    """
     # Move polygon to origin (0,0) to avoid rounding errors
     # This makes a copy of the polygon to avoid destroying it
     input_polygon = ensure_numeric(input_polygon)
-    min_x = min(input_polygon[:,0])
-    min_y = min(input_polygon[:,1])
+    min_x = min(input_polygon[:, 0])
+    min_y = min(input_polygon[:, 1])
     polygon = input_polygon - [min_x, min_y]
 
@@ -742,5 +693,4 @@
     Outputs:
 
-    - list of min and max of x and y coordinates
     - plot of polygons
     """
@@ -754,9 +704,4 @@
     ion()
     hold(True)
-
-    minx = 1e10
-    maxx = 0.0
-    miny = 1e10
-    maxy = 0.0
 
     if label is None:
@@ -772,35 +717,28 @@
             alpha = max(0.0, min(1.0, alpha))
 
-    n = len(polygons_points)
+    num_points = len(polygons_points)
     colour = []
     if style is None:
         style_type = 'line'
         style = []
-        for i in range(n):
+        for i in range(num_points):
             style.append(style_type)
             colour.append('b-')
     else:
-        for s in style:
-            if s == 'line': colour.append('b-')
-            if s == 'outside': colour.append('r.')
-            if s == 'point': colour.append('g.')
-            if s != 'line':
-                if s != 'outside':
-                    if s != 'point':
-                        colour.append(s)
+        for style_name in style:
+            if style_name == 'line':
+                colour.append('b-')
+            if style_name == 'outside':
+                colour.append('r.')
+            if style_name == 'point':
+                colour.append('g.')
+            if style_name not in ['line', 'outside', 'point']:
+                colour.append(style_name)
 
     for i, item in enumerate(polygons_points):
-        x, y = poly_xy(item)
-        if min(x) < minx:
-            minx = min(x)
-        if max(x) > maxx:
-            maxx = max(x)
-        if min(y) < miny:
-            miny = min(y)
-        if max(y) > maxy:
-            maxy = max(y)
-        plot(x, y, colour[i])
+        pt_x, pt_y = _poly_xy(item)
+        plot(pt_x, pt_y, colour[i])
         if alpha:
-            fill(x, y, colour[i], alpha=alpha)
+            fill(pt_x, pt_y, colour[i], alpha=alpha)
         xlabel('x')
         ylabel('y')
@@ -814,9 +752,6 @@
     close('all')
 
-    vec = [minx, maxx, miny, maxy]
-    return vec
-
-
-def poly_xy(polygon):
+
+def _poly_xy(polygon):
     """ this is used within plot_polygons so need to duplicate
         the first point so can have closed polygon in plot
@@ -829,6 +764,6 @@
     try:
         polygon = ensure_numeric(polygon, num.float)
-    except NameError, e:
-        raise NameError, e
+    except NameError, err:
+        raise NameError, err
     except:
         msg = ('Polygon %s could not be converted to numeric array'
@@ -836,135 +771,9 @@
         raise Exception, msg
 
-    x = polygon[:, 0]
-    y = polygon[:, 1]
-    x = num.concatenate((x, [polygon[0, 0]]), axis = 0)
-    y = num.concatenate((y, [polygon[0, 1]]), axis = 0)
-
-    return x, y
-
-
-class Polygon_function:
-    """Create callable object f: x,y -> z, where a,y,z are vectors and
-    where f will return different values depending on whether x,y belongs
-    to specified polygons.
-
-    To instantiate:
-
-       Polygon_function(polygons)
-
-    where polygons is a list of tuples of the form
-
-      [ (P0, v0), (P1, v1), ...]
-
-      with Pi being lists of vertices defining polygons and vi either
-      constants or functions of x,y to be applied to points with the polygon.
-
-    The function takes an optional argument, default which is the value
-    (or function) to used for points not belonging to any polygon.
-    For example:
-
-       Polygon_function(polygons, default = 0.03)
-
-    If omitted the default value will be 0.0
-
-    Note: If two polygons overlap, the one last in the list takes precedence
-
-    Coordinates specified in the call are assumed to be relative to the
-    origin (georeference) e.g. used by domain.
-    By specifying the optional argument georeference,
-    all points are made relative.
-
-    FIXME: This should really work with geo_spatial point sets.
-    """
-
-    def __init__(self, regions, default=0.0, geo_reference=None):
-        """
-        # @brief Create instance of a polygon function.
-        # @param regions A list of (x,y) tuples defining a polygon.
-        # @param default Value or function returning value for points outside poly.
-        # @param geo_reference ??
-        """
-
-        try:
-            len(regions)
-        except:
-            msg = ('Polygon_function takes a list of pairs (polygon, value).'
-                   'Got %s' % str(regions))
-            raise Exception, msg
-
-        T = regions[0]
-
-        if isinstance(T, basestring):
-            msg = ('You passed in a list of text values into polygon_function '
-                   'instead of a list of pairs (polygon, value): "%s"'
-                   % str(T))
-            raise Exception, msg
-
-        try:
-            a = len(T)
-        except:
-            msg = ('Polygon_function takes a list of pairs (polygon, value). '
-                   'Got %s' % str(T))
-            raise Exception, msg
-
-        msg = ('Each entry in regions have two components: (polygon, value). '
-               'I got %s' % str(T))
-        assert a == 2, msg
-
-        if geo_reference is None:
-            from anuga.coordinate_transforms.geo_reference import Geo_reference
-            geo_reference = Geo_reference()
-
-        self.default = default
-
-        # Make points in polygons relative to geo_reference
-        self.regions = []
-        for polygon, value in regions:
-            P = geo_reference.change_points_geo_ref(polygon)
-            self.regions.append((P, value))
-
-    def __call__(self, x, y):
-        """
-        # @brief Implement the 'callable' property of Polygon_function.
-        # @param x List of x coordinates of points ot interest.
-        # @param y List of y coordinates of points ot interest.
-        """
-        x = num.array(x, num.float)
-        y = num.array(y, num.float)
-
-        # x and y must be one-dimensional and same length
-        assert len(x.shape) == 1 and len(y.shape) == 1
-        N = x.shape[0]
-        assert y.shape[0] == N
-
-        points = num.ascontiguousarray(num.concatenate((x[:, num.newaxis],
-                                                        y[:, num.newaxis]),
-                                                       axis = 1 ))
-
-        if callable(self.default):
-            z = self.default(x, y)
-        else:
-            z = num.ones(N, num.float) * self.default
-
-        for polygon, value in self.regions:
-            indices = inside_polygon(points, polygon)
-
-            # FIXME: This needs to be vectorised
-            if callable(value):
-                for i in indices:
-                    xx = num.array([x[i]])
-                    yy = num.array([y[i]])
-                    z[i] = value(xx, yy)[0]
-            else:
-                for i in indices:
-                    z[i] = value
-
-        if len(z) == 0:
-            msg = ('Warning: points provided to Polygon function did not fall '
-                   'within its regions in [%.2f, %.2f], y in [%.2f, %.2f]'
-                   % (min(x), max(x), min(y), max(y)))
-            log.critical(msg)
-
-        return z
+    pts_x = num.concatenate((polygon[:, 0], [polygon[0, 0]]), axis = 0)
+    pts_y = num.concatenate((polygon[:, 1], [polygon[0, 1]]), axis = 0)
+
+    return pts_x, pts_y
+
 
 ################################################################################
@@ -973,13 +782,14 @@
 
 def read_polygon(filename, delimiter=','):
-    """Read points assumed to form a polygon.
-
-# @param filename Path to file containing polygon data.
-# @param delimiter Delimiter to split polygon data with.
-# @return A list of point data from the polygon file.
-
-
-    There must be exactly two numbers in each line separated by the delimiter.
-    No header.
+    """ Read points assumed to form a polygon.
+
+        Also checks to make sure polygon is not complex (self-intersecting).
+
+        filename Path to file containing polygon data.
+        delimiter Delimiter to split polygon data with.
+        A list of point data from the polygon file.
+
+        There must be exactly two numbers in each line separated by the delimiter.
+        No header.
     """
 
@@ -1002,9 +812,5 @@
     return polygon
 
-##
-# @brief Write polygon data to a file.
-# @param polygon Polygon points to write to file.
-# @param filename Path to file to write.
-# @note Delimiter is assumed to be a comma.
+
 def write_polygon(polygon, filename=None):
     """Write polygon to csv file.
@@ -1047,20 +853,12 @@
 
     # Find outer extent of polygon
-    max_x = min_x = polygon[0][0]
-    max_y = min_y = polygon[0][1]
-    for point in polygon[1:]:
-        x = point[0]
-        if x > max_x: max_x = x
-        if x < min_x: min_x = x
-        y = point[1]
-        if y > max_y: max_y = y
-        if y < min_y: min_y = y
-
+    extents = AABB(polygon)
+    
     while len(points) < number_of_points:
-        x = uniform(min_x, max_x)
-        y = uniform(min_y, max_y)
+        rand_x = uniform(extents.xmin, extents.xmax)
+        rand_y = uniform(extents.ymin, extents.ymax)
 
         append = False
-        if is_inside_polygon([x,y], polygon):
+        if is_inside_polygon([rand_x, rand_y], polygon):
             append = True
 
@@ -1068,9 +866,9 @@
             if exclude is not None:
                 for ex_poly in exclude:
-                    if is_inside_polygon([x, y], ex_poly):
+                    if is_inside_polygon([rand_x, rand_y], ex_poly):
                         append = False
 
         if append is True:
-            points.append([x, y])
+            points.append([rand_x, rand_y])
 
     return points
@@ -1091,40 +889,27 @@
     """
 
-    import exceptions
-
-    class Found(exceptions.Exception):
-        pass
-
     polygon = ensure_numeric(polygon)
     
-    point_in = False
-    while not point_in:
-        try:
-            for poly_point in polygon:     # [1:]:
-                for x_mult in range(-1, 2):
-                    for y_mult in range(-1, 2):
-                        x = poly_point[0]
-                        y = poly_point[1]
-
-                        if x == 0:
-                            x_delta = x_mult * delta
-                        else:
-                            x_delta = x + x_mult*x*delta
-
-                        if y == 0:
-                            y_delta = y_mult * delta
-                        else:
-                            y_delta = y + y_mult*y*delta
-
-                        point = [x_delta, y_delta]
-
-                        if is_inside_polygon(point, polygon, closed=False):
-                            raise Found
-        except Found:
-            point_in = True
-        else:
-            delta = delta * 0.1
-
-    return point
+    while True:
+        for poly_point in polygon:
+            for x_mult in range(-1, 2):
+                for y_mult in range(-1, 2):
+                    pt_x, pt_y = poly_point
+
+                    if pt_x == 0:
+                        x_delta = x_mult * delta
+                    else:
+                        x_delta = pt_x + x_mult*pt_x*delta
+
+                    if pt_y == 0:
+                        y_delta = y_mult * delta
+                    else:
+                        y_delta = pt_y + y_mult*pt_y*delta
+
+                    point = [x_delta, y_delta]
+
+                    if is_inside_polygon(point, polygon, closed=False):
+                        return point
+        delta = delta * 0.1
 
 
@@ -1212,7 +997,5 @@
     reduced_polygon = []
     for i, point in enumerate(polygon):
-        x = point[0]
-        y = point[1]
-        if x in [min_x, max_x] and y in [min_y, max_y]:
+        if point[0] in [min_x, max_x] and point[1] in [min_y, max_y]:
             # Keep
             reduced_polygon.append(point)
@@ -1318,8 +1101,8 @@
 
 else:
-    error_msg = 'C implementations could not be accessed by %s.\n ' %__file__
-    error_msg += 'Make sure compile_all.py has been run as described in '
-    error_msg += 'the ANUGA installation guide.'
-    raise Exception(error_msg)
+    ERROR_MSG = 'C implementations could not be accessed by %s.\n ' % __file__
+    ERROR_MSG += 'Make sure compile_all.py has been run as described in '
+    ERROR_MSG += 'the ANUGA installation guide.'
+    raise Exception(ERROR_MSG)
 
 
Index: /trunk/anuga_core/source/anuga/geometry/polygon_function.py
===================================================================
--- /trunk/anuga_core/source/anuga/geometry/polygon_function.py	(revision 7858)
+++ /trunk/anuga_core/source/anuga/geometry/polygon_function.py	(revision 7858)
@@ -0,0 +1,133 @@
+"""
+    Callable function to determine if points lie inside or outside a polygon.
+    
+    As of June 2010 this module has a pylint quality rating of 8.85/10.
+"""
+
+import anuga.utilities.log as log
+import numpy as num
+from polygon import inside_polygon
+
+class Polygon_function:
+    """Create callable object f: x,y -> z, where a,y,z are vectors and
+    where f will return different values depending on whether x,y belongs
+    to specified polygons.
+
+    To instantiate:
+
+       Polygon_function(polygons)
+
+    where polygons is a list of tuples of the form
+
+      [ (P0, v0), (P1, v1), ...]
+
+      with Pi being lists of vertices defining polygons and vi either
+      constants or functions of x,y to be applied to points with the polygon.
+
+    The function takes an optional argument, default which is the value
+    (or function) to used for points not belonging to any polygon.
+    For example:
+
+       Polygon_function(polygons, default = 0.03)
+
+    If omitted the default value will be 0.0
+
+    Note: If two polygons overlap, the one last in the list takes precedence
+
+    Coordinates specified in the call are assumed to be relative to the
+    origin (georeference) e.g. used by domain.
+    By specifying the optional argument georeference,
+    all points are made relative.
+
+    FIXME: This should really work with geo_spatial point sets.
+    """
+
+    def __init__(self, regions, default=0.0, geo_reference=None):
+        """
+        # @brief Create instance of a polygon function.
+        # @param regions A list of (x,y) tuples defining a polygon.
+        # @param default Value or function returning value for points outside poly.
+        # @param geo_reference ??
+        """
+
+        try:
+            len(regions)
+        except:
+            msg = ('Polygon_function takes a list of pairs (polygon, value).'
+                   'Got %s' % str(regions))
+            raise Exception, msg
+
+        first_region = regions[0]
+
+        if isinstance(first_region, basestring):
+            msg = ('You passed in a list of text values into polygon_function '
+                   'instead of a list of pairs (polygon, value): "%s"'
+                   % str(first_region))
+            raise Exception, msg
+
+        try:
+            num_region_components = len(first_region)
+        except:
+            msg = ('Polygon_function takes a list of pairs (polygon, value). '
+                   'Got %s' % str(num_region_components))
+            raise Exception, msg
+
+        msg = ('Each entry in regions have two components: (polygon, value). '
+               'I got %s' % str(num_region_components))
+        assert num_region_components == 2, msg
+
+        if geo_reference is None:
+            from anuga.coordinate_transforms.geo_reference import Geo_reference
+            geo_reference = Geo_reference()
+
+        self.default = default
+
+        # Make points in polygons relative to geo_reference
+        self.regions = []
+        for polygon, value in regions:
+            georeffed_poly = geo_reference.change_points_geo_ref(polygon)
+            self.regions.append((georeffed_poly, value))
+
+    def __call__(self, pts_x, pts_y):
+        """
+        # @brief Implement the 'callable' property of Polygon_function.
+        # @param x List of x coordinates of points ot interest.
+        # @param y List of y coordinates of points ot interest.
+        """
+        pts_x = num.array(pts_x, num.float)
+        pts_y = num.array(pts_y, num.float)
+
+        # x and y must be one-dimensional and same length
+        assert len(pts_x.shape) == 1 and len(pts_y.shape) == 1
+        pts_len = pts_x.shape[0]
+        assert pts_y.shape[0] == pts_len, 'x and y must be same length'
+
+        points = num.ascontiguousarray(num.concatenate((pts_x[:, num.newaxis],
+                                                        pts_y[:, num.newaxis]),
+                                                       axis = 1 ))
+
+        if callable(self.default):
+            result = self.default(pts_x, pts_y)
+        else:
+            result = num.ones(pts_len, num.float) * self.default
+
+        for polygon, value in self.regions:
+            indices = inside_polygon(points, polygon)
+
+            # FIXME: This needs to be vectorised
+            if callable(value):
+                for i in indices:
+                    xx = num.array([pts_x[i]])
+                    yy = num.array([pts_y[i]])
+                    result[i] = value(xx, yy)[0]
+            else:
+                for i in indices:
+                    result[i] = value
+
+        if len(result) == 0:
+            msg = ('Warning: points provided to Polygon function did not fall '
+                   'within its regions in [%.2f, %.2f], y in [%.2f, %.2f]'
+                   % (min(pts_x), max(pts_x), min(pts_y), max(pts_y)))
+            log.critical(msg)
+
+        return result
Index: /trunk/anuga_core/source/anuga/geometry/quad.py
===================================================================
--- /trunk/anuga_core/source/anuga/geometry/quad.py	(revision 7857)
+++ /trunk/anuga_core/source/anuga/geometry/quad.py	(revision 7858)
@@ -8,4 +8,6 @@
 may then be iterated over with a proper intersection test to detect actual
 geometry intersections.
+
+As of June 2010 this module has a pylint quality rating of 10/10.
 
 """
Index: /trunk/anuga_core/source/anuga/geometry/test_polygon.py
===================================================================
--- /trunk/anuga_core/source/anuga/geometry/test_polygon.py	(revision 7857)
+++ /trunk/anuga_core/source/anuga/geometry/test_polygon.py	(revision 7858)
@@ -8,4 +8,6 @@
 
 from polygon import *
+from polygon import _poly_xy
+from polygon_function import Polygon_function
 from anuga.coordinate_transforms.geo_reference import Geo_reference
 from anuga.geospatial_data.geospatial_data import Geospatial_data
@@ -1570,5 +1572,5 @@
         # Simplest case: Polygon is the unit square
         polygon = [[0,0], [1,0], [1,1], [0,1]]
-        x, y = poly_xy(polygon)
+        x, y = _poly_xy(polygon)
         assert len(x) == len(polygon)+1
         assert len(y) == len(polygon)+1
@@ -1584,5 +1586,5 @@
     # Arbitrary polygon
         polygon = [[1,5], [1,1], [100,10], [1,10], [3,6]]
-        x, y = poly_xy(polygon)
+        x, y = _poly_xy(polygon)
         assert len(x) == len(polygon)+1
         assert len(y) == len(polygon)+1
@@ -1605,22 +1607,14 @@
         polygon1 = [[0,0], [1,0], [1,1], [0,1]]
         polygon2 = [[1,1], [2,1], [3,2], [2,2]]
-        v = plot_polygons([polygon1, polygon2], figname='test1')
-        assert len(v) == 4
-        assert v[0] == 0
-        assert v[1] == 3
-        assert v[2] == 0
-        assert v[3] == 2
+        plot_polygons([polygon1, polygon2], figname='test1')
 
         # Another case
         polygon3 = [[1,5], [10,1], [100,10], [50,10], [3,6]]
         v = plot_polygons([polygon2,polygon3], figname='test2')
-        assert len(v) == 4
-        assert v[0] == 1
-        assert v[1] == 100
-        assert v[2] == 1
-        assert v[3] == 10
-
-        os.remove('test1.png')
-        os.remove('test2.png')
+
+        for file in ['test1.png', 'test2.png']:
+            assert os.access(file, os.R_OK)
+            os.remove(file)
+
 
     def test_inside_polygon_geospatial(self):
Index: /trunk/anuga_core/source/anuga/pmesh/view_tsh.py
===================================================================
--- /trunk/anuga_core/source/anuga/pmesh/view_tsh.py	(revision 7857)
+++ /trunk/anuga_core/source/anuga/pmesh/view_tsh.py	(revision 7858)
@@ -10,5 +10,6 @@
 from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
 from anuga.pyvolution.util import file_function
-from anuga.geometry.polygon import Polygon_function, read_polygon
+from anuga.geometry.polygon import read_polygon
+from anuga.geometry.polygon_function import Polygon_function
 from realtime_visualisation_new import *
 
Index: /trunk/anuga_core/source/anuga/shallow_water/sww_interrogate.py
===================================================================
--- /trunk/anuga_core/source/anuga/shallow_water/sww_interrogate.py	(revision 7857)
+++ /trunk/anuga_core/source/anuga/shallow_water/sww_interrogate.py	(revision 7858)
@@ -4,4 +4,5 @@
 
 import os
+import anuga.utilities.log as log
 import numpy as num
 
