Changeset 9062 for trunk/anuga_core/source/anuga/utilities/plot_utils.py

Timestamp:

Mar 14, 2014, 11:31:16 AM (11 years ago)

Author:

davies

Message:

Allowing plot_utils.Make_Geotif to take an xyz array

File:

• trunk/anuga_core/source/anuga/utilities/plot_utils.py (modified) (10 diffs)

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

@@ -583 +583 @@
 ##################################################################################
 
-def Make_Geotif(swwFile, output_quantities=['depth'],
+def Make_Geotif(swwFile=None, 
+             output_quantities=['depth'],
              myTimeStep=1, CellSize=5.0, 
              lower_left=None, upper_right=None,
@@ -594 +595 @@
              verbose=False):
     """
-        Make a georeferenced tif by nearest-neighbour interpolation of sww file outputs.
+        Make a georeferenced tif by nearest-neighbour interpolation of sww file outputs (or a 3-column array with xyz Points)
 
         You must supply projection information as either a proj4string or an integer EPSG_CODE (but not both!)
 
-        INPUTS: swwFile -- name of sww file
+        INPUTS: swwFile -- name of sww file, OR a 3-column array with x/y/z
+                    points. In the latter case x and y are assumed to be in georeferenced
+                    coordinates.  The output raster will contain 'z', and will have a name-tag
+                    based on the name in 'output_quantities'.
                 output_quantities -- list of quantitiies to plot, e.g. ['depth', 'velocity', 'stage','elevation','depthIntegratedVelocity']
                 myTimeStep -- list containing time-index of swwFile to plot (e.g. [1, 10, 32] ) or 'last', or 'max', or 'all'
@@ -614 +618 @@
                 
     """
+
+    #import pdb
+    #pdb.set_trace()
+
     try:
         import gdal
@@ -624 +632 @@
         from matplotlib import nxutils
     except:
-        raise Exception, 'Required modules not installed for mkgeotif'
-
+        raise Exception, 'Required modules not installed for Make_Geotif'
+
+
+    # Check whether swwFile is an array, and if so, redefine various inputs to
+    # make the code work
+    if(type(swwFile)==scipy.ndarray):
+        import copy
+        xyzPoints=copy.copy(swwFile)
+        swwFile=None
 
     if(((EPSG_CODE is None) & (proj4string is None) )|
@@ -638 +653 @@
         pass
 
-    # Read in ANUGA outputs
-    # FIXME: It would be good to support reading of data subsets
-    if(verbose):
-        print 'Reading sww File ...'
-    p=util.get_output(swwFile,min_allowed_height)
-    p2=util.get_centroids(p,velocity_extrapolation)
-    swwIn=scipy.io.netcdf_file(swwFile)
-    xllcorner=swwIn.xllcorner
-    yllcorner=swwIn.yllcorner
-
-    if(myTimeStep=='all'):
-        myTimeStep=range(len(p2.time))
-    # Ensure myTimeStep is a list
-    if type(myTimeStep)!=list:
-        myTimeStep=[myTimeStep]
-
-    if(verbose):
-        print 'Extracting required data ...'
-    # Get ANUGA points
-    swwX=p2.x+xllcorner
-    swwY=p2.y+yllcorner
+    if(swwFile is not None):
+        # Read in ANUGA outputs
+        # FIXME: It would be good to support reading of data subsets
+        if(verbose):
+            print 'Reading sww File ...'
+        p=util.get_output(swwFile,min_allowed_height)
+        p2=util.get_centroids(p,velocity_extrapolation)
+        swwIn=scipy.io.netcdf_file(swwFile)
+        xllcorner=swwIn.xllcorner
+        yllcorner=swwIn.yllcorner
+
+        if(myTimeStep=='all'):
+            myTimeStep=range(len(p2.time))
+        # Ensure myTimeStep is a list
+        if type(myTimeStep)!=list:
+            myTimeStep=[myTimeStep]
+
+        if(verbose):
+            print 'Extracting required data ...'
+        # Get ANUGA points
+        swwX=p2.x+xllcorner
+        swwY=p2.y+yllcorner
+    else:
+        # Get the point data from the 3-column array
+        if(xyzPoints.shape[1]!=3):
+            raise Exception, 'If an array is passed, it must have exactly 3 columns'
+        if(len(output_quantities)!=1):
+            raise Exception, 'Can only have 1 output quantity when passing an array'
+        swwX=xyzPoints[:,0]
+        swwY=xyzPoints[:,1]
+        myTimeStep=['pointData']
 
     # Grid for meshing
@@ -681 +706 @@
     swwXY=scipy.array([swwX[:],swwY[:]]).transpose()
     # Get index of nearest point
-    index_qFun=scipy.interpolate.NearestNDInterpolator(swwXY,scipy.arange(len(swwX),dtype='int32').transpose())
+    index_qFun=scipy.interpolate.NearestNDInterpolator(swwXY,scipy.arange(len(swwX),dtype='int64').transpose())
     gridXY_array=scipy.array([scipy.concatenate(gridX),scipy.concatenate(gridY)]).transpose()
     gridqInd=index_qFun(gridXY_array)
@@ -702 +727 @@
         
 
-            if(myTS!='max'):
+            #if(myTS!='max'):
+            if(type(myTS)=='int'):
                 if(output_quantity=='stage'):
                     gridq=p2.stage[myTS,:][gridqInd]
@@ -715 +741 @@
                     gridq=p2.elev[gridqInd]
                 timestepString=str(round(p2.time[myTS]))
-            else:
+            elif (myTS=='max'):
                 if(output_quantity=='stage'):
                     gridq=p2.stage.max(axis=0)[gridqInd]
@@ -728 +754 @@
                     gridq=p2.elev[gridqInd]
                 timestepString='max'
+            elif(myTS=='pointData'):
+                gridq=xyzPoints[:,2][gridqInd]
+
 
             if(bounding_polygon is not None):
@@ -734 +763 @@
 
             # Make name for output file
-            output_name=output_dir+'/'+os.path.splitext(os.path.basename(swwFile))[0] + '_'+\
-                        output_quantity+'_'+timestepString+\
-                        '_'+str(myTS)+'.tif'
+            if(myTS!='pointData'):
+                output_name=output_dir+'/'+os.path.splitext(os.path.basename(swwFile))[0] + '_'+\
+                            output_quantity+'_'+timestepString+\
+                            '_'+str(myTS)+'.tif'
+            else:
+                output_name=output_dir+'/'+'PointData_'+output_quantity+'.tif'
 
             if(verbose):