Changeset 7688

Timestamp:

Apr 19, 2010, 3:11:12 PM (15 years ago)

Author:

Leharne Fountain

Message:

updates to bunbury storm surge model

File:

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

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

@@ -42 +42 @@
 print "Elevation file name: ", elevation_file_name
 
-grid_file_names = glob.glob(os.path.join(event_folder, 'gcom_txt19780403.18*.gz'))
+grid_file_names = glob.glob(os.path.join(event_folder, 'gcom_txt19780404.0*.gz'))
 grid_file_names.sort()
-number_of_timesteps = len(grid_file_names)
+
+# grid_file_names_temp = []
+# for i in range(0, 120, 4):
+    # grid_file_names_temp.append(grid_file_names[i])
+
+# grid_file_names = grid_file_names_temp
+    
+# number_of_timesteps = len(grid_file_names)
 
 print "Number of timesteps: ", number_of_timesteps
@@ -139 +146 @@
 number_of_points = ncols * nrows
 
-# insert nodata test here
+# ---------------------------------------------
+# Create mask of no_data values across stage, speed and theta 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
+# ---------------------------------------------
 
 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
+    # 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)
+
+# 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
@@ -160 +182 @@
 
 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))
 
     if t > 0.0 and t < 90.0:
-        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))
         assert mx > 0 and my > 0
+        
     elif t > 90.0 and t < 180.0:
-        mx = momentum[i]*math.sin(math.radians(t))
-        my = momentum[i]*math.cos(math.radians(t))
         assert mx > 0 and my < 0
+        
     elif t > 180.0 and t < 270.0:
-        mx = momentum[i]*math.cos(math.radians(t))
-        my = momentum[i]*math.sin(math.radians(t))
         assert mx < 0 and my < 0
+        
     elif t > 270.0 and t < 360.0:
-        mx = momentum[i]*math.sin(math.radians(t))
-        my = momentum[i]*math.cos(math.radians(t))
         assert mx < 0 and my > 0
+        
     elif t == 0.0 or t == 360.0:
-        mx = 0
-        my = momentum[i]
-        assert my >= 0
+        assert my > 0
+        
     elif t == 90.0:
-        mx = momentum[i]
-        my = 0
-        assert mx >= 0
+        assert mx > 0
+        
     elif t == 180.0:
-        mx = 0
-        my = -momentum[i]
-        assert my <= 0
+        assert my < 0
+        
     elif t == 270.0:
-        mx = -momentum[i]
-        my = 0
-        assert mx <= 0
+        assert mx < 0
+        
     elif t == -9999:
         mx = 0
@@ -219 +236 @@
 assert momentum.size == stage.size
 
-#-----------------------------
+# -----------------------------
 # Create the STS file
-#-----------------------------
+# -----------------------------
 print "Creating the STS NetCDF file"