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

Last change on this file since 5580 was 5580, checked in by sexton, 14 years ago

update script to write time series to file from sts file (for Perth)

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