| 1 | import os |
|---|
| 2 | from os import sep |
|---|
| 3 | import project |
|---|
| 4 | try: |
|---|
| 5 | from pylab import plot,show,xlabel,ylabel,legend,title,savefig,hold |
|---|
| 6 | except: |
|---|
| 7 | print 'Cannot import pylab plotting will not work. Csv files are still created' |
|---|
| 8 | from anuga.utilities.polygon import inside_polygon |
|---|
| 9 | from anuga.utilities.numerical_tools import ensure_numeric |
|---|
| 10 | from Scientific.IO.NetCDF import NetCDFFile |
|---|
| 11 | |
|---|
| 12 | def get_sts_gauge_data(filename,verbose=False): |
|---|
| 13 | from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\ |
|---|
| 14 | compress,zeros,fabs,take |
|---|
| 15 | fid = NetCDFFile(filename+'.sts', 'r') #Open existing file for read |
|---|
| 16 | permutation = fid.variables['permutation'][:] |
|---|
| 17 | x = fid.variables['x'][:]+fid.xllcorner #x-coordinates of vertices |
|---|
| 18 | y = fid.variables['y'][:]+fid.yllcorner #y-coordinates of vertices |
|---|
| 19 | points=transpose(asarray([x.tolist(),y.tolist()])) |
|---|
| 20 | time=fid.variables['time'][:]+fid.starttime |
|---|
| 21 | elevation=fid.variables['elevation'][:] |
|---|
| 22 | |
|---|
| 23 | basename='sts_gauge' |
|---|
| 24 | quantity_names=['stage','xmomentum','ymomentum'] |
|---|
| 25 | quantities = {} |
|---|
| 26 | for i, name in enumerate(quantity_names): |
|---|
| 27 | quantities[name] = fid.variables[name][:] |
|---|
| 28 | |
|---|
| 29 | maxname = 'max_sts_stage.csv' |
|---|
| 30 | fid_max = open(project.boundaries_dir+sep+maxname,'w') |
|---|
| 31 | s = 'x, y, max_stage \n' |
|---|
| 32 | fid_max.write(s) |
|---|
| 33 | for j in range(len(x)): |
|---|
| 34 | index= permutation[j] |
|---|
| 35 | stage = quantities['stage'][:,j] |
|---|
| 36 | xmomentum = quantities['xmomentum'][:,j] |
|---|
| 37 | ymomentum = quantities['ymomentum'][:,j] |
|---|
| 38 | |
|---|
| 39 | s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], max(stage)) |
|---|
| 40 | fid_max.write(s) |
|---|
| 41 | |
|---|
| 42 | fid_sts = open(project.boundaries_dir+sep+basename+'_'+ str(index)+'.csv', 'w') |
|---|
| 43 | s = 'time, stage, xmomentum, ymomentum \n' |
|---|
| 44 | fid_sts.write(s) |
|---|
| 45 | |
|---|
| 46 | for k in range(len(time)-1): |
|---|
| 47 | s = '%.6f, %.6f, %.6f, %.6f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k]) |
|---|
| 48 | fid_sts.write(s) |
|---|
| 49 | |
|---|
| 50 | fid_sts.close() |
|---|
| 51 | |
|---|
| 52 | fid.close() |
|---|
| 53 | |
|---|
| 54 | return quantities,elevation,time |
|---|
| 55 | |
|---|
| 56 | quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir,project.scenario_name),verbose=False) |
|---|
| 57 | |
|---|
| 58 | print len(elevation), len(quantities['stage'][0,:]) |
|---|
