Ignore:
Timestamp:
Jul 11, 2008, 4:15:59 PM (14 years ago)
Author:
duncan
Message:

Current Hinwood scenario - added plotting of froude number

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/Hinwood_2008/slope.py

    r5455 r5494  
    55from os import sep
    66import project
     7from copy import deepcopy
    78#from scipy import arange
    8 
    9 from Numeric import arange, array, zeros, Float
     9from csv import writer
     10
     11from Numeric import arange, array, zeros, Float, where, greater, less, \
     12     compress, argmin, choose
    1013
    1114from anuga.fit_interpolate.interpolate import interpolate_sww2csv
    1215from anuga.shallow_water.data_manager import csv2dict
    13 
    14 def load_depths(slope_file):
     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):
    1526    #slope, _ = csv2dict(file_sim)
    1627   
    1728    # Read the depth file
    18     dfid = open(slope_file)
     29    dfid = open(quantity_file)
    1930    lines = dfid.readlines()
    2031    dfid.close()
     
    2637    depths = zeros(n_time, Float)  #
    2738    sensors = zeros((n_time,n_sensors), Float)
    28     depth_locations = title.split(',') #(',')
    29     depth_locations.pop(0) # remove 'time'
    30     # !!!! -3 assumes a
    31     depths = [float(j.split(':')[0]) for j in depth_locations]
     39    quantity_locations = title.split(',') #(',')
     40    quantity_locations.pop(0) # remove 'time'
     41   
     42    locations = [float(j.split(':')[0]) for j in quantity_locations]
    3243   
    3344    for i, line in enumerate(lines):
     
    3849
    3950    #print "dtimes",dtimes
    40     #print "depths", depths
     51    #print "locations", locations
    4152    #print "sensors", sensors
    42     return dtimes, depths, sensors
     53    return dtimes, locations, sensors
    4354   
    44 def load_slopes(slope_file):
    45     times, depths, sensors = load_depths(slope_file)
    46     n_slope_locations = len(depths)-1
     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
    4763    n_time = len(times)
    4864    slope_locations = zeros(n_slope_locations, Float)  #
     
    5066
    5167    # An array of the sensor spacing values
    52     delta_depths = zeros(n_slope_locations, Float)
     68    delta_locations = zeros(n_slope_locations, Float)
    5369   
    5470    for i in arange(n_slope_locations):
    55         slope_locations[i] = (depths[i+1] + depths[i+1])/2.
    56         delta_depths[i] = (depths[i+1] - depths[i])
     71        slope_locations[i] = (locations[i+1] + locations[i+1])/2.
     72        delta_locations[i] = (locations[i+1] - locations[i])
    5773   
    5874    for j in arange(n_time):
    5975        for i in arange(n_slope_locations):
    60             slopes[j,i] = (sensors[j,i+1] - sensors[j,i])/delta_depths[i]
     76            slopes[j,i] = (sensors[j,i+1] - sensors[j,i])/delta_locations[i]
    6177
    6278    return times, slope_locations, slopes
    6379
    64 def graph_slopes(slope_file, break_xs=None, break_times=None):
     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
    6590    from pylab import meshgrid, cm, contourf, contour, ion, plot, xlabel, \
    6691         ylabel, close, legend, savefig, title, figure ,colorbar, show , axis
     92
    6793    origin = 'lower'
    68     times, slope_locations, slopes = load_slopes(slope_file)
    69     #print "times", times
    70     #print "slope_locations", slope_locations
    71     #print "slopes", slopes
    72 
     94
     95    if is_interactive:
     96        ion()
     97       
    7398    # Can't seem to reshape this info once it is in the function
    74     CS = contourf(slope_locations, times, slopes, 10, # [-1, -0.1, 0, 0.1],
    75                   #alpha=0.5,
     99    CS = contourf(x_data, times, z_data, 10,
    76100                  cmap=cm.bone,
    77101                  origin=origin)
    78102   
    79     #CS2 = contour(slope_locations, times, slopes, CS.levels[::2],
    80        #           colors = 'r',
    81         #          origin=origin,
    82          #         hold='on')
    83 
    84     title('slope')
    85     xlabel('x location')
    86     ylabel('Time, seconds')
    87     #axis([5.0, 5.5, 30, 60])
    88 
     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)
    89111
    90112    if break_times is not None and break_xs is not None:
    91113        plot(break_xs, break_times, 'ro')
    92114   
     115   
    93116    # Make a colorbar for the ContourSet returned by the contourf call.
    94117    cbar = colorbar(CS)
    95     cbar.ax.set_ylabel('slope')
     118
    96119    # Add the contour line levels to the colorbar
     120    cbar.ax.set_ylabel('verbosity coefficient')
    97121    #cbar.add_lines(CS2)
    98 
    99     figure()
    100     show()
    101 
    102 def auto_graph_slopes():
    103     from scenarios import scenarios
    104    
    105     outputdir_tag = "_good_tri_area_0.01_A"
    106     outputdir_tag = "_good_tri_area_0.01_old"
    107     #outputdir_tag = "_test"
    108     scenarios = [scenarios[1]] # !!!!!!!!!!!!!!!!!!!!!!
     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"
    109182    for run_data in scenarios:
    110183        id = run_data['scenario_id']
     
    113186                                       outputdir_name=outputdir_name)
    114187        end = id + ".csv"
    115         slope_file = pro_instance.outputdir + "bslope_stage_" + end
    116         graph_slopes(slope_file, run_data['break_xs'],
    117                      run_data['break_times'])
     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       
    118280   
    119281def gauges_for_slope(outputdir_tag, scenarios):
    120    
    121     outputdir_tag = "_good_tri_area_0.01_A"
    122     #outputdir_tag = "_test"
    123     #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!!
     282    """
     283    This is used to create a stage file, using gauges relivent to
     284    finding a slope.
     285    """
    124286    dx = 0.05
    125287    for run_data in scenarios:
     
    138300                                       outputdir_name=outputdir_name)
    139301        end = id + ".csv"
    140         interpolate_sww2csv(pro_instance.outputdir + basename +".sww",
    141                             points,
    142                             pro_instance.outputdir + "fslope_depth_" + end,
    143                             pro_instance.outputdir + "fslope_velocity_x_" + end,
    144                             pro_instance.outputdir + "fslope_velocity_y_" + end,
    145                             pro_instance.outputdir + "fslope_stage_" + end,
    146                             time_thinning=1)
    147  
     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   
    148393   
    149394   
     
    154399    """
    155400    from scenarios import scenarios
    156     scenarios = [scenarios[0]]
    157     #gauges_for_slope("_good_tri_area_0.01_A", scenarios)
    158     #auto_graph_slopes()
    159 
     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 TracChangeset for help on using the changeset viewer.