""" Plot up files from the Hinwood project. """ from os import sep import project from copy import deepcopy #from scipy import arange 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 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(quantity_file) lines = dfid.readlines() dfid.close() title = lines.pop(0) n_time = len(lines) n_sensors = len(lines[0].split(','))-1 # -1 to remove time dtimes = zeros(n_time, Float) #Time depths = zeros(n_time, Float) # sensors = zeros((n_time,n_sensors), Float) 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): fields = line.split(',') #(',') fields = [float(j) for j in fields] dtimes[i] = fields[0] sensors[i] = fields[1:] # 1: to remove time #print "dtimes",dtimes #print "locations", locations #print "sensors", sensors return dtimes, locations, sensors 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) # slopes = zeros((n_time,n_slope_locations), Float) # An array of the sensor spacing values delta_locations = zeros(n_slope_locations, Float) for i in arange(n_slope_locations): 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_locations[i] return times, slope_locations, slopes 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' if is_interactive: ion() # Can't seem to reshape this info once it is in the function CS = contourf(x_data, times, z_data, 10, 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, 'ro') # 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 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'] 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']) plot_title = "Stage slope " + id 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, plot_title=plot_title, 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']) plot_title = "Froude Number" + id 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, plot_title=plot_title, 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): """ This is used to create a stage file, using gauges relivent to finding a slope. """ dx = 0.05 for run_data in scenarios: point_x = arange(run_data['start_slope_x'], run_data['finish_slope_x'], dx).tolist() flume_y_middle = 0.5 points = [] for gauge_x in point_x: points.append([gauge_x, flume_y_middle]) id = run_data['scenario_id'] basename = 'zz_' + run_data['scenario_id'] outputdir_name = id + outputdir_tag pro_instance = project.Project(['data','flumes','Hinwood_2008'], 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, 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 #------------------------------------------------------------- if __name__ == "__main__": """ """ from scenarios import scenarios #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)