Index: trunk/anuga_core/source/anuga/utilities/quantity_setting_functions.py
===================================================================
--- trunk/anuga_core/source/anuga/utilities/quantity_setting_functions.py	(revision 9335)
+++ trunk/anuga_core/source/anuga/utilities/quantity_setting_functions.py	(revision 9338)
@@ -8,5 +8,9 @@
 import anuga.utilities.spatialInputUtil as su
 
-def make_nearestNeighbour_quantity_function(quantity_xyValueIn, domain, threshold_distance = 9.0e+100, background_value = 9.0e+100):
+def make_nearestNeighbour_quantity_function(
+        quantity_xyValueIn, 
+        domain, 
+        threshold_distance = 9.0e+100, 
+        background_value = 9.0e+100):
     """
     Function which makes another function, which can be used in set_quantity 
@@ -14,10 +18,13 @@
     Idea: For every point x,y in the domain, we want to set a quantity based on
           the 'nearest-neighbours' from quantity_xyValue (a 3 column array with
-          x,y,quantity-value), UNLESS the distance from x,y to the nearest-neighbour is > threshold_distance.
-          In the latter case, we want to set the quantity value to 'background_value'
+          x,y,quantity-value), 
+          UNLESS the distance from x,y to the nearest-neighbour is >
+            threshold_distance.
+          In the latter case, we want to set the quantity value to
+            'background_value'
             
           We need a function f(x,y) to do that. This routine makes the
-          function, with the desired quantity_xyValue points, threshold_distance, and
-          background_value
+          function, with the desired quantity_xyValue points, 
+          threshold_distance, and background_value
 
     INPUTS:
@@ -26,5 +33,6 @@
                             defining the points used to set the new quantity values
         domain -- The ANUGA domain
-        threshold_distance -- Points greater than this distance from their nearest quantity_xyValue point are set to background_value
+        threshold_distance -- Points greater than this distance from their 
+            nearest quantity_xyValue point are set to background_value
         background_value -- see 'threshold_distance'
 
@@ -42,5 +50,7 @@
         quantity_xyValue=copy.copy(quantity_xyValueIn.reshape((1,3)))
     # Make a function which gives us the ROW-INDEX of the nearest xy point in quantity_xyValue
-    quantity_xy_interpolator=scipy.interpolate.NearestNDInterpolator(quantity_xyValue[:,0:2], scipy.arange(len(quantity_xyValue[:,2])))
+    quantity_xy_interpolator=scipy.interpolate.NearestNDInterpolator(
+                                quantity_xyValue[:,0:2], 
+                                scipy.arange(len(quantity_xyValue[:,2])))
 
     #
@@ -67,5 +77,7 @@
         q_index=quantity_xy_interpolator(z)
         # Next find indices with distance < threshold_distance
-        dist_lt_thresh=( (z[:,0]-quantity_xyValue[q_index,0])**2 + (z[:,1]-quantity_xyValue[q_index,1])**2 < threshold_distance**2)
+        dist_lt_thresh=( (z[:,0]-quantity_xyValue[q_index,0])**2 + \
+                         (z[:,1]-quantity_xyValue[q_index,1])**2 < \
+                        threshold_distance**2)
         dist_lt_thresh=dist_lt_thresh.nonzero()[0]
         quantity_output[dist_lt_thresh] = quantity_xyValue[q_index[dist_lt_thresh],2] 
@@ -78,5 +90,6 @@
 def composite_quantity_setting_function(poly_fun_pairs, domain):
     """
-        Make a 'composite function' to set quantities -- applies different functions inside different polygon regions.
+        Make a 'composite function' to set quantities -- applies different 
+        functions inside different polygon regions.
              
         poly_fun_pairs = [ [p0, f0], [p1, f1], ...] 
@@ -84,39 +97,60 @@
                     Where:
 
-                      fi is a function, or a constant, or the name of a
-                        gdal-compatible rasterFile, or a numpy array with 3 columns; 
-
-                      pi is a polygon, or a polygon filename (shapefile or a csv format that anuga.read_polygon will read),
-                        or None, or 'All' (or it can be 'Extent' in the case that fi is a
-                        rasterFile name)
+                      fi is a function, 
+                         or a constant, 
+                         or a '.txt' or '.csv' file with comma separated xyz data 
+                            and a single header row, 
+                         or the name of a gdal-compatible rasterFile 
+                            (not ending in .txt or .csv), 
+                         or a numpy array with 3 columns; 
+
+                      pi is a polygon, 
+                        or a polygon filename (shapefile or a csv format that 
+                                                anuga.read_polygon will read),
+                        or None ( equivalent to a polygon with zero area),
+                        or 'All' (equivalent to a polygon covering everything)
+                        or 'Extent' in the case that fi is a rasterFile name 
+                            (equivalent to a polygon with the same extent as the raster)
               
-        IMPORTANT: When polygons overlap, the first elements of the list are given priority. 
+        IMPORTANT: When polygons overlap, the first elements of the list are 
+                   given priority. 
                    The approach is:
-                       First f0 is applied to all points in p0, and we record that these points have been 'set'
-                       Next f1 is applied to all points in p1 which have not been 'set', and then we record those points as being 'set'
-                       Next f2 is applied to all points in p2 which have not been 'set', and then we record those points as being 'set'
+                       First f0 is applied to all points in p0, and we record
+                         that these points have been 'set'
+                       Next f1 is applied to all points in p1 which have not 
+                         been 'set', and then we record those points as being 'set'
+                       Next f2 is applied to all points in p2 which have not 
+                         been 'set', and then we record those points as being 'set'
                        ... etc
 
         INPUT: 
-              poly_fun_pairs = [ [p0, f0], [p1, f1], ...]
-
-                  where fi(x,y) is a function returning quantity values at points, or any of the special cases below
-                  SPECIAL fi CASES:
-                  fi = a constant in which case points in the polygon are set to that value, 
-                  fi = a string rasterFile name which can be passed to quantityRasterFun to make a function,
-                  fi = a numpy array with 3 columns (x,y,Value) in which case nearest-neighbour interpolation is used on the points
-                
-
-                  pi are polygons where we want to use fi inside
-                  SPECIAL pi CASES: 
-                  If pi is a filename ending in .shp or a csv format that anuga.read_polygon can read, we assume it contains a polygon we have to read
-                  If any pi = 'All', then we assume that ALL unset points are set
-                     using the function. This CAN ONLY happen in the last [fi,pi] pair where pi is
-                     not None (since fi will be applied to all remaining points -- so anything else is probably an input mistake)
-                  If any pi = None, then that [fi,pi] pair is skipped
-                  If pi = 'Extent' and fi is the name of a raster file, then the
-                    extent of the raster file is used to define the polygon
-
-              domain = ANUGA domain object
+          poly_fun_pairs = [ [p0, f0], [p1, f1], ...]
+
+              where fi(x,y) is a function returning quantity values at points, 
+                or any of the special cases below
+              SPECIAL fi CASES:
+              fi = a constant in which case points in the polygon are 
+                   set to that value, 
+              fi = a .txt or .csv file name containing x, y, z data,
+                     with comma separators and a single header row
+              fi = a string rasterFile name (not ending in .txt or .csv)
+                    which can be passed to quantityRasterFun to make a function,
+              fi = a numpy array with 3 columns (x,y,Value) in which case 
+                   nearest-neighbour interpolation is used on the points
+            
+
+              pi are polygons where we want to use fi inside
+              SPECIAL pi CASES: 
+              If pi is a filename ending in .shp or a csv format that 
+                anuga.read_polygon can read, we assume it contains a polygon we have to read
+              If any pi = 'All', then we assume that ALL unset points are set
+                 using the function. This CAN ONLY happen in the last [fi,pi] pair where pi is
+                 not None (since fi will be applied to all remaining points 
+                 -- so anything else is probably an input mistake)
+              If any pi = None, then that [fi,pi] pair is skipped
+              If pi = 'Extent' and fi is the name of a raster file, then the
+                extent of the raster file is used to define the polygon
+
+          domain = ANUGA domain object
 
 
@@ -149,5 +183,6 @@
             if (poly_fun_pairs[i][0]=='All'):
                 # This is only ok if all the othe poly_fun_pairs are None
-                remaining_poly_fun_pairs_are_None=[ poly_fun_pairs[j][0] is None for j in range(i+1,lfp)]
+                remaining_poly_fun_pairs_are_None = \
+                    [ poly_fun_pairs[j][0] is None for j in range(i+1,lfp)]
                 if(not all(remaining_poly_fun_pairs_are_None)):
                     raise Exception, 'Can only have the last polygon = All'
@@ -191,7 +226,17 @@
                 elif ( type(fi) is str and os.path.exists(fi)):
                     # fi is a file which is assumed to be 
-                    # a gdal-compatible raster
-                    newfi = quantityRasterFun(domain, fi)
-                    quantityVal[fInds] = newfi(x[fInds], y[fInds])
+                    # a gdal-compatible raster OR an x,y,z elevation file
+                    if os.path.splitext(fi)[1] in ['.txt', '.csv']:
+                        # Treating input file ' + fi + ' as xyz array with 1 header row
+                        fi_array = numpy.genfromtxt(fi, delimiter=",", skip_header=1)
+                        if fi_array.shape[1] is not 3:
+                            print 'Treated input file ' + fi + ' as xyz array with 1 header row'
+                            raise Exception, 'Array should have 3 columns -- x,y,value'
+                        newfi = make_nearestNeighbour_quantity_function(fi_array, domain)
+                        quantityVal[fInds] = newfi(x[fInds], y[fInds])
+                    else:
+                        # Treating input file as a raster
+                        newfi = quantityRasterFun(domain, fi)
+                        quantityVal[fInds] = newfi(x[fInds], y[fInds])
                 elif(type(fi) is numpy.ndarray):
                     if fi.shape[1] is not 3:
Index: trunk/anuga_core/source/anuga/utilities/spatialInputUtil.py
===================================================================
--- trunk/anuga_core/source/anuga/utilities/spatialInputUtil.py	(revision 9335)
+++ trunk/anuga_core/source/anuga/utilities/spatialInputUtil.py	(revision 9338)
@@ -251,5 +251,6 @@
             except:
                 msg = 'Could not read points from ' + filename +\
-                      '. Make sure it has a single header row, with comma separator, and the first 2 columns are x,y'
+                      '. Make sure it has a single header row, ' +\
+                      'with comma separator, and the first 2 columns are x,y'
                 raise Exception, msg
         return points
@@ -408,5 +409,6 @@
                              as joined and so defining line segments
             
-            OUTPUT: The squared distance, and the index i of the segment [x_i,y_i],[x_i+1,y_i+1] in segments which is closest to pt
+            OUTPUT: The squared distance, and the index i of the segment
+                    [x_i,y_i],[x_i+1,y_i+1] in segments which is closest to pt
         """
         ll=len(segments)
@@ -425,5 +427,6 @@
     def shift_point_on_line(pt, lineIn, nearest_segment_index):
         """
-            Support pt is a point, which is near to the 'nearest_segment_index' segment of the line 'lineIn'
+            Support pt is a point, which is near to the 'nearest_segment_index'
+                segment of the line 'lineIn'
     
             This routine moves the nearest end point of that segment on line to pt.
@@ -431,5 +434,6 @@
             INPUTS: pt -- [x,y] point
                     lineIn -- [ [x0, y0], [x1, y1], ..., [xN,yN]] 
-                    nearest_segment_index = index where the distance of pt to the line from [x_i,y_i] to [x_i+1,y_i+1] is minimum
+                    nearest_segment_index = index where the distance of pt to
+                        the line from [x_i,y_i] to [x_i+1,y_i+1] is minimum
     
             OUTPUT: The new line
@@ -473,5 +477,6 @@
         p0 = L1_pts[tmp[1]]
         p1 = L1_pts[tmp[1]+1]
-        endPt_Dist_Sq=min( ( (p0[0]-iP[0])**2 + (p0[1]-iP[1])**2), ( (p1[0]-iP[0])**2 + (p1[1]-iP[1])**2))
+        endPt_Dist_Sq=min( ( (p0[0]-iP[0])**2 + (p0[1]-iP[1])**2), 
+                           ( (p1[0]-iP[0])**2 + (p1[1]-iP[1])**2))
         #
         if(endPt_Dist_Sq>point_movement_threshold**2):
@@ -560,5 +565,6 @@
                      tol2 = [see check_polygon_is_small] Probably doesn't need to change
     
-                     nameFlag = This will be printed if intersection occurs. Convenient way to display intersecting filenames
+                     nameFlag = This will be printed if intersection occurs.
+                            Convenient way to display intersecting filenames
             
             OUTPUTS: L1,L2 with intersection points added in the right places
@@ -567,6 +573,6 @@
         if(L1.Intersects(L2)):
             # Get points on the lines 
-            L1_pts=Wkb2ListPts(L1) #L1.GetPoints()
-            L2_pts=Wkb2ListPts(L2) #L2.GetPoints()
+            L1_pts=Wkb2ListPts(L1) 
+            L2_pts=Wkb2ListPts(L2) 
             
             # Buffer lines by a small amount
@@ -578,6 +584,7 @@
             if(L1_L2_intersect.GetGeometryCount()==1):
                 if(not check_polygon_is_small(L1_L2_intersect, buf, tol2)):
-                    #print L1_L2_intersect.GetEnvelope()
-                    raise Exception, 'line intersection is not allowed. Envelope %s '% str(L1_L2_intersect.GetEnvelope())
+                    msg = 'line intersection is not allowed. ' + \
+                          'Envelope %s '% str(L1_L2_intersect.GetEnvelope())
+                    raise Exception(msg)
                 # Seems to need special treatment with only 1 intersection point
                 intersectionPts=[L1_L2_intersect.Centroid().GetPoint()]
@@ -597,6 +604,8 @@
             # Insert the points into the line segments
             for i in range(len(intersectionPts)):
-                L1_pts = insert_intersection_point(intersectionPts[i], L1_pts, point_movement_threshold, verbose=verbose)
-                L2_pts = insert_intersection_point(intersectionPts[i], L2_pts, point_movement_threshold, verbose=verbose)
+                L1_pts = insert_intersection_point(intersectionPts[i], L1_pts, 
+                            point_movement_threshold, verbose=verbose)
+                L2_pts = insert_intersection_point(intersectionPts[i], L2_pts, 
+                            point_movement_threshold, verbose=verbose)
                 
             # Convert to the input format
@@ -716,5 +725,6 @@
                get a range of raster values inside a polygon
     
-            Approach: A 'trial-grid' of points is created which is 'almost' covering the range of the polygon (xmin-xmax,ymin-ymax). 
+            Approach: A 'trial-grid' of points is created which is 'almost'
+                      covering the range of the polygon (xmin-xmax,ymin-ymax). 
     
                       (Actually it is just inside this region, to avoid polygon-boundary issues, see below)
@@ -763,14 +773,5 @@
         keepers = inside_polygon(Grid, polygon)
         if(len(keepers)==0):
-            import pdb
-            pdb.set_trace()
-        ## Now find the xy values which are actually inside the polygon
-        #polpath=path.Path(polygonArr)
-        #keepers=[]
-        #for i in range(len(Grid[:,0])):
-        #   if(polpath.contains_point(Grid[i,:])):
-        #        keepers=keepers+[i]
-        ##import pdb
-        ##pdb.set_trace()
+            raise Exception('No points extracted from polygon')
         xyInside=Grid[keepers,:]
     
@@ -958,5 +959,6 @@
                         msg='Breakline/riverwall point '+str(blCat[n1][j][0:2]) +' on '+ n1+\
                             ' is outside the bounding polygon.\n'+\
-                            'Check that it exceeds the bounding polygon by a distance < point_movement_threshold \n'+\
+                            'Check that it exceeds the bounding polygon'+\
+                            ' by a distance < point_movement_threshold \n'+\
                             ' so it can be moved back onto the polygon'
                         print 'Polygon\n '
@@ -1013,5 +1015,6 @@
             Read them in
             
-            INPUT: fileList -- a list of shapefile and/or anuga_polygon csv filenames [e.g. from glob.glob('GIS/Breaklines/*.shp')]
+            INPUT: fileList -- a list of shapefile and/or anuga_polygon csv filenames 
+                               [e.g. from glob.glob('GIS/Breaklines/*.shp')]
     
             OUTPUT: dictionary with breaklines [filenames are keys]
@@ -1027,8 +1030,9 @@
     def readListOfRiverWalls(rwfileList):
         """
-            Take a list with the names of riverwall input files [should be comma-separated x,y,elevation files]
-
-            The input file can OPTIONALLY have a first line defining some hydraulic parameters. A valid example 
-            is
+            Take a list with the names of riverwall input files 
+            [should be comma-separated x,y,elevation files]
+
+            The input file can OPTIONALLY have a first line defining some
+            hydraulic parameters. A valid example is
 
             Qfactor: 1.5, s1: 0.94 
@@ -1039,5 +1043,6 @@
             Read their coordinates into a dict with their names, read for use by ANUGA
     
-            INPUT: rwfileList -- a list of riverwall filenames [e.g. from glob.glob('GIS/RiverWalls/*.csv')]
+            INPUT: rwfileList -- a list of riverwall filenames 
+                        [e.g. from glob.glob('GIS/RiverWalls/*.csv')]
     
             OUTPUT: 
@@ -1074,9 +1079,11 @@
     def combine_breakLines_and_riverWalls_for_mesh(breakLines, riverWalls):
         """
-        Combine breaklines and riverwalls for input in mesh generation, ensuring both have 2 coordinates only
+        Combine breaklines and riverwalls for input in mesh generation,
+            ensuring both have 2 coordinates only
         """
         mesh_breakLines=riverWalls.values()+breakLines.values()
         for i in range(len(mesh_breakLines)):
-            mesh_breakLines[i] = [mesh_breakLines[i][j][0:2] for j in range(len(mesh_breakLines[i]))]
+            mesh_breakLines[i] =\
+             [mesh_breakLines[i][j][0:2] for j in range(len(mesh_breakLines[i]))]
         return mesh_breakLines
     
@@ -1089,12 +1096,16 @@
             Can do this with the current function
     
-            INPUTS: pattern == character string containing pattern which is inside exactly 2 keys in breaklines
+            INPUTS: pattern == character string containing pattern which 
+                        is inside exactly 2 keys in breaklines
         
                     breakLines = the breakLinesIn dictionary
                     
-                    reverse2nd = True/False or None. Reverse the order of the 2nd set of edges before making the polygon.
+                    reverse2nd = True/False or None. Reverse the order of the 
+                                 2nd set of edges before making the polygon.
                                  If None, then we compute the distance between the
-                                 first point on breakline1 and the first/last points on breakline2, and reverse2nd if
-                                 the 'distance from the first point' < 'distance from the last point'
+                                  first point on breakline1 and the first/last
+                                  points on breakline2, and reverse2nd if
+                                 the 'distance from the first point' <
+                                      'distance from the last point'
         
             OUTPUT: Polygon made with the 2 breaklines
@@ -1129,7 +1140,7 @@
             breakLines[bk[matchers[1]]].reverse()
         polyOut=breakLines[bk[matchers[0]]] +  breakLines[bk[matchers[1]]]
-        #polyOut=polyOut+[polyOut[0]]
-
-        # If the breakLines values have > 2 entries (i.e. like for riverwalls), remove the third
+
+        # If the breakLines values have > 2 entries (i.e. like for riverwalls),
+        # remove the third
         if(len(polyOut[0])>2):
             polyOut=[polyOut[i][0:2] for i in range(len(polyOut))]
