Changeset 9040 for trunk/anuga_work/development/gareth/experimental

Timestamp:

Dec 10, 2013, 12:31:52 PM (11 years ago)

Author:

davies

Message:

Removing unique vertex storage from DE1 and tsunami. Using centroid storage instead
Also some changes to anuga_work/development/gareth

File:

• trunk/anuga_work/development/gareth/experimental/bal_and/plot_utils.py (modified) (1 diff)

trunk/anuga_work/development/gareth/experimental/bal_and/plot_utils.py

@@ -494 +494 @@
 
 
-
+## FOLLOWING FUNCTION MODIFIED FROM STACK OVERFLOW
+def make_grid(data, lats, lons, fileName, EPSG_CODE=None):
+    """
+        Convert data,lats,lons to a georeferenced raster tif
+        INPUT: data -- array with desired raster cell values
+               lats -- 1d array with 'latitude' or 'y' range
+               lons -- 1D array with 'longitude' or 'x' range
+               fileName -- name of file to write to
+               EPSG_CODE -- Integer code with projection information in EPSG format 
+    """
+    try:
+        import gdal
+        import osr
+    except:
+        raise Exception, 'Cannot find gdal and/or osr python modules'
+
+    #data = scipy.random.rand(5,6)
+    #lats = scipy.array([-5.5, -5.0, -4.5, -4.0, -3.5])
+    #lons = scipy.array([135.0, 135.5, 136.0, 136.5, 137.0, 137.5])
+    #import pdb
+    #pdb.set_trace()
+    xres = lons[1] - lons[0]
+    yres = lats[1] - lats[0]
+
+    ysize = len(lats)
+    xsize = len(lons)
+
+    # Assume data/lats/longs refer to cell centres, as per gdal
+    #ulx = lons[0] - (xres / 2.)
+    #uly = lats[-1] - (yres / 2.)
+    llx = lons[0] - (xres / 2.)
+    lly = lats[0] - (yres / 2.)
+
+    driver = gdal.GetDriverByName('GTiff')
+    ds = driver.Create(fileName, xsize, ysize, 1, gdal.GDT_Float32)
+
+    # this assumes the projection is Geographic lat/lon WGS 84
+    srs = osr.SpatialReference()
+    srs.ImportFromEPSG(EPSG_CODE)
+    ds.SetProjection(srs.ExportToWkt())
+
+    #gt = [ulx, xres, 0, uly, 0, yres ]
+    gt = [llx, xres, 0, lly, 0, yres ]
+    ds.SetGeoTransform(gt)
+
+    outband = ds.GetRasterBand(1)
+    outband.WriteArray(data)
+
+    ds = None
+    return
+
+##################################################################################
+
+def Make_Geotif(swwFile, output_quantity='depth',
+             myTimeStep=1, CellSize=5.0, 
+             lower_left=None, upper_right=None,
+             EPSG_CODE=None, 
+             velocity_extrapolation=True,
+             min_allowed_height=1.0e-05,
+             output_dir='TIFS',
+             verbose=False):
+    """
+        Make a georeferenced tif by nearest-neighbour interpolation of sww file outputs.
+        INPUTS: swwFile -- name of sww file
+                output_quantity -- quantity to plot ('depth' or 'velocity' or 'stage' or 'elevation')
+                myTimeStep -- time-index of swwFile to plot, or 'last', or 'max'
+                CellSize -- approximate pixel size for output raster [adapted to fit lower_left / upper_right]
+                lower_left -- [x0,y0] of lower left corner. If None, use extent of swwFile.
+                upper_right -- [x1,y1] of upper right corner. If None, use extent of swwFile.
+                EPSG_CODE -- Projection information as an integer EPSG code (e.g. 3123 for PRS92 Zone 3). 
+                             Google for info on EPSG Codes
+                velocity_extrapolation -- Compute velocity assuming the code extrapolates with velocity (instead of momentum)?
+                min_allowed_height -- Minimum allowed height from ANUGA
+                output_dir -- Write outputs to this directory
+                
+    """
+    try:
+        import gdal
+        import osr
+        import scipy.io
+        import scipy.interpolate
+        import anuga
+        from anuga.utilities import plot_utils as util
+        import os
+    except:
+        raise Exception, 'Required modules not installed for mkgeotif'
+
+
+    if(EPSG_CODE is None):
+        raise Exception, 'Must specify an integer EPSG_CODE describing the file projection'
+    ###
+    ## INPUTS
+    ###
+    #swwFile='taguig_rainfall_ondoy.sww'
+    #myTimeStep=40 # Index of time-slice in sWW to plot
+    ## Plot extents
+    #lower_left=None #[x0,y0]
+    #upper_right=None #[x1,y1]
+    ## Projection as an epsg code
+    #EPSG_CODE=3123
+    ## CellSize for output raster
+    #CellSize=5.0 # APPROXIMATE CellSize
+    ## ANUGA Parameters
+    #min_allowed_height=1.0e-05
+    #velocity_extrapolation=True
+
+    #output_quantity='depth' # 'velocity' 'depth' 'stage'
+    #output_dir='TEST'
+
+    ##
+    # END INPUTS
+    ##
+
+    # Make output_dir
+    try:
+        os.mkdir(output_dir)
+    except:
+        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=='last'):
+        myTimeStep=len(p.time)-1
+    
+
+    if(verbose):
+        print 'Extracting required data ...'
+    # Get ANUGA points
+    swwX=p2.x+xllcorner
+    swwY=p2.y+yllcorner
+
+    if(myTimeStep!='max'):
+        swwStage=p2.stage[myTimeStep,:]
+        swwHeight=p2.height[myTimeStep,:]
+        swwVel=p2.vel[myTimeStep,:]
+        timestepString=str(round(p2.time[myTimeStep]))
+    else:
+        swwStage=p2.stage.max(axis=0)
+        swwHeight=p2.height.max(axis=0)
+        swwVel=p2.vel.max(axis=0)
+        timestepString='max'
+       
+    #import pdb
+    #pdb.set_trace() 
+
+    # Make name for output file
+    output_name=output_dir+'/'+os.path.splitext(os.path.basename(swwFile))[0] + '_'+\
+                output_quantity+'_'+timestepString+\
+                '_'+str(myTimeStep)+'.tif'
+
+    if(verbose):
+        print 'Computing grid of output locations...'
+    # Get points where we want raster cells
+    if(lower_left is None):
+        lower_left=[swwX.min(),swwY.min()]
+    if(upper_right is None):
+        upper_right=[swwX.max(),swwY.max()]
+
+    nx=round((upper_right[0]-lower_left[0])*1.0/(1.0*CellSize)) + 1
+    xres=(upper_right[0]-lower_left[0])*1.0/(1.0*(nx-1))
+    desiredX=scipy.arange(lower_left[0], upper_right[0],xres )
+    ny=round((upper_right[1]-lower_left[1])*1.0/(1.0*CellSize)) + 1
+    yres=(upper_right[1]-lower_left[1])*1.0/(1.0*(ny-1))
+    desiredY=scipy.arange(lower_left[1], upper_right[1], yres)
+
+    gridX, gridY=scipy.meshgrid(desiredX,desiredY)
+
+    # Make functions to interpolate quantities at points
+    if(verbose):
+        print 'Making interpolation functions...'
+    swwXY=scipy.array([swwX[:],swwY[:]]).transpose()
+    #swwXY=scipy.array([scipy.concatenate(gridX), scipy.concatenate(gridY)]).transpose()
+    if(output_quantity=='stage'):
+        qFun=scipy.interpolate.NearestNDInterpolator(swwXY,swwStage.transpose())
+    elif(output_quantity=='velocity'):
+        qFun=scipy.interpolate.NearestNDInterpolator(swwXY,swwVel.transpose())
+    elif(output_quantity=='elevation'):
+        qFun=scipy.interpolate.NearestNDInterpolator(swwXY,p2.elev.transpose())
+    elif(output_quantity=='depth'):
+        qFun=scipy.interpolate.NearestNDInterpolator(swwXY,swwHeight.transpose())
+    else:
+        raise Exception, 'output_quantity not available -- check if it is spelled correctly'
+
+    if(verbose):
+        print 'Interpolating'
+    gridq=qFun(scipy.array([scipy.concatenate(gridX),scipy.concatenate(gridY)]).transpose()).transpose()
+    gridq.shape=(len(desiredY),len(desiredX))
+    #gridq=gridq.transpose()
+
+    if(verbose):
+        print 'Making raster ...'
+    make_grid(gridq,desiredY,desiredX, output_name,EPSG_CODE)