[5576] | 1 | import os |
---|
| 2 | from os import sep |
---|
| 3 | #import create_boundary |
---|
| 4 | import project |
---|
| 5 | try: |
---|
| 6 | from pylab import plot,show,xlabel,ylabel,legend,title,savefig,hold |
---|
| 7 | except: |
---|
| 8 | print 'Cannot import pylab plotting will not work. Csv files are still created' |
---|
| 9 | from anuga.utilities.polygon import inside_polygon |
---|
| 10 | from anuga.utilities.numerical_tools import ensure_numeric |
---|
| 11 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 12 | |
---|
| 13 | |
---|
| 14 | def get_sts_gauge_data(filename,polygon,verbose=False): |
---|
| 15 | from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\ |
---|
| 16 | compress,zeros,fabs,take |
---|
| 17 | fid = NetCDFFile(filename+'.sts', 'r') #Open existing file for read |
---|
| 18 | x = fid.variables['x'][:]+fid.xllcorner #x-coordinates of vertices |
---|
| 19 | y = fid.variables['y'][:]+fid.yllcorner #y-coordinates of vertices |
---|
| 20 | print 'heloo', len(x), len(y) |
---|
| 21 | points=transpose(asarray([x.tolist(),y.tolist()])) |
---|
| 22 | time=fid.variables['time'][:] |
---|
| 23 | elevation=fid.variables['elevation'][:] |
---|
| 24 | indices=inside_polygon(points, polygon, closed=True, verbose=False) |
---|
| 25 | inside_points=take(points,indices,0) |
---|
| 26 | elevation=take(elevation,indices,0) |
---|
| 27 | if verbose: |
---|
| 28 | inside_points_array=ensure_numeric(inside_points)#for plotting only |
---|
| 29 | #plot(points[:,0],points[:,1],'o') |
---|
| 30 | #plot(inside_points_array[:,0][::5],inside_points_array[:,1][::5],'o') |
---|
| 31 | polygon=ensure_numeric(polygon) |
---|
| 32 | #plot(polygon[:,0],polygon[:,1]) |
---|
| 33 | #show() |
---|
| 34 | |
---|
| 35 | quantity_names=['stage','xmomentum','ymomentum'] |
---|
| 36 | |
---|
| 37 | fid_sts = open(project.boundaries_dir+sep+'sts_gauge.csv', 'w') |
---|
| 38 | s = 'time, stage, xmomentum, ymomentum \n' |
---|
| 39 | fid_sts.write(s) |
---|
| 40 | quantities = {} |
---|
| 41 | for i, name in enumerate(quantity_names): |
---|
| 42 | quantities[name] = fid.variables[name][:] |
---|
| 43 | print len(quantities[name]) |
---|
| 44 | if inside_points is not None: |
---|
| 45 | quantities[name] = take(quantities[name],indices,1) |
---|
| 46 | for j in range(len(x)): |
---|
| 47 | stage = quantities[name][:,j] |
---|
| 48 | xmomentum = quantities[name][:,j] |
---|
| 49 | ymomentum = quantities[name][:,j] |
---|
| 50 | #print 'hello', len(stage) |
---|
| 51 | #s = '%.2f, %.2f, %.2f, %.2f\n' %(time[j], stage, xmomentum, ymomentum) |
---|
| 52 | #fid_sts.write(s) |
---|
| 53 | else: |
---|
| 54 | msg = 'No gauges found inside polygon' |
---|
| 55 | raise msg |
---|
| 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(project.scenario_name,polygon,verbose=False) |
---|
| 62 | from Numeric import size |
---|
| 63 | print 'hello', len(quantities['stage']) |
---|
| 64 | |
---|
| 65 | assert len(points)==len(elevation)==len(quantities['stage'][0,:]) |
---|
| 66 | |
---|
| 67 | ##def plot_gauge(time,quantity,quantity_name,elevation,gauge_id,pos): |
---|
| 68 | ## figname=figname = '%s_%s%s.png' %('tgs/urs_gauge', |
---|
| 69 | ## quantity_name, |
---|
| 70 | ## str(gauge_id)) |
---|
| 71 | ## hold(False) |
---|
| 72 | ## plot(time/60.0,quantity,label='%s (elevation = %sm)'%('urs',str(elevation))) |
---|
| 73 | ## title('Gauge %d at (%f,%f)' %(gauge_id,pos[0],pos[1])) |
---|
| 74 | ## xlabel('Time (minutes)') |
---|
| 75 | ## ylabel(quantity_name) |
---|
| 76 | ## legend() |
---|
| 77 | ## print 'saving figure here %s' %figname |
---|
| 78 | ## savefig(figname) |
---|
| 79 | ## |
---|
| 80 | ##points=points[::5] |
---|
| 81 | ##elevation=elevation[::5] |
---|
| 82 | ##gauge_filename=project.boundaries_dir+sep+'urs_gauges' |
---|
| 83 | ##d=',' |
---|
| 84 | ##if os.path.exists(gauge_filename+ '.txt'): |
---|
| 85 | ## print 'file with list of gauges already created' |
---|
| 86 | ##else: |
---|
| 87 | ## fid = open(gauge_filename+ '.txt', 'w') |
---|
| 88 | ## fid.write('easting,northing,name,elevation\n') |
---|
| 89 | ## for i in range(len(points)): |
---|
| 90 | ## fid.write(str(points[i][0])+d+str(points[i][1])+d+str(i)+d+str(elevation[i])+'\n') |
---|
| 91 | ## fid.close() |
---|
| 92 | |
---|