Changeset 3006


Ignore:
Timestamp:
May 29, 2006, 1:07:46 PM (18 years ago)
Author:
ole
Message:

Work on plot_spread.py

Location:
development/stochastic_study
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed
  • development/stochastic_study/plot_spread.py

    r2986 r3006  
    77# Related major packages
    88from Numeric import zeros, Float, allclose, arange
    9 from caching import cache
    109
    1110# Application specific imports
    1211from utilities.numerical_tools import histogram, create_bins
     12from read_realisations import read_realisations
    1313import project
    1414
    1515# Initialise
    16 gauge = '%s.txt' %project.gauge_names[0]
    17 number_of_realisations = project.number_of_realisations
    18 
    19 time = zeros(project.number_of_timesteps, Float)
    20 data = zeros((project.number_of_timesteps, number_of_realisations), Float)
    21 
    22 
    2316# Read in all realisations-timeseries
    24 
    25 j = 0 # Count realisations
    26 for filename in os.listdir(project.working_dir):
    27     if filename.startswith(project.basename) and filename.endswith(gauge):
    28         if j < data.shape[1]:       
    29             print 'Reading filename %s (column %d)' %(filename, j)
    30             fid = open(filename)
    31             for i, line in enumerate(fid.readlines()):
    32                 if i < data.shape[0]:
    33                     fields = line.strip().split()
    34                     time[i] = float(fields[0])
    35                     data[i,j] = float(fields[1])
    36         else:
    37             print 'Ignored %s (column %d)' %(filename, j)
    38 
    39         fid.close()
    40         j += 1   
     17time, data, filenames = read_realisations('cyclone1',
     18                                          max_realisations = 200,
     19                                          use_cache=True)
     20#time, data, filenames = read_realisations('test', exclude='_0_',
     21#                                          max_realisations = 200,
     22#                                          use_cache=True)
     23#time, data, filenames = read_realisations('nautilus1', exclude='_0_')
     24number_of_realisations = data.shape[1]
    4125
    4226
     
    4529
    4630ion()
    47 #hold(True)
    48 hold(False)
     31hold(True)
     32#hold(False)
    4933
    5034# Simple plot of timeseries for different realisations
    51 for j in range(number_of_realisations):
    52     print j, data[:,j]
     35for j in range(20):
     36    print j, filenames[j]
    5337    plot(time, data[:,j], 'k-')
    5438   
    5539    xlabel('time(s)')
    5640    ylabel('stage (m)')
    57     title('Realisation %d of %d' %(j, number_of_realisations-1))    #
    58 #raw_input('Next')
     41    title('Timeseries for realisation %d of %d' %(j, number_of_realisations-1))
     42    raw_input('Next')
     43 
    5944
    6045
    6146hold(False)
    6247# Plot spread of stage values for each timestep
    63 for i in range(200,project.number_of_timesteps):
    64     # Plot histogram
     48for i in range(300,320): #project.number_of_timesteps):
     49    # Plot histogram#
    6550
    6651    w = data[i,:]
     
    6853    xlabel('realisations')
    6954    ylabel('stage (m)')
    70     title('Timestep %d of %d (t=%.2f)'\
     55    title('Spread at timestep %d of %d (t=%.2f)'\
    7156          %(i, project.number_of_timesteps-1, time[i]))       
    7257    raw_input('Next')
     
    7560hold(False)
    7661# Plot histogram of stage values for each timestep
    77 for i in range(200,project.number_of_timesteps):
     62for i in range(300,320):#project.number_of_timesteps):
    7863    # Plot histogram
    7964
     
    8974    xlabel('stage (m)')
    9075    ylabel('frequency')
    91     title('Timestep %d of %d (t=%.2f)'\
     76    title('Histogram at timestep %d of %d (t=%.2f)'\
    9277          %(i, project.number_of_timesteps-1, time[i]))       
    93     #raw_input('Next')
     78    raw_input('Next')
    9479
    9580
  • development/stochastic_study/project.py

    r2989 r3006  
    3434blocksize = 10 #How many realisations to fit at a time
    3535
    36 number_of_bins = 10
     36number_of_bins = 100
Note: See TracChangeset for help on using the changeset viewer.