Changeset 7794

Timestamp:

Jun 7, 2010, 9:38:23 AM (15 years ago)

Author:

sexton

Message:

updates to bunbury scripts

Location:

branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009

Files:

• build_gems_boundary.py (modified) (6 diffs)
• process_gems_grids.py (modified) (7 diffs)

branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/build_gems_boundary.py

@@ -13 +13 @@
 import sys
 # sys.path.append('W:\\sandpits\\lfountain\\anuga_work\\production\\mandurah_storm_surge_2009\\gems_comparison')
-sys.path.append('/nas/gemd/georisk_models/inundation/sandpits/lfountain/anuga_work/production/mandurah_storm_surge_2009/gems_comparison')
-import project_GEMS_grid
+#sys.path.append('/nas/gemd/georisk_models/inundation/sandpits/lfountain/anuga_work/production/mandurah_storm_surge_2009/gems_comparison')
+import project
 from Scientific.IO.NetCDF import NetCDFFile
 from anuga.coordinate_transforms.geo_reference import Geo_reference, write_NetCDF_georeference
@@ -23 +23 @@
 state = 'western_australia'
 scenario_folder = 'bunbury_storm_surge_scenario_2009'
-event = 'case_1' # Baseline TC Alby event
+event = '20100527_gcom_12min' # Baseline TC Alby event
 
 ENV_INUNDATIONHOME = 'INUNDATIONHOME'
@@ -35 +35 @@
 # Define boundary points to extract from master STS file
 #--------------------------------------------------------
-gems_boundary = numpy.array(read_polygon(project_GEMS_grid.gems_order))
+gems_boundary = numpy.array(read_polygon(project.gems_order))
+#gems_boundary = numpy.array(read_polygon(project.gauges))
 number_of_points = gems_boundary.shape[0]
+print 'hello', number_of_points
+print 'points', gems_boundary
 
 #---------------------------------------------------------
@@ -42 +45 @@
 #---------------------------------------------------------
 print "Reading master NetCDF file"
-infile = NetCDFFile(os.path.join(event_folder, event + '_master.sts'), 'r')
+infile = NetCDFFile(os.path.join(event_folder, event + '_master_2_1.sts'), 'rl') #rl for large netcdf files
 
 grid_size = infile.grid_size
@@ -69 +72 @@
 origin = numpy.array([x_origin[0], y_origin[0]])
 index_array_2D = ((gems_boundary - origin)/grid_size).astype(int)
-assert index_array_2D[:,0].min() >= 0 and index_array_2D[:,0].max() <= ncols
+assert index_array_2D[:,0].min() >= 0 and index_array_2D[:,0].max() <= ncols # These test that the boundary points lie within the data area
 assert index_array_2D[:,1].min() >= 0 and index_array_2D[:,1].max() <= nrows
 
@@ -135 +138 @@
     s[i] = stage[i,:]
     xm[i] = xmomentum[i,:]
-    ym[i] = ymomentum[i,:]
- 
+    ym[i] = ymomentum[i,:] 
+
+#---------------------------------------------------------------------------
+# Get timeseries values for wave height and components of momentum
+#---------------------------------------------------------------------------
+
+basename = 'sts_gauge_tide'
+for i in xrange(number_of_points):
+    out_file = os.path.join(event_folder,
+                            basename+'_'+str(i)+'.csv')
+    fid_sts = open(out_file, 'w')
+    fid_sts.write('time, stage, xmomentum, ymomentum, elevation \n')
+    for j in xrange(number_of_timesteps):
+        fid_sts.write('%.6f, %.6f, %.6f, %.6f, %.6f\n'
+                    % (time[j], stage[j,i], xmomentum[j,i], ymomentum[j,i], elevation[i]))
+
+    fid_sts.close() 
+
 origin = Geo_reference(refzone, x_origin, y_origin)
 geo_ref = write_NetCDF_georeference(origin, fid)
 
 fid.variables['x'][:] = x 
-fid.variables['y'][:] = y 
+fid.variables['y'][:] = y  
 
 fid.close()
+
+sys.exit()
+
+print 'get timeseries'
+#-------------------------------------------------------------------------------
+# Get gauges (timeseries of index points)
+#-------------------------------------------------------------------------------
+
+basename = 'sts_gauge'
+quantity_names = ['stage', 'xmomentum', 'ymomentum']
+quantities = {}
+for i, name in enumerate(quantity_names):
+    quantities[name] = fid.variables[name][:]
+
+
+for j in range(len(x)):
+    index = j # permutation[j]
+    stage = quantities['stage'][:,j]
+    #print 'hello', j, stage
+    xmomentum = quantities['xmomentum'][:,j]
+    ymomentum = quantities['ymomentum'][:,j]
+
+    out_file = os.path.join(event_folder,
+                            basename+'_'+str(index)+'.csv')
+    fid_sts = open(out_file, 'w')
+    fid_sts.write('time, stage, xmomentum, ymomentum \n')
+
+    #-----------------------------------------------------------------------
+    # End of the get gauges
+    #-----------------------------------------------------------------------
+    for k in range(len(time)-1):
+        fid_sts.write('%.6f, %.6f, %.6f, %.6f\n'
+                      % (time[k], stage[k], xmomentum[k], ymomentum[k]))
+
+    

branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/process_gems_grids.py

@@ -22 +22 @@
 state = 'western_australia'
 scenario_folder = 'bunbury_storm_surge_scenario_2009'
-event = 'gcom_60min'
+event = '20100527_gcom_12min'
 
 print "Processing event: ", event
@@ -45 +45 @@
 grid_file_names.sort()
 
-grid_file_names_temp = []
-for i in range(1, 32, 1):
-    grid_file_names_temp.append(grid_file_names[i])
-
-grid_file_names = grid_file_names_temp
+#grid_file_names_temp = []
+# this was to subsample the gz files
+#for i in range(1, 32, 1):
+#    grid_file_names_temp.append(grid_file_names[i])
+
+#grid_file_names = grid_file_names_temp
     
 number_of_timesteps = len(grid_file_names)
@@ -56 +57 @@
 
 start_time = 0      # all times in seconds
-timestep = 3600 #1440 #720     # all times in seconds
+timestep = 720     # all times in seconds
 end_time =  start_time + (timestep*number_of_timesteps) # all times in seconds
 refzone = 50        # UTM zone of model
@@ -93 +94 @@
 
 #------------------------------------------------------
-# Create arrays of elevation, stage, speed and theta
+# Create arrays of elevation, depth, current_x and current_y
 #-------------------------------------------------------
-print 'creating numpy arrays: elevation, stage, speed, theta'
-
-stage = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float)
-speed = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float)
-theta = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float)
+print 'creating numpy arrays: elevation, depth, current_x, current_y'
+
+depth = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float)
+current_x = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float)
+current_y = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float)
 
 for j, f in enumerate(grid_file_names):
@@ -106 +107 @@
     gz_file = gzip.open(f, 'rb')
 
-    st = []
-    sp = []
-    t = []
+    d = []
+    cx = []
+    cy = []
 
     for i, L in enumerate(gz_file):
         
-        if i > 7 and i < (8 + nrows):   #  7 refers to the number of rows of header info that we skip - always CHECK format
-            st += L.strip().split()
+        #if i > 7 and i < (8 + nrows): #  7 refers to the number of rows of header info that we skip - always CHECK format
+        #    print 'first block'
+
+        #if i > (9 + nrows) and i < (10 + 2*nrows):
+        #    print 'second block' 
+
+        #if i > (11 + 2*nrows):
+        #    print 'third block'
+
+        if i > 13 + 3*nrows and i < (14 + 4*nrows):   
+            d += L.strip().split()
             
-        if i > (9 + nrows) and i < (10 + 2*nrows):
-            sp += L.strip().split()
+        if i > (15 + 4*nrows) and i < (16 + 5*nrows):
+            cx += L.strip().split()
             
-        if i > (11 + 2*nrows):
-            t += L.strip().split()
-
-    stage[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(st).astype('d'))
-    speed[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(sp).astype('d'))
-    theta[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(t).astype('d'))
+        if i > (17 + 5*nrows):
+            cy += L.strip().split()
+
+    depth[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(d).astype('d'))
+    current_x[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(cx).astype('d'))
+    current_y[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(cy).astype('d'))
 
     gz_file.close()
@@ -130 +140 @@
 
 print 'Size of elevation array: ', elevation.size
-print 'Size of stage array: ', stage.size
-print 'Size of speed array: ', speed.size
-print 'Size of theta array: ', theta.size
-
-assert stage.size == speed.size == theta.size == ncols * nrows * number_of_timesteps
-assert stage.size == number_of_timesteps * elevation.size
-
-depth = numpy.empty(stage.size, dtype='d')
+print 'Size of depth array: ', depth.size
+print 'Size of current x array: ', current_x.size
+print 'Size of current y array: ', current_y.size
+
+assert depth.size == current_x.size == current_y.size == ncols * nrows * number_of_timesteps
+assert depth.size == number_of_timesteps * elevation.size
+
+stage = numpy.empty(depth.size, dtype='d')
 number_of_points = ncols * nrows
 
 # ---------------------------------------------
-# Create mask of no_data values across stage, speed and theta to ensure all three quantities
+# Create mask of no_data values across depth, current_x and current_y to ensure all three quantities
 # have no_data values in corresponding cells - this is a work-around until GEMS can produce a 
 # a consistant dataset
 # ---------------------------------------------
-import pdb
-pdb.set_trace()
-no_value_index = numpy.where(((stage < -100) + (speed < 0) + (theta < 0)) == True)[0]
-depth = stage - numpy.tile(elevation, number_of_timesteps)
+
+no_value_index = numpy.where(((depth < -9000) + (current_x < -9000) + (current_y < -9000)) == True)[0]
 
 numpy.put(stage, no_value_index, -9999)
-numpy.put(speed, no_value_index, -9999)
-numpy.put(theta, no_value_index, -9999)
+numpy.put(current_x, no_value_index, -9999)
+numpy.put(current_y, no_value_index, -9999)
 numpy.put(depth, no_value_index, 0)
     
-# Taking absolute value is to account for -ve depths obtained when stage-elevation 
+# Taking absolute value is to account for -ve depths obtained when depth-elevation 
 # is slightly -ve  - why I don't know, possbly a rounding error?
-momentum = numpy.absolute(depth * speed)
-
-print 'Calculating xmomentum and ymomentum'
-
-xmomentum = momentum*numpy.sin(numpy.radians(theta))
-ymomentum = momentum*numpy.cos(numpy.radians(theta))
+#momentum = numpy.absolute(depth * current_x)
+
+print 'Calculating stage, xmomentum and ymomentum'
+
+#stage = depth - numpy.tile(numpy.absolute(elevation), number_of_timesteps)
+stage = depth + numpy.tile(elevation, number_of_timesteps)
+xmomentum = current_x*depth #momentum*numpy.sin(numpy.radians(current_y))
+ymomentum = current_y*depth #momentum*numpy.cos(numpy.radians(current_y))
 
 numpy.put(xmomentum, no_value_index, 0)
 numpy.put(ymomentum, no_value_index, 0)
 
-#    if t > 0.0 and t < 90.0:
-#        assert mx >= 0 and my >= 0
-        
-#    elif t > 90.0 and t < 180.0:
-#        assert mx >= 0 and my <= 0
-        
-#    elif t > 180.0 and t < 270.0:
-#        assert mx <= 0  and my <= 0
-        
-#    elif t > 270.0 and t < 360.0:
-#        assert mx <= 0 and my >= 0
-        
-#    elif t == 0.0 or t == 360.0:
-#        assert my >= 0
-        
-#    elif t == 90.0:
-#        assert mx >= 0
-        
-#    elif t == 180.0:
-#        assert my <= 0
-        
-#    elif t == 270.0:
-#        assert mx <= 0
-        
-#    elif t == -9999:
-#        mx = 0
-#        my = 0
-#    else:
-#        print "Unexpected value of theta"
-#        exit()
-        
-#    xmomentum[i] = mx
-#    ymomentum[i] = my
-
-assert len(numpy.where((numpy.sqrt(xmomentum**2 + ymomentum**2) - momentum) > 1.e-06)[0]) == 0
+#assert len(numpy.where((numpy.sqrt(xmomentum**2 + ymomentum**2) - momentum) > 1.e-06)[0]) == 0
     
 x_origin_int = int(10000*x_origin)
@@ -215 +191 @@
 
 assert time.size == number_of_timesteps
-assert momentum.size == stage.size
+assert xmomentum.size == depth.size == ymomentum.size
 
 # -----------------------------