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

Last change on this file since 5607 was 5607, checked in by kristy, 16 years ago
File size: 2.1 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,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'][:]+fid.starttime
20    elevation=fid.variables['elevation'][:]
21       
22    basename='sts_gauge'
23    quantity_names=['stage','xmomentum','ymomentum']
24    quantities = {}
25    for i, name in enumerate(quantity_names):
26        quantities[name] = fid.variables[name][:]
27
28    maxname = 'max_sts_stage.csv'
29    fid_max = open(project.boundaries_dir+sep+maxname,'w')
30    s = 'x, y, max_stage \n'
31    fid_max.write(s)   
32    for j in range(len(x)):
33        locx=int(x[j])
34        locy=int(y[j])
35        stage = quantities['stage'][:,j]
36        xmomentum = quantities['xmomentum'][:,j]
37        ymomentum = quantities['ymomentum'][:,j]
38
39        s = '%.6f, %.6f, %.6f\n' %(x[j], y[j], max(stage))
40        fid_max.write(s)
41       
42        fid_sts = open(project.boundaries_dir+sep+basename+'_'+str(locx)+'_'+str(locy)+'.csv', 'w')
43        s = 'time, stage, xmomentum, ymomentum \n'
44        fid_sts.write(s)
45       
46        for k in range(len(time)-1):
47            s = '%.6f, %.6f, %.6f, %.6f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k])
48            fid_sts.write(s)
49
50        fid_sts.close()     
51
52    fid.close()
53       
54    return quantities,elevation,time
55
56quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir,project.scenario_name),verbose=False)
57
58print len(elevation), len(quantities['stage'][0,:])
Note: See TracBrowser for help on using the repository browser.