source: anuga_work/development/boxingday08/get_gauges.py @ 5599

Last change on this file since 5599 was 5553, checked in by jakeman, 16 years ago

fixed unit tests that check urs2sts for non-zero tide values

File size: 4.8 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
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    quantity_names=['stage','xmomentum','ymomentum']
33    quantities = {}
34    for i, name in enumerate(quantity_names):
35        quantities[name] = fid.variables[name][:]
36        if inside_points is not None:
37            quantities[name] = take(quantities[name],indices,1)
38        else:
39            msg = 'No gauges found inside polygon'
40            raise msg
41    fid.close()
42
43    return inside_points,quantities,elevation,time
44
45from anuga.shallow_water.data_manager import urs2sts
46base_name='tgs/gauges'
47if os.path.exists(base_name+'.sts'):
48    print 'sts file already created'
49else:
50    print 'creating sts file'
51    urs2sts(base_name,
52            basename_out=base_name+'.sts',
53            mean_stage=0.0, verbose=False)
54
55#polygon=create_boundary.bounding_polygon
56polygon=project.bounding_polygon
57points,quantities,elevation,time=get_sts_gauge_data(base_name,polygon,verbose=False)
58assert len(points)==len(elevation)==len(quantities['stage'][0,:])
59
60def plot_gauge(time,quantity,quantity_name,elevation,gauge_id,pos):
61    figname=figname = '%s_%s%s.png' %('tgs/urs_gauge',
62                                      quantity_name,
63                                      str(gauge_id))
64    hold(False)
65    plot(time/60.0,quantity,label='%s (elevation = %sm)'%('urs',str(elevation)))
66    title('Gauge %d at (%f,%f)' %(gauge_id,pos[0],pos[1]))
67    xlabel('Time (minutes)')
68    ylabel(quantity_name)
69    legend()
70    print 'saving figure here %s' %figname
71    savefig(figname)
72   
73points=points[::5]
74elevation=elevation[::5]
75gauge_filename='urs_gauges'
76d=','
77if os.path.exists(gauge_filename+ '.txt'):
78    print 'file with list of gauges already created'
79else:
80    fid = open(gauge_filename+ '.txt', 'w')
81    fid.write('easting,northing,name,elevation\n')
82    for i in range(len(points)):
83        fid.write(str(points[i][0])+d+str(points[i][1])+d+str(i)+d+str(elevation[i])+'\n')
84    fid.close()
85
86def get_sww_gauge_data(sww_filename,gauge_filename):
87    from anuga.abstract_2d_finite_volumes.util import sww2csv_gauges, csv2timeseries_graphs
88    dir_name, base = os.path.split(sww_filename)
89    if os.path.exists(dir_name+sep+'gauge_0.csv'):
90        print 'csv gauge files exist already'
91    else:
92        print 'writing gauge csv files to '+dir_name
93        sww2csv_gauges(sww_file=sww_filename+'.sww',
94                   quantities = ['stage','elevation','xmomentum','ymomentum'],
95                   gauge_file=gauge_filename+'.txt')
96       
97    print 'generating graphs'
98    csv2timeseries_graphs(directories_dic={dir_name:[dir_name,0, 0]},
99                          base_name='gauge_',
100                          plot_numbers='',#['1'],
101                          quantities=['stage','xmomentum','ymomentum'],
102                          extra_plot_name='',
103                          assess_all_csv_files=True,
104                          create_latex=False,
105                          verbose=True)
106
107sww_filename='good_simulation/good_polyline20080721_130608'
108#sww_filename='polyline/good_polyline20080607_205151'
109#sww_filename='most/good_most20080607_205101'
110#sww_filename='cloud/good_urs20080607_205134'
111#get_sww_gauge_data(sww_filename,gauge_filename)
112
113print 'plotting from sww file', sww_filename
114
115for i in range(len(points)):
116    quantity=quantities['stage'][:,i]
117    plot_gauge(time,quantity,'stage',elevation[i],i,points[i])
118    quantity=quantities['xmomentum'][:,i]
119    plot_gauge(time,quantity,'xmomentum',elevation[i],i,points[i])
120    #quantity=quantities['ymomentum'][:,i]
121    #plot_gauge(time,quantity,'ymomentum',elevation[i],i,points[i])
Note: See TracBrowser for help on using the repository browser.