source: anuga_work/production/dampier_2006/get_timeseries_and.py @ 4430

Last change on this file since 4430 was 4357, checked in by nick, 17 years ago

update to dampier

File size: 6.1 KB
RevLine 
[4308]1"""
2Generate images of "gauges" for production run
3
4Inputs:
5
6production dirs: dictionary of production directories with a
7                 association to that simulation run, eg high tide,
8                 low tide and MSL.
9                   
10Outputs:
11
12* figures used for report stored in the report_figure directory
13
14"""
15
16from os import sep, getcwd, access, F_OK, mkdir, getenv
17#from anuga.abstract_2d_finite_volumes.util import get_data_from_file
18import os
19from Numeric import zeros, array, allclose
20from os import getcwd, sep, altsep, mkdir, access, F_OK
21import project
22from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, sww2timeseries, get_data_from_file
23
24#makes the csv files from the evolved model
25
26
27home = getenv('INUNDATIONHOME') #Sandpit's parent dir   
28#user = get_user_name()
29data = 'data'
30state = 'western_australia'
31scenario_name = 'dampier.sww'
32scenario = 'dampier_tsunami_scenario_2006'
33#scenario = 'test_dampier'
34an = 'anuga'
35bo = 'boundaries'
36
[4357]37run_time = '20070311_225041_run'
38production_dirs = {run_time: 'URS evolved data'#,
39                   #'boundaries': 'URS boundary condition'
[4308]40                   }
41
42topo = 'topographies'
43out = 'outputs'
44urs = 'urs'
45gridded = '1_10000'
46
47#gauge_filename = os.path.join(home,data,state,scenario,an,out,scenario_name,'.sww')
48#sww_filename ={}
49#swwfiles = {}
50
51#sww_filename[os.path.join(home,data,state,scenario,an,out,run_time,scenario_name)] = run_time
52
53gauge_boundary_filename = 'gauges_time_series_near_top.csv'
54gauge_evolved_filename = 'gauges_time_series_near_top.csv'
55#gauge_boundary_filename = 'gauges_time_series_middle.csv'
56#gauge_evolved_filename = 'gauges_time_series_middle.csv'
57
58#gauge_boundary_filename = 'gauges_time_series_first.csv'
59#gauge_evolved_filename = 'gauges_time_series_first.csv'
60
61boundary_dir_filename = os.path.join(home,data,state,scenario,an,bo,gauge_boundary_filename)
62print'boundary_dir_filename',boundary_dir_filename
63out_dir=os.path.join(home,data,state,scenario,an,out,run_time)
64out_dir_name= out_dir+sep
65
66evolved_dir_filename= os.path.join(home,data,state,scenario,an,out,run_time,gauge_evolved_filename)
67
68#start_screen_catcher(out_dir_name, print_to_screen=True)
69start_screen_catcher(out_dir_name, extra_info='get_time', print_to_screen=True)
70print'boundary_dir_filename',boundary_dir_filename
71print'evolved_dir_filename',evolved_dir_filename
72
73#file_loc = project.output_dir + label_id + sep
74#swwfile = file_loc + project.scenario_name + '.sww'
75#swwfiles[swwfile] = label_id
76print "hella what the!"
77'''
78#get the timeseries for the evolved sww file
79texname, elev_output = sww2timeseries(sww_filename,
80                                      project.gauges_dir_name_simple,
81                                      production_dirs,
82                                      report = False,
83                                      #reportname = report_name,
84                                      plot_quantity = ['stage', 'momentum'],
85                                      surface = False,
86                                      time_min = None,
87                                      time_max = None,
88                                      title_on = False,
89                                      verbose = True)
90                                     
91'''                                     
92swwfiles = {}
93                                     
94for label_id in production_dirs.keys():
95       
96    if label_id == 'boundaries':
97        print 'boundaries'
98#        file_loc = project.boundaries_in_dir
[4311]99        swwfile = project.boundaries_dir_name3 + '.sww'
100#        swwfile = project.boundaries_dir_name6+'.sww'
[4308]101    else:
102#        file_loc = project.output_dir + label_id + sep
103        file_loc = out_dir_name
104        swwfile = file_loc + project.scenario_name + '.sww'
105    swwfiles[swwfile] = label_id
106    print"swwfiles",swwfiles
107       
108texname, elev_output = sww2timeseries(swwfiles,
109                                      project.gauges_dir_name_simple,
110                                      production_dirs,
111                                      report = False,
112                                      #reportname = report_name,
113                                      plot_quantity = ['stage', 'xmomentum', 'ymomentum'],
114#                                      plot_quantity = ['stage', 'momentum'],
115                                      surface = False,
116                                      time_min = None,
117                                      time_max = None,
118                                      title_on = False,
119#                                      use_cache = True,
120                                      verbose = True)
121
122#makes the csv files from the evolved model
123
124
125
[4311]126e_header, e_data = get_data_from_file(evolved_dir_filename)
[4308]127
[4311]128e_time = e_data[:,0]
129e_stage = e_data[:,1]
130e_momentum = e_data[:,2]
131e_speed = e_data[:,3]
132e_elevation = e_data[:,4]
133
[4308]134print'boundary_dir_filename',boundary_dir_filename
[4311]135b_header, b_data = get_data_from_file(boundary_dir_filename)
[4308]136
[4311]137b_time = b_data[:,0]
138b_stage = b_data[:,1]
139b_momentum = b_data[:,2]
140b_speed = b_data[:,3]
141b_elevation = b_data[:,4]
142
143
[4308]144# compares the 2 models
145j=0
146k=2
147b_sample = []
148e_sample = []
149for i in range(len(b_time)):
150#    if j<(len(e_time)) and b_time[i] == e_time[j]:
151    if j<(len(e_time)) and k<(len(e_time)) and b_time[i] == e_time[j]:
152        b_sample.append(float(b_stage[i]))
153        e_sample.append(float(e_stage[k]))
[4357]154#        print 'out:', b_time[i], b_stage[i],i, j, b_sample[j],e_stage[j]
[4308]155        if k <len(e_time): print 'time e equal b:', b_time[i], b_stage[i],i, j, b_sample[j],e_stage[j], e_stage[k],(len(e_time)-1)
156        j = j +1
157        k = k +1
158
159
[4357]160print len(b_sample), len(e_stage),e_stage,type(e_stage)
161e_stage = e_stage.tolist()
162e_stage.pop()
163#e_stage = e_stage.tolist().pop()
164e_stage.pop()
165#e_stage.pop()
[4308]166 
[4357]167print len(b_sample), len(e_stage),e_stage
[4308]168#assert allclose (b_sample, e_sample, 0.5, 0.5)
169
[4357]170assert allclose (b_sample, e_stage, 0.2, 0.2) 
[4308]171print "test successful" 
172
173
174
175
176                             
[4311]177                             
[4308]178                                         
Note: See TracBrowser for help on using the repository browser.