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

Last change on this file since 5426 was 5426, checked in by duncan, 16 years ago

Current Hinwood scenario

File size: 4.7 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
10from pylab import meshgrid, cm, contourf, contour, ion, plot, xlabel, \
11     ylabel, close, legend, savefig, title, figure ,colorbar, show #, cbar
12
13from anuga.fit_interpolate.interpolate import interpolate_sww2csv
14from anuga.shallow_water.data_manager import csv2dict
15
16def load_depths(slope_file):
17    #slope, _ = csv2dict(file_sim)
18   
19    # Read the depth file
20    dfid = open(slope_file)
21    lines = dfid.readlines()
22    dfid.close()
23
24    title = lines.pop(0)
25    n_time = len(lines)
26    n_sensors = len(lines[0].split(','))-1  # -1 to remove time
27    dtimes = zeros(n_time, Float)  #Time
28    depths = zeros(n_time, Float)  #
29    sensors = zeros((n_time,n_sensors), Float)
30    depth_locations = title.split(',') #(',')
31    depth_locations.pop(0) # remove 'time'
32    # !!!! -3 assumes a
33    depths = [float(j.split(':')[0]) for j in depth_locations]
34   
35    for i, line in enumerate(lines):
36        fields = line.split(',') #(',')
37        fields = [float(j) for j in fields]
38        dtimes[i] = fields[0]
39        sensors[i] = fields[1:] # 1: to remove time
40
41    #print "dtimes",dtimes
42    #print "depths", depths
43    #print "sensors", sensors
44    return dtimes, depths, sensors
45   
46def load_slopes(slope_file):
47    times, depths, sensors = load_depths(slope_file)
48    n_slope_locations = len(depths)-1
49    n_time = len(times)
50    slope_locations = zeros(n_slope_locations, Float)  #
51    slopes = zeros((n_time,n_slope_locations), Float)
52
53    # An array of the sensor spacing values
54    delta_depths = zeros(n_slope_locations, Float)
55   
56    for i in arange(n_slope_locations):
57        slope_locations[i] = (depths[i+1] + depths[i+1])/2.
58        delta_depths[i] = (depths[i+1] - depths[i])
59   
60    for j in arange(n_time):
61        for i in arange(n_slope_locations):
62            slopes[j,i] = (sensors[j,i+1] - sensors[j,i])/delta_depths[i]
63
64    return times, slope_locations, slopes
65
66def graph_slopes(slope_file):
67    origin = 'lower'
68    times, slope_locations, slopes = load_slopes(slope_file)
69    X, Y = meshgrid(times, slope_locations)
70    #print "X", times
71    #print "Y", slope_locations
72    #print "slopes", slopes
73    CS = contourf(slope_locations, times, slopes, 10, # [-1, -0.1, 0, 0.1],
74                        #alpha=0.5,
75                        cmap=cm.bone,
76                        origin=origin)
77   
78    #CS2 = contour(slope_locations, times, slopes, CS.levels[::2],
79       #           colors = 'r',
80        #          origin=origin,
81         #         hold='on')
82
83    title('slope')
84    xlabel('x location')
85    ylabel('Time, seconds')
86
87    # Make a colorbar for the ContourSet returned by the contourf call.
88    cbar = colorbar(CS)
89    cbar.ax.set_ylabel('verbosity coefficient')
90    # Add the contour line levels to the colorbar
91    #cbar.add_lines(CS2)
92
93    figure()
94    show()
95
96def auto_graph_slopes():
97    from scenarios import scenarios
98   
99    outputdir_tag = "_good_tri_area_0.01_A"
100    scenarios = [scenarios[1]] # !!!!!!!!!!!!!!!!!!!!!!
101    for run_data in scenarios:
102        id = run_data['scenario_id']
103        outputdir_name = id + outputdir_tag
104        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
105                                       outputdir_name=outputdir_name)
106        end = id + ".csv"
107        slope_file = pro_instance.outputdir + "slope_depth_" + end
108        graph_slopes(slope_file) 
109   
110def gauges_for_slope():
111    from scenarios import scenarios
112   
113    outputdir_tag = "_good_tri_area_0.01_A"
114    scenarios = [scenarios[1]] # !!!!!!!!!!!!!!!!!!!!!!
115    dx = 0.001
116    for run_data in scenarios:
117        point_x = arange(run_data['start_slope_x'],
118                        run_data['finish_slope_x'],
119                        dx).tolist()
120        flume_y_middle = 0.5
121        points = []
122        for gauge_x in point_x:
123            points.append([gauge_x, flume_y_middle])
124        id = run_data['scenario_id']
125   
126        basename = 'zz_' + run_data['scenario_id']
127        outputdir_name = id + outputdir_tag
128        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
129                                       outputdir_name=outputdir_name)
130        end = id + ".csv"
131        interpolate_sww2csv(pro_instance.outputdir + basename +".sww",
132                            points,
133                            pro_instance.outputdir + "slope_depth_" + end,
134                            pro_instance.outputdir + "slope_velocity_x_" + end,
135                            pro_instance.outputdir + "slope_velocity_y_" + end)
136 
137   
138   
139
140#-------------------------------------------------------------
141if __name__ == "__main__":
142    """
143    """
144    auto_graph_slopes() 
145
Note: See TracBrowser for help on using the repository browser.