source: development/stochastic_study/plot_spread.py @ 3006

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

# Standard modules
import sys, os

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

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

# Initialise
# Read in all realisations-timeseries
time, data, filenames = read_realisations('cyclone1', 
                                          max_realisations = 200,
                                          use_cache=True)
#time, data, filenames = read_realisations('test', exclude='_0_',
#                                          max_realisations = 200,
#                                          use_cache=True)
#time, data, filenames = read_realisations('nautilus1', exclude='_0_')
number_of_realisations = data.shape[1]


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

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

# Simple plot of timeseries for different realisations 
for j in range(20):
    print j, filenames[j]
    plot(time, data[:,j], 'k-')
    
    xlabel('time(s)')
    ylabel('stage (m)')
    title('Timeseries for realisation %d of %d' %(j, number_of_realisations-1))
    raw_input('Next')
 


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

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


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('Histogram at timestep %d of %d (t=%.2f)'\
          %(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()