| 1 | """ |
|---|
| 2 | Script for building boundary to run tsunami inundation scenario for sydney, |
|---|
| 3 | WA, Australia. The boundary is based on the National Hazard Map. |
|---|
| 4 | |
|---|
| 5 | Input: order_filename from project.py |
|---|
| 6 | event_number needs to be reflected in the project file |
|---|
| 7 | Output: creates a sts file and csv files stored in project.event_sts. |
|---|
| 8 | The run_sydney.py is reliant on the output of this script. |
|---|
| 9 | """ |
|---|
| 10 | import os |
|---|
| 11 | from os import sep |
|---|
| 12 | from os.path import join |
|---|
| 13 | import project |
|---|
| 14 | |
|---|
| 15 | from anuga.utilities.numerical_tools import ensure_numeric |
|---|
| 16 | from Scientific.IO.NetCDF import NetCDFFile |
|---|
| 17 | from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\ |
|---|
| 18 | compress,zeros,fabs,allclose,ones |
|---|
| 19 | from anuga.utilities.polygon import inside_polygon,read_polygon |
|---|
| 20 | from time import localtime, strftime, gmtime |
|---|
| 21 | from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary |
|---|
| 22 | |
|---|
| 23 | try: |
|---|
| 24 | from pylab import plot,show,xlabel,ylabel,legend,title,savefig,hold |
|---|
| 25 | except: |
|---|
| 26 | print 'Cannot import pylab plotting will not work. Csv files are still created' |
|---|
| 27 | |
|---|
| 28 | #-------------------------------------------------------------------------- |
|---|
| 29 | # Create sts boundary from mux2files |
|---|
| 30 | #-------------------------------------------------------------------------- |
|---|
| 31 | |
|---|
| 32 | dir=os.path.join(project.muxhome,'mux') |
|---|
| 33 | |
|---|
| 34 | # Refer to event_07875.list in home+state+sep+event08 for event details |
|---|
| 35 | # taken from David's event list for 7875 |
|---|
| 36 | urs_filenames = [ |
|---|
| 37 | "%s" % os.path.join(dir, 'Altiplano-0003-z.grd'), |
|---|
| 38 | "%s" % os.path.join(dir, 'Altiplano-0004-z.grd'), |
|---|
| 39 | "%s" % os.path.join(dir, 'Altiplano-0005-z.grd'), |
|---|
| 40 | "%s" % os.path.join(dir, 'Altiplano-0006-z.grd'), |
|---|
| 41 | "%s" % os.path.join(dir, 'Altiplano-0007-z.grd'), |
|---|
| 42 | "%s" % os.path.join(dir, 'Altiplano-0008-z.grd'), |
|---|
| 43 | "%s" % os.path.join(dir, 'Altiplano-0009-z.grd'), |
|---|
| 44 | "%s" % os.path.join(dir, 'Altiplano-0010-z.grd'), |
|---|
| 45 | "%s" % os.path.join(dir, 'Altiplano-0011-z.grd'), |
|---|
| 46 | "%s" % os.path.join(dir, 'Altiplano-0012-z.grd'), |
|---|
| 47 | "%s" % os.path.join(dir, 'Altiplano-0013-z.grd'), |
|---|
| 48 | "%s" % os.path.join(dir, 'Altiplano-0014-z.grd'), |
|---|
| 49 | "%s" % os.path.join(dir, 'Altiplano-0015-z.grd'), |
|---|
| 50 | "%s" % os.path.join(dir, 'Altiplano-0016-z.grd'), |
|---|
| 51 | "%s" % os.path.join(dir, 'Altiplano-0017-z.grd'), |
|---|
| 52 | "%s" % os.path.join(dir, 'Altiplano-0018-z.grd'), |
|---|
| 53 | "%s" % os.path.join(dir, 'Altiplano-0019-z.grd'), |
|---|
| 54 | "%s" % os.path.join(dir, 'Altiplano-0020-z.grd'), |
|---|
| 55 | "%s" % os.path.join(dir, 'Altiplano-0021-z.grd'), |
|---|
| 56 | "%s" % os.path.join(dir, 'Altiplano-0022-z.grd'), |
|---|
| 57 | "%s" % os.path.join(dir, 'Altiplano-0023-z.grd')] |
|---|
| 58 | |
|---|
| 59 | print 'number of sources', len(urs_filenames) |
|---|
| 60 | |
|---|
| 61 | if len(urs_filenames) == 21: #change per event |
|---|
| 62 | |
|---|
| 63 | weight_factor = 31.215400 #change per event |
|---|
| 64 | weights=weight_factor*ones(len(urs_filenames),Float) |
|---|
| 65 | |
|---|
| 66 | scenario_name=project.scenario_name |
|---|
| 67 | urs_order = join(project.urs_order) |
|---|
| 68 | |
|---|
| 69 | print 'reading', urs_order |
|---|
| 70 | |
|---|
| 71 | urs2sts(urs_filenames, |
|---|
| 72 | basename_out = join(project.event_sts), |
|---|
| 73 | ordering_filename = urs_order, |
|---|
| 74 | weights = weights, |
|---|
| 75 | mean_stage = project.tide, |
|---|
| 76 | verbose = True) |
|---|
| 77 | else: |
|---|
| 78 | print 'number of sources do not match event.list' |
|---|
| 79 | |
|---|
| 80 | #------------------------------------------------------------------------------ |
|---|
| 81 | # Get gauges (timeseries of index points) |
|---|
| 82 | #------------------------------------------------------------------------------ |
|---|
| 83 | def get_sts_gauge_data(filename,verbose=False): |
|---|
| 84 | from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\ |
|---|
| 85 | compress,zeros,fabs,take |
|---|
| 86 | fid = NetCDFFile(filename+'.sts', 'r') #Open existing file for read |
|---|
| 87 | permutation = fid.variables['permutation'][:] |
|---|
| 88 | x = fid.variables['x'][:]+fid.xllcorner #x-coordinates of vertices |
|---|
| 89 | y = fid.variables['y'][:]+fid.yllcorner #y-coordinates of vertices |
|---|
| 90 | points=transpose(asarray([x.tolist(),y.tolist()])) |
|---|
| 91 | time=fid.variables['time'][:]+fid.starttime |
|---|
| 92 | elevation=fid.variables['elevation'][:] |
|---|
| 93 | |
|---|
| 94 | basename='sts_gauge' |
|---|
| 95 | quantity_names=['stage','xmomentum','ymomentum'] |
|---|
| 96 | quantities = {} |
|---|
| 97 | for i, name in enumerate(quantity_names): |
|---|
| 98 | quantities[name] = fid.variables[name][:] |
|---|
| 99 | |
|---|
| 100 | #------------------------------------------------------------------------------ |
|---|
| 101 | # Get Maxium wave height throughout timeseries at each index point |
|---|
| 102 | #------------------------------------------------------------------------------ |
|---|
| 103 | |
|---|
| 104 | maxname = 'max_sts_stage.csv' |
|---|
| 105 | fid_max = open(project.event_sts+sep+maxname,'w') |
|---|
| 106 | s = 'index, x, y, max_stage \n' |
|---|
| 107 | fid_max.write(s) |
|---|
| 108 | for j in range(len(x)): |
|---|
| 109 | index= permutation[j] |
|---|
| 110 | stage = quantities['stage'][:,j] |
|---|
| 111 | xmomentum = quantities['xmomentum'][:,j] |
|---|
| 112 | ymomentum = quantities['ymomentum'][:,j] |
|---|
| 113 | |
|---|
| 114 | s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], max(stage)) |
|---|
| 115 | fid_max.write(s) |
|---|
| 116 | |
|---|
| 117 | #------------------------------------------------------------------------------ |
|---|
| 118 | # Get Minium wave height throughout timeseries at each index point |
|---|
| 119 | #------------------------------------------------------------------------------ |
|---|
| 120 | |
|---|
| 121 | minname = 'min_sts_stage.csv' |
|---|
| 122 | fid_min = open(project.event_sts+sep+minname,'w') |
|---|
| 123 | s = 'index, x, y, max_stage \n' |
|---|
| 124 | fid_min.write(s) |
|---|
| 125 | for j in range(len(x)): |
|---|
| 126 | index= permutation[j] |
|---|
| 127 | stage = quantities['stage'][:,j] |
|---|
| 128 | xmomentum = quantities['xmomentum'][:,j] |
|---|
| 129 | ymomentum = quantities['ymomentum'][:,j] |
|---|
| 130 | |
|---|
| 131 | s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], min(stage)) |
|---|
| 132 | fid_min.write(s) |
|---|
| 133 | |
|---|
| 134 | |
|---|
| 135 | |
|---|
| 136 | fid_sts = open(project.event_sts+sep+basename+'_'+ str(index)+'.csv', 'w') |
|---|
| 137 | s = 'time, stage, xmomentum, ymomentum \n' |
|---|
| 138 | fid_sts.write(s) |
|---|
| 139 | |
|---|
| 140 | #------------------------------------------------------------------------------ |
|---|
| 141 | # End of the get gauges |
|---|
| 142 | #------------------------------------------------------------------------------ |
|---|
| 143 | for k in range(len(time)-1): |
|---|
| 144 | s = '%.6f, %.6f, %.6f, %.6f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k]) |
|---|
| 145 | fid_sts.write(s) |
|---|
| 146 | |
|---|
| 147 | fid_sts.close() |
|---|
| 148 | |
|---|
| 149 | fid.close() |
|---|
| 150 | |
|---|
| 151 | |
|---|
| 152 | return quantities,elevation,time |
|---|
| 153 | |
|---|
| 154 | quantities,elevation,time=get_sts_gauge_data(os.path.join(project.event_sts),verbose=False) |
|---|
| 155 | |
|---|
| 156 | print len(elevation), len(quantities['stage'][0,:]) |
|---|
| 157 | |
|---|
