Changeset 7683
 Timestamp:
 Apr 13, 2010, 2:34:11 PM (15 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/process_gems_grids.py
r7681 r7683 42 42 print "Elevation file name: ", elevation_file_name 43 43 44 grid_file_names = glob.glob(os.path.join(event_folder, 'gcom *.gz'))44 grid_file_names = glob.glob(os.path.join(event_folder, 'gcom_txt19780403.18*.gz')) 45 45 grid_file_names.sort() 46 46 number_of_timesteps = len(grid_file_names) 47 grid_size = 20 # cellsize from the GEMSURGE output grids 47 48 print "Number of timesteps: ", number_of_timesteps 49 48 50 start_time = 0 # all times in seconds 49 end_time = 61200 # end_time = start_time + timestep*number of grid_file_names. all times in seconds50 51 timestep = 720 # all times in seconds 52 end_time = start_time + (timestep*number_of_timesteps) # all times in seconds 51 53 refzone = 50 # UTM zone of model 52 x_origin = 368600 #364221.03 # xllcorner from GEMSURGE output grids53 y_origin = 6304600 #6371062.71 # yllcorner from GEMSURGE output grids54 54 event_sts = os.path.join(event_folder, event) 55 55 … … 64 64 65 65 # Strip off the ASCII header and also read ncols, nrows, 66 # no_data value, and read in elevation grid:66 # x_origin, y_origin, grid_size, no_data value, and read in elevation grid: 67 67 for i, L in enumerate(lines): 68 68 if i == 0: … … 70 70 if i == 1: 71 71 nrows = int(L.strip().split()[1]) 72 if i == 2: 73 x_origin = int(L.strip().split()[1]) 74 if i == 3: 75 y_origin = int(L.strip().split()[1]) 76 if i == 4: 77 grid_size = int(L.strip().split()[1]) 72 78 if i == 5: 73 79 no_data = int(float(L.strip().split()[1])) … … 77 83 print 'Number or columns: ', ncols 78 84 print 'Number of rows: ', nrows 85 print 'X origin: ', x_origin 86 print 'Y origin: ', y_origin 87 print 'Grid size: ', grid_size 88 print 'No data value: ', no_data 79 89 80 90 for f in grid_file_names: … … 133 143 k = j + number_of_points 134 144 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 stage145 depth[j:k] = [0 if x<500 else x for x in d] # sets stage to zero for nodata values of stage 136 146 137 147 # 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 153 163 my = momentum[i]*math.cos(math.radians(t)) 154 assert mx > 0 and my >0164 assert mx >= 0 and my >= 0 155 165 if t > 90.0 and t < 180.0: 156 166 mx = momentum[i]*math.sin(math.radians(t)) 157 167 my = momentum[i]*math.cos(math.radians(t)) 158 assert mx > 0 and my <0168 assert mx >= 0 and my <= 0 159 169 if t > 180.0 and t < 270.0: 160 170 mx = momentum[i]*math.cos(math.radians(t)) 161 171 my = momentum[i]*math.sin(math.radians(t)) 162 assert mx < 0 and my <0172 assert mx <= 0 and my <= 0 163 173 if t > 270.0 and t < 360.0: 164 174 mx = momentum[i]*math.sin(math.radians(t)) 165 175 my = momentum[i]*math.cos(math.radians(t)) 166 assert mx < 0 and my >0176 assert mx <= 0 and my >= 0 167 177 if t == 0.0 or t == 360.0: 168 178 mx = 0 … … 181 191 my = 0 182 192 assert mx <= 0 193 if t == 9999: 194 mx = 0 195 my = 0 183 196 184 197 xmomentum[i] = mx
Note: See TracChangeset
for help on using the changeset viewer.