source: anuga_work/development/Hinwood_2008/slope.py @ 5494

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

Current Hinwood scenario - added plotting of froude number

File size: 14.1 KB
Line 
1
2"""
3Plot up files from the Hinwood project.
4"""
5from os import sep
6import project
7from copy import deepcopy
8#from scipy import arange
9from csv import writer
10
11from Numeric import arange, array, zeros, Float, where, greater, less, \
12     compress, argmin, choose
13
14from anuga.fit_interpolate.interpolate import interpolate_sww2csv
15from anuga.shallow_water.data_manager import csv2dict
16from anuga.utilities.numerical_tools import ensure_numeric
17
18
19SLOPE_STR = 'stage_slopes'
20TIME_STR = 'times'
21
22TIME_BORDER = 5
23LOCATION_BORDER = .5
24
25def load_sensors(quantity_file):
26    #slope, _ = csv2dict(file_sim)
27   
28    # Read the depth file
29    dfid = open(quantity_file)
30    lines = dfid.readlines()
31    dfid.close()
32
33    title = lines.pop(0)
34    n_time = len(lines)
35    n_sensors = len(lines[0].split(','))-1  # -1 to remove time
36    dtimes = zeros(n_time, Float)  #Time
37    depths = zeros(n_time, Float)  #
38    sensors = zeros((n_time,n_sensors), Float)
39    quantity_locations = title.split(',') #(',')
40    quantity_locations.pop(0) # remove 'time'
41   
42    locations = [float(j.split(':')[0]) for j in quantity_locations]
43   
44    for i, line in enumerate(lines):
45        fields = line.split(',') #(',')
46        fields = [float(j) for j in fields]
47        dtimes[i] = fields[0]
48        sensors[i] = fields[1:] # 1: to remove time
49
50    #print "dtimes",dtimes
51    #print "locations", locations
52    #print "sensors", sensors
53    return dtimes, locations, sensors
54   
55def load_slopes(stage_file):
56    """
57    Finds the slope, wrt distance of a distance, time, quantity csv file.
58
59    returns the times and slope_locations vectors and the slopes array.
60    """
61    times, locations, sensors = load_sensors(stage_file)
62    n_slope_locations = len(locations)-1
63    n_time = len(times)
64    slope_locations = zeros(n_slope_locations, Float)  #
65    slopes = zeros((n_time,n_slope_locations), Float)
66
67    # An array of the sensor spacing values
68    delta_locations = zeros(n_slope_locations, Float)
69   
70    for i in arange(n_slope_locations):
71        slope_locations[i] = (locations[i+1] + locations[i+1])/2.
72        delta_locations[i] = (locations[i+1] - locations[i])
73   
74    for j in arange(n_time):
75        for i in arange(n_slope_locations):
76            slopes[j,i] = (sensors[j,i+1] - sensors[j,i])/delta_locations[i]
77
78    return times, slope_locations, slopes
79
80
81def graph_contours(times, x_data, z_data,
82                 y_label='Time, seconds',
83                 plot_title="slope",
84                 x_label='x location, m',
85                 save_as=None,
86                 is_interactive=False,
87                 break_xs=None,
88                 break_times=None):
89    # Do not move these imports.  Tornado doesn't have pylab
90    from pylab import meshgrid, cm, contourf, contour, ion, plot, xlabel, \
91         ylabel, close, legend, savefig, title, figure ,colorbar, show , axis
92
93    origin = 'lower'
94
95    if is_interactive:
96        ion()
97       
98    # Can't seem to reshape this info once it is in the function
99    CS = contourf(x_data, times, z_data, 10,
100                  cmap=cm.bone,
101                  origin=origin)
102   
103    #CS2 = contour(x_data, times, z_data, CS.levels[::1],
104     #             colors = 'r',
105      #            origin=origin,
106       #           hold='on')
107   
108    title(plot_title)
109    xlabel(x_label)
110    ylabel(y_label)
111
112    if break_times is not None and break_xs is not None:
113        plot(break_xs, break_times, 'ro')
114   
115   
116    # Make a colorbar for the ContourSet returned by the contourf call.
117    cbar = colorbar(CS)
118
119    # Add the contour line levels to the colorbar
120    cbar.ax.set_ylabel('verbosity coefficient')
121    #cbar.add_lines(CS2)
122         
123    if is_interactive:       
124        raw_input() # Wait for enter pressed
125       
126    if save_as is not None:
127        savefig(save_as)
128    close()  #Need to close this plot
129
130def graph_froude(times, x_data, z_data,
131                 y_label='Time, seconds',
132                 plot_title="Froude Number",
133                 x_label='x location, m',
134                 save_as=None,
135                 is_interactive=False,
136                 break_xs=None,
137                 break_times=None):
138    # Do not move these imports.  Tornado doesn't have pylab
139    from pylab import meshgrid, cm, contourf, contour, ion, plot, xlabel, \
140         ylabel, close, legend, savefig, title, figure ,colorbar, show , axis
141
142    origin = 'lower'
143
144    if is_interactive:
145        ion()
146       
147    # Can't seem to reshape this info once it is in the function
148    CS = contourf(x_data, times, z_data, [-1,0.6,0.8,1,2,4],
149                  colors = ('black', 'r', 'g', 'b','r'),
150                  #cmap=cm.bone,
151                  origin=origin)
152   
153    #CS2 = contour(x_data, times, z_data, CS.levels[::1],
154     #             colors = 'r',
155      #            origin=origin,
156       #           hold='on')
157   
158    title(plot_title)
159    xlabel(x_label)
160    ylabel(y_label)
161
162    if break_times is not None and break_xs is not None:
163        plot(break_xs, break_times, 'yo')
164   
165   
166    # Make a colorbar for the ContourSet returned by the contourf call.
167    cbar = colorbar(CS)
168
169    # Add the contour line levels to the colorbar
170    cbar.ax.set_ylabel('verbosity coefficient')
171    #cbar.add_lines(CS2)
172         
173    if is_interactive:       
174        raw_input() # Wait for enter pressed
175       
176    if save_as is not None:
177        savefig(save_as)
178    close()  #Need to close this plot
179   
180def auto_graph_slopes(outputdir_tag, scenarios, is_interactive=False):
181    plot_type = ".pdf"
182    for run_data in scenarios:
183        id = run_data['scenario_id']
184        outputdir_name = id + outputdir_tag
185        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
186                                       outputdir_name=outputdir_name)
187        end = id + ".csv"
188        anuga_break_times = []
189        for break_time in run_data['break_times']:
190            anuga_break_times.append( \
191                break_time -  run_data['ANUGA_start_time'])
192           
193        stage_file = pro_instance.outputdir + "fslope_stage_" + end
194        save_as = pro_instance.plots_dir + sep + \
195                            outputdir_name + "_slope_stage" + plot_type
196        times, locations, slopes = load_slopes(stage_file)
197        times, slopes = get_band(anuga_break_times[0]-TIME_BORDER,
198                                   100, times, slopes, 0)
199        locations, slopes = get_band(
200            run_data['break_xs'][0]- 2*LOCATION_BORDER,
201            100, locations, slopes, -1)
202        graph_contours(times, locations, slopes,
203                     break_xs=run_data['break_xs'],
204                     break_times=anuga_break_times,
205                     save_as=save_as,
206                       is_interactive=is_interactive)
207
208def auto_graph_froudes(outputdir_tag, scenarios, is_interactive=False):
209   
210    plot_type = ".pdf"
211   
212    for run_data in scenarios:
213        id = run_data['scenario_id']
214        outputdir_name = id + outputdir_tag
215        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
216                                       outputdir_name=outputdir_name)
217        end = id + ".csv"
218        anuga_break_times = []
219        for break_time in run_data['break_times']:
220            anuga_break_times.append( \
221                break_time -  run_data['ANUGA_start_time'])
222           
223        froude_file = pro_instance.outputdir + "fslope_froude_" + end
224        save_as = pro_instance.plots_dir + sep + \
225                            outputdir_name + "_froude" + plot_type
226        dtimes, locations, sensors = load_sensors(froude_file)
227        dtimes, sensors = get_band(anuga_break_times[0]-TIME_BORDER,
228                                   100, dtimes, sensors, 0)
229        locations, sensors = get_band(run_data['break_xs'][0]-LOCATION_BORDER,
230                                      100, locations, sensors, -1)
231        #print "dtimes", dtimes
232        #print "sensors", sensors
233        #times, slope_locations, slopes = load_slopes(stage_file)
234        graph_froude(dtimes, locations, sensors,
235                     break_xs=run_data['break_xs'],
236                     break_times=anuga_break_times,
237                     save_as=save_as,
238                     is_interactive=is_interactive)
239       
240 
241def auto_find_min_slopes(outputdir_tag, scenarios):
242    """
243   
244    """
245   
246    for run_data in scenarios:
247        id = run_data['scenario_id']
248        outputdir_name = id + outputdir_tag
249        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
250                                       outputdir_name=outputdir_name)
251        end = id + ".csv"
252        anuga_break_times = []
253        for break_time in run_data['break_times']:
254            anuga_break_times.append( \
255                break_time -  run_data['ANUGA_start_time'])
256           
257        stage_file = pro_instance.outputdir + "fslope_stage_" + end
258       
259        times, slope_locations, slopes = load_slopes(stage_file)
260        waves = find_min_slopes(times, slope_locations, slopes,
261                                run_data['break_times'],
262                                anuga_break_times,
263                                run_data['band_offset'])
264        # write the wave info here
265        # return the keys actually.
266        for wave_key in waves.iterkeys():
267            wave_file = stage_file[:-3] + '_'+ wave_key + ".csv"
268            print "wave_file", wave_file
269            wave_writer = writer(file(wave_file, "wb"))
270            wave_writer.writerow(["x location", "min slope", "Time"])
271            wave_writer.writerows(map(None,
272                                      slope_locations,
273                                      waves[wave_key][SLOPE_STR],
274                                      waves[wave_key][TIME_STR]))
275           
276            #wave_writer.close()
277           
278           
279       
280   
281def gauges_for_slope(outputdir_tag, scenarios):
282    """
283    This is used to create a stage file, using gauges relivent to
284    finding a slope.
285    """
286    dx = 0.05
287    for run_data in scenarios:
288        point_x = arange(run_data['start_slope_x'],
289                        run_data['finish_slope_x'],
290                        dx).tolist()
291        flume_y_middle = 0.5
292        points = []
293        for gauge_x in point_x:
294            points.append([gauge_x, flume_y_middle])
295        id = run_data['scenario_id']
296   
297        basename = 'zz_' + run_data['scenario_id']
298        outputdir_name = id + outputdir_tag
299        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
300                                       outputdir_name=outputdir_name)
301        end = id + ".csv"
302        interpolate_sww2csv( \
303            pro_instance.outputdir + basename +".sww",
304            points,
305            pro_instance.outputdir + "fslope_depth_" + end,
306            pro_instance.outputdir + "fslope_velocity_x_" + end,
307            pro_instance.outputdir + "fslope_velocity_y_" + end,
308            pro_instance.outputdir + "fslope_stage_" + end,
309            pro_instance.outputdir + "fslope_froude_" + end,
310            time_thinning=1)
311
312def find_min_slopes(times, slope_locations, slopes,
313                    break_times, anuga_break_times, band_offset):
314    bands = break_times2bands(break_times, band_offset)
315
316    waves = {}
317    for i,_ in enumerate(bands[0:-1]):
318        id = "wave " + str(i)
319        max_q, max_q_times = get_min_in_band(bands[i], bands[i+1],
320                                             times, slopes)
321        waves[id] = {SLOPE_STR:max_q,
322                     TIME_STR:max_q_times}
323    return waves
324
325   
326def get_band(min, max, vector, quantity_array, axis):
327    """
328    Return a band of vector and quantity, within (not including) the
329    min, max.
330
331    For a time band, set the axis to 0.
332    For a location band, set the axis to -1.
333   
334    """
335
336    SMALL_MIN = -1e10  # Not that small, but small enough
337    vector = ensure_numeric(vector)
338    quantity_array = ensure_numeric(quantity_array)
339
340    assert min > SMALL_MIN
341    no_maxs = where(less(vector,max), vector, SMALL_MIN)
342    #print "no_maxs", no_maxs
343    band_condition = greater(no_maxs, min)
344    band_vector = compress(band_condition, vector, axis=axis)
345    #print "band_time", band_time
346    #print "quantity_array", quantity_array.shape
347    band_quantity = compress(band_condition, quantity_array, axis=axis)
348    return band_vector, band_quantity
349
350def get_min_in_band(min_time, max_time, time_vector, quantity_array):
351    """
352    given a quantity array, with the 2nd axis being
353    time, represented by the time_vector, find the minimum within
354    the time band.
355
356    Assumes times are positive
357    """
358   
359    time_vector = ensure_numeric(time_vector)
360    quantity_array = ensure_numeric(quantity_array)
361
362    band_time, band_quantity  = get_band(min_time, max_time,
363                                         time_vector, quantity_array, 0)
364    #print "band_quantity",band_quantity
365    max_quantity_indices = argmin(band_quantity, axis=0)
366    #print "max_quantity_indices", max_quantity_indices
367    max_quantity_times = choose(max_quantity_indices, band_time)
368    #print "max_quantity_times", max_quantity_times
369    max_quantities = choose(max_quantity_indices, band_quantity)
370    #print "max_quantities", max_quantities
371
372    return max_quantities, max_quantity_times
373
374def break_times2bands(break_times, band_offset):
375    """
376    Break_times is a list of times, ascending.
377    bands is a list of times, being the midpoints of break_times, with a min
378    and max band added.
379    """
380    assert len(break_times)>2
381   
382    bands = [] #deepcopy(break_times)
383    bands.append(break_times[0]-0.5*(break_times[1]-break_times[0]))
384
385    for i,break_x in enumerate(break_times[0:-1]):
386        bands.append(break_times[i]+0.5*(break_times[i+1]-break_times[i]))
387       
388    bands.append(break_times[-1]+0.5*(break_times[-1]-break_times[-2]))
389    bands = ensure_numeric(bands)
390    bands += band_offset
391    return bands
392   
393   
394   
395
396#-------------------------------------------------------------
397if __name__ == "__main__":
398    """
399    """
400    from scenarios import scenarios
401    #scenarios = [scenarios[0]]
402    outputdir_tag = "_good_tri_area_0.01_limiterC"
403    #outputdir_tag = "_test_limiterC"
404    #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!!
405    #gauges_for_slope(outputdir_tag, scenarios)
406    auto_graph_slopes(outputdir_tag, scenarios) #, is_interactive=True)
407    #auto_find_min_slopes(outputdir_tag, scenarios)
408    auto_graph_froudes(outputdir_tag, scenarios) 
Note: See TracBrowser for help on using the repository browser.