[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] | 5 | import sys, os, string |
---|
[2972] | 6 | |
---|
| 7 | # Related major packages |
---|
| 8 | from Numeric import zeros, Float, allclose, arange |
---|
| 9 | |
---|
| 10 | # Application specific imports |
---|
| 11 | from utilities.numerical_tools import histogram, create_bins |
---|
[3006] | 12 | from read_realisations import read_realisations |
---|
[2967] | 13 | import 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 |
---|
| 20 | study = 'cyclone1' #~4000 realisations, blocks of 100, parallel, stddev= 0.0006 |
---|
| 21 | #study = 'cyclone3' #~120 realisations, blocks of 10, parallel, stddev= 0.0013 |
---|
| 22 | |
---|
[3010] | 23 | time, 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] | 30 | number_of_realisations = data.shape[1] |
---|
[3035] | 31 | number_of_timesteps = data.shape[0] |
---|
[3011] | 32 | print 'Read %d realisations' %number_of_realisations |
---|
[2972] | 33 | |
---|
[3035] | 34 | |
---|
[2967] | 35 | # Plot |
---|
| 36 | from pylab import ion, hold, plot, title, xlabel, ylabel, legend, savefig, show |
---|
| 37 | |
---|
| 38 | ion() |
---|
[3006] | 39 | hold(True) |
---|
| 40 | #hold(False) |
---|
[2967] | 41 | |
---|
[3010] | 42 | |
---|
[3035] | 43 | for 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 | |
---|
| 57 | raw_input('Stats') |
---|
[3006] | 58 | |
---|
[2972] | 59 | |
---|
[2986] | 60 | |
---|
| 61 | hold(False) |
---|
| 62 | # Plot spread of stage values for each timestep |
---|
[3035] | 63 | for 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] | 72 | raw_input('Next') |
---|
[2972] | 73 | |
---|
[3035] | 74 | show() |
---|
| 75 | import sys; sys.exit() |
---|
[2975] | 76 | |
---|
[2986] | 77 | hold(False) |
---|
[2972] | 78 | # Plot histogram of stage values for each timestep |
---|
[3006] | 79 | for 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] | 109 | show() |
---|
| 110 | |
---|
| 111 | |
---|
| 112 | |
---|
| 113 | #from pylab import * |
---|
| 114 | #plot(time, stage) |
---|
| 115 | #show() |
---|