source: branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/build_gems_boundary.py @ 7794

# ---------------------------------------------------------------------------
# This script reads in both the master STS (output from process_gems_grids.py) 
# and the gems boundary order file (defined as gems_order in project.py)
# to generate a boundary STS file which contains data just at the boundary 
# points for all timesteps, which is input into run_model.py in ANUGA.
# Written by Nariman Habili and Leharne Fountain, 2010
# ---------------------------------------------------------------------------

import os
import numpy
from Scientific.IO.NetCDF import NetCDFFile
from anuga.utilities.polygon import read_polygon
import sys
# sys.path.append('W:\\sandpits\\lfountain\\anuga_work\\production\\mandurah_storm_surge_2009\\gems_comparison')
#sys.path.append('/nas/gemd/georisk_models/inundation/sandpits/lfountain/anuga_work/production/mandurah_storm_surge_2009/gems_comparison')
import project
from Scientific.IO.NetCDF import NetCDFFile
from anuga.coordinate_transforms.geo_reference import Geo_reference, write_NetCDF_georeference

#-----------------
# Directory setup
#-------------------
state = 'western_australia'
scenario_folder = 'bunbury_storm_surge_scenario_2009'
event = '20100527_gcom_12min' # Baseline TC Alby event

ENV_INUNDATIONHOME = 'INUNDATIONHOME'
home = os.path.join(os.getenv(ENV_INUNDATIONHOME), 'data') # Absolute path for data folder
# GEMS_folder = os.path.join(home, state, scenario_folder, 'GEMS')
anuga_folder = os.path.join(home, state, scenario_folder, 'anuga')
boundaries_folder = os.path.join(anuga_folder, 'boundaries')
event_folder = os.path.join(boundaries_folder, event)

#--------------------------------------------------------
# Define boundary points to extract from master STS file
#--------------------------------------------------------
gems_boundary = numpy.array(read_polygon(project.gems_order))
#gems_boundary = numpy.array(read_polygon(project.gauges))
number_of_points = gems_boundary.shape[0]
print 'hello', number_of_points
print 'points', gems_boundary

#---------------------------------------------------------
# Read master sts file generated from GEMS ascii grids
#---------------------------------------------------------
print "Reading master NetCDF file"
infile = NetCDFFile(os.path.join(event_folder, event + '_master_2_1.sts'), 'rl') #rl for large netcdf files

grid_size = infile.grid_size
ncols = infile.ncols
nrows = infile.nrows
no_data = infile.no_data
refzone = infile.zone
x_origin = infile.xllcorner
y_origin = infile.yllcorner
x = infile.variables['x'][:]
y = infile.variables['y'][:]
time = infile.variables['time'][:]
elevation = infile.variables['elevation'][:]
stage_array = infile.variables['stage'][:]
xmomentum_array = infile.variables['xmomentum'][:]
ymomentum_array = infile.variables['ymomentum'][:]
number_of_timesteps = stage_array.shape[0]

infile.close()

#--------------------------------------------
# Creating index from boundary points
#--------------------------------------------
print "Creating index from boundary points" # finding location of boundary points in master file

origin = numpy.array([x_origin[0], y_origin[0]])
index_array_2D = ((gems_boundary - origin)/grid_size).astype(int)
assert index_array_2D[:,0].min() >= 0 and index_array_2D[:,0].max() <= ncols # These test that the boundary points lie within the data area
assert index_array_2D[:,1].min() >= 0 and index_array_2D[:,1].max() <= nrows

index_array_1D = numpy.empty(index_array_2D.shape[0], int)
index_array_1D = index_array_2D[:,1]*ncols + index_array_2D[:,0]

#--------------------------------------------------
# Extract values from master sts at index points
#--------------------------------------------------
print "Extract values from master sts at index points"

x = x.take(index_array_1D)
y = y.take(index_array_1D)
elevation = elevation.take(index_array_1D)
stage = numpy.empty([number_of_timesteps, number_of_points], float)
xmomentum = numpy.empty([number_of_timesteps, number_of_points], float)
ymomentum = numpy.empty([number_of_timesteps, number_of_points], float)

for i, k in enumerate(index_array_1D):
    stage[:,i] = stage_array[:,k] 
    xmomentum[:,i] = xmomentum_array[:,k] 
    ymomentum[:,i] = ymomentum_array[:,k] 

#-----------------------------
# Create the boundary STS file
#-----------------------------
print "Creating the boundary STS NetCDF file"

fid = NetCDFFile(os.path.join(event_folder, event + '_boundary.sts'), 'w')
fid.institution = 'Geoscience Australia'
fid.description = "description"
fid.starttime = 0.0
fid.ncols = ncols
fid.nrows = nrows
fid.grid_size = grid_size
fid.no_data = no_data
fid.createDimension('number_of_points', number_of_points)
fid.createDimension('number_of_timesteps', number_of_timesteps)
fid.createDimension('numbers_in_range', 2)

fid.createVariable('x', 'd', ('number_of_points',))
fid.createVariable('y', 'd', ('number_of_points',))
fid.createVariable('elevation', 'd', ('number_of_points',))
fid.createVariable('elevation_range', 'd', ('numbers_in_range',))
fid.createVariable('time', 'i', ('number_of_timesteps',))
fid.createVariable('stage', 'd', ('number_of_timesteps', 'number_of_points'))
fid.createVariable('stage_range', 'd', ('numbers_in_range', ))
fid.createVariable('xmomentum', 'd', ('number_of_timesteps', 'number_of_points'))
fid.createVariable('xmomentum_range', 'd', ('numbers_in_range',))
fid.createVariable('ymomentum', 'd', ('number_of_timesteps', 'number_of_points'))
fid.createVariable('ymomentum_range', 'd', ('numbers_in_range',))

fid.variables['elevation_range'][:] = numpy.array([1e+036, -1e+036])
fid.variables['stage_range'][:] = numpy.array([1e+036, -1e+036])
fid.variables['xmomentum_range'][:] = numpy.array([1e+036, -1e+036])
fid.variables['ymomentum_range'][:] = numpy.array([1e+036, -1e+036])
fid.variables['elevation'][:] = elevation
fid.variables['time'][:] = time

s = fid.variables['stage']
xm = fid.variables['xmomentum']
ym = fid.variables['ymomentum']

for i in xrange(number_of_timesteps):
    s[i] = stage[i,:]
    xm[i] = xmomentum[i,:]
    ym[i] = ymomentum[i,:] 

#---------------------------------------------------------------------------
# Get timeseries values for wave height and components of momentum
#---------------------------------------------------------------------------

basename = 'sts_gauge_tide'
for i in xrange(number_of_points):
    out_file = os.path.join(event_folder,
                            basename+'_'+str(i)+'.csv')
    fid_sts = open(out_file, 'w')
    fid_sts.write('time, stage, xmomentum, ymomentum, elevation \n')
    for j in xrange(number_of_timesteps):
        fid_sts.write('%.6f, %.6f, %.6f, %.6f, %.6f\n'
                    % (time[j], stage[j,i], xmomentum[j,i], ymomentum[j,i], elevation[i]))

    fid_sts.close() 

origin = Geo_reference(refzone, x_origin, y_origin)
geo_ref = write_NetCDF_georeference(origin, fid)

fid.variables['x'][:] = x 
fid.variables['y'][:] = y  

fid.close()

sys.exit()

print 'get timeseries'
#-------------------------------------------------------------------------------
# Get gauges (timeseries of index points)
#-------------------------------------------------------------------------------

basename = 'sts_gauge'
quantity_names = ['stage', 'xmomentum', 'ymomentum']
quantities = {}
for i, name in enumerate(quantity_names):
    quantities[name] = fid.variables[name][:]


for j in range(len(x)):
    index = j # permutation[j]
    stage = quantities['stage'][:,j]
    #print 'hello', j, stage
    xmomentum = quantities['xmomentum'][:,j]
    ymomentum = quantities['ymomentum'][:,j]

    out_file = os.path.join(event_folder,
                            basename+'_'+str(index)+'.csv')
    fid_sts = open(out_file, 'w')
    fid_sts.write('time, stage, xmomentum, ymomentum \n')

    #-----------------------------------------------------------------------
    # End of the get gauges
    #-----------------------------------------------------------------------
    for k in range(len(time)-1):
        fid_sts.write('%.6f, %.6f, %.6f, %.6f\n'
                      % (time[k], stage[k], xmomentum[k], ymomentum[k]))