Changeset 7681 for branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/process_gems_grids.py
 Timestamp:
 Apr 12, 2010, 4:19:19 PM (13 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/process_gems_grids.py
r7678 r7681 22 22 state = 'western_australia' 23 23 scenario_folder = 'bunbury_storm_surge_scenario_2009' 24 event = 'case_1' # Baseline TC Alby event 24 event = 'case_1' 25 26 print "Processing event: ", event 25 27 26 28 ENV_INUNDATIONHOME = 'INUNDATIONHOME' … … 36 38 # Input files from GEMS 37 39 # 38 elevation_file_name = os.path.join(event_folder, ***'topog_flip.asc'***) # Name of the vertically flipped elevation grid 39 grid_file_names = glob.glob(os.path.join(event_folder, '*.gz')) 40 elevation_file_name = os.path.join(event_folder, 'buntopog_20m_flip.asc') # Name of the vertically flipped elevation grid 41 42 print "Elevation file name: ", elevation_file_name 43 44 grid_file_names = glob.glob(os.path.join(event_folder, 'gcom*.gz')) 40 45 grid_file_names.sort() 41 46 number_of_timesteps = len(grid_file_names) 42 grid_size = 20 #0.0005# cellsize from the GEMSURGE output grids43 start_time = 0 # all times in seconds44 end_time = 142560 #88200 #all times in seconds45 timestep = 720 #1800# all times in seconds46 refzone = 50 # UTM zone47 x_origin = ***#115.550 #364221.03 # xllcorner from GEMSURGE output grids48 y_origin = ***#32.790 #6371062.71 # yllcorner from GEMSURGE output grids47 grid_size = 20 # cellsize from the GEMSURGE output grids 48 start_time = 0 # all times in seconds 49 end_time = 61200 # end_time = start_time + timestep*number of grid_file_names. all times in seconds 50 timestep = 720 # all times in seconds 51 refzone = 50 # UTM zone of model 52 x_origin = 368600 #364221.03 # xllcorner from GEMSURGE output grids 53 y_origin = 6304600 #6371062.71 # yllcorner from GEMSURGE output grids 49 54 event_sts = os.path.join(event_folder, event) 50 51 assert number_of_timesteps * timestep == end_time  start_time52 55 53 56 elevation = [] … … 61 64 62 65 # Strip off the ASCII header and also read ncols, nrows, 63 # x_origin, y_origin, grid_size,no_data value, and read in elevation grid:66 # no_data value, and read in elevation grid: 64 67 for i, L in enumerate(lines): 65 68 if i == 0: … … 67 70 if i == 1: 68 71 nrows = int(L.strip().split()[1]) 69 # if i == 2: 70 # x_origin = int(float(L.strip().split()[1])) # this assumes an integer xllcorner 71 # if i == 3: 72 # y_origin = int(float(L.strip().split()[1])) # this assumes an integer yllcorner 73 # if i == 4: 74 # grid_size = int(float(L.strip().split()[1])) 75 # if i == 5: 76 # no_data = int(float(L.strip().split()[1])) 72 if i == 5: 73 no_data = int(float(L.strip().split()[1])) 77 74 if i > 5: 78 75 elevation+= L.strip().split() … … 80 77 print 'Number or columns: ', ncols 81 78 print 'Number of rows: ', nrows 82 print 'Grid size: ', grid_size 83 print 'X origin: ', x_origin 84 print 'Y origin: ', y_origin 85 79 86 80 for f in grid_file_names: 87 81 print 'file: ', f … … 93 87 for i, L in enumerate(lines): 94 88 95 if i > 7 and i < (8 + nrows): # Number 7 refers to the number of rows of header info that we skip CHECK format89 if i > 7 and i < (8 + nrows): # 7 refers to the number of rows of header info that we skip  always CHECK format 96 90 stage += [L.strip().split()] 97 91 … … 131 125 assert stage.size == number_of_timesteps * elevation.size 132 126 127 d = numpy.empty(elevation.size, dtype='d') 133 128 depth = numpy.empty(stage.size, dtype='d') 134 129 number_of_points = ncols * nrows … … 137 132 j = number_of_points * i 138 133 k = j + number_of_points 139 depth[j:k] = stage[j:k]  elevation # Check GEMS elevations below sea level are negative 140 134 d = stage[j:k]  elevation # Assumes elevations below sea level are negative 135 depth[j:k] = [0 if x<9900 else x for x in d] # sets stage to zero for nodata values of stage 136 137 # 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 138 # for q in xrange(number_of_points): 139 # if elevation[q] >= 0 and stage[j+q]==0: 140 # depth[j+q] = 0 141 141 142 momentum = depth * speed 142 143 … … 153 154 assert mx > 0 and my > 0 154 155 if t > 90.0 and t < 180.0: 155 mx = momentum[i]*math.sin(math.radians( 180 t))156 my = momentum[i]*math.cos(math.radians( 180 t))156 mx = momentum[i]*math.sin(math.radians(t)) 157 my = momentum[i]*math.cos(math.radians(t)) 157 158 assert mx > 0 and my < 0 158 159 if t > 180.0 and t < 270.0: 159 mx = momentum[i]*math.cos(math.radians( 270 t))160 my = momentum[i]*math.sin(math.radians( 270 t))160 mx = momentum[i]*math.cos(math.radians(t)) 161 my = momentum[i]*math.sin(math.radians(t)) 161 162 assert mx < 0 and my < 0 162 163 if t > 270.0 and t < 360.0: 163 mx = momentum[i]*math.sin(math.radians( 360 t))164 my = momentum[i]*math.cos(math.radians( 360 t))164 mx = momentum[i]*math.sin(math.radians(t)) 165 my = momentum[i]*math.cos(math.radians(t)) 165 166 assert mx < 0 and my > 0 166 167 if t == 0.0 or t == 360.0: 167 168 mx = 0 168 169 my = momentum[i] 169 assert my > 0170 assert my >= 0 170 171 if t == 90.0: 171 172 mx = momentum[i] 172 173 my = 0 173 assert mx > 0174 assert mx >= 0 174 175 if t == 180.0: 175 176 mx = 0 176 my = momentum[i]177 assert my < 0177 my = momentum[i] 178 assert my <= 0 178 179 if t == 270.0: 179 mx = momentum[i]180 mx = momentum[i] 180 181 my = 0 181 assert mx < 0182 assert mx <= 0 182 183 183 184 xmomentum[i] = mx 184 185 ymomentum[i] = my 185 186 186 assert m x[i]^2 + my[i]^2 == momentum[i]187 assert math.sqrt(mx**2 + my**2)  momentum[i] < 1.e06 187 188 188 189 x_origin_int = int(10000*x_origin) … … 205 206 print "Creating the STS NetCDF file" 206 207 207 fid = NetCDFFile(os.path.join(event_folder, event + '_ master.sts'), 'w')208 fid = NetCDFFile(os.path.join(event_folder, event + '_test.sts'), 'w') 208 209 fid.institution = 'Geoscience Australia' 209 210 fid.description = "description"
Note: See TracChangeset
for help on using the changeset viewer.