""" Script for building boundary to run tsunami inundation scenario for sydney, WA, Australia. The boundary is based on the National Hazard Map. Input: order_filename from project.py event_number needs to be reflected in the project file Output: creates a sts file and csv files stored in project.event_sts. The run_sydney.py is reliant on the output of this script. """ import os from os import sep from os.path import join 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.muxhome,'mux') # Refer to event_07875.list in home+state+sep+event08 for event details # taken from David's event list for 7875 urs_filenames = [ "%s" % os.path.join(dir, 'Altiplano-0003-z.grd'), "%s" % os.path.join(dir, 'Altiplano-0004-z.grd'), "%s" % os.path.join(dir, 'Altiplano-0005-z.grd'), "%s" % os.path.join(dir, 'Altiplano-0006-z.grd'), "%s" % os.path.join(dir, 'Altiplano-0007-z.grd'), "%s" % os.path.join(dir, 'Altiplano-0008-z.grd'), "%s" % os.path.join(dir, 'Altiplano-0009-z.grd'), "%s" % os.path.join(dir, 'Altiplano-0010-z.grd'), "%s" % os.path.join(dir, 'Altiplano-0011-z.grd'), "%s" % os.path.join(dir, 'Altiplano-0012-z.grd'), "%s" % os.path.join(dir, 'Altiplano-0013-z.grd'), "%s" % os.path.join(dir, 'Altiplano-0014-z.grd'), "%s" % os.path.join(dir, 'Altiplano-0015-z.grd'), "%s" % os.path.join(dir, 'Altiplano-0016-z.grd'), "%s" % os.path.join(dir, 'Altiplano-0017-z.grd'), "%s" % os.path.join(dir, 'Altiplano-0018-z.grd'), "%s" % os.path.join(dir, 'Altiplano-0019-z.grd'), "%s" % os.path.join(dir, 'Altiplano-0020-z.grd'), "%s" % os.path.join(dir, 'Altiplano-0021-z.grd'), "%s" % os.path.join(dir, 'Altiplano-0022-z.grd'), "%s" % os.path.join(dir, 'Altiplano-0023-z.grd')] print 'number of sources', len(urs_filenames) if len(urs_filenames) == 21: #change per event weight_factor = 31.215400 #change per event weights=weight_factor*ones(len(urs_filenames),Float) scenario_name=project.scenario_name urs_order = join(project.urs_order) print 'reading', urs_order urs2sts(urs_filenames, basename_out = join(project.event_sts), ordering_filename = urs_order, 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.event_sts+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.event_sts+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.event_sts+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.event_sts),verbose=False) print len(elevation), len(quantities['stage'][0,:])