[6291] | 1 | """ |
---|
| 2 | Script for building boundary to run a tsunami inundation scenario for busselton, |
---|
| 3 | WA, Australia. The boundary is sourced from 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.boundaries_dir_event. |
---|
| 8 | The run_busselton.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 | # location of mux files |
---|
| 32 | # Note, this may change when all mux files are on the GA network |
---|
| 33 | dir=join(project.muxhome,'mux') |
---|
| 34 | |
---|
| 35 | # Refer to event_027255.list in home+state+sep+event08 for event details |
---|
| 36 | # taken from David's event list for 27255 |
---|
| 37 | # Note, this is specific to the event and will be supplied by David. |
---|
| 38 | # In time, this will become more automated. |
---|
| 39 | urs_filenames = [ |
---|
| 40 | os.path.join(dir,'Java-0017-z.grd'), |
---|
| 41 | os.path.join(dir,'Java-0018-z.grd'), |
---|
| 42 | os.path.join(dir,'Java-0019-z.grd'), |
---|
| 43 | os.path.join(dir,'Java-0022-z.grd'), |
---|
| 44 | os.path.join(dir,'Java-0023-z.grd'), |
---|
| 45 | os.path.join(dir,'Java-0024-z.grd'), |
---|
| 46 | os.path.join(dir,'Java-0027-z.grd'), |
---|
| 47 | os.path.join(dir,'Java-0028-z.grd'), |
---|
| 48 | os.path.join(dir,'Java-0031-z.grd'), |
---|
| 49 | os.path.join(dir,'Java-0032-z.grd')] |
---|
| 50 | |
---|
| 51 | print 'number of sources', len(urs_filenames) |
---|
| 52 | |
---|
| 53 | if len(urs_filenames) == 10: #change per event |
---|
| 54 | |
---|
| 55 | weight_factor = 30.84050 #change per event |
---|
| 56 | weights=weight_factor*ones(len(urs_filenames),Float) |
---|
| 57 | |
---|
| 58 | scenario_name = project.scenario_name |
---|
| 59 | urs_order = join(project.urs_order) |
---|
| 60 | |
---|
| 61 | print 'reading', urs_order |
---|
| 62 | |
---|
| 63 | urs2sts(urs_filenames, |
---|
| 64 | basename_out = join(project.event_sts), |
---|
| 65 | ordering_filename = urs_order, |
---|
| 66 | weights = weights, |
---|
| 67 | mean_stage = project.tide, |
---|
| 68 | verbose = True) |
---|
| 69 | else: |
---|
| 70 | print 'number of sources do not match event.list' |
---|
| 71 | |
---|
| 72 | #------------------------------------------------------------------------------ |
---|
| 73 | # Get gauges (timeseries of index points) |
---|
| 74 | #------------------------------------------------------------------------------ |
---|
| 75 | def get_sts_gauge_data(filename,verbose=False): |
---|
| 76 | from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\ |
---|
| 77 | compress,zeros,fabs,take |
---|
| 78 | fid = NetCDFFile(filename+'.sts', 'r') #Open existing file for read |
---|
| 79 | permutation = fid.variables['permutation'][:] |
---|
| 80 | x = fid.variables['x'][:]+fid.xllcorner #x-coordinates of vertices |
---|
| 81 | y = fid.variables['y'][:]+fid.yllcorner #y-coordinates of vertices |
---|
| 82 | points=transpose(asarray([x.tolist(),y.tolist()])) |
---|
| 83 | time=fid.variables['time'][:]+fid.starttime |
---|
| 84 | elevation=fid.variables['elevation'][:] |
---|
| 85 | |
---|
| 86 | basename='sts_gauge' |
---|
| 87 | quantity_names=['stage','xmomentum','ymomentum'] |
---|
| 88 | quantities = {} |
---|
| 89 | for i, name in enumerate(quantity_names): |
---|
| 90 | quantities[name] = fid.variables[name][:] |
---|
| 91 | |
---|
| 92 | #------------------------------------------------------------------------------ |
---|
| 93 | # Get maxium wave height throughout timeseries at each index point |
---|
| 94 | #------------------------------------------------------------------------------ |
---|
| 95 | |
---|
| 96 | maxname = 'max_sts_stage.csv' |
---|
| 97 | fid_max = open(project.boundaries_dir_event+sep+maxname,'w') |
---|
| 98 | s = 'index, x, y, max_stage \n' |
---|
| 99 | fid_max.write(s) |
---|
| 100 | for j in range(len(x)): |
---|
| 101 | index= permutation[j] |
---|
| 102 | stage = quantities['stage'][:,j] |
---|
| 103 | xmomentum = quantities['xmomentum'][:,j] |
---|
| 104 | ymomentum = quantities['ymomentum'][:,j] |
---|
| 105 | |
---|
| 106 | s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], max(stage)) |
---|
| 107 | fid_max.write(s) |
---|
| 108 | |
---|
| 109 | #------------------------------------------------------------------------------ |
---|
| 110 | # Get minium wave height throughout timeseries at each index point |
---|
| 111 | #------------------------------------------------------------------------------ |
---|
| 112 | |
---|
| 113 | minname = 'min_sts_stage.csv' |
---|
| 114 | fid_min = open(project.boundaries_dir_event+sep+minname,'w') |
---|
| 115 | s = 'index, x, y, max_stage \n' |
---|
| 116 | fid_min.write(s) |
---|
| 117 | for j in range(len(x)): |
---|
| 118 | index= permutation[j] |
---|
| 119 | stage = quantities['stage'][:,j] |
---|
| 120 | xmomentum = quantities['xmomentum'][:,j] |
---|
| 121 | ymomentum = quantities['ymomentum'][:,j] |
---|
| 122 | |
---|
| 123 | s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], min(stage)) |
---|
| 124 | fid_min.write(s) |
---|
| 125 | |
---|
| 126 | |
---|
| 127 | |
---|
| 128 | fid_sts = open(project.boundaries_dir_event+sep+basename+'_'+ str(index)+'.csv', 'w') |
---|
| 129 | s = 'time, stage, xmomentum, ymomentum \n' |
---|
| 130 | fid_sts.write(s) |
---|
| 131 | |
---|
| 132 | #------------------------------------------------------------------------------ |
---|
| 133 | # End of the get gauges |
---|
| 134 | #------------------------------------------------------------------------------ |
---|
| 135 | for k in range(len(time)-1): |
---|
| 136 | s = '%.6f, %.6f, %.6f, %.6f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k]) |
---|
| 137 | fid_sts.write(s) |
---|
| 138 | |
---|
| 139 | fid_sts.close() |
---|
| 140 | |
---|
| 141 | fid.close() |
---|
| 142 | |
---|
| 143 | |
---|
| 144 | return quantities,elevation,time |
---|
| 145 | |
---|
| 146 | quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir_event,project.scenario_name),verbose=False) |
---|
| 147 | |
---|
| 148 | print len(elevation), len(quantities['stage'][0,:]) |
---|
| 149 | |
---|