source: anuga_work/production/perth/get_gauges.py @ 5576

Last change on this file since 5576 was 5576, checked in by kristy, 16 years ago

updated scripts for polyline boundary

File size: 3.6 KB
Line 
1import os
2from os import sep
3#import create_boundary
4import project
5try:
6    from pylab import plot,show,xlabel,ylabel,legend,title,savefig,hold
7except:
8    print 'Cannot import pylab plotting will not work. Csv files are still created'
9from anuga.utilities.polygon import inside_polygon
10from anuga.utilities.numerical_tools import ensure_numeric
11from Scientific.IO.NetCDF import NetCDFFile
12
13
14def 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
60polygon=project.poly_all
61points,quantities,elevation,time=get_sts_gauge_data(project.scenario_name,polygon,verbose=False)
62from Numeric import size
63print 'hello', len(quantities['stage'])
64
65assert 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
Note: See TracBrowser for help on using the repository browser.