Changeset 7723

Timestamp:

May 13, 2010, 9:51:30 AM (15 years ago)

Author:

sexton

Message:

updates

File:

• branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/process_gems_grids.py (modified) (11 diffs)

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 = 'case_1'
+event = 'gcom_60min'
 
 print "Processing event: ", event
@@ -46 +46 @@
 
 grid_file_names_temp = []
-for i in range(0, 240, 2):
+for i in range(1, 32, 1):
     grid_file_names_temp.append(grid_file_names[i])
 
@@ -56 +56 @@
 
 start_time = 0      # all times in seconds
-timestep = 1440 #720     # all times in seconds
+timestep = 3600 #1440 #720     # all times in seconds
 end_time =  start_time + (timestep*number_of_timesteps) # all times in seconds
 refzone = 50        # UTM zone of model
@@ -62 +62 @@
 
 elevation = []
-#stage = []
-#speed = []
-#theta = []
 
 elevation_file = open(elevation_file_name, 'rb')
@@ -95 +92 @@
 print 'No data value: ', no_data
 
+#------------------------------------------------------
+# Create arrays of elevation, stage, speed and theta
+#-------------------------------------------------------
+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)
 
-import pdb
-
 for j, f in enumerate(grid_file_names):
     print 'file: ', f
     
     gz_file = gzip.open(f, 'rb')
-    #lines = gz_file.readlines()
-    #gz_file.close()
-
-    #pdb.set_trace() 
 
     st = []
     sp = []
-    th = []
+    t = []
 
     for i, L in enumerate(gz_file):
@@ -123 +119 @@
             
         if i > (11 + 2*nrows):
-            th += L.strip().split()
+            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(th).astype('d'))
+    theta[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(t).astype('d'))
 
     gz_file.close()
-    
-#pdb.set_trace()
-
-#------------------------------------------------------
-# Create arrays of elevation, stage, speed and theta
-#-------------------------------------------------------
-print 'creating numpy arrays: elevation, stage, speed, theta'
-            
+
 elevation = numpy.array(elevation).astype('d')
 
-#stage_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float)
-#speed_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float)
-#theta_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float)
-
-#for i in range(number_of_timesteps):
-#    ii = nrows*i
-#    stage_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(stage[ii : ii + nrows]).astype('d'))
-#    speed_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(speed[ii:ii + nrows]).astype('d'))
-#    theta_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(theta[ii:ii + nrows]).astype('d'))
-
-#stage = stage_array#.ravel()
-#speed = speed_array#.ravel()
-#theta = theta_array#.ravel()
-    
 print 'Size of elevation array: ', elevation.size
 print 'Size of stage array: ', stage.size
@@ -162 +137 @@
 assert stage.size == number_of_timesteps * elevation.size
 
-d = numpy.empty(elevation.size, dtype='d')
 depth = numpy.empty(stage.size, dtype='d')
 number_of_points = ncols * nrows
@@ -171 +145 @@
 # a consistant dataset
 # ---------------------------------------------
-
-#for i in xrange(number_of_timesteps):
-#    j = number_of_points * i
-#    k = j + number_of_points
-    # no_value_mask = (stage[j:k] < -9000) + (speed[j:k] < -9000) + (theta[j:k] < -9000)
-#    no_value_index = numpy.where(((stage[j:k] < -100) + (speed[j:k] < 0) + (theta[j:k] < 0)) == True)[0]
-#    depth[j:k] = stage[j:k] - elevation
-#    numpy.put(stage[j:k], no_value_index, -9999)
-#    numpy.put(speed[j:k], no_value_index, -9999)
-#    numpy.put(theta[j:k], no_value_index, -9999)
-#    numpy.put(depth[j:k], no_value_index, 0)
-
+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)
+
 numpy.put(stage, no_value_index, -9999)
 numpy.put(speed, no_value_index, -9999)
 numpy.put(theta, no_value_index, -9999)
 numpy.put(depth, no_value_index, 0)
-
-# for i in xrange(number_of_timesteps):
-    # j = number_of_points * i
-    # k = j + number_of_points
-    # d = stage[j:k] - elevation     # Assumes elevations below sea level are negative
-    # depth[j:k] = [0 if x<-9000 else x for x in d]   # sets stage to zero for no-data values of stage
-    
-    # depth[j:k] = [x if x>=0 and elevation>0 else 0 for x in d] # This is to correct for GEMS data which assumes stage=0 onshore when dry but should = elevation
-    # for q in xrange(number_of_points):
-        # if elevation[q] >= 0 and stage[j+q]==0:
-            # depth[j+q] = 0
     
 # Taking absolute value is to account for -ve depths obtained when stage-elevation 
@@ -205 +159 @@
 momentum = numpy.absolute(depth * speed)
 
-#xmomentum = numpy.empty(momentum.size, dtype='d')
-#ymomentum = numpy.empty(momentum.size, dtype='d')
-
 print 'Calculating xmomentum and ymomentum'
-
-#for i, t in enumerate(theta):
-    
-    #mx = momentum[i]*math.sin(math.radians(t)) #Assuming t is in the "to" direction and 0 degrees is north
-    #my = momentum[i]*math.cos(math.radians(t))
 
 xmomentum = momentum*numpy.sin(numpy.radians(theta))
@@ -256 +202 @@
 
 assert len(numpy.where((numpy.sqrt(xmomentum**2 + ymomentum**2) - momentum) > 1.e-06)[0]) == 0
-
-#assert numpy.sqrt(xmomentum**2 + ymomentum**2) - momentum < 1.e-06
     
 x_origin_int = int(10000*x_origin)
@@ -280 +224 @@
 #for j in range(2):
 
-fid = NetCDFFile(os.path.join(event_folder, event + '_master_2_.sts'), 'wl')
+fid = NetCDFFile(os.path.join(event_folder, event + '_master_2_1.sts'), 'wl')
 fid.institution = 'Geoscience Australia'
 fid.description = "description"