Changeset 7723


Ignore:
Timestamp:
May 13, 2010, 9:51:30 AM (14 years ago)
Author:
sexton
Message:

updates

File:
1 edited

Legend:

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

    r7692 r7723  
    2222state = 'western_australia'
    2323scenario_folder = 'bunbury_storm_surge_scenario_2009'
    24 event = 'case_1'
     24event = 'gcom_60min'
    2525
    2626print "Processing event: ", event
     
    4646
    4747grid_file_names_temp = []
    48 for i in range(0, 240, 2):
     48for i in range(1, 32, 1):
    4949    grid_file_names_temp.append(grid_file_names[i])
    5050
     
    5656
    5757start_time = 0      # all times in seconds
    58 timestep = 1440 #720     # all times in seconds
     58timestep = 3600 #1440 #720     # all times in seconds
    5959end_time =  start_time + (timestep*number_of_timesteps) # all times in seconds
    6060refzone = 50        # UTM zone of model
     
    6262
    6363elevation = []
    64 #stage = []
    65 #speed = []
    66 #theta = []
    6764
    6865elevation_file = open(elevation_file_name, 'rb')
     
    9592print 'No data value: ', no_data
    9693
     94#------------------------------------------------------
     95# Create arrays of elevation, stage, speed and theta
     96#-------------------------------------------------------
     97print 'creating numpy arrays: elevation, stage, speed, theta'
     98
    9799stage = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float)
    98100speed = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float)
    99101theta = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float)
    100102
    101 import pdb
    102 
    103103for j, f in enumerate(grid_file_names):
    104104    print 'file: ', f
    105105   
    106106    gz_file = gzip.open(f, 'rb')
    107     #lines = gz_file.readlines()
    108     #gz_file.close()
    109 
    110     #pdb.set_trace()
    111107
    112108    st = []
    113109    sp = []
    114     th = []
     110    t = []
    115111
    116112    for i, L in enumerate(gz_file):
     
    123119           
    124120        if i > (11 + 2*nrows):
    125             th += L.strip().split()
     121            t += L.strip().split()
    126122
    127123    stage[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(st).astype('d'))
    128124    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'))
     125    theta[j*nrows*ncols : (j+1)*nrows*ncols] = numpy.flipud(numpy.array(t).astype('d'))
    130126
    131127    gz_file.close()
    132    
    133 #pdb.set_trace()
    134 
    135 #------------------------------------------------------
    136 # Create arrays of elevation, stage, speed and theta
    137 #-------------------------------------------------------
    138 print 'creating numpy arrays: elevation, stage, speed, theta'
    139            
     128
    140129elevation = numpy.array(elevation).astype('d')
    141130
    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()
    155    
    156131print 'Size of elevation array: ', elevation.size
    157132print 'Size of stage array: ', stage.size
     
    162137assert stage.size == number_of_timesteps * elevation.size
    163138
    164 d = numpy.empty(elevation.size, dtype='d')
    165139depth = numpy.empty(stage.size, dtype='d')
    166140number_of_points = ncols * nrows
     
    171145# a consistant dataset
    172146# ---------------------------------------------
    173 
    174 #for i in xrange(number_of_timesteps):
    175 #    j = number_of_points * i
    176 #    k = j + number_of_points
    177     # no_value_mask = (stage[j:k] < -9000) + (speed[j:k] < -9000) + (theta[j:k] < -9000)
    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 
     147import pdb
     148pdb.set_trace()
    185149no_value_index = numpy.where(((stage < -100) + (speed < 0) + (theta < 0)) == True)[0]
    186150depth = stage - numpy.tile(elevation, number_of_timesteps)
     151
    187152numpy.put(stage, no_value_index, -9999)
    188153numpy.put(speed, no_value_index, -9999)
    189154numpy.put(theta, no_value_index, -9999)
    190155numpy.put(depth, no_value_index, 0)
    191 
    192 # for i in xrange(number_of_timesteps):
    193     # j = number_of_points * i
    194     # k = j + number_of_points
    195     # d = stage[j:k] - elevation     # Assumes elevations below sea level are negative
    196     # depth[j:k] = [0 if x<-9000 else x for x in d]   # sets stage to zero for no-data values of stage
    197    
    198     # 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
    199     # for q in xrange(number_of_points):
    200         # if elevation[q] >= 0 and stage[j+q]==0:
    201             # depth[j+q] = 0
    202156   
    203157# Taking absolute value is to account for -ve depths obtained when stage-elevation
     
    205159momentum = numpy.absolute(depth * speed)
    206160
    207 #xmomentum = numpy.empty(momentum.size, dtype='d')
    208 #ymomentum = numpy.empty(momentum.size, dtype='d')
    209 
    210161print 'Calculating xmomentum and ymomentum'
    211 
    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))
    216162
    217163xmomentum = momentum*numpy.sin(numpy.radians(theta))
     
    256202
    257203assert len(numpy.where((numpy.sqrt(xmomentum**2 + ymomentum**2) - momentum) > 1.e-06)[0]) == 0
    258 
    259 #assert numpy.sqrt(xmomentum**2 + ymomentum**2) - momentum < 1.e-06
    260204   
    261205x_origin_int = int(10000*x_origin)
     
    280224#for j in range(2):
    281225
    282 fid = NetCDFFile(os.path.join(event_folder, event + '_master_2_.sts'), 'wl')
     226fid = NetCDFFile(os.path.join(event_folder, event + '_master_2_1.sts'), 'wl')
    283227fid.institution = 'Geoscience Australia'
    284228fid.description = "description"
Note: See TracChangeset for help on using the changeset viewer.