source: development/stochastic_study/plot_spread.py @ 3011

"""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 utilities.numerical_tools import histogram, create_bins
from read_realisations import read_realisations
import project

# Initialise
# Read in all realisations-timeseries
study = 'cyclone1'
study = 'cyclone3'
#study = 'nautilus3'

time, data, filenames = read_realisations(study, 
                                          #max_realisations = 200,
                                          gauge_number=2,
                                          use_cache=False)
number_of_realisations = data.shape[1]
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)


for j in range(20):
    #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))
    
    #title('Study %s: timeseries for %s (realisation %d of %d)'\
    #      %(study, filenames[j], 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('Study %s: spread at timestep %d of %d (t=%.2f)'\
          %(study, 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('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()