source: anuga_work/production/patong/build_boundary.py @ 6816

"""
Script for building boundary to run tsunami inundation scenario for onslow,
WA, Australia.  The boundary is based on the National Hazard Map.

Input: order_filename from project.py
Output: creates a sts file and csv files stored in project.boundaries_dir
The run_patong.py is reliant on the output of this script.
"""
import os
from os import sep
import project

from anuga.utilities.numerical_tools import ensure_numeric
from Scientific.IO.NetCDF import NetCDFFile
from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
    compress,zeros,fabs,allclose,ones
from anuga.utilities.polygon import inside_polygon,read_polygon
from time import localtime, strftime, gmtime
from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary

try:
    from pylab import plot,show,xlabel,ylabel,legend,title,savefig,hold
except:
    print 'Cannot import pylab plotting will not work. Csv files are still created'

#--------------------------------------------------------------------------
# Create sts boundary from mux2files 
#--------------------------------------------------------------------------

dir=os.path.join(project.boundaries_dir,'mux')

# Refer to event_
urs_filenames = [os.path.join(dir,'tgs.out')]

if len(urs_filenames) == 1: #change per event    
    weight_factor = 1
    weights=weight_factor*ones(len(urs_filenames),Float)
        
    scenario_name=project.scenario_name
    order_filename=os.path.join(project.order_filename_dir)

    print 'reading', order_filename
    # Create ordered sts file
    print 'creating sts file'

    urs2sts(urs_filenames,
            basename_out=os.path.join(project.boundaries_dir,scenario_name),
            ordering_filename=order_filename,
            weights=weights,
            mean_stage=project.tide,
            verbose=True)
else:
    print 'number of sources do not match event.list'

#------------------------------------------------------------------------------
# Get gauges (timeseries of index points)
#------------------------------------------------------------------------------
def get_sts_gauge_data(filename,verbose=False):
    from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
        compress,zeros,fabs,take
    fid = NetCDFFile(filename+'.sts', 'r')    #Open existing file for read
    permutation = fid.variables['permutation'][:]
    x = fid.variables['x'][:]+fid.xllcorner   #x-coordinates of vertices
    y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
    points=transpose(asarray([x.tolist(),y.tolist()]))
    time=fid.variables['time'][:]+fid.starttime
    elevation=fid.variables['elevation'][:]
       
    basename='sts_gauge'
    quantity_names=['stage','xmomentum','ymomentum']
    quantities = {}
    for i, name in enumerate(quantity_names):
        quantities[name] = fid.variables[name][:]

#------------------------------------------------------------------------------
# Get Maxium wave height throughout timeseries at each index point
#------------------------------------------------------------------------------

    maxname = 'max_sts_stage.csv'
    fid_max = open(project.boundaries_dir+sep+maxname,'w')
    s = 'index, x, y, max_stage \n'
    fid_max.write(s)    
    for j in range(len(x)):
        index= permutation[j]
        stage = quantities['stage'][:,j]
        xmomentum = quantities['xmomentum'][:,j]
        ymomentum = quantities['ymomentum'][:,j]

        s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], max(stage))
        fid_max.write(s)
      
#------------------------------------------------------------------------------
# Get Minium wave height throughout timeseries at each index point
#------------------------------------------------------------------------------

    minname = 'min_sts_stage.csv'
    fid_min = open(project.boundaries_dir+sep+minname,'w')
    s = 'index, x, y, max_stage \n'
    fid_min.write(s)    
    for j in range(len(x)):
        index= permutation[j]
        stage = quantities['stage'][:,j]
        xmomentum = quantities['xmomentum'][:,j]
        ymomentum = quantities['ymomentum'][:,j]

        s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], min(stage))
        fid_min.write(s)
        


        fid_sts = open(project.boundaries_dir+sep+basename+'_'+ str(index)+'.csv', 'w')
        s = 'time, stage, xmomentum, ymomentum \n'
        fid_sts.write(s)

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

        fid_sts.close()      

    fid.close()


    return quantities,elevation,time

quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir,project.scenario_name),verbose=False)

print len(elevation), len(quantities['stage'][0,:])