source: anuga_work/production/busselton/ph2_compare_table.py @ 6874

"""
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

def calc_M(mean_stages,boundary_mean):
    #M_value=zeros(no_models,Float)
    M_value=mean_stages[1,:]/boundary_mean
    return M_value

def calc_errors(mean_stages,no_models):
    abs_error = zeros(len(mean_stages[1,:]),Float)
    rel_error = zeros(len(mean_stages[1,:]),Float)
    for i in range(len(mean_stages[1,:])):
        abs_error[i] = (mean_stages[1,i]-mean_stages[1,no_models-1])
        rel_error[i]= (abs_error[i]/mean_stages[1,no_models-1])*100
    return abs_error,rel_error

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

event_number = 68693
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 = 2
#plot?
plot_flag = False
#Table??
table_flag = True
filename = directory + sep + 'Busselton_Ph2_compare_'+str(event_number)+'.csv'

mean_stages, max_stages = get_stage_data(filename, depth, no_models)
print mean_stages
###################################################
# 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)

###################################################
#calculate M values (scaling of 100 m depth mean to 10 m)
###################################################
M_value = calc_M(mean_stages,boundary_mean)


###################################################
#calculate absolute and relative errors 
###################################################
abs_error,rel_error=calc_errors(mean_stages,no_models)


##if plot_flag==True:
##    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 Busselton')
##    legend(['original','Green\'s law'])
##    #axis([5,105,min(min(stages))*0.9,max(max(stages))*1.1])
##    grid(True)
##    figname = 'ph2compare_Busselton_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 Busselton')
##    legend(['250m poly','original','Green\'s law'])
##    #axis([5,105,min(min(stages))*0.9,max(max(stages))*1.1])
##    grid(True)
##    figname = 'ph2compare_Busselton_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 Busselton')
##    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_Busselton_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 Busselton')
##    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_Busselton_max_ALL_' + str(event_number) + '_mean'
##    savefig(figname)
##    close('all')

if table_flag==True:
    import csv
    csvfile = directory+sep+'ph2compare_Busselton_max_ALL_' + str(event_number) + '_mean.csv'
    fid1 = open(csvfile,'w')
    writer = csv.writer(fid1,delimiter=',',)

    if no_models==1:
        writer.writerow(('depth (m)','250m no poly'))
        for i in range(len(depth)):
            writer.writerow((depth[i],mean_stages[i,0]))
        writer.writerow(('100',boundary_mean,))
        writer.writerow(('M-value',M_value[0]))
        #writer.writerow(('Absolute Error', abs_error[0],abs_error[1]))
        #writer.writerow(('Relative Error %', rel_error[0],rel_error[1]))

    if no_models==2:
        writer.writerow(('depth (m)','250m no poly','original'))
        for i in range(len(depth)):
            writer.writerow((depth[i],mean_stages[i,0],mean_stages[i,1]))
        writer.writerow(('100',boundary_mean,boundary_mean))
        writer.writerow(('M-value',M_value[0],M_value[1]))
        writer.writerow(('Absolute Error', abs_error[0],abs_error[1]))
        writer.writerow(('Relative Error %', rel_error[0],rel_error[1]))
        
    if no_models==3:
        writer.writerow(('depth (m)','250m no poly','250m poly','original'))
        for i in range(len(depth)):
            writer.writerow((depth[i],mean_stages[i,0],mean_stages[i,1],mean_stages[i,2]))
        writer.writerow(('100',boundary_mean,boundary_mean,boundary_mean))
        writer.writerow(('M-value',M_value[0],M_value[1],M_value[2]))
        writer.writerow(('Absolute Error', abs_error[0],abs_error[1],abs_error[2]))
        writer.writerow(('Relative Error %', rel_error[0],rel_error[1],rel_error[2]))

    fid1.close()