Changeset 9299

Timestamp:

Aug 13, 2014, 5:11:19 PM (11 years ago)

Author:

davies

Message:

Smoother interpolation for makeGeotif + speedups in quantity_setting_functions and spatialInputUtils

Location:

trunk/anuga_core/source/anuga/utilities

Files:

• plot_utils.py (modified) (6 diffs)
• quantity_setting_functions.py (modified) (3 diffs)
• spatialInputUtil.py (modified) (4 diffs)
• test_quantity_setting_functions.py (modified) (2 diffs)

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

@@ -828 +828 @@
              output_dir='TIFS',
              bounding_polygon=None,
-             verbose=False):
+             verbose=False,
+             k_nearest_neighbours=3):
     """
         Make a georeferenced tif by nearest-neighbour interpolation of sww file outputs (or a 3-column array with xyz Points)
@@ -852 +853 @@
                 output_dir -- Write outputs to this directory
                 bounding_polygon -- polygon (e.g. from read_polygon) If present, only set values of raster cells inside the bounding_polygon
+                k_nearest_neighbours -- how many neighbours to use in interpolation. If k>1, inverse-distance-weighted interpolation is used
                 
     """
@@ -860 +862 @@
     import scipy.io
     import scipy.interpolate
+    import scipy.spatial
     import anuga
     from anuga.utilities import plot_utils as util
@@ -956 +959 @@
         print 'Making interpolation functions...'
     swwXY=scipy.array([swwX[:],swwY[:]]).transpose()
-    # Get index of nearest point
-    index_qFun=scipy.interpolate.NearestNDInterpolator(swwXY,scipy.arange(len(swwX),dtype='int64').transpose())
-    #index_qFun=scipy.interpolate.LinearNDInterpolator(swwXY,scipy.arange(len(swwX),dtype='int64').transpose())
-
-    #print 'index_qFun', index_qFun
-    
+
+    # Get function to interpolate quantity onto gridXY_array
     gridXY_array=scipy.array([scipy.concatenate(gridX),scipy.concatenate(gridY)]).transpose()
-    #gridXY_array=scipy.concatenate((gridX,gridY)).transpose()
-    
     gridXY_array=scipy.ascontiguousarray(gridXY_array)
-    #print gridXY_array.flags['C_CONTIGUOUS']
-    #print 'GRIDXY_ARRAY', gridXY_array.shape
-
-    gridqInd=index_qFun(gridXY_array)
+
+    # Create Interpolation function
+    #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)
+        # Function to do the interpolation
+        def myInterpFun(quantity):
+            return quantity[gridqInd]
+    else:
+        # Combined nearest neighbours and inverse-distance interpolation
+        index_qFun=scipy.spatial.cKDTree(swwXY)
+        NNInfo=index_qFun.query(gridXY_array,k=3)
+        # 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]
+            return (num/denom)
+
+
 
     if(bounding_polygon is not None):
@@ -995 +1012 @@
             if(type(myTS)==int):
                 if(output_quantity=='stage'):
-                    gridq=p2.stage[myTS,:][gridqInd]
+                    gridq=myInterpFun(p2.stage[myTS,:])
                 if(output_quantity=='depth'):
-                    gridq=p2.height[myTS,:][gridqInd]
+                    gridq=myInterpFun(p2.height[myTS,:])
                     gridq=gridq*(gridq>=0.) # Force positive depth (tsunami alg)
                 if(output_quantity=='velocity'):
-                    gridq=p2.vel[myTS,:][gridqInd]
+                    gridq=myInterpFun(p2.vel[myTS,:])
                 if(output_quantity=='friction'):
-                    gridq=p2.friction[gridqInd]
+                    gridq=myInterpFun(p2.friction)
                 if(output_quantity=='depthIntegratedVelocity'):
                     swwDIVel=(p2.xmom[myTS,:]**2+p2.ymom[myTS,:]**2)**0.5
-                    gridq=swwDIVel[gridqInd]
+                    gridq=myInterpFun(swwDIVel)
                 if(output_quantity=='elevation'):
-                    gridq=p2.elev[gridqInd]
+                    gridq=myInterpFun(p2.elev)
     
                 if(myTSi is 'max'):
@@ -1014 +1031 @@
                     timestepString=str(round(p2.time[myTS]))
             elif(myTS=='pointData'):
-                gridq=xyzPoints[:,2][gridqInd]
+                gridq=myInterpFun(xyzPoints[:,2])
 
 

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

@@ -116 +116 @@
     import os
     import numpy
-    from matplotlib import path
+    #from matplotlib import path
+    from anuga.geometry.polygon import inside_polygon
 
     def F(x,y):
@@ -151 +152 @@
                 continue
 
-            # Get a 'true/false' flag for whether the point is in pi, and not set
+            # Get indices fInds of points in pi which are not set
             if(pi is 'All'):
                 fInside=(1-isSet)
+                fInds=(fInside==1).nonzero()[0]
             else:
                 if(pi is 'Extent' and type(fi) is str and os.path.exists(fi)):
@@ -159 +161 @@
                     # Then we get the extent from the raster itself
                     import anuga.utilities.spatialInputUtil as su
-                    newpi=su.getRasterExtent(fi,asPolygon=True)
-                    pi_path=path.Path(numpy.array(newpi))
+                    pi_path=su.getRasterExtent(fi,asPolygon=True)
                 else:
-                    # Try matplotlib -- make the path a matplotlib path object
-                    pi_path=path.Path(numpy.array(pi))
-
-                fInside=xy_array_trans[:,0]*0.
-                for j in range(len(fInside)):
-                    fInside[j] = pi_path.contains_point(xy_array_trans[j,:])
-                
-                fInside=fInside*(1-isSet)
-            #
-            # Find indices of points we need to adjust 
-            fInds=(fInside==1.).nonzero()[0]
+                    pi_path=pi
+                notSet=(isSet==0.).nonzero()[0]
+                fInds = inside_polygon(xy_array_trans[notSet,:], pi_path)
+                fInds = notSet[fInds]
+             
             if(len(fInds)>0):
                 if(hasattr(fi,'__call__')):

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

@@ -48 +48 @@
 import struct
 import numpy
-from matplotlib import path
+#from matplotlib import path
+from anuga.geometry.polygon import inside_polygon
 
 try:
@@ -627 +628 @@
     
                       The x/y point spacing on the trial-grid will be close to
-                      approx_grid_spacing, but we ensure there are at least 3x3 points on the trial grid.
+                      approx_grid_spacing, but we ensure there are at least 4x4 points on the trial grid.
                       Also, we reduce the spacing so that the min_x+R and max_x-R
                         are both similarly close to the polygon extents [see definition of R below]
@@ -655 +656 @@
     
         # Make a 'grid' of points which covers the polygon
-        xGridCount=max( numpy.ceil( (poly_xmax-poly_xmin)/approx_grid_spacing[0]+1. ).astype(int), 3)
+        xGridCount=max( numpy.ceil( (poly_xmax-poly_xmin)/approx_grid_spacing[0]+1. ).astype(int), 4)
         R=(poly_xmax-poly_xmin)*eps
         Xvals=numpy.linspace(poly_xmin+R,poly_xmax-R, xGridCount)
-        yGridCount=max( numpy.ceil( (poly_ymax-poly_ymin)/approx_grid_spacing[1]+1. ).astype(int), 3)
+        yGridCount=max( numpy.ceil( (poly_ymax-poly_ymin)/approx_grid_spacing[1]+1. ).astype(int), 4)
         R=(poly_ymax-poly_ymin)*eps
         Yvals=numpy.linspace(poly_ymin+R,poly_ymax-R, yGridCount)
@@ -665 +666 @@
         Grid=numpy.vstack([xGrid.flatten(),yGrid.flatten()]).transpose()
     
-        # 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()
+        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()
         xyInside=Grid[keepers,:]
     

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

@@ -153 +153 @@
         allDat=numpy.vstack([xs,ys,elev]).transpose()
         util.Make_Geotif(allDat, output_quantities=['ElevTest'], EPSG_CODE=32756, 
-                        output_dir='.', CellSize=1.)
+                        output_dir='.', CellSize=1.,k_nearest_neighbours=1)
 
         # Make a polygon-point pair which we use to set elevation in a 'channel'
@@ -232 +232 @@
         allDat=numpy.vstack([xs,ys,elev]).transpose()
         util.Make_Geotif(allDat, output_quantities=['ElevTest'], EPSG_CODE=32756, 
-                        output_dir='.', CellSize=1.)
+                        output_dir='.', CellSize=1.,k_nearest_neighbours=1)
 
         # Make a polygon-point pair which we use to set elevation in a 'channel'