source: anuga_work/production/perth/ph2_compare.py @ 6032

Last change on this file since 6032 was 5992, checked in by Leharne, 15 years ago

Phase 2 hazard map comparison plot

File size: 4.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
9from 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)
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
27def 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###################################################
59# Determine max stage at boundary points
60###################################################
61
62event_number = 27283
63boundary_mean = get_max_boundary_data(project.boundaries_dir+str(event_number)+sep+'max_sts_stage.csv')
64
65###################################################
66# Read in max data from all models
67###################################################
68
69from Numeric import zeros, Float
70directory = project.home+project.state+sep+project.scenario+sep+'map_work'
71depth = [5.0,10.,20.,50.]
72no_models = 3
73filename = directory + sep + 'perth_ph2_compare_v5.csv'
74mean_stages, max_stages = get_stage_data(filename, depth, no_models)
75
76###################################################
77# Compare with Green's function and plot
78###################################################
79
80from anuga.abstract_2d_finite_volumes.util import greens_law
81from Numeric import arange
82d1 = 100.
83d2 = arange(d1,1,-0.1)
84h1 = boundary_mean
85green = []
86for d in d2:
87    h2 = greens_law(d1,d,h1)
88    green.append(h2)
89
90ion()
91figure(1)
92plot(depth,mean_stages[:,2],'>g',d2,green,'-g')
93xlabel('depth (m)') 
94ylabel('stage (m)')
95title('ANUGA outputs (average stage) versus Green\'s approximation \n \
96for event 27283 at Perth')
97legend(['original','Green\'s law'])
98#axis([5,105,min(min(stages))*0.9,max(max(stages))*1.1])
99grid(True)
100figname = 'ph2compare_perth_mean_ORIG_' + str(event_number) + '_mean'
101savefig(figname)
102
103figure(2)
104plot(depth,mean_stages[:,1],'+r',depth,mean_stages[:,2],'>g',d2,green,'-g')
105xlabel('depth (m)') 
106ylabel('stage (m)')
107title('ANUGA outputs (average stage) versus Green\'s approximation \n \
108for event 27283 at Perth')
109legend(['250m poly','original','Green\'s law'])
110#axis([5,105,min(min(stages))*0.9,max(max(stages))*1.1])
111grid(True)
112figname = 'ph2compare_perth_mean_250AP_' + str(event_number) + '_mean'
113savefig(figname)
114
115figure(3)
116plot(depth,mean_stages[:,0],'ob',depth,mean_stages[:,1],'+r',depth,mean_stages[:,2],'>g',d2,green,'-g')
117xlabel('depth (m)') 
118ylabel('stage (m)')
119title('ANUGA outputs (average stage) versus Green\'s approximation \n \
120for event 27283 at Perth')
121legend(['250m no poly','250m poly','original','Green\'s law'])
122#axis([5,105,min(min(stages))*0.9,max(max(stages))*1.1])
123grid(True)
124figname = 'ph2compare_perth_mean_ALL_' + str(event_number) + '_mean'
125savefig(figname)
126
127figure(4)
128plot(depth,max_stages[:,0],'ob',depth,max_stages[:,1],'+r',depth,max_stages[:,2],'>g',d2,green,'-g')
129xlabel('depth (m)') 
130ylabel('stage (m)')
131title('ANUGA outputs (max stage) versus Green\'s approximation \n \
132for event 27283 at Perth')
133legend(['250m no poly','250m poly','original','Green\'s law'])
134#axis([5,105,min(min(stages))*0.9,max(max(stages))*1.1])
135grid(True)
136figname = 'ph2compare_perth_max_ALL_' + str(event_number) + '_mean'
137savefig(figname)
138close('all')
139
Note: See TracBrowser for help on using the repository browser.