source: development/stochastic_study/plot_spread.py @ 2986

"""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
from caching import cache

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

# Initialise
gauge = '%s.txt' %project.gauge_names[0]
number_of_realisations = project.number_of_realisations

time = zeros(project.number_of_timesteps, Float)
data = zeros((project.number_of_timesteps, number_of_realisations), Float)


# Read in all realisations-timeseries

j = 0 # Count realisations
for filename in os.listdir(project.working_dir):
    if filename.startswith(project.basename) and filename.endswith(gauge):
        if j < data.shape[1]:        
            print 'Reading filename %s (column %d)' %(filename, j)
            fid = open(filename)
            for i, line in enumerate(fid.readlines()):
                if i < data.shape[0]:
                    fields = line.strip().split()
                    time[i] = float(fields[0])
                    data[i,j] = float(fields[1])
        else:
            print 'Ignored %s (column %d)' %(filename, j)

        fid.close()
        j += 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(number_of_realisations):
    print j, data[:,j]
    plot(time, data[:,j], 'k-')
    
    xlabel('time(s)')
    ylabel('stage (m)')
    title('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(200,project.number_of_timesteps):
    # Plot histogram

    w = data[i,:]
    plot(w, 'k.')
    xlabel('realisations')
    ylabel('stage (m)')
    title('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(200,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('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()