1 | """ |
---|
2 | Read in event time series and determine max stage for Ph2 comparison |
---|
3 | Compare with Green's function (no focussing) |
---|
4 | Compare with ANUGA outputs (full model, 250m no polys, 250m all polys) |
---|
5 | Leharne Fountain and Jane Sexton, 2008 |
---|
6 | """ |
---|
7 | |
---|
8 | import project |
---|
9 | #from pylab import plot, xlabel, ylabel, savefig, ion, close, axis, title, legend, grid, figure |
---|
10 | from os import sep |
---|
11 | |
---|
12 | ################################################### |
---|
13 | # Relevant definitions |
---|
14 | ################################################### |
---|
15 | |
---|
16 | def get_max_boundary_data(filename): |
---|
17 | from anuga.utilities.numerical_tools import mean |
---|
18 | fid = open(filename) |
---|
19 | lines = fid.readlines() |
---|
20 | fid.close() |
---|
21 | stage = [] |
---|
22 | for line in lines[1:]: |
---|
23 | fields = line.split(',') |
---|
24 | stage.append(float(fields[3])) |
---|
25 | return mean(stage) |
---|
26 | |
---|
27 | def get_stage_data(filename,depth,no_models): |
---|
28 | from anuga.utilities.numerical_tools import mean |
---|
29 | fid = open(filename) |
---|
30 | lines = fid.readlines() |
---|
31 | fid.close() |
---|
32 | |
---|
33 | mean_stages = zeros((len(depth),no_models), Float) |
---|
34 | max_stages = zeros((len(depth),no_models), Float) |
---|
35 | for i in range(no_models): |
---|
36 | stage5 = [] |
---|
37 | stage10 = [] |
---|
38 | stage20 = [] |
---|
39 | stage50 = [] |
---|
40 | for line in lines[1:]: |
---|
41 | fields = line.split(',') |
---|
42 | x = float(fields[0]) |
---|
43 | # based on csv file output from ArcGIS: depth, x,y, no poly, all poly, and orig |
---|
44 | y = float(fields[i+3]) |
---|
45 | if x == -depth[0]: |
---|
46 | stage5.append(y) |
---|
47 | if x == -depth[1]: |
---|
48 | stage10.append(y) |
---|
49 | if x == -depth[2]: |
---|
50 | stage20.append(y) |
---|
51 | if x == -depth[3]: |
---|
52 | stage50.append(y) |
---|
53 | |
---|
54 | mean_stages[:,i] = [mean(stage5), mean(stage10), mean(stage20), mean(stage50)] |
---|
55 | max_stages[:,i] = [max(stage5), max(stage10), max(stage20), max(stage50)] |
---|
56 | return mean_stages, max_stages |
---|
57 | |
---|
58 | def calc_M(mean_stages,boundary_mean): |
---|
59 | #M_value=zeros(no_models,Float) |
---|
60 | M_value=mean_stages[1,:]/boundary_mean |
---|
61 | return M_value |
---|
62 | |
---|
63 | def calc_errors(mean_stages,no_models): |
---|
64 | abs_error = zeros(len(mean_stages[1,:]),Float) |
---|
65 | rel_error = zeros(len(mean_stages[1,:]),Float) |
---|
66 | for i in range(len(mean_stages[1,:])): |
---|
67 | abs_error[i] = (mean_stages[1,i]-mean_stages[1,no_models-1]) |
---|
68 | rel_error[i]= (abs_error[i]/mean_stages[1,no_models-1])*100 |
---|
69 | return abs_error,rel_error |
---|
70 | |
---|
71 | ################################################### |
---|
72 | # Determine max stage at boundary points |
---|
73 | ################################################### |
---|
74 | |
---|
75 | event_number = 68693 |
---|
76 | boundary_mean = get_max_boundary_data(project.boundaries_dir+str(event_number)+sep+'max_sts_stage.csv') |
---|
77 | |
---|
78 | ################################################### |
---|
79 | # Read in max data from all models |
---|
80 | ################################################### |
---|
81 | |
---|
82 | from Numeric import zeros, Float |
---|
83 | directory = project.home+project.state+sep+project.scenario+sep+'map_work' |
---|
84 | depth = [5.0,10.,20.,50.] |
---|
85 | no_models = 2 |
---|
86 | #plot? |
---|
87 | plot_flag = False |
---|
88 | #Table?? |
---|
89 | table_flag = True |
---|
90 | filename = directory + sep + 'Busselton_Ph2_compare_'+str(event_number)+'.csv' |
---|
91 | |
---|
92 | mean_stages, max_stages = get_stage_data(filename, depth, no_models) |
---|
93 | print mean_stages |
---|
94 | ################################################### |
---|
95 | # Compare with Green's function and plot |
---|
96 | ################################################### |
---|
97 | |
---|
98 | from anuga.abstract_2d_finite_volumes.util import greens_law |
---|
99 | from Numeric import arange |
---|
100 | d1 = 100. |
---|
101 | d2 = arange(d1,1,-0.1) |
---|
102 | h1 = boundary_mean |
---|
103 | green = [] |
---|
104 | for d in d2: |
---|
105 | h2 = greens_law(d1,d,h1) |
---|
106 | green.append(h2) |
---|
107 | |
---|
108 | ################################################### |
---|
109 | #calculate M values (scaling of 100 m depth mean to 10 m) |
---|
110 | ################################################### |
---|
111 | M_value = calc_M(mean_stages,boundary_mean) |
---|
112 | |
---|
113 | |
---|
114 | ################################################### |
---|
115 | #calculate absolute and relative errors |
---|
116 | ################################################### |
---|
117 | abs_error,rel_error=calc_errors(mean_stages,no_models) |
---|
118 | |
---|
119 | |
---|
120 | ##if plot_flag==True: |
---|
121 | ## ion() |
---|
122 | ## figure(1) |
---|
123 | ## plot(depth,mean_stages[:,2],'>g',d2,green,'-g') |
---|
124 | ## xlabel('depth (m)') |
---|
125 | ## ylabel('stage (m)') |
---|
126 | ## title('ANUGA outputs (average stage) versus Green\'s approximation \n \ |
---|
127 | ## for event 27283 at Busselton') |
---|
128 | ## legend(['original','Green\'s law']) |
---|
129 | ## #axis([5,105,min(min(stages))*0.9,max(max(stages))*1.1]) |
---|
130 | ## grid(True) |
---|
131 | ## figname = 'ph2compare_Busselton_mean_ORIG_' + str(event_number) + '_mean' |
---|
132 | ## savefig(figname) |
---|
133 | ## |
---|
134 | ## figure(2) |
---|
135 | ## plot(depth,mean_stages[:,1],'+r',depth,mean_stages[:,2],'>g',d2,green,'-g') |
---|
136 | ## xlabel('depth (m)') |
---|
137 | ## ylabel('stage (m)') |
---|
138 | ## title('ANUGA outputs (average stage) versus Green\'s approximation \n \ |
---|
139 | ## for event 27283 at Busselton') |
---|
140 | ## legend(['250m poly','original','Green\'s law']) |
---|
141 | ## #axis([5,105,min(min(stages))*0.9,max(max(stages))*1.1]) |
---|
142 | ## grid(True) |
---|
143 | ## figname = 'ph2compare_Busselton_mean_250AP_' + str(event_number) + '_mean' |
---|
144 | ## savefig(figname) |
---|
145 | ## |
---|
146 | ## figure(3) |
---|
147 | ## plot(depth,mean_stages[:,0],'ob',depth,mean_stages[:,1],'+r',depth,mean_stages[:,2],'>g',d2,green,'-g') |
---|
148 | ## xlabel('depth (m)') |
---|
149 | ## ylabel('stage (m)') |
---|
150 | ## title('ANUGA outputs (average stage) versus Green\'s approximation \n \ |
---|
151 | ## for event 27283 at Busselton') |
---|
152 | ## legend(['250m no poly','250m poly','original','Green\'s law']) |
---|
153 | ## #axis([5,105,min(min(stages))*0.9,max(max(stages))*1.1]) |
---|
154 | ## grid(True) |
---|
155 | ## figname = 'ph2compare_Busselton_mean_ALL_' + str(event_number) + '_mean' |
---|
156 | ## savefig(figname) |
---|
157 | ## |
---|
158 | ## figure(4) |
---|
159 | ## plot(depth,max_stages[:,0],'ob',depth,max_stages[:,1],'+r',depth,max_stages[:,2],'>g',d2,green,'-g') |
---|
160 | ## xlabel('depth (m)') |
---|
161 | ## ylabel('stage (m)') |
---|
162 | ## title('ANUGA outputs (max stage) versus Green\'s approximation \n \ |
---|
163 | ## for event 27283 at Busselton') |
---|
164 | ## legend(['250m no poly','250m poly','original','Green\'s law']) |
---|
165 | ## #axis([5,105,min(min(stages))*0.9,max(max(stages))*1.1]) |
---|
166 | ## grid(True) |
---|
167 | ## figname = 'ph2compare_Busselton_max_ALL_' + str(event_number) + '_mean' |
---|
168 | ## savefig(figname) |
---|
169 | ## close('all') |
---|
170 | |
---|
171 | if table_flag==True: |
---|
172 | import csv |
---|
173 | csvfile = directory+sep+'ph2compare_Busselton_max_ALL_' + str(event_number) + '_mean.csv' |
---|
174 | fid1 = open(csvfile,'w') |
---|
175 | writer = csv.writer(fid1,delimiter=',',) |
---|
176 | |
---|
177 | if no_models==1: |
---|
178 | writer.writerow(('depth (m)','250m no poly')) |
---|
179 | for i in range(len(depth)): |
---|
180 | writer.writerow((depth[i],mean_stages[i,0])) |
---|
181 | writer.writerow(('100',boundary_mean,)) |
---|
182 | writer.writerow(('M-value',M_value[0])) |
---|
183 | #writer.writerow(('Absolute Error', abs_error[0],abs_error[1])) |
---|
184 | #writer.writerow(('Relative Error %', rel_error[0],rel_error[1])) |
---|
185 | |
---|
186 | if no_models==2: |
---|
187 | writer.writerow(('depth (m)','250m no poly','original')) |
---|
188 | for i in range(len(depth)): |
---|
189 | writer.writerow((depth[i],mean_stages[i,0],mean_stages[i,1])) |
---|
190 | writer.writerow(('100',boundary_mean,boundary_mean)) |
---|
191 | writer.writerow(('M-value',M_value[0],M_value[1])) |
---|
192 | writer.writerow(('Absolute Error', abs_error[0],abs_error[1])) |
---|
193 | writer.writerow(('Relative Error %', rel_error[0],rel_error[1])) |
---|
194 | |
---|
195 | if no_models==3: |
---|
196 | writer.writerow(('depth (m)','250m no poly','250m poly','original')) |
---|
197 | for i in range(len(depth)): |
---|
198 | writer.writerow((depth[i],mean_stages[i,0],mean_stages[i,1],mean_stages[i,2])) |
---|
199 | writer.writerow(('100',boundary_mean,boundary_mean,boundary_mean)) |
---|
200 | writer.writerow(('M-value',M_value[0],M_value[1],M_value[2])) |
---|
201 | writer.writerow(('Absolute Error', abs_error[0],abs_error[1],abs_error[2])) |
---|
202 | writer.writerow(('Relative Error %', rel_error[0],rel_error[1],rel_error[2])) |
---|
203 | |
---|
204 | fid1.close() |
---|
205 | |
---|
206 | |
---|