Changeset 9051

Timestamp:

Feb 11, 2014, 9:09:53 PM (11 years ago)

Author:

davies

Message:

Fix makeGeoTif + minor tweaks

Location:

trunk

Files:

• anuga_core/source/anuga/shallow_water/swDE1_domain_ext.c (modified) (1 diff)
• anuga_core/source/anuga/utilities/plot_utils.py (modified) (6 diffs)
• anuga_core/source/anuga_validation_tests/case_studies/towradgi/Compare_results_with_fieldObs.py (modified) (2 diffs)
• anuga_core/source/anuga_validation_tests/case_studies/towradgi/Towradgi_historic_flood.py (modified) (1 diff)
• anuga_work/development/gareth/tests/Petar/Towradgi/Aug98_hort_discontinous.py (modified) (2 diffs)
• anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_sparse.py (modified) (2 diffs)
• anuga_work/development/gareth/tests/taguig_parallel/runThis.sh (modified) (1 diff)
• anuga_work/development/gareth/tests/taguig_parallel/runtaguig.py (modified) (5 diffs)

trunk/anuga_core/source/anuga/shallow_water/swDE1_domain_ext.c

@@ -558 +558 @@
 
 
-            // If neighbour is a boundary condition, add the flux to the boundary_flux_integral
-            if(n<0){ 
+            // If this cell is not a ghost, and the neighbour is a boundary
+            // condition OR a ghost cell, then add the flux to the
+            // boundary_flux_integral
+            if( (n<0 & tri_full_flag[k]==1) | ( n>=0 && (tri_full_flag[k]==1 & tri_full_flag[n]==0)) ){ 
                 boundary_flux_sum[0] += edgeflux_store[ki3];
             }

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

@@ -550 +550 @@
     xsize = len(lons)
 
-    # Assume data/lats/longs refer to cell centres, and compute lower left coordinate
-    llx = lons[0] - (xres / 2.)
-    lly = lats[0] - (yres / 2.)
+    # Assume data/lats/longs refer to cell centres, and compute upper left coordinate
+    ulx = lons[0] - (xres / 2.)
+    uly = lats[lats.shape[0]-1] + (yres / 2.)
 
     # GDAL magic to make the tif
@@ -568 +568 @@
     ds.SetProjection(srs.ExportToWkt())
 
-    gt = [llx, xres, 0, lly, 0, yres ]
+    gt = [ulx, xres, 0, uly, 0, -yres ]
+    #gt = [llx, xres, 0, lly, yres,0 ]
     ds.SetGeoTransform(gt)
+
+    #import pdb
+    #pdb.set_trace()
 
     outband = ds.GetRasterBand(1)
@@ -595 +599 @@
         INPUTS: swwFile -- name of sww file
                 output_quantities -- list of quantitiies to plot, e.g. ['depth', 'velocity', 'stage','elevation','depthIntegratedVelocity']
-                myTimeStep -- time-index of swwFile to plot, or 'last', or 'max'
+                myTimeStep -- list containing 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.
@@ -624 +628 @@
         raise Exception, 'Must specify EITHER an integer EPSG_CODE describing the file projection, OR a proj4string'
 
+    # Ensure myTimeStep is a list
+    if type(myTimeStep)!=list:
+        myTimeStep=[myTimeStep]
+
     # Make output_dir
     try:
@@ -640 +648 @@
     yllcorner=swwIn.yllcorner
 
-    if(myTimeStep=='last'):
-        myTimeStep=len(p.time)-1
-    
 
     if(verbose):
@@ -650 +655 @@
     swwY=p2.y+yllcorner
 
+    # Grid for meshing
+    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)
+
+    if(verbose):
+        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='int32').transpose())
+    gridqInd=index_qFun(scipy.array([scipy.concatenate(gridX),scipy.concatenate(gridY)]).transpose())
+
     # Loop over all output quantities and produce the output
-    for output_quantity in output_quantities:
-
-        if(myTimeStep!='max'):
-            if(output_quantity=='stage'):
-                swwStage=p2.stage[myTimeStep,:]
-            if(output_quantity=='depth'):
-                swwHeight=p2.height[myTimeStep,:]
-            if(output_quantity=='velocity'):
-                swwVel=p2.vel[myTimeStep,:]
-            if(output_quantity=='depthIntegratedVelocity'):
-                swwDIVel=(p2.xmom[myTimeStep,:]**2+p2.ymom[myTimeStep,:]**2)**0.5
-            timestepString=str(round(p2.time[myTimeStep]))
-        else:
-            if(output_quantity=='stage'):
-                swwStage=p2.stage.max(axis=0)
-            if(output_quantity=='depth'):
-                swwHeight=p2.height.max(axis=0)
-            if(output_quantity=='velocity'):
-                swwVel=p2.vel.max(axis=0)
-            if(output_quantity=='depthIntegratedVelocity'):
-                swwDIVel=((p2.xmom**2+p2.ymom**2).max(axis=0))**0.5
-            timestepString='max'
-           
-        # Make name for output file
-        output_name=output_dir+'/'+os.path.splitext(os.path.basename(swwFile))[0] + '_'+\
-                    output_quantity+'_'+timestepString+\
-                    '_'+str(myTimeStep)+'.tif'
-
+    for myTS in myTimeStep:
         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()
-        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())
-        elif(output_quantity=='depthIntegratedVelocity'):
-            qFun=scipy.interpolate.NearestNDInterpolator(swwXY,swwDIVel.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))
-
-        if(verbose):
-            print 'Making raster ...'
-        make_grid(gridq,desiredY,desiredX, output_name,EPSG_CODE=EPSG_CODE, proj4string=proj4string)
+            print myTS
+        for output_quantity in output_quantities:
+
+            if(myTS=='last'):
+                myTS=len(p.time)-1
+        
+
+            if(myTS!='max'):
+                if(output_quantity=='stage'):
+                    gridq=p2.stage[myTS,:][gridqInd]
+                if(output_quantity=='depth'):
+                    gridq=p2.height[myTS,:][gridqInd]
+                if(output_quantity=='velocity'):
+                    gridq=p2.vel[myTS,:][gridqInd]
+                if(output_quantity=='depthIntegratedVelocity'):
+                    swwDIVel=(p2.xmom[myTS,:]**2+p2.ymom[myTS,:]**2)**0.5
+                    gridq=swwDIVel[gridqInd]
+                timestepString=str(round(p2.time[myTS]))
+            else:
+                if(output_quantity=='stage'):
+                    gridq=p2.stage.max(axis=0)[gridqInd]
+                if(output_quantity=='depth'):
+                    gridq=p2.height.max(axis=0)[gridqInd]
+                if(output_quantity=='velocity'):
+                    gridq=p2.vel.max(axis=0)[gridqInd]
+                if(output_quantity=='depthIntegratedVelocity'):
+                    swwDIVel=((p2.xmom**2+p2.ymom**2).max(axis=0))**0.5
+                    gridq=swwDIVel[gridqInd]
+                timestepString='max'
+               
+            # Make name for output file
+            output_name=output_dir+'/'+os.path.splitext(os.path.basename(swwFile))[0] + '_'+\
+                        output_quantity+'_'+timestepString+\
+                        '_'+str(myTS)+'.tif'
+
+            if(verbose):
+                print 'Making raster ...'
+            gridq.shape=(len(desiredY),len(desiredX))
+            make_grid(scipy.flipud(gridq),desiredY,desiredX, output_name,EPSG_CODE=EPSG_CODE, proj4string=proj4string)
 
     return

trunk/anuga_core/source/anuga_validation_tests/case_studies/towradgi/Compare_results_with_fieldObs.py

@@ -8 +8 @@
 from anuga.utilities import plot_utils as util
 
-p=util.get_output('MODEL_OUTPUTS/Towradgi_historic_flood.sww')
+swwdir='MODEL_OUTPUTS/'
+#swwdir='MODEL_bridges_n0p04/'
+
+p=util.get_output(swwdir+'Towradgi_historic_flood.sww')
 p2=util.get_centroids(p,velocity_extrapolation=True)
 floodLevels=scipy.genfromtxt('Validation/historic_1998_flood_levels_towradgi_ck.csv', delimiter=',',skip_header=1)
@@ -50 +53 @@
     tif_outdir='OUTPUT_TIFS'
     CellSize=5.0
-    util.Make_Geotif('MODEL_OUTPUTS/Towradgi_historic_flood.sww',
+    util.Make_Geotif(swwdir+'Towradgi_historic_flood.sww',
                       ['depth','velocity','depthIntegratedVelocity','elevation'],'max',
                       CellSize=CellSize,EPSG_CODE=32756,output_dir=tif_outdir)

trunk/anuga_core/source/anuga_validation_tests/case_studies/towradgi/Towradgi_historic_flood.py

@@ -983 +983 @@
         domain.write_time()
 
+barrier()
 if myid == 0:
     print 'Number of processors %g ' %numprocs

trunk/anuga_work/development/gareth/tests/Petar/Towradgi/Aug98_hort_discontinous.py

@@ -249 +249 @@
  
 if myid == 0 and verbose: print 'CREATING INLETS'
-"""
+
 #------------------------------------------------------------------------------
 # ENTER CULVERT DATA
@@ -670 +670 @@
                                     label=join('Runs', '7H'), # this puts the culvert hydraulics output into the same dir as the sww's
                                     verbose=False) 
-"""
+
 
 #------------------------------------------------------------------------------

trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_sparse.py

@@ -21 +21 @@
 from anuga.structures.inlet_operator import Inlet_operator
 from anuga.shallow_water.shallow_water_domain import Domain as Domain
-from bal_and import *
+#from bal_and import *
 #from balanced_dev import Domain as Domain
 #from anuga_tsunami import *
@@ -32 +32 @@
 domain = Domain(points, vertices, boundary) # Create domain
 domain.set_name('channel_SU_2_v2') # Output name
-#domain.set_flow_algorithm('DE1')
+domain.set_flow_algorithm('DE1')
 #------------------------------------------------------------------------------
 # Setup initial conditions
 #------------------------------------------------------------------------------
 def topography(x, y):
-        return -x/10.  #+ 1.*(numpy.sin(x/10.) +abs(y-50.)/10.) -0.*(x>80.) # linear bed slope
+        return -x/10.  + 1.*(numpy.sin(x/10.) +abs(y-50.)/10.) -0.*(x>80.) # linear bed slope
 
 def stagetopo(x,y):

trunk/anuga_work/development/gareth/tests/taguig_parallel/runThis.sh

@@ -2 +2 @@
 rm -r MODEL_OUTPUTS*
 # Run in parallel
-mpirun -np 12 python runtaguig.py 
+mpirun -np 40 python runtaguig.py 
 # Run in serial
-mpirun -np 5 python runtaguig.py 
+mpirun -np 80 python runtaguig.py 
 # Compare serial and parallel
 ipython compare_parallel_serial.py

trunk/anuga_work/development/gareth/tests/taguig_parallel/runtaguig.py

@@ -59 +59 @@
     #domain.ghost_layer_width=4
     domain.set_store_centroids(True)
-    #domain.set_flow_algorithm('DE1')
+    domain.set_flow_algorithm('DE1')
 
     # Define function to compute elevation
@@ -69 +69 @@
                             alpha=0.0001)
 
+
 else:
     domain=None
@@ -77 +78 @@
 barrier()
 
+domain.set_quantity('stage', 1.0)
+domain.set_quantity('friction', 0.03) 
 ## Force first order
 #domain.beta_w=0.0
@@ -85 +88 @@
 #domain.beta_vh_dry=0.0
 ###
+#print domain.quantities['stage'].centroid_values.max()
+#print domain.quantities['stage'].vertex_values.max()
 
-domain.set_quantity('stage', 1.0)
-domain.set_quantity('friction', 0.03) 
-
+#barrier()
+#raise Exception, 'Finish here'
 # Rain
-#rain_mps=90.0/(1000*3600) # 90mm/hr is meters-per-second
-#def myrain(t):
-#    return(rain_mps)
-#rainin_term=anuga.operators.rate_operators.Rate_operator(domain,rate=myrain)
+rain_mps=90.0/(1000.*3600.) # 90mm/hr is meters-per-second
+def myrain(t):
+    return(rain_mps)
+#myindices=range(len(domain.centroid_coordinates[:,0]))
+#rainin_term=anuga.operators.rate_operators.Rate_operator(domain,rate=myrain,indices=myindices)
+rainin_term=anuga.operators.rate_operators.Rate_operator(domain,rate=myrain,polygon=project.bounding_polygon)
 
 ## Boundary conditions
@@ -119 +125 @@
 
 barrier()
+finalize()