source: development/stochastic_study/plot_spread.py @ 3322

Last change on this file since 3322 was 3035, checked in by ole, 19 years ago

Working with Suresh at ACFR

File size: 2.9 KB
RevLine 
[2967]1"""Read in timeseries from N simulations and plot spread for each timestep at specified locations (ch 5,7,9)
2"""
3
[2972]4# Standard modules
[3010]5import sys, os, string
[2972]6
7# Related major packages
8from Numeric import zeros, Float, allclose, arange
9
10# Application specific imports
11from utilities.numerical_tools import histogram, create_bins
[3006]12from read_realisations import read_realisations
[2967]13import project
14
[2972]15# Initialise
16# Read in all realisations-timeseries
[3010]17
[3035]18#study = 'nautilus1' #~70 realisations, blocks of 100, sequential, stddev= 0.0006
19#study = 'nautilus3' #~100 realisations, blocks of 10, sequential, stddev= 0.0013
20study = 'cyclone1' #~4000 realisations, blocks of 100, parallel, stddev= 0.0006
21#study = 'cyclone3' #~120 realisations, blocks of 10, parallel, stddev= 0.0013
22
[3010]23time, data, filenames = read_realisations(study, 
[3011]24                                          #max_realisations = 200,
[3035]25                                          gauge_number=0,
26                                          sorting='randomised',
27                                          #sorting='numerical',
28                                          verbose=True,
29                                          use_cache=True)
[3006]30number_of_realisations = data.shape[1]
[3035]31number_of_timesteps = data.shape[0]
[3011]32print 'Read %d realisations' %number_of_realisations
[2972]33
[3035]34
[2967]35# Plot
36from pylab import ion, hold, plot, title, xlabel, ylabel, legend, savefig, show
37
38ion()
[3006]39hold(True)
40#hold(False)
[2967]41
[3010]42
[3035]43for j in range(1):
[3010]44    #print j, filenames[j]
[2975]45    plot(time, data[:,j], 'k-')
[2967]46   
[2972]47    xlabel('time(s)')
48    ylabel('stage (m)')
[3010]49
50    shortname = string.join(os.path.split(filenames[j])[-1].split('_')[-2:], '_')
51    title('Study %s: timeseries for %s (realisation %d of %d)'\
52          %(study, shortname, j, number_of_realisations-1))
53   
[3035]54    #raw_input('Next')
55
56
57raw_input('Stats')   
[3006]58 
[2972]59
[2986]60
61hold(False)
62# Plot spread of stage values for each timestep
[3035]63for i in range(0,project.number_of_timesteps):
[3006]64    # Plot histogram#
[2986]65
66    w = data[i,:]
67    plot(w, 'k.')
68    xlabel('realisations')
69    ylabel('stage (m)')
[3010]70    title('Study %s: spread at timestep %d of %d (t=%.2f)'\
71          %(study, i, project.number_of_timesteps-1, time[i]))       
[3035]72raw_input('Next')
[2972]73
[3035]74show()
75import sys; sys.exit()
[2975]76
[2986]77hold(False)
[2972]78# Plot histogram of stage values for each timestep
[3006]79for i in range(300,320):#project.number_of_timesteps):
[2972]80    # Plot histogram
81
82    w = data[i,:]
83     
84    bins = create_bins(w, project.number_of_bins)
85    print 'bins', bins
[2975]86    hist = histogram(w, bins, relative=True)
87    #print 'hist', hist
[2972]88   
89    #plot(w, 'k.')
[2975]90    plot(bins, hist, 'k.')
[2972]91    xlabel('stage (m)')
[2975]92    ylabel('frequency')
[3010]93    title('Study %s: histogram at timestep %d of %d (t=%.2f)'\
94          %(study, i, project.number_of_timesteps-1, time[i]))       
[3006]95    raw_input('Next')
[2972]96
97
[2967]98#title('Gauge %s' %name)
99#xlabel('time(s)')
100#ylabel('stage (m)')   
101#legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
102#savefig(name, dpi = 300)
103
104
[2972]105
106
107
108
[2967]109show()
110
111
112
113#from pylab import *
114#plot(time, stage)
115#show()
Note: See TracBrowser for help on using the repository browser.