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

"""
Read in event time series and determine max stage for Ph2 comparison
Compare with Green's function (no focussing)
Compare with ANUGA outputs (full model, 250m no polys, 250m all polys)
Leharne Fountain and Jane Sexton, 2008
"""

import project
from pylab import plot, xlabel, ylabel, savefig, ion, close, axis, title, legend, grid, figure
from os import sep

###################################################
# Relevant definitions
###################################################

def get_max_boundary_data(filename):
    from anuga.utilities.numerical_tools import mean
    fid = open(filename)
    lines = fid.readlines()
    fid.close()
    stage = []
    for line in lines[1:]:
        fields = line.split(',')
        stage.append(float(fields[3]))
    return mean(stage)

def get_stage_data(filename,depth,no_models):
    from anuga.utilities.numerical_tools import mean
    fid = open(filename)
    lines = fid.readlines()
    fid.close()
    
    mean_stages = zeros((len(depth),no_models), Float)
    max_stages = zeros((len(depth),no_models), Float)
    for i in range(no_models):
        stage5 = []
        stage10 = []
        stage20 = []
        stage50 = []
        for line in lines[1:]:
            fields = line.split(',')
            x = float(fields[0])
            # based on csv file output from ArcGIS: depth, x,y, no poly, all poly, and orig
            y = float(fields[i+3]) 
            if x == -depth[0]:
                stage5.append(y)
            if x == -depth[1]:
                stage10.append(y)
            if x == -depth[2]:
                stage20.append(y)
            if x == -depth[3]:
                stage50.append(y)
            
        mean_stages[:,i] = [mean(stage5), mean(stage10), mean(stage20), mean(stage50)]
        max_stages[:,i] = [max(stage5), max(stage10), max(stage20), max(stage50)]
    return mean_stages, max_stages

###################################################
# Determine max stage at boundary points
###################################################

event_number = 27283
boundary_mean = get_max_boundary_data(project.boundaries_dir+str(event_number)+sep+'max_sts_stage.csv')

###################################################
# Read in max data from all models
###################################################

from Numeric import zeros, Float
directory = project.home+project.state+sep+project.scenario+sep+'map_work'
depth = [5.0,10.,20.,50.]
no_models = 3
filename = directory + sep + 'perth_ph2_compare_v5.csv'
mean_stages, max_stages = get_stage_data(filename, depth, no_models)

###################################################
# Compare with Green's function and plot
###################################################

from anuga.abstract_2d_finite_volumes.util import greens_law
from Numeric import arange
d1 = 100.
d2 = arange(d1,1,-0.1)
h1 = boundary_mean
green = []
for d in d2:
    h2 = greens_law(d1,d,h1)
    green.append(h2)

ion()
figure(1)
plot(depth,mean_stages[:,2],'>g',d2,green,'-g')
xlabel('depth (m)') 
ylabel('stage (m)')
title('ANUGA outputs (average stage) versus Green\'s approximation \n \
for event 27283 at Perth')
legend(['original','Green\'s law'])
#axis([5,105,min(min(stages))*0.9,max(max(stages))*1.1])
grid(True)
figname = 'ph2compare_perth_mean_ORIG_' + str(event_number) + '_mean'
savefig(figname)

figure(2)
plot(depth,mean_stages[:,1],'+r',depth,mean_stages[:,2],'>g',d2,green,'-g')
xlabel('depth (m)') 
ylabel('stage (m)')
title('ANUGA outputs (average stage) versus Green\'s approximation \n \
for event 27283 at Perth')
legend(['250m poly','original','Green\'s law'])
#axis([5,105,min(min(stages))*0.9,max(max(stages))*1.1])
grid(True)
figname = 'ph2compare_perth_mean_250AP_' + str(event_number) + '_mean'
savefig(figname)

figure(3)
plot(depth,mean_stages[:,0],'ob',depth,mean_stages[:,1],'+r',depth,mean_stages[:,2],'>g',d2,green,'-g')
xlabel('depth (m)') 
ylabel('stage (m)')
title('ANUGA outputs (average stage) versus Green\'s approximation \n \
for event 27283 at Perth')
legend(['250m no poly','250m poly','original','Green\'s law'])
#axis([5,105,min(min(stages))*0.9,max(max(stages))*1.1])
grid(True)
figname = 'ph2compare_perth_mean_ALL_' + str(event_number) + '_mean'
savefig(figname)

figure(4)
plot(depth,max_stages[:,0],'ob',depth,max_stages[:,1],'+r',depth,max_stages[:,2],'>g',d2,green,'-g')
xlabel('depth (m)') 
ylabel('stage (m)')
title('ANUGA outputs (max stage) versus Green\'s approximation \n \
for event 27283 at Perth')
legend(['250m no poly','250m poly','original','Green\'s law'])
#axis([5,105,min(min(stages))*0.9,max(max(stages))*1.1])
grid(True)
figname = 'ph2compare_perth_max_ALL_' + str(event_number) + '_mean'
savefig(figname)
close('all')