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

Last change on this file since 5592 was 5592, checked in by kristy, 14 years ago
File size: 2.7 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
12def get_sts_gauge_data(filename,polygon,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    x = fid.variables['x'][:]+fid.xllcorner   #x-coordinates of vertices
17    y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
18    points=transpose(asarray([x.tolist(),y.tolist()]))
19    time=fid.variables['time'][:]
20    elevation=fid.variables['elevation'][:]
21    indices=inside_polygon(points, polygon, closed=True, verbose=False)
22    inside_points=take(points,indices,0)
23    elevation=take(elevation,indices,0)
24    if verbose:
25        inside_points_array=ensure_numeric(inside_points)#for plotting only
26        #plot(points[:,0],points[:,1],'o')
27        #plot(inside_points_array[:,0][::5],inside_points_array[:,1][::5],'o')
28        polygon=ensure_numeric(polygon)
29        #plot(polygon[:,0],polygon[:,1])
30        #show()
31
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##        else:
41##            msg = 'No gauges found inside polygon'
42##            raise msg
43##           
44    for j in range(len(x)-1):
45        locx=int(x[j])
46        locy=int(y[j])
47        fid_sts = open(project.boundaries_dir+sep+basename+'_'+str(locx)+'_'+str(locy)+'.csv', 'w')
48        s = 'time, stage, xmomentum, ymomentum \n'
49        fid_sts.write(s)
50        stage = quantities['stage'][:,j]
51        xmomentum = quantities['xmomentum'][:,j]
52        ymomentum = quantities['ymomentum'][:,j]
53        for k in range(len(time)-1):
54            s = '%.2f, %.2f, %.2f, %.2f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k])
55            fid_sts.write(s)
56
57        fid_sts.close()     
58
59    fid.close()
60       
61    return inside_points,quantities,elevation,time
62
63polygon=project.poly_all
64points,quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir,project.scenario_name),polygon,verbose=False)
65print quantities['stage'][0,:]
66print len(points), len(elevation), len(quantities['stage'][0,:])
67assert len(points)==len(elevation)==len(quantities['stage'][0,:])
Note: See TracBrowser for help on using the repository browser.