Changeset 7683


Ignore:
Timestamp:
Apr 13, 2010, 2:34:11 PM (10 years ago)
Author:
fountain
Message:
 
File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/process_gems_grids.py

    r7681 r7683  
    4242print "Elevation file name: ", elevation_file_name
    4343
    44 grid_file_names = glob.glob(os.path.join(event_folder, 'gcom*.gz'))
     44grid_file_names = glob.glob(os.path.join(event_folder, 'gcom_txt19780403.18*.gz'))
    4545grid_file_names.sort()
    4646number_of_timesteps = len(grid_file_names)
    47 grid_size = 20  # cellsize from the GEMSURGE output grids
     47
     48print "Number of timesteps: ", number_of_timesteps
     49
    4850start_time = 0      # all times in seconds
    49 end_time = 61200    # end_time = start_time + timestep*number of grid_file_names. all times in seconds
    5051timestep = 720     # all times in seconds
     52end_time =  start_time + (timestep*number_of_timesteps) # all times in seconds
    5153refzone = 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
    5454event_sts = os.path.join(event_folder, event)
    5555
     
    6464
    6565# 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:
    6767for i, L in enumerate(lines):
    6868    if i == 0:
     
    7070    if i == 1:
    7171        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])
    7278    if i == 5:
    7379        no_data = int(float(L.strip().split()[1]))
     
    7783print 'Number or columns: ', ncols
    7884print 'Number of rows: ', nrows
     85print 'X origin: ', x_origin
     86print 'Y origin: ', y_origin
     87print 'Grid size: ', grid_size
     88print 'No data value: ', no_data
    7989       
    8090for f in grid_file_names:
     
    133143    k = j + number_of_points
    134144    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 no-data values of stage
     145    depth[j:k] = [0 if x<-500 else x for x in d]   # sets stage to zero for no-data values of stage
    136146   
    137147    # 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
     
    152162        mx = momentum[i]*math.sin(math.radians(t)) #Assuming t is in the "to" direction and 0 degrees is north
    153163        my = momentum[i]*math.cos(math.radians(t))
    154         assert mx > 0 and my > 0
     164        assert mx >= 0 and my >= 0
    155165    if t > 90.0 and t < 180.0:
    156166        mx = momentum[i]*math.sin(math.radians(t))
    157167        my = momentum[i]*math.cos(math.radians(t))
    158         assert mx > 0 and my < 0
     168        assert mx >= 0 and my <= 0
    159169    if t > 180.0 and t < 270.0:
    160170        mx = momentum[i]*math.cos(math.radians(t))
    161171        my = momentum[i]*math.sin(math.radians(t))
    162         assert mx < 0 and my < 0
     172        assert mx <= 0 and my <= 0
    163173    if t > 270.0 and t < 360.0:
    164174        mx = momentum[i]*math.sin(math.radians(t))
    165175        my = momentum[i]*math.cos(math.radians(t))
    166         assert mx < 0 and my > 0
     176        assert mx <= 0 and my >= 0
    167177    if t == 0.0 or t == 360.0:
    168178        mx = 0
     
    181191        my = 0
    182192        assert mx <= 0
     193    if t == -9999:
     194        mx = 0
     195        my = 0
    183196       
    184197    xmomentum[i] = mx
Note: See TracChangeset for help on using the changeset viewer.