Changeset 7692
 Timestamp:
 Apr 22, 2010, 3:07:19 PM (12 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/process_gems_grids.py
r7688 r7692 42 42 print "Elevation file name: ", elevation_file_name 43 43 44 grid_file_names = glob.glob(os.path.join(event_folder, 'gcom _txt19780404.0*.gz'))44 grid_file_names = glob.glob(os.path.join(event_folder, 'gcom*.gz')) 45 45 grid_file_names.sort() 46 46 47 #grid_file_names_temp = []48 # for i in range(0, 120, 4):49 #grid_file_names_temp.append(grid_file_names[i])50 51 #grid_file_names = grid_file_names_temp52 53 #number_of_timesteps = len(grid_file_names)47 grid_file_names_temp = [] 48 for i in range(0, 240, 2): 49 grid_file_names_temp.append(grid_file_names[i]) 50 51 grid_file_names = grid_file_names_temp 52 53 number_of_timesteps = len(grid_file_names) 54 54 55 55 print "Number of timesteps: ", number_of_timesteps 56 56 57 57 start_time = 0 # all times in seconds 58 timestep = 720 # all times in seconds58 timestep = 1440 #720 # all times in seconds 59 59 end_time = start_time + (timestep*number_of_timesteps) # all times in seconds 60 60 refzone = 50 # UTM zone of model … … 62 62 63 63 elevation = [] 64 stage = []65 speed = []66 theta = []64 #stage = [] 65 #speed = [] 66 #theta = [] 67 67 68 68 elevation_file = open(elevation_file_name, 'rb') … … 94 94 print 'Grid size: ', grid_size 95 95 print 'No data value: ', no_data 96 97 for f in grid_file_names: 96 97 stage = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float) 98 speed = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float) 99 theta = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float) 100 101 import pdb 102 103 for j, f in enumerate(grid_file_names): 98 104 print 'file: ', f 99 105 100 106 gz_file = gzip.open(f, 'rb') 101 lines = gz_file.readlines() 102 gz_file.close() 103 104 for i, L in enumerate(lines): 107 #lines = gz_file.readlines() 108 #gz_file.close() 109 110 #pdb.set_trace() 111 112 st = [] 113 sp = [] 114 th = [] 115 116 for i, L in enumerate(gz_file): 105 117 106 118 if i > 7 and i < (8 + nrows): # 7 refers to the number of rows of header info that we skip  always CHECK format 107 st age += [L.strip().split()]119 st += L.strip().split() 108 120 109 121 if i > (9 + nrows) and i < (10 + 2*nrows): 110 sp eed += [L.strip().split()]122 sp += L.strip().split() 111 123 112 124 if i > (11 + 2*nrows): 113 theta += [L.strip().split()] 125 th += L.strip().split() 126 127 stage[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(st).astype('d')) 128 speed[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(sp).astype('d')) 129 theta[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(th).astype('d')) 130 131 gz_file.close() 132 133 #pdb.set_trace() 114 134 115 135 # … … 120 140 elevation = numpy.array(elevation).astype('d') 121 141 122 stage_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float)123 speed_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float)124 theta_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float)125 126 for i in range(number_of_timesteps):127 ii = nrows*i128 stage_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(stage[ii : ii + nrows]).astype('d'))129 speed_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(speed[ii:ii + nrows]).astype('d'))130 theta_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(theta[ii:ii + nrows]).astype('d'))131 132 stage = stage_array.ravel()133 speed = speed_array.ravel()134 theta = theta_array.ravel()142 #stage_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float) 143 #speed_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float) 144 #theta_array = numpy.empty([number_of_timesteps*nrows, ncols], dtype=float) 145 146 #for i in range(number_of_timesteps): 147 # ii = nrows*i 148 # stage_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(stage[ii : ii + nrows]).astype('d')) 149 # speed_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(speed[ii:ii + nrows]).astype('d')) 150 # theta_array[ii:ii+nrows, :] = numpy.flipud(numpy.array(theta[ii:ii + nrows]).astype('d')) 151 152 #stage = stage_array#.ravel() 153 #speed = speed_array#.ravel() 154 #theta = theta_array#.ravel() 135 155 136 156 print 'Size of elevation array: ', elevation.size … … 152 172 #  153 173 154 for i in xrange(number_of_timesteps):155 j = number_of_points * i156 k = j + number_of_points174 #for i in xrange(number_of_timesteps): 175 # j = number_of_points * i 176 # k = j + number_of_points 157 177 # no_value_mask = (stage[j:k] < 9000) + (speed[j:k] < 9000) + (theta[j:k] < 9000) 158 no_value_index = numpy.where(((stage[j:k] < 100) + (speed[j:k] < 0) + (theta[j:k] < 0)) == True)[0] 159 depth[j:k] = stage[j:k]  elevation 160 numpy.put(stage[j:k], no_value_index, 9999) 161 numpy.put(speed[j:k], no_value_index, 9999) 162 numpy.put(theta[j:k], no_value_index, 9999) 163 numpy.put(depth[j:k], no_value_index, 0) 178 # no_value_index = numpy.where(((stage[j:k] < 100) + (speed[j:k] < 0) + (theta[j:k] < 0)) == True)[0] 179 # depth[j:k] = stage[j:k]  elevation 180 # numpy.put(stage[j:k], no_value_index, 9999) 181 # numpy.put(speed[j:k], no_value_index, 9999) 182 # numpy.put(theta[j:k], no_value_index, 9999) 183 # numpy.put(depth[j:k], no_value_index, 0) 184 185 no_value_index = numpy.where(((stage < 100) + (speed < 0) + (theta < 0)) == True)[0] 186 depth = stage  numpy.tile(elevation, number_of_timesteps) 187 numpy.put(stage, no_value_index, 9999) 188 numpy.put(speed, no_value_index, 9999) 189 numpy.put(theta, no_value_index, 9999) 190 numpy.put(depth, no_value_index, 0) 164 191 165 192 # for i in xrange(number_of_timesteps): … … 174 201 # depth[j+q] = 0 175 202 176 momentum = depth * speed 177 178 xmomentum = numpy.empty(momentum.size, dtype='d') 179 ymomentum = numpy.empty(momentum.size, dtype='d') 203 # Taking absolute value is to account for ve depths obtained when stageelevation 204 # is slightly ve  why I don't know, possbly a rounding error? 205 momentum = numpy.absolute(depth * speed) 206 207 #xmomentum = numpy.empty(momentum.size, dtype='d') 208 #ymomentum = numpy.empty(momentum.size, dtype='d') 180 209 181 210 print 'Calculating xmomentum and ymomentum' 182 211 183 for i, t in enumerate(theta): 184 185 mx = momentum[i]*math.sin(math.radians(t)) #Assuming t is in the "to" direction and 0 degrees is north 186 my = momentum[i]*math.cos(math.radians(t)) 187 188 if t > 0.0 and t < 90.0: 189 assert mx > 0 and my > 0 190 191 elif t > 90.0 and t < 180.0: 192 assert mx > 0 and my < 0 193 194 elif t > 180.0 and t < 270.0: 195 assert mx < 0 and my < 0 196 197 elif t > 270.0 and t < 360.0: 198 assert mx < 0 and my > 0 199 200 elif t == 0.0 or t == 360.0: 201 assert my > 0 202 203 elif t == 90.0: 204 assert mx > 0 205 206 elif t == 180.0: 207 assert my < 0 208 209 elif t == 270.0: 210 assert mx < 0 211 212 elif t == 9999: 213 mx = 0 214 my = 0 215 else: 216 print "Unexpected value of theta" 217 exit() 218 219 xmomentum[i] = mx 220 ymomentum[i] = my 221 222 assert math.sqrt(mx**2 + my**2)  momentum[i] < 1.e06 212 #for i, t in enumerate(theta): 213 214 #mx = momentum[i]*math.sin(math.radians(t)) #Assuming t is in the "to" direction and 0 degrees is north 215 #my = momentum[i]*math.cos(math.radians(t)) 216 217 xmomentum = momentum*numpy.sin(numpy.radians(theta)) 218 ymomentum = momentum*numpy.cos(numpy.radians(theta)) 219 220 numpy.put(xmomentum, no_value_index, 0) 221 numpy.put(ymomentum, no_value_index, 0) 222 223 # if t > 0.0 and t < 90.0: 224 # assert mx >= 0 and my >= 0 225 226 # elif t > 90.0 and t < 180.0: 227 # assert mx >= 0 and my <= 0 228 229 # elif t > 180.0 and t < 270.0: 230 # assert mx <= 0 and my <= 0 231 232 # elif t > 270.0 and t < 360.0: 233 # assert mx <= 0 and my >= 0 234 235 # elif t == 0.0 or t == 360.0: 236 # assert my >= 0 237 238 # elif t == 90.0: 239 # assert mx >= 0 240 241 # elif t == 180.0: 242 # assert my <= 0 243 244 # elif t == 270.0: 245 # assert mx <= 0 246 247 # elif t == 9999: 248 # mx = 0 249 # my = 0 250 # else: 251 # print "Unexpected value of theta" 252 # exit() 253 254 # xmomentum[i] = mx 255 # ymomentum[i] = my 256 257 assert len(numpy.where((numpy.sqrt(xmomentum**2 + ymomentum**2)  momentum) > 1.e06)[0]) == 0 258 259 #assert numpy.sqrt(xmomentum**2 + ymomentum**2)  momentum < 1.e06 223 260 224 261 x_origin_int = int(10000*x_origin) … … 241 278 print "Creating the STS NetCDF file" 242 279 243 fid = NetCDFFile(os.path.join(event_folder, event + '_test.sts'), 'w') 280 #for j in range(2): 281 282 fid = NetCDFFile(os.path.join(event_folder, event + '_master_2_.sts'), 'wl') 244 283 fid.institution = 'Geoscience Australia' 245 284 fid.description = "description" … … 270 309 fid.variables['ymomentum_range'][:] = numpy.array([1e+036, 1e+036]) 271 310 fid.variables['elevation'][:] = elevation 272 fid.variables['time'][:] = time 311 fid.variables['time'][:] = time#[j*60 : j*60 + 60] 273 312 274 313 s = fid.variables['stage'] … … 276 315 ym = fid.variables['ymomentum'] 277 316 317 #jj = j*number_of_points*(number_of_timesteps/2) 318 278 319 for i in xrange(number_of_timesteps): 279 ii = i*number_of_points 320 ii = i*number_of_points# + jj 280 321 s[i] = stage[ii : ii + number_of_points] 281 322 xm[i] = xmomentum[ii : ii + number_of_points] 282 323 ym[i] = ymomentum[ii : ii + number_of_points] 283 324 284 325 origin = Geo_reference(refzone, min(x), min(y)) # Check this section for inputs in eastings and northings  it works for long and lat 285 326 geo_ref = write_NetCDF_georeference(origin, fid) … … 289 330 290 331 fid.close() 332
Note: See TracChangeset
for help on using the changeset viewer.