| 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 | |
|---|
| 13 | def get_sts_gauge_data(filename,polygon,verbose=False): |
|---|
| 14 | from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\ |
|---|
| 15 | compress,zeros,fabs,take |
|---|
| 16 | fid = NetCDFFile(filename+'.sts', 'r') #Open existing file for read |
|---|
| 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'][:] |
|---|
| 21 | elevation=fid.variables['elevation'][:] |
|---|
| 22 | indices=inside_polygon(points, polygon, closed=True, verbose=False) |
|---|
| 23 | inside_points=take(points,indices,0) |
|---|
| 24 | elevation=take(elevation,indices,0) |
|---|
| 25 | if verbose: |
|---|
| 26 | inside_points_array=ensure_numeric(inside_points)#for plotting only |
|---|
| 27 | #plot(points[:,0],points[:,1],'o') |
|---|
| 28 | #plot(inside_points_array[:,0][::5],inside_points_array[:,1][::5],'o') |
|---|
| 29 | polygon=ensure_numeric(polygon) |
|---|
| 30 | #plot(polygon[:,0],polygon[:,1]) |
|---|
| 31 | #show() |
|---|
| 32 | |
|---|
| 33 | basename='sts_gauge' |
|---|
| 34 | quantity_names=['stage','xmomentum','ymomentum'] |
|---|
| 35 | quantities = {} |
|---|
| 36 | for i, name in enumerate(quantity_names): |
|---|
| 37 | quantities[name] = fid.variables[name][:] |
|---|
| 38 | if inside_points is not None: |
|---|
| 39 | quantities[name] = take(quantities[name],indices,1) |
|---|
| 40 | for j in range(len(x)-1): |
|---|
| 41 | locx=int(x[j]) |
|---|
| 42 | locy=int(y[j]) |
|---|
| 43 | fid_sts = open(project.boundaries_dir+sep+basename+'_'+str(locx)+'_'+str(locy)+'.csv', 'w') |
|---|
| 44 | s = 'time, stage, xmomentum, ymomentum \n' |
|---|
| 45 | fid_sts.write(s) |
|---|
| 46 | stage = quantities[name][:,j] |
|---|
| 47 | xmomentum = quantities[name][:,j] |
|---|
| 48 | ymomentum = quantities[name][:,j] |
|---|
| 49 | for k in range(len(time)-1): |
|---|
| 50 | s = '%.2f, %.2f, %.2f, %.2f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k]) |
|---|
| 51 | fid_sts.write(s) |
|---|
| 52 | else: |
|---|
| 53 | msg = 'No gauges found inside polygon' |
|---|
| 54 | raise msg |
|---|
| 55 | fid_sts.close() |
|---|
| 56 | fid.close() |
|---|
| 57 | |
|---|
| 58 | return inside_points,quantities,elevation,time |
|---|
| 59 | |
|---|
| 60 | polygon=project.poly_all |
|---|
| 61 | points,quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir,project.scenario_name),polygon,verbose=False) |
|---|
| 62 | |
|---|
| 63 | |
|---|
| 64 | assert len(points)==len(elevation)==len(quantities['stage'][0,:]) |
|---|
