source: anuga_work/development/stochastic_study/plot_spread.py @ 5599

"""Read in timeseries from N simulations and plot spread for each timestep at specified locations (ch 5,7,9)
"""

# Standard modules
import sys, os, string

# Related major packages
from Numeric import zeros, Float, allclose, arange

# Application specific imports
from anuga.utilities.numerical_tools import histogram, create_bins
from read_realisations import read_realisations
import project

# Initialise
# Read in all realisations-timeseries

#study = 'nautilus1' #~70 realisations, blocks of 100, sequential, stddev= 0.0006
#study = 'nautilus3' #~100 realisations, blocks of 10, sequential, stddev= 0.0013 
study = 'cyclone1' #~4000 realisations, blocks of 100, parallel, stddev= 0.0006
#study = 'cyclone3' #~120 realisations, blocks of 10, parallel, stddev= 0.0013

time, data, filenames = read_realisations(study, 
                                          #max_realisations = 200,
                                          gauge_number=0,
                                          sorting='randomised',
                                          #sorting='numerical',
                                          verbose=True,
                                          use_cache=True)
number_of_realisations = data.shape[1]
number_of_timesteps = data.shape[0]
print 'Read %d realisations' %number_of_realisations


# Plot
from pylab import ion, hold, plot, title, xlabel, ylabel, legend, savefig, show

ion()
hold(True)
#hold(False)

N=250
for j in range(N):
    #print j, filenames[j]
    plot(time, data[:,j], 'k-')
    
    xlabel('time(s)')
    ylabel('stage (m)')

    shortname = string.join(os.path.split(filenames[j])[-1].split('_')[-2:], '_')
    title('Study %s: timeseries for %s (realisation %d of %d)'\
          %(study, shortname, j, number_of_realisations-1))
    
    #raw_input('Next')

title('%d realisations' %N)
savefig('realisations', dpi = 300)
#raw_input('Stats')    
 


hold(False)
# Plot spread of stage values for each timestep
for i in range(0, project.number_of_timesteps, 10):
    # Plot histogram#

    w = data[i,:]
    plot(w, 'k.')
    xlabel('realisations')
    ylabel('stage (m)')
    #title('Study %s: spread at timestep %d of %d (t=%.2f)'\
    #      %(study, i, project.number_of_timesteps-1, time[i]))
    title('Spread at timestep %d of %d (t=%.2f)'\
          %(i, project.number_of_timesteps-1, time[i]))

    savefig('spread_%d' %i, dpi = 300)    
    #raw_input('Next')

show()
import sys; sys.exit()

hold(False)
# Plot histogram of stage values for each timestep
for i in range(300,320):#project.number_of_timesteps):
    # Plot histogram

    w = data[i,:]
     
    bins = create_bins(w, project.number_of_bins)
    print 'bins', bins
    hist = histogram(w, bins, relative=True)
    #print 'hist', hist
    
    #plot(w, 'k.')
    plot(bins, hist, 'k.')
    xlabel('stage (m)')
    ylabel('frequency')
    title('Study %s: histogram at timestep %d of %d (t=%.2f)'\
          %(study, i, project.number_of_timesteps-1, time[i]))        
    raw_input('Next')


#title('Gauge %s' %name)
#xlabel('time(s)')
#ylabel('stage (m)')    
#legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
#savefig(name, dpi = 300)






show()



#from pylab import *
#plot(time, stage)
#show()