Changeset 7692


Ignore:
Timestamp:
Apr 22, 2010, 3:07:19 PM (10 years ago)
Author:
habili
Message:

Changes to make code more memory efficient

File:
1 edited

Legend:

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

    r7688 r7692  
    4242print "Elevation file name: ", elevation_file_name
    4343
    44 grid_file_names = glob.glob(os.path.join(event_folder, 'gcom_txt19780404.0*.gz'))
     44grid_file_names = glob.glob(os.path.join(event_folder, 'gcom*.gz'))
    4545grid_file_names.sort()
    4646
    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_temp
    52    
    53 # number_of_timesteps = len(grid_file_names)
     47grid_file_names_temp = []
     48for i in range(0, 240, 2):
     49    grid_file_names_temp.append(grid_file_names[i])
     50
     51grid_file_names = grid_file_names_temp
     52   
     53number_of_timesteps = len(grid_file_names)
    5454
    5555print "Number of timesteps: ", number_of_timesteps
    5656
    5757start_time = 0      # all times in seconds
    58 timestep = 720     # all times in seconds
     58timestep = 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 = []
     64#stage = []
     65#speed = []
     66#theta = []
    6767
    6868elevation_file = open(elevation_file_name, 'rb')
     
    9494print 'Grid size: ', grid_size
    9595print 'No data value: ', no_data
    96        
    97 for f in grid_file_names:
     96
     97stage = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float)
     98speed = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float)
     99theta = numpy.empty(number_of_timesteps*nrows*ncols, dtype=float)
     100
     101import pdb
     102
     103for j, f in enumerate(grid_file_names):
    98104    print 'file: ', f
    99105   
    100106    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):
    105117       
    106118        if i > 7 and i < (8 + nrows):   #  7 refers to the number of rows of header info that we skip - always CHECK format
    107             stage += [L.strip().split()]
     119            st += L.strip().split()
    108120           
    109121        if i > (9 + nrows) and i < (10 + 2*nrows):
    110             speed += [L.strip().split()]
     122            sp += L.strip().split()
    111123           
    112124        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()
    114134
    115135#------------------------------------------------------
     
    120140elevation = numpy.array(elevation).astype('d')
    121141
    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*i
    128     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()
    135155   
    136156print 'Size of elevation array: ', elevation.size
     
    152172# ---------------------------------------------
    153173
    154 for i in xrange(number_of_timesteps):
    155     j = number_of_points * i
    156     k = j + number_of_points
     174#for i in xrange(number_of_timesteps):
     175#    j = number_of_points * i
     176#    k = j + number_of_points
    157177    # 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
     185no_value_index = numpy.where(((stage < -100) + (speed < 0) + (theta < 0)) == True)[0]
     186depth = stage - numpy.tile(elevation, number_of_timesteps)
     187numpy.put(stage, no_value_index, -9999)
     188numpy.put(speed, no_value_index, -9999)
     189numpy.put(theta, no_value_index, -9999)
     190numpy.put(depth, no_value_index, 0)
    164191
    165192# for i in xrange(number_of_timesteps):
     
    174201            # depth[j+q] = 0
    175202   
    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 stage-elevation
     204# is slightly -ve  - why I don't know, possbly a rounding error?
     205momentum = numpy.absolute(depth * speed)
     206
     207#xmomentum = numpy.empty(momentum.size, dtype='d')
     208#ymomentum = numpy.empty(momentum.size, dtype='d')
    180209
    181210print 'Calculating xmomentum and ymomentum'
    182211
    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.e-06
     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
     217xmomentum = momentum*numpy.sin(numpy.radians(theta))
     218ymomentum = momentum*numpy.cos(numpy.radians(theta))
     219
     220numpy.put(xmomentum, no_value_index, 0)
     221numpy.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
     257assert 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
    223260   
    224261x_origin_int = int(10000*x_origin)
     
    241278print "Creating the STS NetCDF file"
    242279
    243 fid = NetCDFFile(os.path.join(event_folder, event + '_test.sts'), 'w')
     280#for j in range(2):
     281
     282fid = NetCDFFile(os.path.join(event_folder, event + '_master_2_.sts'), 'wl')
    244283fid.institution = 'Geoscience Australia'
    245284fid.description = "description"
     
    270309fid.variables['ymomentum_range'][:] = numpy.array([1e+036, -1e+036])
    271310fid.variables['elevation'][:] = elevation
    272 fid.variables['time'][:] = time
     311fid.variables['time'][:] = time#[j*60 : j*60 + 60]
    273312
    274313s = fid.variables['stage']
     
    276315ym = fid.variables['ymomentum']
    277316
     317#jj = j*number_of_points*(number_of_timesteps/2)
     318
    278319for i in xrange(number_of_timesteps):
    279     ii = i*number_of_points
     320    ii = i*number_of_points# + jj
    280321    s[i] = stage[ii : ii + number_of_points]
    281322    xm[i] = xmomentum[ii : ii + number_of_points]
    282323    ym[i] = ymomentum[ii : ii + number_of_points]
    283        
     324
    284325origin = Geo_reference(refzone, min(x), min(y)) # Check this section for inputs in eastings and northings - it works for long and lat
    285326geo_ref = write_NetCDF_georeference(origin, fid)
     
    289330
    290331fid.close()
     332
Note: See TracChangeset for help on using the changeset viewer.