Changeset 7681


Ignore:
Timestamp:
Apr 12, 2010, 4:19:19 PM (10 years ago)
Author:
fountain
Message:

updates to bunbury storm surge model

File:
1 edited

Legend:

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

    r7678 r7681  
    2222state = 'western_australia'
    2323scenario_folder = 'bunbury_storm_surge_scenario_2009'
    24 event = 'case_1' # Baseline TC Alby event
     24event = 'case_1'
     25
     26print "Processing event: ", event
    2527
    2628ENV_INUNDATIONHOME = 'INUNDATIONHOME'
     
    3638# Input files from GEMS
    3739#-----------------------
    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'))
     40elevation_file_name = os.path.join(event_folder, 'buntopog_20m_flip.asc')  # Name of the vertically flipped elevation grid
     41
     42print "Elevation file name: ", elevation_file_name
     43
     44grid_file_names = glob.glob(os.path.join(event_folder, 'gcom*.gz'))
    4045grid_file_names.sort()
    4146number_of_timesteps = len(grid_file_names)
    42 grid_size = 20 #0.0005 # cellsize from the GEMSURGE output grids
    43 start_time = 0  # all times in seconds
    44 end_time = 142560 #88200 # all times in seconds
    45 timestep = 720 #1800 # all times in seconds
    46 refzone = 50 # UTM zone
    47 x_origin = ***#115.550 #364221.03   # xllcorner from GEMSURGE output grids
    48 y_origin = ***#-32.790 #6371062.71  # yllcorner from GEMSURGE output grids
     47grid_size = 20 # cellsize from the GEMSURGE output grids
     48start_time = 0      # all times in seconds
     49end_time = 61200    # end_time = start_time + timestep*number of grid_file_names. all times in seconds
     50timestep = 720     # all times in seconds
     51refzone = 50        # UTM zone of model
     52x_origin = 368600 #364221.03   # xllcorner from GEMSURGE output grids
     53y_origin = 6304600 #6371062.71  # yllcorner from GEMSURGE output grids
    4954event_sts = os.path.join(event_folder, event)
    50 
    51 assert number_of_timesteps * timestep == end_time - start_time
    5255
    5356elevation = []
     
    6164
    6265# 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:
    6467for i, L in enumerate(lines):
    6568    if i == 0:
     
    6770    if i == 1:
    6871        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]))
    7774    if i > 5:
    7875        elevation+= L.strip().split()
     
    8077print 'Number or columns: ', ncols
    8178print 'Number of rows: ', nrows
    82 print 'Grid size: ', grid_size
    83 print 'X origin: ', x_origin
    84 print 'Y origin: ', y_origin
    85  
     79       
    8680for f in grid_file_names:
    8781    print 'file: ', f
     
    9387    for i, L in enumerate(lines):
    9488       
    95         if i > 7 and i < (8 + nrows):        # Number 7 refers to the number of rows of header info that we skip - CHECK format
     89        if i > 7 and i < (8 + nrows):   #  7 refers to the number of rows of header info that we skip - always CHECK format
    9690            stage += [L.strip().split()]
    9791           
     
    131125assert stage.size == number_of_timesteps * elevation.size
    132126
     127d = numpy.empty(elevation.size, dtype='d')
    133128depth = numpy.empty(stage.size, dtype='d')
    134129number_of_points = ncols * nrows
     
    137132    j = number_of_points * i
    138133    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 no-data 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   
    141142momentum = depth * speed
    142143
     
    153154        assert mx > 0 and my > 0
    154155    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))
    157158        assert mx > 0 and my < 0
    158159    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))
    161162        assert mx < 0 and my < 0
    162163    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))
    165166        assert mx < 0 and my > 0
    166167    if t == 0.0 or t == 360.0:
    167168        mx = 0
    168169        my = momentum[i]
    169         assert my > 0
     170        assert my >= 0
    170171    if t == 90.0:
    171172        mx = momentum[i]
    172173        my = 0
    173         assert mx > 0
     174        assert mx >= 0
    174175    if t == 180.0:
    175176        mx = 0
    176         my = momentum[i]
    177         assert my < 0
     177        my = -momentum[i]
     178        assert my <= 0
    178179    if t == 270.0:
    179         mx = momentum[i]
     180        mx = -momentum[i]
    180181        my = 0
    181         assert mx < 0
     182        assert mx <= 0
    182183       
    183184    xmomentum[i] = mx
    184185    ymomentum[i] = my
    185186
    186     assert mx[i]^2 + my[i]^2 == momentum[i]
     187    assert math.sqrt(mx**2 + my**2) - momentum[i] < 1.e-06
    187188   
    188189x_origin_int = int(10000*x_origin)
     
    205206print "Creating the STS NetCDF file"
    206207
    207 fid = NetCDFFile(os.path.join(event_folder, event + '_master.sts'), 'w')
     208fid = NetCDFFile(os.path.join(event_folder, event + '_test.sts'), 'w')
    208209fid.institution = 'Geoscience Australia'
    209210fid.description = "description"
Note: See TracChangeset for help on using the changeset viewer.