Changeset 5494 for anuga_work/development/Hinwood_2008/slope.py

Timestamp:

Jul 11, 2008, 4:15:59 PM (16 years ago)

Author:

duncan

Message:

Current Hinwood scenario - added plotting of froude number

File:

• anuga_work/development/Hinwood_2008/slope.py (modified) (7 diffs)

anuga_work/development/Hinwood_2008/slope.py

@@ -5 +5 @@
 from os import sep
 import project
+from copy import deepcopy
 #from scipy import arange
-
-from Numeric import arange, array, zeros, Float
+from csv import writer
+
+from Numeric import arange, array, zeros, Float, where, greater, less, \
+     compress, argmin, choose
 
 from anuga.fit_interpolate.interpolate import interpolate_sww2csv
 from anuga.shallow_water.data_manager import csv2dict
-
-def load_depths(slope_file):
+from anuga.utilities.numerical_tools import ensure_numeric
+
+
+SLOPE_STR = 'stage_slopes'
+TIME_STR = 'times'
+
+TIME_BORDER = 5
+LOCATION_BORDER = .5
+
+def load_sensors(quantity_file):
     #slope, _ = csv2dict(file_sim)
     
     # Read the depth file
-    dfid = open(slope_file)
+    dfid = open(quantity_file)
     lines = dfid.readlines()
     dfid.close()
@@ -26 +37 @@
     depths = zeros(n_time, Float)  #
     sensors = zeros((n_time,n_sensors), Float)
-    depth_locations = title.split(',') #(',')
-    depth_locations.pop(0) # remove 'time'
-    # !!!! -3 assumes a 
-    depths = [float(j.split(':')[0]) for j in depth_locations]
+    quantity_locations = title.split(',') #(',')
+    quantity_locations.pop(0) # remove 'time'
+    
+    locations = [float(j.split(':')[0]) for j in quantity_locations]
     
     for i, line in enumerate(lines):
@@ -38 +49 @@
 
     #print "dtimes",dtimes
-    #print "depths", depths
+    #print "locations", locations
     #print "sensors", sensors
-    return dtimes, depths, sensors
+    return dtimes, locations, sensors
    
-def load_slopes(slope_file):
-    times, depths, sensors = load_depths(slope_file)
-    n_slope_locations = len(depths)-1
+def load_slopes(stage_file):
+    """
+    Finds the slope, wrt distance of a distance, time, quantity csv file.
+
+    returns the times and slope_locations vectors and the slopes array.
+    """
+    times, locations, sensors = load_sensors(stage_file)
+    n_slope_locations = len(locations)-1
     n_time = len(times)
     slope_locations = zeros(n_slope_locations, Float)  #
@@ -50 +66 @@
 
     # An array of the sensor spacing values
-    delta_depths = zeros(n_slope_locations, Float)
+    delta_locations = zeros(n_slope_locations, Float)
     
     for i in arange(n_slope_locations):
-        slope_locations[i] = (depths[i+1] + depths[i+1])/2.
-        delta_depths[i] = (depths[i+1] - depths[i])
+        slope_locations[i] = (locations[i+1] + locations[i+1])/2.
+        delta_locations[i] = (locations[i+1] - locations[i])
     
     for j in arange(n_time):
         for i in arange(n_slope_locations):
-            slopes[j,i] = (sensors[j,i+1] - sensors[j,i])/delta_depths[i]
+            slopes[j,i] = (sensors[j,i+1] - sensors[j,i])/delta_locations[i]
 
     return times, slope_locations, slopes
 
-def graph_slopes(slope_file, break_xs=None, break_times=None):
+
+def graph_contours(times, x_data, z_data,
+                 y_label='Time, seconds',
+                 plot_title="slope",
+                 x_label='x location, m',
+                 save_as=None,
+                 is_interactive=False,
+                 break_xs=None,
+                 break_times=None):
+    # Do not move these imports.  Tornado doesn't have pylab
     from pylab import meshgrid, cm, contourf, contour, ion, plot, xlabel, \
          ylabel, close, legend, savefig, title, figure ,colorbar, show , axis
+
     origin = 'lower'
-    times, slope_locations, slopes = load_slopes(slope_file)
-    #print "times", times
-    #print "slope_locations", slope_locations
-    #print "slopes", slopes
-
+
+    if is_interactive:
+        ion()
+        
     # Can't seem to reshape this info once it is in the function
-    CS = contourf(slope_locations, times, slopes, 10, # [-1, -0.1, 0, 0.1],
-                  #alpha=0.5,
+    CS = contourf(x_data, times, z_data, 10,
                   cmap=cm.bone,
                   origin=origin)
     
-    #CS2 = contour(slope_locations, times, slopes, CS.levels[::2],
-       #           colors = 'r',
-        #          origin=origin,
-         #         hold='on')
-
-    title('slope')
-    xlabel('x location')
-    ylabel('Time, seconds')
-    #axis([5.0, 5.5, 30, 60])
-
+    #CS2 = contour(x_data, times, z_data, CS.levels[::1],
+     #             colors = 'r',
+      #            origin=origin,
+       #           hold='on')
+    
+    title(plot_title)
+    xlabel(x_label)
+    ylabel(y_label)
 
     if break_times is not None and break_xs is not None:
         plot(break_xs, break_times, 'ro')
     
+    
     # Make a colorbar for the ContourSet returned by the contourf call.
     cbar = colorbar(CS)
-    cbar.ax.set_ylabel('slope')
+
     # Add the contour line levels to the colorbar
+    cbar.ax.set_ylabel('verbosity coefficient')
     #cbar.add_lines(CS2)
-
-    figure()
-    show()
-
-def auto_graph_slopes():
-    from scenarios import scenarios
-    
-    outputdir_tag = "_good_tri_area_0.01_A"
-    outputdir_tag = "_good_tri_area_0.01_old"
-    #outputdir_tag = "_test"
-    scenarios = [scenarios[1]] # !!!!!!!!!!!!!!!!!!!!!!
+         
+    if is_interactive:       
+        raw_input() # Wait for enter pressed
+        
+    if save_as is not None:
+        savefig(save_as)
+    close()  #Need to close this plot
+
+def graph_froude(times, x_data, z_data,
+                 y_label='Time, seconds',
+                 plot_title="Froude Number",
+                 x_label='x location, m',
+                 save_as=None,
+                 is_interactive=False,
+                 break_xs=None,
+                 break_times=None):
+    # Do not move these imports.  Tornado doesn't have pylab
+    from pylab import meshgrid, cm, contourf, contour, ion, plot, xlabel, \
+         ylabel, close, legend, savefig, title, figure ,colorbar, show , axis
+
+    origin = 'lower'
+
+    if is_interactive:
+        ion()
+        
+    # Can't seem to reshape this info once it is in the function
+    CS = contourf(x_data, times, z_data, [-1,0.6,0.8,1,2,4],
+                  colors = ('black', 'r', 'g', 'b','r'),
+                  #cmap=cm.bone,
+                  origin=origin)
+    
+    #CS2 = contour(x_data, times, z_data, CS.levels[::1],
+     #             colors = 'r',
+      #            origin=origin,
+       #           hold='on')
+    
+    title(plot_title)
+    xlabel(x_label)
+    ylabel(y_label)
+
+    if break_times is not None and break_xs is not None:
+        plot(break_xs, break_times, 'yo')
+    
+    
+    # Make a colorbar for the ContourSet returned by the contourf call.
+    cbar = colorbar(CS)
+
+    # Add the contour line levels to the colorbar
+    cbar.ax.set_ylabel('verbosity coefficient')
+    #cbar.add_lines(CS2)
+         
+    if is_interactive:       
+        raw_input() # Wait for enter pressed
+        
+    if save_as is not None:
+        savefig(save_as)
+    close()  #Need to close this plot
+    
+def auto_graph_slopes(outputdir_tag, scenarios, is_interactive=False):
+    plot_type = ".pdf"
     for run_data in scenarios:
         id = run_data['scenario_id']
@@ -113 +186 @@
                                        outputdir_name=outputdir_name)
         end = id + ".csv"
-        slope_file = pro_instance.outputdir + "bslope_stage_" + end
-        graph_slopes(slope_file, run_data['break_xs'],
-                     run_data['break_times']) 
+        anuga_break_times = []
+        for break_time in run_data['break_times']:
+            anuga_break_times.append( \
+                break_time -  run_data['ANUGA_start_time'])
+            
+        stage_file = pro_instance.outputdir + "fslope_stage_" + end
+        save_as = pro_instance.plots_dir + sep + \
+                            outputdir_name + "_slope_stage" + plot_type
+        times, locations, slopes = load_slopes(stage_file)
+        times, slopes = get_band(anuga_break_times[0]-TIME_BORDER,
+                                   100, times, slopes, 0)
+        locations, slopes = get_band(
+            run_data['break_xs'][0]- 2*LOCATION_BORDER,
+            100, locations, slopes, -1)
+        graph_contours(times, locations, slopes,
+                     break_xs=run_data['break_xs'],
+                     break_times=anuga_break_times,
+                     save_as=save_as,
+                       is_interactive=is_interactive)
+
+def auto_graph_froudes(outputdir_tag, scenarios, is_interactive=False):
+    
+    plot_type = ".pdf"
+    
+    for run_data in scenarios:
+        id = run_data['scenario_id']
+        outputdir_name = id + outputdir_tag
+        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
+                                       outputdir_name=outputdir_name)
+        end = id + ".csv"
+        anuga_break_times = []
+        for break_time in run_data['break_times']:
+            anuga_break_times.append( \
+                break_time -  run_data['ANUGA_start_time'])
+            
+        froude_file = pro_instance.outputdir + "fslope_froude_" + end
+        save_as = pro_instance.plots_dir + sep + \
+                            outputdir_name + "_froude" + plot_type
+        dtimes, locations, sensors = load_sensors(froude_file)
+        dtimes, sensors = get_band(anuga_break_times[0]-TIME_BORDER,
+                                   100, dtimes, sensors, 0)
+        locations, sensors = get_band(run_data['break_xs'][0]-LOCATION_BORDER,
+                                      100, locations, sensors, -1)
+        #print "dtimes", dtimes
+        #print "sensors", sensors
+        #times, slope_locations, slopes = load_slopes(stage_file)
+        graph_froude(dtimes, locations, sensors,
+                     break_xs=run_data['break_xs'],
+                     break_times=anuga_break_times,
+                     save_as=save_as,
+                     is_interactive=is_interactive)
+        
+ 
+def auto_find_min_slopes(outputdir_tag, scenarios):
+    """
+    
+    """
+    
+    for run_data in scenarios:
+        id = run_data['scenario_id']
+        outputdir_name = id + outputdir_tag
+        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
+                                       outputdir_name=outputdir_name)
+        end = id + ".csv"
+        anuga_break_times = []
+        for break_time in run_data['break_times']:
+            anuga_break_times.append( \
+                break_time -  run_data['ANUGA_start_time'])
+            
+        stage_file = pro_instance.outputdir + "fslope_stage_" + end
+        
+        times, slope_locations, slopes = load_slopes(stage_file)
+        waves = find_min_slopes(times, slope_locations, slopes,
+                                run_data['break_times'],
+                                anuga_break_times,
+                                run_data['band_offset'])
+        # write the wave info here
+        # return the keys actually.
+        for wave_key in waves.iterkeys():
+            wave_file = stage_file[:-3] + '_'+ wave_key + ".csv"
+            print "wave_file", wave_file
+            wave_writer = writer(file(wave_file, "wb"))
+            wave_writer.writerow(["x location", "min slope", "Time"])
+            wave_writer.writerows(map(None,
+                                      slope_locations,
+                                      waves[wave_key][SLOPE_STR],
+                                      waves[wave_key][TIME_STR]))
+            
+            #wave_writer.close()
+            
+            
+        
     
 def gauges_for_slope(outputdir_tag, scenarios):
-    
-    outputdir_tag = "_good_tri_area_0.01_A"
-    #outputdir_tag = "_test"
-    #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!!
+    """
+    This is used to create a stage file, using gauges relivent to
+    finding a slope.
+    """
     dx = 0.05
     for run_data in scenarios:
@@ -138 +300 @@
                                        outputdir_name=outputdir_name)
         end = id + ".csv"
-        interpolate_sww2csv(pro_instance.outputdir + basename +".sww",
-                            points,
-                            pro_instance.outputdir + "fslope_depth_" + end,
-                            pro_instance.outputdir + "fslope_velocity_x_" + end,
-                            pro_instance.outputdir + "fslope_velocity_y_" + end,
-                            pro_instance.outputdir + "fslope_stage_" + end,
-                            time_thinning=1)
-  
+        interpolate_sww2csv( \
+            pro_instance.outputdir + basename +".sww",
+            points,
+            pro_instance.outputdir + "fslope_depth_" + end,
+            pro_instance.outputdir + "fslope_velocity_x_" + end,
+            pro_instance.outputdir + "fslope_velocity_y_" + end,
+            pro_instance.outputdir + "fslope_stage_" + end,
+            pro_instance.outputdir + "fslope_froude_" + end,
+            time_thinning=1)
+
+def find_min_slopes(times, slope_locations, slopes,
+                    break_times, anuga_break_times, band_offset):
+    bands = break_times2bands(break_times, band_offset)
+
+    waves = {}
+    for i,_ in enumerate(bands[0:-1]):
+        id = "wave " + str(i)
+        max_q, max_q_times = get_min_in_band(bands[i], bands[i+1],
+                                             times, slopes)
+        waves[id] = {SLOPE_STR:max_q,
+                     TIME_STR:max_q_times}
+    return waves
+
+    
+def get_band(min, max, vector, quantity_array, axis):
+    """
+    Return a band of vector and quantity, within (not including) the
+    min, max.
+
+    For a time band, set the axis to 0.
+    For a location band, set the axis to -1.
+    
+    """
+
+    SMALL_MIN = -1e10  # Not that small, but small enough
+    vector = ensure_numeric(vector)
+    quantity_array = ensure_numeric(quantity_array)
+
+    assert min > SMALL_MIN
+    no_maxs = where(less(vector,max), vector, SMALL_MIN)
+    #print "no_maxs", no_maxs
+    band_condition = greater(no_maxs, min)
+    band_vector = compress(band_condition, vector, axis=axis)
+    #print "band_time", band_time
+    #print "quantity_array", quantity_array.shape
+    band_quantity = compress(band_condition, quantity_array, axis=axis)
+    return band_vector, band_quantity
+
+def get_min_in_band(min_time, max_time, time_vector, quantity_array):
+    """
+    given a quantity array, with the 2nd axis being
+    time, represented by the time_vector, find the minimum within
+    the time band.
+
+    Assumes times are positive
+    """
+    
+    time_vector = ensure_numeric(time_vector)
+    quantity_array = ensure_numeric(quantity_array)
+
+    band_time, band_quantity  = get_band(min_time, max_time,
+                                         time_vector, quantity_array, 0)
+    #print "band_quantity",band_quantity
+    max_quantity_indices = argmin(band_quantity, axis=0)
+    #print "max_quantity_indices", max_quantity_indices
+    max_quantity_times = choose(max_quantity_indices, band_time)
+    #print "max_quantity_times", max_quantity_times
+    max_quantities = choose(max_quantity_indices, band_quantity)
+    #print "max_quantities", max_quantities
+
+    return max_quantities, max_quantity_times
+
+def break_times2bands(break_times, band_offset):
+    """
+    Break_times is a list of times, ascending.
+    bands is a list of times, being the midpoints of break_times, with a min
+    and max band added.
+    """
+    assert len(break_times)>2
+    
+    bands = [] #deepcopy(break_times)
+    bands.append(break_times[0]-0.5*(break_times[1]-break_times[0]))
+
+    for i,break_x in enumerate(break_times[0:-1]):
+        bands.append(break_times[i]+0.5*(break_times[i+1]-break_times[i]))
+        
+    bands.append(break_times[-1]+0.5*(break_times[-1]-break_times[-2]))
+    bands = ensure_numeric(bands)
+    bands += band_offset
+    return bands
+    
     
     
@@ -154 +399 @@
     """
     from scenarios import scenarios
-    scenarios = [scenarios[0]]
-    #gauges_for_slope("_good_tri_area_0.01_A", scenarios)
-    #auto_graph_slopes() 
-
+    #scenarios = [scenarios[0]]
+    outputdir_tag = "_good_tri_area_0.01_limiterC"
+    #outputdir_tag = "_test_limiterC"
+    #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!!
+    #gauges_for_slope(outputdir_tag, scenarios)
+    auto_graph_slopes(outputdir_tag, scenarios) #, is_interactive=True) 
+    #auto_find_min_slopes(outputdir_tag, scenarios)
+    auto_graph_froudes(outputdir_tag, scenarios)