Index: /trunk/anuga_core/source/anuga/utilities/plot_utils.py
===================================================================
--- /trunk/anuga_core/source/anuga/utilities/plot_utils.py	(revision 9510)
+++ /trunk/anuga_core/source/anuga/utilities/plot_utils.py	(revision 9511)
@@ -736,6 +736,5 @@
     return near_points[keepers], local_coord[keepers]
 
-########################
-# TRIANGLE AREAS, WATER VOLUME
+
 def triangle_areas(p, subset=None):
     # Compute areas of triangles in p -- assumes p contains vertex information
@@ -764,5 +763,4 @@
     return area
 
-###
 
 def water_volume(p, p2, per_unit_area=False, subset=None):
@@ -779,6 +777,6 @@
    
     # This accounts for how volume is measured in ANUGA 
-    # Compute in 2 steps to reduce precision error (important sometimes)
-    # Is this really needed?
+    # Compute in 2 steps to reduce precision error from limited SWW precision
+    # FIXME: Is this really needed?
     for i in range(l):
         #volume[i]=((p2.stage[i,subset]-p2.elev[subset])*(p2.stage[i,subset]>p2.elev[subset])*area).sum()
@@ -1050,5 +1048,6 @@
 
     # Get function to interpolate quantity onto gridXY_array
-    gridXY_array=scipy.array([scipy.concatenate(gridX),scipy.concatenate(gridY)]).transpose()
+    gridXY_array=scipy.array([scipy.concatenate(gridX),
+                              scipy.concatenate(gridY)]).transpose()
     gridXY_array=scipy.ascontiguousarray(gridXY_array)
 
@@ -1056,6 +1055,8 @@
     #basic_nearest_neighbour=False
     if(k_nearest_neighbours==1):
-        index_qFun=scipy.interpolate.NearestNDInterpolator(swwXY,scipy.arange(len(swwX),dtype='int64').transpose())
-        gridqInd=index_qFun(gridXY_array)
+        index_qFun = scipy.interpolate.NearestNDInterpolator(
+            swwXY,
+            scipy.arange(len(swwX),dtype='int64').transpose())
+        gridqInd = index_qFun(gridXY_array)
         # Function to do the interpolation
         def myInterpFun(quantity):
@@ -1063,18 +1064,16 @@
     else:
         # Combined nearest neighbours and inverse-distance interpolation
-        index_qFun=scipy.spatial.cKDTree(swwXY)
-        NNInfo=index_qFun.query(gridXY_array,k=k_nearest_neighbours)
+        index_qFun = scipy.spatial.cKDTree(swwXY)
+        NNInfo = index_qFun.query(gridXY_array, k=k_nearest_neighbours)
         # Weights for interpolation
-        nn_wts=1./(NNInfo[0]+1.0e-100)
-        nn_inds=NNInfo[1]
+        nn_wts = 1./(NNInfo[0]+1.0e-100)
+        nn_inds = NNInfo[1]
         def myInterpFun(quantity):
-            denom=0.
-            num=0.
+            denom = 0.
+            num = 0.
             for i in range(k_nearest_neighbours):
-                denom+=nn_wts[:,i]
-                num+= quantity[nn_inds[:,i]]*nn_wts[:,i]
+                denom += nn_wts[:,i]
+                num += quantity[nn_inds[:,i]]*nn_wts[:,i]
             return (num/denom)
-
-
 
     if(bounding_polygon is not None):
Index: /trunk/anuga_core/source/anuga/utilities/quantity_setting_functions.py
===================================================================
--- /trunk/anuga_core/source/anuga/utilities/quantity_setting_functions.py	(revision 9510)
+++ /trunk/anuga_core/source/anuga/utilities/quantity_setting_functions.py	(revision 9511)
@@ -12,5 +12,7 @@
         domain, 
         threshold_distance = 9.0e+100, 
-        background_value = 9.0e+100):
+        background_value = 9.0e+100,
+        k_nearest_neighbours = 1,
+    ):
     """
     Function which makes another function, which can be used in set_quantity 
@@ -32,4 +34,6 @@
             defining the points used to set the new quantity values
         @param domain -- The ANUGA domain
+        @param k_nearest_neighbors If > 1, then an inverse-distance-weighted average
+               of the k nearest neighbours is used
         @param threshold_distance -- Points greater than this distance from 
             their nearest quantity_xyValue point are set to background_value
@@ -38,50 +42,85 @@
     OUTPUTS: 
         A function f which can be passed to domain.set_quantity('myQuantity', f)
-    """
+
+    """
+
     import scipy
     import scipy.interpolate
-
-
-    if(len(quantity_xyValueIn.shape)>1):
-        quantity_xyValue=copy.copy(quantity_xyValueIn) # Pointer, no copy
+    import scipy.spatial
+
+    if(len(quantity_xyValueIn.shape) > 1):
+        quantity_xyValue = quantity_xyValueIn
     else:
         # Treat the single-point case
-        quantity_xyValue=copy.copy(quantity_xyValueIn.reshape((1,3)))
+        quantity_xyValue = 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])))
+
+    # Make a function which returns k-nearest-neighbour indices + distances
+    quantity_xy_interpolator = scipy.spatial.cKDTree(quantity_xyValue[:,0:2])
+
     # Make a function of x,y which we can pass to domain.set_quantity
     def quant_NN_fun(x,y):
         """
         Function to assign quantity from the nearest point in quantity_xyValue,
-        UNLESS the point is more than 'threshold_distance' away from the nearest point,
-        in which case the background friction value is used
-       
+        UNLESS the point is more than 'threshold_distance' away from the
+        nearest point, in which case the background value is used
+
         """
+
         import scipy
         import scipy.interpolate
+        import scipy.spatial
+
         # Since ANUGA stores x,y internally in non-georeferenced coordinates,
         # we adjust them here
-        xll=domain.geo_reference.xllcorner
-        yll=domain.geo_reference.yllcorner
-        z=scipy.zeros(shape=(len(x), 2))
-        z[:,0]=x+xll
-        z[:,1]=y+yll
+        xll = domain.geo_reference.xllcorner
+        yll = domain.geo_reference.yllcorner
+        z = scipy.zeros(shape=(len(x), 2))
+        z[:,0] = x+xll
+        z[:,1] = y+yll
+
         # This will hold the quantity values
-        quantity_output=x*0. +background_value
+        quantity_output = x*0. + background_value
         # Compute the index of the nearest-neighbour in quantity_xyValue
-        q_index=quantity_xy_interpolator(z)
+        neighbour_data = quantity_xy_interpolator.query(z, 
+            k=k_nearest_neighbours)
+
         # 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=dist_lt_thresh.nonzero()[0]
-        quantity_output[dist_lt_thresh] =\
-            quantity_xyValue[q_index[dist_lt_thresh],2]
+        if(k_nearest_neighbours==1):
+            dist_lt_thresh = neighbour_data[0] < threshold_distance
+        else:
+            dist_lt_thresh = neighbour_data[0][:,0] < threshold_distance
+
+        dist_lt_thresh = dist_lt_thresh.nonzero()[0]
+
+        # Initialise output
+        quantity_output = x*0 + background_value
+
+        # Interpolate
+        if len(dist_lt_thresh)>0:
+            numerator = 0
+            denominator = 0
+            for i in range(k_nearest_neighbours):
+                if(k_nearest_neighbours==1):
+                    distances = neighbour_data[0][dist_lt_thresh]
+                    indices = neighbour_data[1][dist_lt_thresh]
+                else:
+                    distances = neighbour_data[0][dist_lt_thresh,i]
+                    indices = neighbour_data[1][dist_lt_thresh,i]
+                
+                inverse_distance = 1.0/(distances+1.0e-100)
+                values = quantity_xyValue[indices,2]
+                numerator += values*inverse_distance
+                denominator += inverse_distance
+           
+            quantity_output[dist_lt_thresh] = numerator/denominator
+
         return quantity_output
+
     # Return the quantity function
     return quant_NN_fun
Index: /trunk/anuga_core/source/anuga/utilities/test/test_quantity_setting_functions.py
===================================================================
--- /trunk/anuga_core/source/anuga/utilities/test/test_quantity_setting_functions.py	(revision 9510)
+++ /trunk/anuga_core/source/anuga/utilities/test/test_quantity_setting_functions.py	(revision 9511)
@@ -43,15 +43,17 @@
         """
 
-        boundaryPolygon=[ [minX, minY], [minX, minY+100.], [minX+100., minY+100.], [minX+100., minY]]
-        anuga.create_mesh_from_regions(boundaryPolygon, 
-                                   boundary_tags={'left': [0],
-                                                'top': [1],
-                                                'right': [2],
-                                                'bottom': [3]},
-                                   maximum_triangle_area = 1.,
-                                   minimum_triangle_angle = 28.0,
-                                   filename = 'test_quantity_setting_functions.msh',
-                                   interior_regions =[ ],
-                                   verbose=False)
+        boundaryPolygon=[ [minX, minY], [minX, minY+100.], 
+            [minX+100., minY+100.], [minX+100., minY]]
+        anuga.create_mesh_from_regions(
+            boundaryPolygon, 
+            boundary_tags={'left': [0],
+                         'top': [1],
+                         'right': [2],
+                         'bottom': [3]},
+            maximum_triangle_area = 1.,
+            minimum_triangle_angle = 28.0,
+            filename = 'test_quantity_setting_functions.msh',
+            interior_regions =[ ],
+            verbose=False)
 
         domain=anuga.create_domain_from_file('test_quantity_setting_functions.msh')
@@ -103,6 +105,6 @@
         inPts=numpy.vstack([xR,yR,zR]).transpose()
 
-        F=qs.make_nearestNeighbour_quantity_function(inPts, domain, 
-                            threshold_distance = 9.0e+100, background_value = 9.0e+100)
+        F = qs.make_nearestNeighbour_quantity_function(inPts, domain, 
+                threshold_distance = 9.0e+100, background_value = 9.0e+100)
 
         # Test that F evaluated in 'ANUGA coordinates' [lower-left = 0,0] is correct
@@ -125,6 +127,6 @@
         inPts=numpy.vstack([xR,yR,zR]).transpose()
 
-        F=qs.make_nearestNeighbour_quantity_function(inPts, domain, 
-                            threshold_distance = 6., background_value = 9.0e+100)
+        F = qs.make_nearestNeighbour_quantity_function(inPts, domain, 
+                threshold_distance = 6., background_value = 9.0e+100)
 
         # Test that F evaluated in 'ANUGA coordinates' [lower-left = 0,0] is correct
@@ -158,11 +160,13 @@
 
         # Make a polygon-point pair which we use to set elevation in a 'channel'
-        trenchPoly=[[minX+40., minY], [minX+40., minY+100.], [minX+60., minY+100.], [minX+60., minY]]
+        trenchPoly = [[minX+40., minY], [minX+40., minY+100.], 
+            [minX+60., minY+100.], [minX+60., minY]]
 
         #################################################################
  
         # This example uses a constant, and a raster, to set the quantity           
-        F=qs.composite_quantity_setting_function([[trenchPoly, -1000.], ['Extent', 'PointData_ElevTest.tif']],\
-                                                domain) 
+        F=qs.composite_quantity_setting_function(
+            [[trenchPoly, -1000.], ['Extent', 'PointData_ElevTest.tif']],
+            domain)
 
         # Points where we test the function
@@ -176,5 +180,6 @@
         # Find the nearest domain point to the second test point
         # This will have been used in constructing the elevation raster
-        nearest=((domain.centroid_coordinates[:,0]-3.)**2 + (domain.centroid_coordinates[:,1]-20.)**2).argmin()
+        nearest=((domain.centroid_coordinates[:,0]-3.)**2 + 
+                 (domain.centroid_coordinates[:,1]-20.)**2).argmin()
         nearest_x=domain.centroid_coordinates[nearest,0]
         assert(numpy.allclose(fitted[1],-nearest_x/150.))
@@ -182,8 +187,10 @@
         #################################################################
  
-        # This example uses a constant, and a raster, to set the quantity, and applies the min/max bound           
-        F=qs.composite_quantity_setting_function([[trenchPoly, -1000.], ['Extent', 'PointData_ElevTest.tif']],\
-                                                domain,
-                                                clip_range = [[-500., 1.0e+100], [-1.0e+100, 1.0e+100]]) 
+        # This example uses a constant, and a raster, to set the quantity, and
+        # applies the min/max bound           
+        F=qs.composite_quantity_setting_function(
+            [[trenchPoly, -1000.], ['Extent', 'PointData_ElevTest.tif']],
+            domain,
+            clip_range = [[-500., 1.0e+100], [-1.0e+100, 1.0e+100]]) 
 
         # Points where we test the function
@@ -197,6 +204,8 @@
         # Find the nearest domain point to the second test point
         # This will have been used in constructing the elevation raster
-        nearest=((domain.centroid_coordinates[:,0]-3.)**2 + (domain.centroid_coordinates[:,1]-20.)**2).argmin()
-        nearest_x=domain.centroid_coordinates[nearest,0]
+        nearest = ((domain.centroid_coordinates[:,0]-3.)**2 + 
+                   (domain.centroid_coordinates[:,1]-20.)**2).argmin()
+        nearest_x = domain.centroid_coordinates[nearest,0]
+
         assert(numpy.allclose(fitted[1],-nearest_x/150.))
 
@@ -206,12 +215,14 @@
         def f0(x,y):
             return x/10.
-        F=qs.composite_quantity_setting_function([[trenchPoly, f0], ['Extent', 'PointData_ElevTest.tif']],\
-                                                domain) 
-        fitted=F(testPts_X,testPts_Y)
+        F = qs.composite_quantity_setting_function(
+            [[trenchPoly, f0], ['Extent', 'PointData_ElevTest.tif']],
+            domain) 
+        fitted = F(testPts_X,testPts_Y)
         # Now the fitted value in the trench should be determined by f0
         assert(numpy.allclose(fitted[0],50./10.))
         # The second test point should be as before
-        nearest=((domain.centroid_coordinates[:,0]-3.)**2 + (domain.centroid_coordinates[:,1]-20.)**2).argmin()
-        nearest_x=domain.centroid_coordinates[nearest,0]
+        nearest = ((domain.centroid_coordinates[:,0]-3.)**2 + 
+                   (domain.centroid_coordinates[:,1]-20.)**2).argmin()
+        nearest_x = domain.centroid_coordinates[nearest,0]
         assert(numpy.allclose(fitted[1],-nearest_x/150.))
 
@@ -219,6 +230,7 @@
 
         # This example uses 'All' as a polygon
-        F=qs.composite_quantity_setting_function([['All', f0 ], [None, 'PointData_ElevTest.tif']],\
-                                                domain) 
+        F = qs.composite_quantity_setting_function(
+            [['All', f0 ], [None, 'PointData_ElevTest.tif']],
+            domain) 
         fitted=F(testPts_X,testPts_Y)
         # Now the fitted value in the trench should be determined by f0
@@ -229,6 +241,7 @@
         # This example should fail
         def should_fail():
-            F=qs.composite_quantity_setting_function([['All', f0], ['All', 'PointData_ElevTest.tif']],\
-                                                    domain) 
+            F=qs.composite_quantity_setting_function(
+                [['All', f0], ['All', 'PointData_ElevTest.tif']],
+                domain)
             # Need to call it to get the error
             F(numpy.array([3.]), numpy.array([3.]))
@@ -239,7 +252,8 @@
         # This example should fail (since the clip_range minimum >= maximum)
         def should_fail_1():
-            F=qs.composite_quantity_setting_function([[trenchPoly, f0], ['All', 'PointData_ElevTest.tif']],\
-                                                    domain,
-                                                    clip_range = [[ -500., -1000.], [-1.0e+100, 1.0e+100]]) 
+            F=qs.composite_quantity_setting_function(
+                [[trenchPoly, f0], ['All', 'PointData_ElevTest.tif']],
+                domain,
+                clip_range = [[ -500., -1000.], [-1.0e+100, 1.0e+100]]) 
             return
         self.assertRaises(Exception, lambda: should_fail_1())
@@ -251,5 +265,5 @@
         # 
         #
-        domain=self.create_domain(1.0, 0.0)
+        domain = self.create_domain(1.0, 0.0)
 
         # Evolve the model
@@ -259,30 +273,36 @@
         # Make a raster from the elevation data
         from anuga.utilities import plot_utils as util
-        xs=domain.centroid_coordinates[:,0]+domain.geo_reference.xllcorner
-        ys=domain.centroid_coordinates[:,1]+domain.geo_reference.yllcorner
-        elev=domain.quantities['elevation'].centroid_values
-
-        allDat=numpy.vstack([xs,ys,elev]).transpose()
-        util.Make_Geotif(allDat, output_quantities=['ElevTest'], EPSG_CODE=32756, 
-                        output_dir='.', CellSize=1.,k_nearest_neighbours=1)
+        xs = domain.centroid_coordinates[:,0]+domain.geo_reference.xllcorner
+        ys = domain.centroid_coordinates[:,1]+domain.geo_reference.yllcorner
+        elev = domain.quantities['elevation'].centroid_values
+
+        allDat = numpy.vstack([xs,ys,elev]).transpose()
+        util.Make_Geotif(allDat, output_quantities=['ElevTest'], 
+            EPSG_CODE=32756, output_dir='.', CellSize=1.,k_nearest_neighbours=1)
 
         # Make a polygon-point pair which we use to set elevation in a 'channel'
-        trenchPoly=[[minX+40., minY], [minX+40., minY+100.], [minX+60., minY+100.], [minX+60., minY]]
-        trenchPts=numpy.array([minX+50., minY+50., -1000.])
+        trenchPoly = [[minX+40., minY], [minX+40., minY+100.], 
+            [minX+60., minY+100.], [minX+60., minY]]
+        trenchPts = numpy.array([minX+50., minY+50., -1000.])
         #
-        PtPolData=[[trenchPoly, trenchPts]]
-        F=qs.quantity_from_Pt_Pol_Data_and_Raster(PtPolData, 'PointData_ElevTest.tif', domain) 
-
-        testPts_X=numpy.array([50., 3.])
-        testPts_Y=numpy.array([1., 20.])
-        fitted=F(testPts_X,testPts_Y)
+        PtPolData = [[trenchPoly, trenchPts]]
+        F = qs.quantity_from_Pt_Pol_Data_and_Raster(PtPolData, 
+            'PointData_ElevTest.tif', domain) 
+
+        testPts_X = numpy.array([50., 3.])
+        testPts_Y = numpy.array([1., 20.])
+        fitted = F(testPts_X,testPts_Y)
+
+        import pdb
+        pdb.set_trace()
     
         # The fitted value in the trench should be -1000.
-        assert(fitted[0]==-1000.)
+        assert(numpy.allclose(fitted[0],-1000.))
 
         # Find the nearest domain point to the second test point
         # This will have been used in constructing the elevation raster
-        nearest=((domain.centroid_coordinates[:,0]-3.)**2 + (domain.centroid_coordinates[:,1]-20.)**2).argmin()
-        nearest_x=domain.centroid_coordinates[nearest,0]
+        nearest = ((domain.centroid_coordinates[:,0]-3.)**2 +\
+                   (domain.centroid_coordinates[:,1]-20.)**2).argmin()
+        nearest_x = domain.centroid_coordinates[nearest,0]
         assert(numpy.allclose(fitted[1],-nearest_x/150.))
 
