Changeset 7681

Timestamp:

Apr 12, 2010, 4:19:19 PM (13 years ago)

Author:

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) (10 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' # Baseline TC Alby event
+event = 'case_1'
+
+print "Processing event: ", event
 
 ENV_INUNDATIONHOME = 'INUNDATIONHOME'
@@ -36 +38 @@
 # Input files from GEMS
 #-----------------------
-elevation_file_name = os.path.join(event_folder, ***'topog_flip.asc'***) # Name of the vertically flipped elevation grid
-grid_file_names = glob.glob(os.path.join(event_folder, '*.gz'))
+elevation_file_name = os.path.join(event_folder, 'buntopog_20m_flip.asc')  # Name of the vertically flipped elevation grid
+
+print "Elevation file name: ", elevation_file_name
+
+grid_file_names = glob.glob(os.path.join(event_folder, 'gcom*.gz'))
 grid_file_names.sort()
 number_of_timesteps = len(grid_file_names)
-grid_size = 20 #0.0005 # cellsize from the GEMSURGE output grids
-start_time = 0  # all times in seconds
-end_time = 142560 #88200 # all times in seconds
-timestep = 720 #1800 # all times in seconds
-refzone = 50 # UTM zone
-x_origin = ***#115.550 #364221.03   # xllcorner from GEMSURGE output grids
-y_origin = ***#-32.790 #6371062.71  # yllcorner from GEMSURGE output grids
+grid_size = 20  # cellsize from the GEMSURGE output grids
+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
+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)
-
-assert number_of_timesteps * timestep == end_time - start_time
 
 elevation = []
@@ -61 +64 @@
 
 # Strip off the ASCII header and also read ncols, nrows, 
-# x_origin, y_origin, grid_size, no_data value, and read in elevation grid:
+# no_data value, and read in elevation grid:
 for i, L in enumerate(lines):
     if i == 0:
@@ -67 +70 @@
     if i == 1:
         nrows = int(L.strip().split()[1])
-    # if i == 2:
-        # x_origin = int(float(L.strip().split()[1])) # this assumes an integer xllcorner
-    # if i == 3:
-        # y_origin = int(float(L.strip().split()[1])) # this assumes an integer yllcorner
-    # if i == 4:
-        # grid_size = int(float(L.strip().split()[1]))
-    # if i == 5:
-        # no_data = int(float(L.strip().split()[1]))
+    if i == 5:
+        no_data = int(float(L.strip().split()[1]))
     if i > 5:
         elevation+= L.strip().split()
@@ -80 +77 @@
 print 'Number or columns: ', ncols
 print 'Number of rows: ', nrows
-print 'Grid size: ', grid_size
-print 'X origin: ', x_origin
-print 'Y origin: ', y_origin
-  
+        
 for f in grid_file_names:
     print 'file: ', f
@@ -93 +87 @@
     for i, L in enumerate(lines):
         
-        if i > 7 and i < (8 + nrows):        # Number 7 refers to the number of rows of header info that we skip - CHECK format
+        if i > 7 and i < (8 + nrows):   #  7 refers to the number of rows of header info that we skip - always CHECK format
             stage += [L.strip().split()]
             
@@ -131 +125 @@
 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
@@ -137 +132 @@
     j = number_of_points * i
     k = j + number_of_points
-    depth[j:k] = stage[j:k] - elevation # Check GEMS elevations below sea level are negative
-
+    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] = [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
+    
 momentum = depth * speed
 
@@ -153 +154 @@
         assert mx > 0 and my > 0
     if t > 90.0 and t < 180.0:
-        mx = momentum[i]*math.sin(math.radians(180 - t))
-        my = momentum[i]*math.cos(math.radians(180 - t))
+        mx = momentum[i]*math.sin(math.radians(t))
+        my = momentum[i]*math.cos(math.radians(t))
         assert mx > 0 and my < 0
     if t > 180.0 and t < 270.0:
-        mx = momentum[i]*math.cos(math.radians(270 - t))
-        my = momentum[i]*math.sin(math.radians(270 - t))
+        mx = momentum[i]*math.cos(math.radians(t))
+        my = momentum[i]*math.sin(math.radians(t))
         assert mx < 0 and my < 0
     if t > 270.0 and t < 360.0:
-        mx = momentum[i]*math.sin(math.radians(360 - t))
-        my = momentum[i]*math.cos(math.radians(360 - t))
+        mx = momentum[i]*math.sin(math.radians(t))
+        my = momentum[i]*math.cos(math.radians(t))
         assert mx < 0 and my > 0
     if t == 0.0 or t == 360.0:
         mx = 0
         my = momentum[i]
-        assert my > 0
+        assert my >= 0
     if t == 90.0:
         mx = momentum[i]
         my = 0
-        assert mx > 0
+        assert mx >= 0
     if t == 180.0:
         mx = 0
-        my = momentum[i]
-        assert my < 0
+        my = -momentum[i]
+        assert my <= 0
     if t == 270.0:
-        mx = momentum[i]
+        mx = -momentum[i]
         my = 0
-        assert mx < 0
+        assert mx <= 0
         
     xmomentum[i] = mx
     ymomentum[i] = my
 
-    assert mx[i]^2 + my[i]^2 == momentum[i]
+    assert math.sqrt(mx**2 + my**2) - momentum[i] < 1.e-06
     
 x_origin_int = int(10000*x_origin)
@@ -205 +206 @@
 print "Creating the STS NetCDF file"
 
-fid = NetCDFFile(os.path.join(event_folder, event + '_master.sts'), 'w')
+fid = NetCDFFile(os.path.join(event_folder, event + '_test.sts'), 'w')
 fid.institution = 'Geoscience Australia'
 fid.description = "description"