source: anuga_work/development/Hinwood_2008/plot.py @ 5698

Last change on this file since 5698 was 5698, checked in by duncan, 16 years ago

Hinwood - continuing fix for line overrunning graph in anuga report, plus words for the report

File size: 13.8 KB
Line 
1
2"""
3Plot up files from the Hinwood project.
4"""
5from os import sep
6import project
7from time import localtime, strftime
8
9from Numeric import compress
10
11from calc_norm import get_max_min_condition_array
12
13def plot_compare_csv(location_sim, file_sim, location_exp, file_exp,
14                     y_label,
15                     save_as=None,
16                     plot_title="",
17                     x_label='Time (s)',
18                     legend_sim='ANUGA simulation',
19                     legend_exp='Measured flume result',
20                     is_interactive=False,
21                     x_axis=None,
22                     y_axis=None):
23    """
24    """
25    from pylab import ion, plot, xlabel, ylabel, close, legend, \
26         savefig, title, axis, setp
27    from anuga.shallow_water.data_manager import csv2dict
28
29
30
31    # Load in the csv files and convert info from strings to floats
32    simulation, _ = csv2dict(file_sim)
33    experiment, _ = csv2dict(file_exp)
34    time_sim = [float(x) for x in simulation['time']]
35    quantity_sim = [float(x) for x in simulation[location_sim]]
36    #print "quantity_sim", quantity_sim
37    time_exp = [float(x) for x in experiment['Time']]
38    quantity_exp = [float(x) for x in experiment[location_exp]]
39
40    if is_interactive:
41        ion()
42   
43    l_sim, l_exp = plot(time_sim, quantity_sim, time_exp, quantity_exp)
44    setp(l_sim, color='r')
45    setp(l_exp, color='b')
46
47    # Add axis stuff and legend
48    xlabel(x_label)
49    ylabel(y_label)
50    # The order defines the label
51    #legend((legend_exp, legend_sim),'upper left')
52    legend((legend_sim, legend_exp),'upper left')
53    title(plot_title)
54    if x_axis is not None:
55        axis(xmin=x_axis[0], xmax=x_axis[1])
56   
57    if y_axis is not None:
58        axis(ymin=y_axis[0], ymax=y_axis[1])
59       
60    if is_interactive:
61        # Wait for enter pressed
62        raw_input()
63
64    if save_as is not None:
65        savefig(save_as)
66   
67    #Need to close this plot
68    close()
69   
70def Hinwood_files_locations(run_data, outputdir_tag, plot_type,
71                            quantity = "depth",
72                            y_location_tag=':0.0'):
73    """
74    run_data is a dictionary of data describing a Hinwood experiment
75    outputdir_tag a string at the end of an output dir; '_good_tri_area_0.01_A'
76    plot_type the file extension of the plot, eg '.pdf'
77    """
78    id = run_data['scenario_id']
79    outputdir_name = id + outputdir_tag
80    pro_instance = project.Project(['data','flumes','Hinwood_2008'],
81                                   outputdir_name=outputdir_name)
82   
83   
84    file_sim = pro_instance.outputdir + quantity + "_" + id + ".csv"
85    #print "file_exp",file_exp
86    file_exp = pro_instance.raw_data_dir + sep + id + 'pressfilt_exp_' \
87               + quantity + '.csv'
88    #print "file_sim", file_sim
89    location_sims = []
90    location_exps = []
91    save_as_list = []
92    for gauge_x in run_data['gauge_x']:
93        gauge_x = str(gauge_x)
94        location_sims.append(gauge_x + y_location_tag)
95        location_exps.append(gauge_x)
96        save_as_list.append(pro_instance.plots_dir + sep + \
97                            outputdir_name + "_" + quantity + "_" + \
98                            gauge_x + plot_type)
99    return file_exp, file_sim, location_sims, location_exps, outputdir_name, \
100           save_as_list
101
102def plot_exp_sim_comparision(scenarios, outputdir_tag,
103                             quantity = "stage",is_interactive=False,
104         y_location_tag=':0.0'):
105    plot_type = ".pdf"
106    plot_files = {} # Index scenario (T1R3) value, list of files
107    for run_data in scenarios:
108       
109        temp = Hinwood_files_locations(run_data, outputdir_tag,
110                                       plot_type, quantity,
111                                       y_location_tag=y_location_tag)
112                                   
113        file_exp, file_sim, location_sims, location_exps, outputdir_name, \
114                  save_as_list = temp
115        plot_files[run_data['scenario_id']] = save_as_list
116        print "file_exp",file_exp
117        print "run_data['scenario_id']", run_data['scenario_id']
118        #location_sims = [location_sims[0]]
119        #location_exps = [location_exps[0]]
120        #save_as_list = [save_as_list[0]]
121        for loc_sim, loc_exp, save_as, gauge, gauge_elev in \
122                map(None, location_sims,
123                    location_exps,
124                    save_as_list,
125                    run_data['gauge_x'],
126                    run_data['gauge_bed_elevation']):
127            time_date = strftime('plot date: %d/%m/%Y Time: %H:%M:%S',
128                                      localtime())
129            plot_title = "Scenario: " + outputdir_name + "\n" + \
130                         "X Gauge (m):" + str(gauge) + "    " + time_date
131
132            # When the shore gauges out of the are used, don't
133            # use the standard y axis scale
134            x_axis = run_data['wave_times']
135            if gauge < run_data['axis_maximum_x']:
136                y_axis = [run_data['axis'][2],run_data['axis'][3]]
137       
138            else:
139                y_axis = [run_data['axis'][2] + gauge_elev,
140                          run_data['axis'][3] + gauge_elev]
141            print "Doing ", plot_title
142            plot_compare_csv(location_sim=loc_sim,
143                             file_sim=file_sim,
144                             location_exp=loc_exp,
145                             file_exp=file_exp,
146                             plot_title=plot_title,
147                             y_label='Water '+ quantity +' (m)',
148                             is_interactive=is_interactive,
149                             save_as=save_as,
150                             x_axis=x_axis,
151                             y_axis=y_axis)
152            if is_interactive is True:
153                break
154    return plot_files
155
156
157def auto_plot_compare_csv_subfigures(scenarios, outputdir_tag,
158                                     to_publish_indexes,
159                                     quantity = "stage",is_interactive=False,
160                                     y_location_tag=':0.0',
161                                     add_run_info=False):
162    """
163    Used by validation graphs to produce the final figures
164    """
165    plot_type = ".pdf"
166    save_as_list = []
167    for run_data in scenarios:
168       
169        id = run_data['scenario_id']
170        if not to_publish_indexes.has_key(id):
171            continue
172           
173        outputdir_name = id + outputdir_tag
174        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
175                                       outputdir_name=outputdir_name)
176   
177   
178        file_sim = pro_instance.outputdir + quantity + "_" + id + ".csv"
179        file_exp = pro_instance.raw_data_dir + sep + id + 'pressfilt_exp_' \
180                   + quantity + '.csv'
181
182        if add_run_info is True:
183            plot_title = "Title will not be in final draft \n" + \
184                         id + outputdir_tag
185        else:
186            plot_title = ""
187        #save_as = pro_instance.plots_dir + sep + \
188        #          outputdir_name + "_" + quantity + "_" + plot_type
189        save_as = pro_instance.plots_dir + sep + 'S' + id[1:2] + \
190                  '-' + quantity + '-compare' + plot_type
191        save_as_list.append(save_as)
192
193               
194        plot_compare_csv_subfigures(file_sim,
195                                    file_exp,
196                                    to_publish_indexes[id],
197                                    run_data,
198                                    plot_title=plot_title,
199                                    save_as=save_as,
200                                    y_label='Water '+ quantity +' (m)',
201                                    is_interactive=is_interactive,
202                                    y_location_tag=y_location_tag)
203    return save_as_list
204
205def plot_compare_csv_subfigures(file_sim,
206                                file_exp,
207                                gauge_indexs,
208                                run_data,
209                                save_as=None,
210                                plot_title="",
211                                x_label='Time (s)',
212                                y_label=None,
213                                legend_sim='ANUGA simulation',
214                                legend_exp='Measured flume result',
215                                is_interactive=False,
216                                x_axis=None,
217                                y_axis=None,
218                                y_location_tag=':0.0'):
219    """
220    Used by validation graphs to produce the final figures
221    """
222    from pylab import ion, plot, xlabel, ylabel, close, legend, \
223         savefig, title, axis, setp, subplot, grid, figlegend, gca, \
224         text, sqrt
225    import pylab
226   
227    from anuga.shallow_water.data_manager import csv2dict
228   
229    print "file_sim", file_sim
230
231    if False:
232        fig_width_pt = 246.0  # Get this from LaTeX using \showthe\columnwidth
233        inches_per_pt = 1.0/72.27               # Convert pt to inch
234        golden_mean = (sqrt(5)-1.0)/2.0         # Aesthetic ratio
235        fig_width = fig_width_pt*inches_per_pt  # width in inches
236        fig_height = fig_width*golden_mean      # height in inches
237        fig_size =  [fig_width,fig_height]
238        print "fig_size", fig_size
239        params = {'backend': 'ps',
240                  'axes.labelsize': 10,
241                  'text.fontsize': 10,
242                  'legend.fontsize': 10,
243                  'xtick.labelsize': 8,
244                  'ytick.labelsize': 8,
245                  'text.usetex': True,
246                  'figure.figsize': [3.4, 2.1]}
247    if False:
248        params = {'backend': 'ps',
249                  'axes.labelsize': 10,
250                  'text.fontsize': 10,
251                  'legend.fontsize': 10,
252                  'xtick.labelsize': 8,
253                  'ytick.labelsize': 8,
254                  'figure.figsize': [3.4, 3.1]}
255        pylab.rcParams.update(params)
256
257
258   
259    # Load in the csv files and convert info from strings to floats
260    simulation, _ = csv2dict(file_sim)
261    experiment, _ = csv2dict(file_exp)
262    time_sim = [float(x) for x in simulation['time']]
263    time_exp = [float(x) for x in experiment['Time']]
264
265
266    # Trim the simulation data, due to strange anuga paper bug
267    condition_sim = get_max_min_condition_array(run_data['wave_times'][0],
268                                                run_data['wave_times'][1],
269                                                time_sim)
270    time_sim = compress(condition_sim, time_sim)
271    condition_exp = get_max_min_condition_array(run_data['wave_times'][0],
272                                                run_data['wave_times'][1],
273                                                time_exp)
274    time_exp = compress(condition_exp, time_exp)
275   
276   
277    if is_interactive:
278        ion()
279       
280    for position, i in enumerate(gauge_indexs):
281        gauge_x = run_data['gauge_x'][i]
282       
283        grid_position = (len(gauge_indexs))*100 + 10 + position +1
284        subplot(grid_position)
285        location_sim = str(gauge_x) + y_location_tag
286        location_exp = str(gauge_x)
287       
288        quantity_sim = [float(x) for x in simulation[location_sim]]
289        quantity_exp = [float(x) for x in experiment[location_exp]]
290
291        # Trim the simulation data, due to strange anuga paper bug
292        quantity_sim = compress(condition_sim, quantity_sim)
293        quantity_exp = compress(condition_exp, quantity_exp)
294
295   
296        l_sim, l_exp = plot(time_sim, quantity_sim, time_exp, quantity_exp)
297        setp(l_sim, color='k', linestyle='--')
298        setp(l_exp, color='k')
299        grid(True)
300       
301        # When the shore gauges are used, don't
302        # use the standard y axis scale
303        gauge_elev = run_data['gauge_bed_elevation'][i]
304        x_axis = run_data['wave_times']
305        if gauge_x < run_data['axis_maximum_x']:
306            y_axis = [run_data['axis'][2],run_data['axis'][3]]   
307            setp(gca(), yticks=[-0.04,-0.02, 0.0, 0.02, 0.04])
308            y_text = -0.035
309        else:
310            y_axis = [run_data['axis'][2] + gauge_elev,
311                      run_data['axis'][3] + gauge_elev]   
312            setp(gca(), yticks=[0.02, 0.04, 0.06, 0.08, 0.1])
313            y_text = 0.025
314        text(x_axis[0]+5,y_text, str(round(gauge_x,1)) + " m gauge")   
315        if position == 0:
316            title(plot_title)
317        if position == 1:
318            ylabel(y_label)
319        #legend((legend_sim, legend_exp),'lower center')
320        if y_axis is not None:
321            axis(ymin=y_axis[0]-0.001, ymax=y_axis[1]-0.001)
322            # the fudge -0.001 is to reduce the ytick's
323
324        axis(xmin=x_axis[0], xmax=x_axis[1])
325
326        # Only do tick marks on the final graph
327        if not position == 2:
328            setp(gca(), xticklabels=[])
329
330    # Add axis stuff and legend
331    xlabel(x_label)
332       
333    #ylabel(y_label)
334    # The order defines the label
335    #legend((legend_exp, legend_sim),'upper left')
336
337    # I couldn't get this working
338    #figlegend((l_sim, l_exp),
339     #         (legend_sim, legend_exp),
340      #        'lower center', axespad = -0.05)
341   
342    if is_interactive:
343        # Wait for enter pressed
344        raw_input()
345
346    if save_as is not None:
347        savefig(save_as)
348   
349    #Need to close this plot
350    close()
351   
352#-------------------------------------------------------------
353if __name__ == "__main__":
354    """ Plot the stage graph for the ANUGA validation papar
355    """
356    from scenarios import scenarios
357    outputdir_tag = "_good_tri_area_0.01_limiterE"
358    outputdir_tag = "_nolmts_wdth_0.1_z_0.012_ys_0.01_mta_0.01_H"
359    outputdir_tag = "_good_lmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01_F"
360    outputdir_tag = "_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01_G"
361    #outputdir_tag = "_nolmts_wdth_0.01_z_0.012_ys_0.01_mta_1e-05_G"
362    #outputdir_tag = "_test_C"
363    #scenarios = scenarios[1:]
364    #scenarios = [scenarios[7]]
365    is_interactive = False
366    #is_interactive = True
367    plot_exp_sim_comparision(scenarios, outputdir_tag,
368                             is_interactive=is_interactive,
369                             y_location_tag=':0.0')
370
Note: See TracBrowser for help on using the repository browser.