Ignore:
Timestamp:
Jun 7, 2010, 9:38:23 AM (13 years ago)
Author:
sexton
Message:

updates to bunbury scripts

File:
1 edited

Legend:

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

    r7723 r7794  
    2222state = 'western_australia'
    2323scenario_folder = 'bunbury_storm_surge_scenario_2009'
    24 event = 'gcom_60min'
     24event = '20100527_gcom_12min'
    2525
    2626print "Processing event: ", event
     
    4545grid_file_names.sort()
    4646
    47 grid_file_names_temp = []
    48 for i in range(1, 32, 1):
    49     grid_file_names_temp.append(grid_file_names[i])
    50 
    51 grid_file_names = grid_file_names_temp
     47#grid_file_names_temp = []
     48# this was to subsample the gz files
     49#for i in range(1, 32, 1):
     50#    grid_file_names_temp.append(grid_file_names[i])
     51
     52#grid_file_names = grid_file_names_temp
    5253   
    5354number_of_timesteps = len(grid_file_names)
     
    5657
    5758start_time = 0      # all times in seconds
    58 timestep = 3600 #1440 #720     # all times in seconds
     59timestep = 720     # all times in seconds
    5960end_time =  start_time + (timestep*number_of_timesteps) # all times in seconds
    6061refzone = 50        # UTM zone of model
     
    9394
    9495#------------------------------------------------------
    95 # Create arrays of elevation, stage, speed and theta
     96# Create arrays of elevation, depth, current_x and current_y
    9697#-------------------------------------------------------
    97 print 'creating numpy arrays: elevation, stage, speed, theta'
    98 
    99 stage = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float)
    100 speed = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float)
    101 theta = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float)
     98print 'creating numpy arrays: elevation, depth, current_x, current_y'
     99
     100depth = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float)
     101current_x = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float)
     102current_y = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float)
    102103
    103104for j, f in enumerate(grid_file_names):
     
    106107    gz_file = gzip.open(f, 'rb')
    107108
    108     st = []
    109     sp = []
    110     t = []
     109    d = []
     110    cx = []
     111    cy = []
    111112
    112113    for i, L in enumerate(gz_file):
    113114       
    114         if i > 7 and i < (8 + nrows):   #  7 refers to the number of rows of header info that we skip - always CHECK format
    115             st += L.strip().split()
     115        #if i > 7 and i < (8 + nrows): #  7 refers to the number of rows of header info that we skip - always CHECK format
     116        #    print 'first block'
     117
     118        #if i > (9 + nrows) and i < (10 + 2*nrows):
     119        #    print 'second block'
     120
     121        #if i > (11 + 2*nrows):
     122        #    print 'third block'
     123
     124        if i > 13 + 3*nrows and i < (14 + 4*nrows):   
     125            d += L.strip().split()
    116126           
    117         if i > (9 + nrows) and i < (10 + 2*nrows):
    118             sp += L.strip().split()
     127        if i > (15 + 4*nrows) and i < (16 + 5*nrows):
     128            cx += L.strip().split()
    119129           
    120         if i > (11 + 2*nrows):
    121             t += L.strip().split()
    122 
    123     stage[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(st).astype('d'))
    124     speed[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(sp).astype('d'))
    125     theta[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(t).astype('d'))
     130        if i > (17 + 5*nrows):
     131            cy += L.strip().split()
     132
     133    depth[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(d).astype('d'))
     134    current_x[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(cx).astype('d'))
     135    current_y[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(cy).astype('d'))
    126136
    127137    gz_file.close()
     
    130140
    131141print 'Size of elevation array: ', elevation.size
    132 print 'Size of stage array: ', stage.size
    133 print 'Size of speed array: ', speed.size
    134 print 'Size of theta array: ', theta.size
    135 
    136 assert stage.size == speed.size == theta.size == ncols * nrows * number_of_timesteps
    137 assert stage.size == number_of_timesteps * elevation.size
    138 
    139 depth = numpy.empty(stage.size, dtype='d')
     142print 'Size of depth array: ', depth.size
     143print 'Size of current x array: ', current_x.size
     144print 'Size of current y array: ', current_y.size
     145
     146assert depth.size == current_x.size == current_y.size == ncols * nrows * number_of_timesteps
     147assert depth.size == number_of_timesteps * elevation.size
     148
     149stage = numpy.empty(depth.size, dtype='d')
    140150number_of_points = ncols * nrows
    141151
    142152# ---------------------------------------------
    143 # Create mask of no_data values across stage, speed and theta to ensure all three quantities
     153# Create mask of no_data values across depth, current_x and current_y to ensure all three quantities
    144154# have no_data values in corresponding cells - this is a work-around until GEMS can produce a
    145155# a consistant dataset
    146156# ---------------------------------------------
    147 import pdb
    148 pdb.set_trace()
    149 no_value_index = numpy.where(((stage < -100) + (speed < 0) + (theta < 0)) == True)[0]
    150 depth = stage - numpy.tile(elevation, number_of_timesteps)
     157
     158no_value_index = numpy.where(((depth < -9000) + (current_x < -9000) + (current_y < -9000)) == True)[0]
    151159
    152160numpy.put(stage, no_value_index, -9999)
    153 numpy.put(speed, no_value_index, -9999)
    154 numpy.put(theta, no_value_index, -9999)
     161numpy.put(current_x, no_value_index, -9999)
     162numpy.put(current_y, no_value_index, -9999)
    155163numpy.put(depth, no_value_index, 0)
    156164   
    157 # Taking absolute value is to account for -ve depths obtained when stage-elevation
     165# Taking absolute value is to account for -ve depths obtained when depth-elevation
    158166# is slightly -ve  - why I don't know, possbly a rounding error?
    159 momentum = numpy.absolute(depth * speed)
    160 
    161 print 'Calculating xmomentum and ymomentum'
    162 
    163 xmomentum = momentum*numpy.sin(numpy.radians(theta))
    164 ymomentum = momentum*numpy.cos(numpy.radians(theta))
     167#momentum = numpy.absolute(depth * current_x)
     168
     169print 'Calculating stage, xmomentum and ymomentum'
     170
     171#stage = depth - numpy.tile(numpy.absolute(elevation), number_of_timesteps)
     172stage = depth + numpy.tile(elevation, number_of_timesteps)
     173xmomentum = current_x*depth #momentum*numpy.sin(numpy.radians(current_y))
     174ymomentum = current_y*depth #momentum*numpy.cos(numpy.radians(current_y))
    165175
    166176numpy.put(xmomentum, no_value_index, 0)
    167177numpy.put(ymomentum, no_value_index, 0)
    168178
    169 #    if t > 0.0 and t < 90.0:
    170 #        assert mx >= 0 and my >= 0
    171        
    172 #    elif t > 90.0 and t < 180.0:
    173 #        assert mx >= 0 and my <= 0
    174        
    175 #    elif t > 180.0 and t < 270.0:
    176 #        assert mx <= 0  and my <= 0
    177        
    178 #    elif t > 270.0 and t < 360.0:
    179 #        assert mx <= 0 and my >= 0
    180        
    181 #    elif t == 0.0 or t == 360.0:
    182 #        assert my >= 0
    183        
    184 #    elif t == 90.0:
    185 #        assert mx >= 0
    186        
    187 #    elif t == 180.0:
    188 #        assert my <= 0
    189        
    190 #    elif t == 270.0:
    191 #        assert mx <= 0
    192        
    193 #    elif t == -9999:
    194 #        mx = 0
    195 #        my = 0
    196 #    else:
    197 #        print "Unexpected value of theta"
    198 #        exit()
    199        
    200 #    xmomentum[i] = mx
    201 #    ymomentum[i] = my
    202 
    203 assert len(numpy.where((numpy.sqrt(xmomentum**2 + ymomentum**2) - momentum) > 1.e-06)[0]) == 0
     179#assert len(numpy.where((numpy.sqrt(xmomentum**2 + ymomentum**2) - momentum) > 1.e-06)[0]) == 0
    204180   
    205181x_origin_int = int(10000*x_origin)
     
    215191
    216192assert time.size == number_of_timesteps
    217 assert momentum.size == stage.size
     193assert xmomentum.size == depth.size == ymomentum.size
    218194
    219195# -----------------------------
Note: See TracChangeset for help on using the changeset viewer.