| 1 | import os |
|---|
| 2 | from os import sep |
|---|
| 3 | #import create_boundary |
|---|
| 4 | import project |
|---|
| 5 | try: |
|---|
| 6 | from pylab import plot,show,xlabel,ylabel,legend,title,savefig,hold |
|---|
| 7 | except: |
|---|
| 8 | print 'Cannot import pylab plotting will not work. Csv files are still created' |
|---|
| 9 | from anuga.utilities.polygon import inside_polygon |
|---|
| 10 | from anuga.utilities.numerical_tools import ensure_numeric |
|---|
| 11 | from Scientific.IO.NetCDF import NetCDFFile |
|---|
| 12 | def 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 | |
|---|
| 45 | from anuga.shallow_water.data_manager import urs2sts |
|---|
| 46 | base_name='tgs/gauges' |
|---|
| 47 | if os.path.exists(base_name+'.sts'): |
|---|
| 48 | print 'sts file already created' |
|---|
| 49 | else: |
|---|
| 50 | print 'creating sts file' |
|---|
| 51 | urs2sts(base_name,mean_stage=0.0,verbose=False) |
|---|
| 52 | |
|---|
| 53 | #polygon=create_boundary.bounding_polygon |
|---|
| 54 | polygon=project.bounding_polygon |
|---|
| 55 | points,quantities,elevation,time=get_sts_gauge_data(base_name,polygon,verbose=False) |
|---|
| 56 | assert len(points)==len(elevation)==len(quantities['stage'][0,:]) |
|---|
| 57 | |
|---|
| 58 | def 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 | |
|---|
| 71 | points=points[::5] |
|---|
| 72 | elevation=elevation[::5] |
|---|
| 73 | gauge_filename='urs_gauges' |
|---|
| 74 | d=',' |
|---|
| 75 | if os.path.exists(gauge_filename+ '.txt'): |
|---|
| 76 | print 'file with list of gauges already created' |
|---|
| 77 | else: |
|---|
| 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 | |
|---|
| 84 | def 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 | |
|---|
| 105 | sww_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) |
|---|
| 109 | for 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]) |
|---|
