Changeset 5494


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

Current Hinwood scenario - added plotting of froude number

Location:
anuga_work/development/Hinwood_2008
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/development/Hinwood_2008/plot.py

    r5461 r5494  
    1919    """
    2020    from pylab import ion, plot, xlabel, ylabel, close, legend, \
    21          savefig, title, axis
     21         savefig, title, axis, setp
    2222    from anuga.shallow_water.data_manager import csv2dict
    2323
    … …  
    3636        ion()
    3737   
    38     plot(time_sim, quantity_sim, time_exp, quantity_exp)
     38    l_sim, l_exp = plot(time_sim, quantity_sim, time_exp, quantity_exp)
     39    setp(l_sim, color='r')
     40    setp(l_exp, color='b')
    3941
    4042    # Add axis stuff and legend
    … …  
    7072                                   outputdir_name=outputdir_name)
    7173   
     74   
    7275    file_sim = pro_instance.outputdir + quantity + "_" + id + ".csv"
    7376    #print "file_exp",file_exp
    … …  
    8790    return file_exp, file_sim, location_sims, location_exps, outputdir_name, \
    8891           save_as_list
    89 def plot(scenarios, outputdir_tag, quantity = "stage"):
     92def plot(scenarios, outputdir_tag, quantity = "stage",is_interactive=False):
    9093    plot_type = ".pdf"
    9194   
    … …  
    120123                             plot_title=plot_title,
    121124                             y_label='Water '+ quantity +' (m)',
    122                              is_interactive=False,
     125                             is_interactive=is_interactive,
    123126                             save_as=save_as,
    124127                             use_axis=use_axis)
     128            if is_interactive is True:
     129                break
    125130   
    126131#-------------------------------------------------------------
    … …  
    129134    """
    130135    from scenarios import scenarios
    131     outputdir_tag = "_good_tri_area_0.01_C"
    132     outputdir_tag = "_test_C"
     136    outputdir_tag = "_good_tri_area_0.01_D"
     137    #outputdir_tag = "_test_C"
    133138    #scenarios = scenarios[1:]
    134     #scenarios = [scenarios[0]]
    135     plot(scenarios, outputdir_tag)
     139    scenarios = [scenarios[0]]
     140    is_interactive = False
     141    #is_interactive = True
     142    plot(scenarios, outputdir_tag,is_interactive=is_interactive)
    136143

  • anuga_work/development/Hinwood_2008/run_dam.py

    r5461 r5494  
    5959         metadata_dic,
    6060         boundary_path=None,
    61          friction=0.01,
     61         friction=0.012,  # planed wood. http://www.lmnoeng.com/manningn.htm
    6262         outputdir_name=None,
    6363         run_type=0):
    … …  
    6565   
    6666    basename = 'zz_' + metadata_dic['scenario_id']
     67    end_tag = '_limiterD'
    6768    if run_type == 1:
    68         outputdir_name += '_test_D'
     69        outputdir_name += '_test' + end_tag
    6970        yieldstep = 1.0
    7071        finaltime = 15.
    … …  
    7273       
    7374    elif run_type == 2:
    74         outputdir_name += '_test_long_time'
     75        outputdir_name += '_test_long_time' + end_tag
    7576        yieldstep = 0.5
    7677        finaltime = None
    … …  
    7879       
    7980    elif run_type == 3:
    80         outputdir_name += '_yieldstep_0.1_tri_area_0.01_D'
     81        outputdir_name += '_yieldstep_0.1_tri_area_0.01' + end_tag
    8182        yieldstep = 0.1
    8283        finaltime = None       
    8384        maximum_triangle_area=0.01
    8485    elif run_type == 4:
    85         outputdir_name += '_good_tri_area_0.01_D'
     86        outputdir_name += '_good_tri_area_0.01' + end_tag
    8687        # this is not a test
    8788        # Output will go to a file
    … …  
    9192        maximum_triangle_area=0.01
    9293    elif run_type == 5:
    93         outputdir_name += '_good_tri_area_0.001_D'
     94        outputdir_name += '_good_tri_area_0.001' + end_tag
    9495        # this is not a test
    9596        # Output will go to a file
    … …  
    117118    #return pro_instance
    118119    if finaltime is None:
    119         finaltime = boundary_final_time
     120        finaltime = boundary_final_time - 0.1 # Edge boundary problems
    120121    # Boundary file manipulation
    121122    if boundary_path is None:
    … …  
    168169    domain.set_datadir(pro_instance.outputdir)
    169170    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
    170     domain.set_minimum_storable_height(0.001)
    171     #domain.set_store_vertices_uniquely(True)  # for writting to sww
     171    domain.set_minimum_storable_height(0.0001)
     172
     173    domain.set_default_order(2) # Use second order spatial scheme
     174    domain.set_timestepping_method('rk2')
     175    domain.use_edge_limiter = True
     176    domain.tight_slope_limiters = True
     177   
     178    domain.beta_w      = 0.6
     179    domain.beta_uh     = 0.6
     180    domain.beta_vh     = 0.6
     181   
    172182
    173183    #-------------------------------------------------------------------------
    … …  
    206216
    207217    # It seems that ANUGA can't handle a starttime that is >0.
    208     domain.starttime = 1.0
     218    #domain.starttime = 1.0 #!!! what was this doing?
    209219    for t in domain.evolve(yieldstep, finaltime):   
    210220        domain.write_time()
    … …  
    230240                            pro_instance.outputdir + "velocity_x_" + id,
    231241                            pro_instance.outputdir + "velocity_y_" + id,
    232                             pro_instance.outputdir + "stage_" + id)
     242                            pro_instance.outputdir + "stage_" + id,
     243                            pro_instance.outputdir + "froude_" + id)
    233244 
    234245    return pro_instance
    … …  
    259270                             run_type = run_type,
    260271                             outputdir_name=run_data['scenario_id'])
    261         gauges_for_slope(pro_instance.outputdir,[run_data])
     272        #gauges_for_slope(pro_instance.outputdir,[run_data])

  • 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)

  • anuga_work/development/Hinwood_2008/test_slope.py

    r5447 r5494  
    1010import tempfile
    1111
    12 from scipy import array, allclose
     12from Numeric import array, allclose
     13from anuga.utilities.numerical_tools import ensure_numeric
     14
    1315
    1416from slope import *
    … …  
    3335        handle.write(sample)
    3436        handle.close()
    35         dtimes, depths, sensors = load_depths(file_name)
     37        dtimes, depths, sensors = load_sensors(file_name)
    3638        actual = array([0,0.01,0.02,0.03])
    3739        self.assert_ (allclose(dtimes, actual, 0.001))
    … …  
    7375        handle.write(sample)
    7476        handle.close()
    75         graph_slopes(file_name)
    7677
     78        # since it currently graphs
     79        #graph_slopes(file_name)
    7780       
    7881        os.remove(file_name)
    7982             
     83    def test_get_min_in_band(self):
     84        time = array([0,1,2,3,4,5])
     85        slopes = array([[10,12],
     86                        [9,10],
     87                        [8,9],
     88                        [9,7],
     89                        [10,6],
     90                        [7,2]])
     91
     92        max_q, max_q_times = get_min_in_band(0,5,time,slopes)
     93        actual = array([8,6])
     94        self.assert_ (allclose(max_q, actual, 0.001))
     95        actual = array([2,4])
     96        self.assert_ (allclose(max_q_times, actual, 0.001))
     97
     98       
     99        max_q, max_q_times = get_min_in_band(1,4,time,slopes)
     100        actual = array([8,7])
     101        self.assert_ (allclose(max_q, actual, 0.001))
     102        actual = array([2,3])
     103        self.assert_ (allclose(max_q_times, actual, 0.001))
     104       
     105    def test_break_times2bands(self):
     106        break_times = [1,3,4,8]
     107        band_offset = 0.0
     108        bands  = break_times2bands(break_times, band_offset)
     109        #print "bands", bands
     110        actual = array([0,2,3.5,6,10 ])
     111        self.assert_ (allclose(ensure_numeric(bands), actual, 0.001))
     112       
     113    def test_get_time_band(self):
     114        time = array([0,1,2,3,4,5])
     115        slopes = array([[10,12],
     116                        [9,10],
     117                        [8,9],
     118                        [9,7],
     119                        [10,6],
     120                        [7,2]])
     121
     122        band_time, band_quantity = get_time_band(0,5,time,slopes)
     123        actual = array([1,2,3,4])
     124        self.assert_ (allclose(band_time, actual, 0.001))
     125        actual = array([[9,10],
     126                        [8,9],
     127                        [9,7],
     128                        [10,6]])
     129        self.assert_ (allclose(band_quantity, actual, 0.001))
     130
     131   
     132    def test_get_band(self):
     133        time = array([0,1,2,3,4,5])
     134        slopes = array([[10,12],
     135                        [9,10],
     136                        [8,9],
     137                        [9,7],
     138                        [10,6],
     139                        [7,2]])
     140
     141        band_time, band_quantity = get_band(0,5,time,slopes,0)
     142        actual = array([1,2,3,4])
     143        self.assert_ (allclose(band_time, actual, 0.001))
     144        actual = array([[9,10],
     145                        [8,9],
     146                        [9,7],
     147                        [10,6]])
     148        self.assert_ (allclose(band_quantity, actual, 0.001))
     149       
     150        location = array([0,1,2,3,4,5])
     151        slopes = array([[0,1,2,3,4,5],
     152                        [0,10,20,30,40,50]])
     153
     154        band_time, band_quantity = get_band(0,5,location,slopes,-1)
     155        actual = array([1,2,3,4])
     156        self.assert_ (allclose(band_time, actual, 0.001))
     157        actual = array([[1,2,3,4],
     158                        [10,20,30,40]])
     159        self.assert_ (allclose(band_quantity, actual, 0.001))     
     160       
    80161if __name__ == "__main__":
    81162
    82163    suite = unittest.makeSuite(slopeTestCase,'test')
     164    suite = unittest.makeSuite(slopeTestCase,'test_get_band')
     165
    83166    runner = unittest.TextTestRunner()  #verbosity=2)
    84167    runner.run(suite)
Note: See TracChangeset for help on using the changeset viewer.