Changeset 9044


Ignore:
Timestamp:
Dec 13, 2013, 12:01:48 PM (11 years ago)
Author:
davies
Message:

Adding Towradgi validation test

Location:
trunk/anuga_core/source
Files:
285 added
1 edited

Legend:

Unmodified
Added
Removed

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

    r9042 r9044  
    142142        self.x, self.y, self.time, self.vols, self.stage, \
    143143                self.height, self.elev, self.xmom, self.ymom, \
    144                 self.xvel, self.yvel, self.vel, self.minimum_allowed_height = \
     144                self.xvel, self.yvel, self.vel, self.minimum_allowed_height,\
     145                self.xllcorner, self.yllcorner = \
    145146                read_output(filename, minimum_allowed_height)
    146147        self.filename=filename
    … …  
    165166    # Open ncdf connection
    166167    fid=NetCDFFile(filename)
     168   
     169    # Get lower-left
     170    xllcorner=fid.xllcorner
     171    yllcorner=fid.yllcorner
    167172
    168173
    … …  
    214219    vel = (xvel**2+yvel**2)**0.5
    215220
    216     return x, y, time, vols, stage, height, elev, xmom, ymom, xvel, yvel, vel, minimum_allowed_height
     221    return x, y, time, vols, stage, height, elev, xmom, ymom, xvel, yvel, vel, minimum_allowed_height, xllcorner,yllcorner
    217222
    218223##############
    … …  
    574579##################################################################################
    575580
    576 def Make_Geotif(swwFile, output_quantity='depth',
     581def Make_Geotif(swwFile, output_quantities=['depth'],
    577582             myTimeStep=1, CellSize=5.0,
    578583             lower_left=None, upper_right=None,
    … …  
    589594
    590595        INPUTS: swwFile -- name of sww file
    591                 output_quantity -- quantity to plot ('depth' or 'velocity' or 'stage' or 'elevation' or 'depthIntegratedVelocity')
     596                output_quantities -- list of quantitiies to plot, e.g. ['depth', 'velocity', 'stage','elevation','depthIntegratedVelocity']
    592597                myTimeStep -- time-index of swwFile to plot, or 'last', or 'max'
    593598                CellSize -- approximate pixel size for output raster [adapted to fit lower_left / upper_right]
    … …  
    645650    swwY=p2.y+yllcorner
    646651
    647     if(myTimeStep!='max'):
    648         swwStage=p2.stage[myTimeStep,:]
    649         swwHeight=p2.height[myTimeStep,:]
    650         swwVel=p2.vel[myTimeStep,:]
    651         # Use if-statement here to reduce computation
    652         if(output_quantity=='depthIntegratedVelocity'):
    653             swwDIVel=(p2.xmom[myTimeStep,:]**2+p2.ymom[myTimeStep,:]**2)**0.5
    654         timestepString=str(round(p2.time[myTimeStep]))
    655     else:
    656         swwStage=p2.stage.max(axis=0)
    657         swwHeight=p2.height.max(axis=0)
    658         swwVel=p2.vel.max(axis=0)
    659         # Use if-statement here to reduce computation
    660         if(output_quantity=='depthIntegratedVelocity'):
    661             swwDIVel=((p2.xmom**2+p2.ymom**2).max(axis=0))**0.5
    662         timestepString='max'
    663        
    664     # Make name for output file
    665     output_name=output_dir+'/'+os.path.splitext(os.path.basename(swwFile))[0] + '_'+\
    666                 output_quantity+'_'+timestepString+\
    667                 '_'+str(myTimeStep)+'.tif'
    668 
    669     if(verbose):
    670         print 'Computing grid of output locations...'
    671     # Get points where we want raster cells
    672     if(lower_left is None):
    673         lower_left=[swwX.min(),swwY.min()]
    674     if(upper_right is None):
    675         upper_right=[swwX.max(),swwY.max()]
    676 
    677     nx=round((upper_right[0]-lower_left[0])*1.0/(1.0*CellSize)) + 1
    678     xres=(upper_right[0]-lower_left[0])*1.0/(1.0*(nx-1))
    679     desiredX=scipy.arange(lower_left[0], upper_right[0],xres )
    680     ny=round((upper_right[1]-lower_left[1])*1.0/(1.0*CellSize)) + 1
    681     yres=(upper_right[1]-lower_left[1])*1.0/(1.0*(ny-1))
    682     desiredY=scipy.arange(lower_left[1], upper_right[1], yres)
    683 
    684     gridX, gridY=scipy.meshgrid(desiredX,desiredY)
    685 
    686     # Make functions to interpolate quantities at points
    687     if(verbose):
    688         print 'Making interpolation functions...'
    689     swwXY=scipy.array([swwX[:],swwY[:]]).transpose()
    690     if(output_quantity=='stage'):
    691         qFun=scipy.interpolate.NearestNDInterpolator(swwXY,swwStage.transpose())
    692     elif(output_quantity=='velocity'):
    693         qFun=scipy.interpolate.NearestNDInterpolator(swwXY,swwVel.transpose())
    694     elif(output_quantity=='elevation'):
    695         qFun=scipy.interpolate.NearestNDInterpolator(swwXY,p2.elev.transpose())
    696     elif(output_quantity=='depth'):
    697         qFun=scipy.interpolate.NearestNDInterpolator(swwXY,swwHeight.transpose())
    698     elif(output_quantity=='depthIntegratedVelocity'):
    699         qFun=scipy.interpolate.NearestNDInterpolator(swwXY,swwDIVel.transpose())
    700     else:
    701         raise Exception, 'output_quantity not available -- check if it is spelled correctly'
    702 
    703     if(verbose):
    704         print 'Interpolating'
    705     gridq=qFun(scipy.array([scipy.concatenate(gridX),scipy.concatenate(gridY)]).transpose()).transpose()
    706     gridq.shape=(len(desiredY),len(desiredX))
    707 
    708     if(verbose):
    709         print 'Making raster ...'
    710     make_grid(gridq,desiredY,desiredX, output_name,EPSG_CODE=EPSG_CODE, proj4string=proj4string)
     652    # Loop over all output quantities and produce the output
     653    for output_quantity in output_quantities:
     654
     655        if(myTimeStep!='max'):
     656            if(output_quantity=='stage'):
     657                swwStage=p2.stage[myTimeStep,:]
     658            if(output_quantity=='depth'):
     659                swwHeight=p2.height[myTimeStep,:]
     660            if(output_quantity=='velocity'):
     661                swwVel=p2.vel[myTimeStep,:]
     662            if(output_quantity=='depthIntegratedVelocity'):
     663                swwDIVel=(p2.xmom[myTimeStep,:]**2+p2.ymom[myTimeStep,:]**2)**0.5
     664            timestepString=str(round(p2.time[myTimeStep]))
     665        else:
     666            if(output_quantity=='stage'):
     667                swwStage=p2.stage.max(axis=0)
     668            if(output_quantity=='depth'):
     669                swwHeight=p2.height.max(axis=0)
     670            if(output_quantity=='velocity'):
     671                swwVel=p2.vel.max(axis=0)
     672            if(output_quantity=='depthIntegratedVelocity'):
     673                swwDIVel=((p2.xmom**2+p2.ymom**2).max(axis=0))**0.5
     674            timestepString='max'
     675           
     676        # Make name for output file
     677        output_name=output_dir+'/'+os.path.splitext(os.path.basename(swwFile))[0] + '_'+\
     678                    output_quantity+'_'+timestepString+\
     679                    '_'+str(myTimeStep)+'.tif'
     680
     681        if(verbose):
     682            print 'Computing grid of output locations...'
     683        # Get points where we want raster cells
     684        if(lower_left is None):
     685            lower_left=[swwX.min(),swwY.min()]
     686        if(upper_right is None):
     687            upper_right=[swwX.max(),swwY.max()]
     688
     689        nx=round((upper_right[0]-lower_left[0])*1.0/(1.0*CellSize)) + 1
     690        xres=(upper_right[0]-lower_left[0])*1.0/(1.0*(nx-1))
     691        desiredX=scipy.arange(lower_left[0], upper_right[0],xres )
     692        ny=round((upper_right[1]-lower_left[1])*1.0/(1.0*CellSize)) + 1
     693        yres=(upper_right[1]-lower_left[1])*1.0/(1.0*(ny-1))
     694        desiredY=scipy.arange(lower_left[1], upper_right[1], yres)
     695
     696        gridX, gridY=scipy.meshgrid(desiredX,desiredY)
     697
     698        # Make functions to interpolate quantities at points
     699        if(verbose):
     700            print 'Making interpolation functions...'
     701        swwXY=scipy.array([swwX[:],swwY[:]]).transpose()
     702        if(output_quantity=='stage'):
     703            qFun=scipy.interpolate.NearestNDInterpolator(swwXY,swwStage.transpose())
     704        elif(output_quantity=='velocity'):
     705            qFun=scipy.interpolate.NearestNDInterpolator(swwXY,swwVel.transpose())
     706        elif(output_quantity=='elevation'):
     707            qFun=scipy.interpolate.NearestNDInterpolator(swwXY,p2.elev.transpose())
     708        elif(output_quantity=='depth'):
     709            qFun=scipy.interpolate.NearestNDInterpolator(swwXY,swwHeight.transpose())
     710        elif(output_quantity=='depthIntegratedVelocity'):
     711            qFun=scipy.interpolate.NearestNDInterpolator(swwXY,swwDIVel.transpose())
     712        else:
     713            raise Exception, 'output_quantity not available -- check if it is spelled correctly'
     714
     715        if(verbose):
     716            print 'Interpolating'
     717        gridq=qFun(scipy.array([scipy.concatenate(gridX),scipy.concatenate(gridY)]).transpose()).transpose()
     718        gridq.shape=(len(desiredY),len(desiredX))
     719
     720        if(verbose):
     721            print 'Making raster ...'
     722        make_grid(gridq,desiredY,desiredX, output_name,EPSG_CODE=EPSG_CODE, proj4string=proj4string)
    711723
    712724    return
     725
     726def plot_triangles(p):
     727    """ Add mesh triangles to a pyplot plot
     728    """
     729    for i in range(len(p.vols)):
     730        k1=p.vols[i][0]
     731        k2=p.vols[i][1]
     732        k3=p.vols[i][2]
     733        pyplot.plot([p.x[k1], p.x[k2], p.x[k3], p.x[k1]], [p.y[k1], p.y[k2], p.y[k3], p.y[k1]],'-',color='black')
     734        #pyplot.plot([p.x[k3], p.x[k2]], [p.y[k3], p.y[k2]],'-',color='black')
     735        #pyplot.plot([p.x[k3], p.x[k1]], [p.y[k3], p.y[k1]],'-',color='black')
     736
     737def find_neighbours(p,ind):
     738    """
     739        Find the triangles neighbouring triangle 'ind'
     740        p is an object from get_output containing mesh vertices
     741    """
     742    ind_nei=p.vols[ind]
     743   
     744    shared_nei0=p.vols[:,1]*0.0
     745    shared_nei1=p.vols[:,1]*0.0
     746    shared_nei2=p.vols[:,1]*0.0
     747    # Compute indices that match one of the vertices of triangle ind
     748    # Note: Each triangle can only match a vertex, at most, once
     749    for i in range(3):
     750        shared_nei0+=1*(p.x[p.vols[:,i]]==p.x[ind_nei[0]])*\
     751            1*(p.y[p.vols[:,i]]==p.y[ind_nei[0]])
     752       
     753        shared_nei1+=1*(p.x[p.vols[:,i]]==p.x[ind_nei[1]])*\
     754            1*(p.y[p.vols[:,i]]==p.y[ind_nei[1]])
     755       
     756        shared_nei2+=1*(p.x[p.vols[:,i]]==p.x[ind_nei[2]])*\
     757            1*(p.y[p.vols[:,i]]==p.y[ind_nei[2]])
     758   
     759    out=(shared_nei2 + shared_nei1 + shared_nei0)
     760    return((out==2).nonzero())
     761
     762def calc_edge_elevations(p):
     763    """
     764        Compute the triangle edge elevations on p
     765        Return x,y,elev for edges
     766    """
     767    pe_x=p.x*0.
     768    pe_y=p.y*0.
     769    pe_el=p.elev*0.
     770
     771   
     772    # Compute coordinates + elevations
     773    pe_x[p.vols[:,0]] = 0.5*(p.x[p.vols[:,1]] + p.x[p.vols[:,2]])
     774    pe_y[p.vols[:,0]] = 0.5*(p.y[p.vols[:,1]] + p.y[p.vols[:,2]])
     775    pe_el[p.vols[:,0]] = 0.5*(p.elev[p.vols[:,1]] + p.elev[p.vols[:,2]])
     776   
     777    pe_x[p.vols[:,1]] = 0.5*(p.x[p.vols[:,0]] + p.x[p.vols[:,2]])
     778    pe_y[p.vols[:,1]] = 0.5*(p.y[p.vols[:,0]] + p.y[p.vols[:,2]])
     779    pe_el[p.vols[:,1]] = 0.5*(p.elev[p.vols[:,0]] + p.elev[p.vols[:,2]])
     780
     781    pe_x[p.vols[:,2]] = 0.5*(p.x[p.vols[:,0]] + p.x[p.vols[:,1]])
     782    pe_y[p.vols[:,2]] = 0.5*(p.y[p.vols[:,0]] + p.y[p.vols[:,1]])
     783    pe_el[p.vols[:,2]] = 0.5*(p.elev[p.vols[:,0]] + p.elev[p.vols[:,1]])
     784
     785    return [pe_x, pe_y, pe_el]
     786
     787
Note: See TracChangeset for help on using the changeset viewer.