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

Last change on this file since 5495 was 5495, checked in by duncan, 14 years ago

Current Hinwood scenario

File size: 14.3 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        plot_title = "Stage slope " + id
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                       plot_title=plot_title,
204                       break_xs=run_data['break_xs'],
205                       break_times=anuga_break_times,
206                       save_as=save_as,
207                       is_interactive=is_interactive)
208
209def auto_graph_froudes(outputdir_tag, scenarios, is_interactive=False):
210   
211    plot_type = ".pdf"
212   
213    for run_data in scenarios:
214        id = run_data['scenario_id']
215        outputdir_name = id + outputdir_tag
216        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
217                                       outputdir_name=outputdir_name)
218        end = id + ".csv"
219        anuga_break_times = []
220        for break_time in run_data['break_times']:
221            anuga_break_times.append( \
222                break_time -  run_data['ANUGA_start_time'])
223        plot_title = "Froude Number" + id
224           
225        froude_file = pro_instance.outputdir + "fslope_froude_" + end
226        save_as = pro_instance.plots_dir + sep + \
227                            outputdir_name + "_froude" + plot_type
228        dtimes, locations, sensors = load_sensors(froude_file)
229        dtimes, sensors = get_band(anuga_break_times[0]-TIME_BORDER,
230                                   100, dtimes, sensors, 0)
231        locations, sensors = get_band(run_data['break_xs'][0]-LOCATION_BORDER,
232                                      100, locations, sensors, -1)
233        #print "dtimes", dtimes
234        #print "sensors", sensors
235        #times, slope_locations, slopes = load_slopes(stage_file)
236        graph_froude(dtimes, locations, sensors,
237                     plot_title=plot_title,
238                     break_xs=run_data['break_xs'],
239                     break_times=anuga_break_times,
240                     save_as=save_as,
241                     is_interactive=is_interactive)
242       
243 
244def auto_find_min_slopes(outputdir_tag, scenarios):
245    """
246   
247    """
248   
249    for run_data in scenarios:
250        id = run_data['scenario_id']
251        outputdir_name = id + outputdir_tag
252        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
253                                       outputdir_name=outputdir_name)
254        end = id + ".csv"
255        anuga_break_times = []
256        for break_time in run_data['break_times']:
257            anuga_break_times.append( \
258                break_time -  run_data['ANUGA_start_time'])
259           
260        stage_file = pro_instance.outputdir + "fslope_stage_" + end
261       
262        times, slope_locations, slopes = load_slopes(stage_file)
263        waves = find_min_slopes(times, slope_locations, slopes,
264                                run_data['break_times'],
265                                anuga_break_times,
266                                run_data['band_offset'])
267        # write the wave info here
268        # return the keys actually.
269        for wave_key in waves.iterkeys():
270            wave_file = stage_file[:-3] + '_'+ wave_key + ".csv"
271            print "wave_file", wave_file
272            wave_writer = writer(file(wave_file, "wb"))
273            wave_writer.writerow(["x location", "min slope", "Time"])
274            wave_writer.writerows(map(None,
275                                      slope_locations,
276                                      waves[wave_key][SLOPE_STR],
277                                      waves[wave_key][TIME_STR]))
278           
279            #wave_writer.close()
280           
281           
282       
283   
284def gauges_for_slope(outputdir_tag, scenarios):
285    """
286    This is used to create a stage file, using gauges relivent to
287    finding a slope.
288    """
289    dx = 0.05
290    for run_data in scenarios:
291        point_x = arange(run_data['start_slope_x'],
292                        run_data['finish_slope_x'],
293                        dx).tolist()
294        flume_y_middle = 0.5
295        points = []
296        for gauge_x in point_x:
297            points.append([gauge_x, flume_y_middle])
298        id = run_data['scenario_id']
299   
300        basename = 'zz_' + run_data['scenario_id']
301        outputdir_name = id + outputdir_tag
302        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
303                                       outputdir_name=outputdir_name)
304        end = id + ".csv"
305        interpolate_sww2csv( \
306            pro_instance.outputdir + basename +".sww",
307            points,
308            pro_instance.outputdir + "fslope_depth_" + end,
309            pro_instance.outputdir + "fslope_velocity_x_" + end,
310            pro_instance.outputdir + "fslope_velocity_y_" + end,
311            pro_instance.outputdir + "fslope_stage_" + end,
312            pro_instance.outputdir + "fslope_froude_" + end,
313            time_thinning=1)
314
315def find_min_slopes(times, slope_locations, slopes,
316                    break_times, anuga_break_times, band_offset):
317    bands = break_times2bands(break_times, band_offset)
318
319    waves = {}
320    for i,_ in enumerate(bands[0:-1]):
321        id = "wave " + str(i)
322        max_q, max_q_times = get_min_in_band(bands[i], bands[i+1],
323                                             times, slopes)
324        waves[id] = {SLOPE_STR:max_q,
325                     TIME_STR:max_q_times}
326    return waves
327
328   
329def get_band(min, max, vector, quantity_array, axis):
330    """
331    Return a band of vector and quantity, within (not including) the
332    min, max.
333
334    For a time band, set the axis to 0.
335    For a location band, set the axis to -1.
336   
337    """
338
339    SMALL_MIN = -1e10  # Not that small, but small enough
340    vector = ensure_numeric(vector)
341    quantity_array = ensure_numeric(quantity_array)
342
343    assert min > SMALL_MIN
344    no_maxs = where(less(vector,max), vector, SMALL_MIN)
345    #print "no_maxs", no_maxs
346    band_condition = greater(no_maxs, min)
347    band_vector = compress(band_condition, vector, axis=axis)
348    #print "band_time", band_time
349    #print "quantity_array", quantity_array.shape
350    band_quantity = compress(band_condition, quantity_array, axis=axis)
351    return band_vector, band_quantity
352
353def get_min_in_band(min_time, max_time, time_vector, quantity_array):
354    """
355    given a quantity array, with the 2nd axis being
356    time, represented by the time_vector, find the minimum within
357    the time band.
358
359    Assumes times are positive
360    """
361   
362    time_vector = ensure_numeric(time_vector)
363    quantity_array = ensure_numeric(quantity_array)
364
365    band_time, band_quantity  = get_band(min_time, max_time,
366                                         time_vector, quantity_array, 0)
367    #print "band_quantity",band_quantity
368    max_quantity_indices = argmin(band_quantity, axis=0)
369    #print "max_quantity_indices", max_quantity_indices
370    max_quantity_times = choose(max_quantity_indices, band_time)
371    #print "max_quantity_times", max_quantity_times
372    max_quantities = choose(max_quantity_indices, band_quantity)
373    #print "max_quantities", max_quantities
374
375    return max_quantities, max_quantity_times
376
377def break_times2bands(break_times, band_offset):
378    """
379    Break_times is a list of times, ascending.
380    bands is a list of times, being the midpoints of break_times, with a min
381    and max band added.
382    """
383    assert len(break_times)>2
384   
385    bands = [] #deepcopy(break_times)
386    bands.append(break_times[0]-0.5*(break_times[1]-break_times[0]))
387
388    for i,break_x in enumerate(break_times[0:-1]):
389        bands.append(break_times[i]+0.5*(break_times[i+1]-break_times[i]))
390       
391    bands.append(break_times[-1]+0.5*(break_times[-1]-break_times[-2]))
392    bands = ensure_numeric(bands)
393    bands += band_offset
394    return bands
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.