Changeset 7683

Timestamp:

Apr 13, 2010, 2:34:11 PM (15 years ago)

Author:

fountain

Message:

 

File:

• branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/process_gems_grids.py (modified) (7 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*.gz'))
+grid_file_names = glob.glob(os.path.join(event_folder, 'gcom_txt19780403.18*.gz'))
 grid_file_names.sort()
 number_of_timesteps = len(grid_file_names)
-grid_size = 20  # cellsize from the GEMSURGE output grids
+
+print "Number of timesteps: ", number_of_timesteps
+
 start_time = 0      # all times in seconds
-end_time = 61200    # end_time = start_time + timestep*number of grid_file_names. 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
-x_origin = 368600 #364221.03   # xllcorner from GEMSURGE output grids
-y_origin = 6304600 #6371062.71  # yllcorner from GEMSURGE output grids
 event_sts = os.path.join(event_folder, event)
 
@@ -64 +64 @@
 
 # Strip off the ASCII header and also read ncols, nrows, 
-# no_data value, and read in elevation grid:
+# x_origin, y_origin, grid_size, no_data value, and read in elevation grid:
 for i, L in enumerate(lines):
     if i == 0:
@@ -70 +70 @@
     if i == 1:
         nrows = int(L.strip().split()[1])
+    if i == 2:
+        x_origin = int(L.strip().split()[1])
+    if i == 3:
+        y_origin = int(L.strip().split()[1])
+    if i == 4:
+        grid_size = int(L.strip().split()[1])
     if i == 5:
         no_data = int(float(L.strip().split()[1]))
@@ -77 +83 @@
 print 'Number or columns: ', ncols
 print 'Number of rows: ', nrows
+print 'X origin: ', x_origin
+print 'Y origin: ', y_origin
+print 'Grid size: ', grid_size
+print 'No data value: ', no_data
         
 for f in grid_file_names:
@@ -133 +143 @@
     k = j + number_of_points
     d = stage[j:k] - elevation     # Assumes elevations below sea level are negative
-    depth[j:k] = [0 if x<-9900 else x for x in d]   # sets stage to zero for no-data values of stage
+    depth[j:k] = [0 if x<-500 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
@@ -152 +162 @@
         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
+        assert mx >= 0 and my >= 0
     if 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
+        assert mx >= 0 and my <= 0
     if 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
+        assert mx <= 0 and my <= 0
     if 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
+        assert mx <= 0 and my >= 0
     if t == 0.0 or t == 360.0:
         mx = 0
@@ -181 +191 @@
         my = 0
         assert mx <= 0
+    if t == -9999:
+        mx = 0
+        my = 0
         
     xmomentum[i] = mx