Changeset 9303

Timestamp:

Aug 14, 2014, 10:33:25 AM (11 years ago)

Author:

davies

Message:

Fixing up recent plot_utils improvements

Location:

trunk/anuga_core/source/anuga/utilities

Files:

• plot_utils.py (modified) (1 diff)
• test_plot_utils.py (modified) (1 diff)

trunk/anuga_core/source/anuga/utilities/plot_utils.py

@@ -975 +975 @@
         # Combined nearest neighbours and inverse-distance interpolation
         index_qFun=scipy.spatial.cKDTree(swwXY)
-        NNInfo=index_qFun.query(gridXY_array,k=3)
+        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]
         def myInterpFun(quantity):
-            denom=nn_wts[:,0]+nn_wts[:,1]+nn_wts[:,2]
-            num = quantity[nn_inds[:,0]]*nn_wts[:,0]+\
-                  quantity[nn_inds[:,1]]*nn_wts[:,1]+\
-                  quantity[nn_inds[:,2]]*nn_wts[:,2]
+            denom=0.
+            num=0.
+            for i in range(k_nearest_neighbours):
+                denom+=nn_wts[:,i]
+                num+= quantity[nn_inds[:,i]]*nn_wts[:,i]
             return (num/denom)
 

trunk/anuga_core/source/anuga/utilities/test_plot_utils.py

@@ -397 +397 @@
         # Delete tif made with Make_Geotif
         os.remove('PointData_TestData.tif')
+
+    def test_Make_Geotif_with_knn(self):
+        # VERY BASIC TEST using knn+inverse distance interpolation to make the grid
+        #
+        # Simply create some data and grid it
+        #
+        # If nothing fails, that's good
+        #
+        # Pick a domain that makes sense in EPSG:32756
+        x=np.linspace(307000., 308000., 100)     
+        y=np.linspace(6193000., 6194000., 100)
+        myCellSize=5.0
+        xG,yG=np.meshgrid(x, y)
+        xG=xG.flatten()
+        yG=yG.flatten()
+        # Surface is z=x+y
+        fakeZ=xG-min(xG)+yG -min(yG)
+        dataToGrid=np.vstack([xG,yG,fakeZ]).transpose()
+        util.Make_Geotif(dataToGrid, output_quantities=['TestData'],
+                         EPSG_CODE=32756, output_dir='.', CellSize=myCellSize,
+                         k_nearest_neighbours=4)
+
+       
+        # Use gdal to check that at least the data extent is ok
+        import gdal
+        raster=gdal.Open('PointData_TestData.tif')
+        rasterGeoTrans=raster.GetGeoTransform()
+        assert(np.allclose(x.min()-myCellSize/2.0, rasterGeoTrans[0]))
+        assert(np.allclose(y.max()+myCellSize/2.0, rasterGeoTrans[3]))
+        #release data file
+        raster = None
+        # Delete tif made with Make_Geotif
+        os.remove('PointData_TestData.tif')