# source:anuga_work/production/perth/ph2_compare_table.py@6540

Last change on this file since 6540 was 6093, checked in by jgriffin, 15 years ago
File size: 7.7 KB
Line
1"""
2Read in event time series and determine max stage for Ph2 comparison
3Compare with Green's function (no focussing)
4Compare with ANUGA outputs (full model, 250m no polys, 250m all polys)
5Leharne Fountain and Jane Sexton, 2008
6"""
7
8import project
9#from pylab import plot, xlabel, ylabel, savefig, ion, close, axis, title, legend, grid, figure
10from os import sep
11
12###################################################
13# Relevant definitions
14###################################################
15
16def get_max_boundary_data(filename):
17    from anuga.utilities.numerical_tools import mean
18    fid = open(filename)
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
27def get_stage_data(filename,depth,no_models):
28    from anuga.utilities.numerical_tools import mean
29    fid = open(filename)
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
58def 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
63def 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
75event_number = 27283
76boundary_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
82from Numeric import zeros, Float
83directory = project.home+project.state+sep+project.scenario+sep+'map_work'
84depth = [5.0,10.,20.,50.]
85no_models = 2
86#plot?
87plot_flag = False
88#Table??
89table_flag = True
90filename = directory + sep + 'perth_ph2_compare_'+str(event_number)+'.csv'
91
92mean_stages, max_stages = get_stage_data(filename, depth, no_models)
93print mean_stages
94###################################################
95# Compare with Green's function and plot
96###################################################
97
98from anuga.abstract_2d_finite_volumes.util import greens_law
99from Numeric import arange
100d1 = 100.
101d2 = arange(d1,1,-0.1)
102h1 = boundary_mean
103green = []
104for 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###################################################
111M_value = calc_M(mean_stages,boundary_mean)
112
113
114###################################################
115#calculate absolute and relative errors
116###################################################
117abs_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 perth')
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_perth_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 perth')
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_perth_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 perth')
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_perth_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 perth')
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_perth_max_ALL_' + str(event_number) + '_mean'
168##    savefig(figname)
169##    close('all')
170
171if table_flag==True:
172    import csv
173    csvfile = directory+sep+'ph2compare_perth_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
Note: See TracBrowser for help on using the repository browser.