| 1 | """ |
|---|
| 2 | Script for building boundary to run tsunami inundation scenario for perth, |
|---|
| 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.boundaries_dir_event. |
|---|
| 8 | The run_perth.py is reliant on the output of this script. |
|---|
| 9 | """ |
|---|
| 10 | import os |
|---|
| 11 | from os import sep |
|---|
| 12 | import project |
|---|
| 13 | |
|---|
| 14 | from anuga.utilities.numerical_tools import ensure_numeric |
|---|
| 15 | from Scientific.IO.NetCDF import NetCDFFile |
|---|
| 16 | from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\ |
|---|
| 17 | compress,zeros,fabs,allclose,ones |
|---|
| 18 | from anuga.utilities.polygon import inside_polygon,read_polygon |
|---|
| 19 | from time import localtime, strftime, gmtime |
|---|
| 20 | from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary |
|---|
| 21 | |
|---|
| 22 | try: |
|---|
| 23 | from pylab import plot,show,xlabel,ylabel,legend,title,savefig,hold |
|---|
| 24 | except: |
|---|
| 25 | print 'Cannot import pylab plotting will not work. Csv files are still created' |
|---|
| 26 | |
|---|
| 27 | #-------------------------------------------------------------------------- |
|---|
| 28 | # Create sts boundary from mux2files |
|---|
| 29 | #-------------------------------------------------------------------------- |
|---|
| 30 | |
|---|
| 31 | dir=os.path.join(project.muxhome,'mux') |
|---|
| 32 | |
|---|
| 33 | # Refer to event_027255.list in home+state+sep+event08 for event details |
|---|
| 34 | # taken from David's event list for 27255 |
|---|
| 35 | urs_filenames = [ |
|---|
| 36 | os.path.join(dir,'Java-0000-z.grd'), |
|---|
| 37 | os.path.join(dir,'Java-0001-z.grd'), |
|---|
| 38 | os.path.join(dir,'Java-0002-z.grd'), |
|---|
| 39 | os.path.join(dir,'Java-0003-z.grd'), |
|---|
| 40 | os.path.join(dir,'Java-0005-z.grd'), |
|---|
| 41 | os.path.join(dir,'Java-0006-z.grd'), |
|---|
| 42 | os.path.join(dir,'Java-0007-z.grd'), |
|---|
| 43 | os.path.join(dir,'Java-0008-z.grd'), |
|---|
| 44 | os.path.join(dir,'Java-0010-z.grd'), |
|---|
| 45 | os.path.join(dir,'Java-0011-z.grd'), |
|---|
| 46 | os.path.join(dir,'Java-0012-z.grd'), |
|---|
| 47 | os.path.join(dir,'Java-0013-z.grd'), |
|---|
| 48 | os.path.join(dir,'Java-0015-z.grd'), |
|---|
| 49 | os.path.join(dir,'Java-0016-z.grd'), |
|---|
| 50 | os.path.join(dir,'Java-0017-z.grd'), |
|---|
| 51 | os.path.join(dir,'Java-0018-z.grd'), |
|---|
| 52 | os.path.join(dir,'Java-0020-z.grd'), |
|---|
| 53 | os.path.join(dir,'Java-0021-z.grd'), |
|---|
| 54 | os.path.join(dir,'Java-0022-z.grd'), |
|---|
| 55 | os.path.join(dir,'Java-0023-z.grd'), |
|---|
| 56 | os.path.join(dir,'Java-0025-z.grd'), |
|---|
| 57 | os.path.join(dir,'Java-0026-z.grd'), |
|---|
| 58 | os.path.join(dir,'Java-0027-z.grd'), |
|---|
| 59 | os.path.join(dir,'Java-0028-z.grd'), |
|---|
| 60 | os.path.join(dir,'Java-0029-z.grd'), |
|---|
| 61 | os.path.join(dir,'Java-0030-z.grd'), |
|---|
| 62 | os.path.join(dir,'Java-0031-z.grd'), |
|---|
| 63 | os.path.join(dir,'Java-0032-z.grd'), |
|---|
| 64 | os.path.join(dir,'Java-0033-z.grd'), |
|---|
| 65 | os.path.join(dir,'Java-0034-z.grd'), |
|---|
| 66 | os.path.join(dir,'Java-0035-z.grd'), |
|---|
| 67 | os.path.join(dir,'Java-0036-z.grd'), |
|---|
| 68 | os.path.join(dir,'Java-0037-z.grd'), |
|---|
| 69 | os.path.join(dir,'Java-0038-z.grd'), |
|---|
| 70 | os.path.join(dir,'Java-0039-z.grd'), |
|---|
| 71 | os.path.join(dir,'Java-0040-z.grd'), |
|---|
| 72 | os.path.join(dir,'Java-0041-z.grd'), |
|---|
| 73 | os.path.join(dir,'Java-0042-z.grd'), |
|---|
| 74 | os.path.join(dir,'Java-0043-z.grd'), |
|---|
| 75 | os.path.join(dir,'Java-0044-z.grd'), |
|---|
| 76 | os.path.join(dir,'Java-0045-z.grd'), |
|---|
| 77 | os.path.join(dir,'Java-0046-z.grd'), |
|---|
| 78 | os.path.join(dir,'Java-0047-z.grd'), |
|---|
| 79 | os.path.join(dir,'Java-0048-z.grd')] |
|---|
| 80 | |
|---|
| 81 | print 'number of sources', len(urs_filenames) |
|---|
| 82 | |
|---|
| 83 | if len(urs_filenames) == 44: #change per event |
|---|
| 84 | |
|---|
| 85 | # AS per David Burbidge email on Friday 4th July the mag 9.3 event |
|---|
| 86 | # has 1m worth of slip on each sub fault therefore mutliple each unit |
|---|
| 87 | # source by the slip (10.4544) and sum the 44 time series together to |
|---|
| 88 | # get the time series for this event at the points on your boundary. |
|---|
| 89 | |
|---|
| 90 | weight_factor = 10.4544 #change per event |
|---|
| 91 | weights=weight_factor*ones(len(urs_filenames),Float) |
|---|
| 92 | |
|---|
| 93 | scenario_name=project.scenario_name |
|---|
| 94 | order_filename=os.path.join(project.order_filename_dir) |
|---|
| 95 | |
|---|
| 96 | print 'reading', order_filename |
|---|
| 97 | # Create ordered sts file |
|---|
| 98 | print 'creating sts file' |
|---|
| 99 | |
|---|
| 100 | urs2sts(urs_filenames, |
|---|
| 101 | basename_out=os.path.join(project.boundaries_dir_event,scenario_name), |
|---|
| 102 | ordering_filename=order_filename, |
|---|
| 103 | weights=weights, |
|---|
| 104 | mean_stage=project.tide, |
|---|
| 105 | verbose=True) |
|---|
| 106 | else: |
|---|
| 107 | print 'number of sources do not match event.list' |
|---|
| 108 | |
|---|
| 109 | #------------------------------------------------------------------------------ |
|---|
| 110 | # Get gauges (timeseries of index points) |
|---|
| 111 | #------------------------------------------------------------------------------ |
|---|
| 112 | def get_sts_gauge_data(filename,verbose=False): |
|---|
| 113 | from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\ |
|---|
| 114 | compress,zeros,fabs,take |
|---|
| 115 | fid = NetCDFFile(filename+'.sts', 'r') #Open existing file for read |
|---|
| 116 | permutation = fid.variables['permutation'][:] |
|---|
| 117 | x = fid.variables['x'][:]+fid.xllcorner #x-coordinates of vertices |
|---|
| 118 | y = fid.variables['y'][:]+fid.yllcorner #y-coordinates of vertices |
|---|
| 119 | points=transpose(asarray([x.tolist(),y.tolist()])) |
|---|
| 120 | time=fid.variables['time'][:]+fid.starttime |
|---|
| 121 | elevation=fid.variables['elevation'][:] |
|---|
| 122 | |
|---|
| 123 | basename='sts_gauge' |
|---|
| 124 | quantity_names=['stage','xmomentum','ymomentum'] |
|---|
| 125 | quantities = {} |
|---|
| 126 | for i, name in enumerate(quantity_names): |
|---|
| 127 | quantities[name] = fid.variables[name][:] |
|---|
| 128 | |
|---|
| 129 | #------------------------------------------------------------------------------ |
|---|
| 130 | # Get Maxium wave height throughout timeseries at each index point |
|---|
| 131 | #------------------------------------------------------------------------------ |
|---|
| 132 | |
|---|
| 133 | maxname = 'max_sts_stage.csv' |
|---|
| 134 | fid_max = open(project.boundaries_dir_event+sep+maxname,'w') |
|---|
| 135 | s = 'index, x, y, max_stage \n' |
|---|
| 136 | fid_max.write(s) |
|---|
| 137 | for j in range(len(x)): |
|---|
| 138 | index= permutation[j] |
|---|
| 139 | stage = quantities['stage'][:,j] |
|---|
| 140 | xmomentum = quantities['xmomentum'][:,j] |
|---|
| 141 | ymomentum = quantities['ymomentum'][:,j] |
|---|
| 142 | |
|---|
| 143 | s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], max(stage)) |
|---|
| 144 | fid_max.write(s) |
|---|
| 145 | |
|---|
| 146 | #------------------------------------------------------------------------------ |
|---|
| 147 | # Get Minium wave height throughout timeseries at each index point |
|---|
| 148 | #------------------------------------------------------------------------------ |
|---|
| 149 | |
|---|
| 150 | minname = 'min_sts_stage.csv' |
|---|
| 151 | fid_min = open(project.boundaries_dir_event+sep+minname,'w') |
|---|
| 152 | s = 'index, x, y, max_stage \n' |
|---|
| 153 | fid_min.write(s) |
|---|
| 154 | for j in range(len(x)): |
|---|
| 155 | index= permutation[j] |
|---|
| 156 | stage = quantities['stage'][:,j] |
|---|
| 157 | xmomentum = quantities['xmomentum'][:,j] |
|---|
| 158 | ymomentum = quantities['ymomentum'][:,j] |
|---|
| 159 | |
|---|
| 160 | s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], min(stage)) |
|---|
| 161 | fid_min.write(s) |
|---|
| 162 | |
|---|
| 163 | |
|---|
| 164 | |
|---|
| 165 | fid_sts = open(project.boundaries_dir_event+sep+basename+'_'+ str(index)+'.csv', 'w') |
|---|
| 166 | s = 'time, stage, xmomentum, ymomentum \n' |
|---|
| 167 | fid_sts.write(s) |
|---|
| 168 | |
|---|
| 169 | #------------------------------------------------------------------------------ |
|---|
| 170 | # End of the get gauges |
|---|
| 171 | #------------------------------------------------------------------------------ |
|---|
| 172 | for k in range(len(time)-1): |
|---|
| 173 | s = '%.6f, %.6f, %.6f, %.6f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k]) |
|---|
| 174 | fid_sts.write(s) |
|---|
| 175 | |
|---|
| 176 | fid_sts.close() |
|---|
| 177 | |
|---|
| 178 | fid.close() |
|---|
| 179 | |
|---|
| 180 | |
|---|
| 181 | return quantities,elevation,time |
|---|
| 182 | |
|---|
| 183 | quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir_event,project.scenario_name),verbose=False) |
|---|
| 184 | |
|---|
| 185 | print len(elevation), len(quantities['stage'][0,:]) |
|---|
| 186 | |
|---|
