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

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

Current Hinwood scenario - cropping the flume length

File size: 5.3 KB
Line 
1
2"""
3Plot up files from the Hinwood project.
4"""
5from os import sep
6import project
7#from scipy import arange
8
9from Numeric import arange, array, zeros, Float
10
11from anuga.fit_interpolate.interpolate import interpolate_sww2csv
12from anuga.shallow_water.data_manager import csv2dict
13
14def load_depths(slope_file):
15    #slope, _ = csv2dict(file_sim)
16   
17    # Read the depth file
18    dfid = open(slope_file)
19    lines = dfid.readlines()
20    dfid.close()
21
22    title = lines.pop(0)
23    n_time = len(lines)
24    n_sensors = len(lines[0].split(','))-1  # -1 to remove time
25    dtimes = zeros(n_time, Float)  #Time
26    depths = zeros(n_time, Float)  #
27    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]
32   
33    for i, line in enumerate(lines):
34        fields = line.split(',') #(',')
35        fields = [float(j) for j in fields]
36        dtimes[i] = fields[0]
37        sensors[i] = fields[1:] # 1: to remove time
38
39    #print "dtimes",dtimes
40    #print "depths", depths
41    #print "sensors", sensors
42    return dtimes, depths, sensors
43   
44def load_slopes(slope_file):
45    times, depths, sensors = load_depths(slope_file)
46    n_slope_locations = len(depths)-1
47    n_time = len(times)
48    slope_locations = zeros(n_slope_locations, Float)  #
49    slopes = zeros((n_time,n_slope_locations), Float)
50
51    # An array of the sensor spacing values
52    delta_depths = zeros(n_slope_locations, Float)
53   
54    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])
57   
58    for j in arange(n_time):
59        for i in arange(n_slope_locations):
60            slopes[j,i] = (sensors[j,i+1] - sensors[j,i])/delta_depths[i]
61
62    return times, slope_locations, slopes
63
64def graph_slopes(slope_file, break_xs=None, break_times=None):
65    from pylab import meshgrid, cm, contourf, contour, ion, plot, xlabel, \
66         ylabel, close, legend, savefig, title, figure ,colorbar, show , axis
67    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
73    # 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,
76                  cmap=cm.bone,
77                  origin=origin)
78   
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
89
90    if break_times is not None and break_xs is not None:
91        plot(break_xs, break_times, 'ro')
92   
93    # Make a colorbar for the ContourSet returned by the contourf call.
94    cbar = colorbar(CS)
95    cbar.ax.set_ylabel('slope')
96    # Add the contour line levels to the colorbar
97    #cbar.add_lines(CS2)
98
99    figure()
100    show()
101
102def 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]] # !!!!!!!!!!!!!!!!!!!!!!
109    for run_data in scenarios:
110        id = run_data['scenario_id']
111        outputdir_name = id + outputdir_tag
112        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
113                                       outputdir_name=outputdir_name)
114        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']) 
118   
119def gauges_for_slope(outputdir_tag, scenarios):
120   
121    outputdir_tag = "_good_tri_area_0.01_A"
122    #outputdir_tag = "_test"
123    #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!!
124    dx = 0.05
125    for run_data in scenarios:
126        point_x = arange(run_data['start_slope_x'],
127                        run_data['finish_slope_x'],
128                        dx).tolist()
129        flume_y_middle = 0.5
130        points = []
131        for gauge_x in point_x:
132            points.append([gauge_x, flume_y_middle])
133        id = run_data['scenario_id']
134   
135        basename = 'zz_' + run_data['scenario_id']
136        outputdir_name = id + outputdir_tag
137        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
138                                       outputdir_name=outputdir_name)
139        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 
148   
149   
150
151#-------------------------------------------------------------
152if __name__ == "__main__":
153    """
154    """
155    from scenarios import scenarios
156    scenarios = [scenarios[0]]
157    #gauges_for_slope("_good_tri_area_0.01_A", scenarios)
158    #auto_graph_slopes()
159
Note: See TracBrowser for help on using the repository browser.