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

Last change on this file since 5402 was 5402, checked in by ole, 14 years ago

Project file and get_gauges to support boxing day validation scripts committed in changeset:5401

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