Changeset 1921


Ignore:
Timestamp:
Oct 14, 2005, 11:49:24 AM (18 years ago)
Author:
ole
Message:

Get timeseries - work for Karratha

File:
1 edited

Legend:

Unmodified
Added
Removed
  • production/karratha_2005/get_timeseries.py

    r1901 r1921  
    77
    88#swwfile = project.outputdir + project.basename + '.sww'
    9 swwfile = project.outputname + '_0.0tide' + '.sww'
     9#swwfile = project.outputname + '_0.0tide' + '.sww'
     10swwfile = project.outputname + '.sww'
     11#swwfile = project.outputname + '_0.0tide_notsunami.sww'
    1012#swwfile = project.outputdir + 'karratha_100m_12m.sww'
    1113
    1214#MOST Boundary directly north of Legendre island
    13 y = project.north
    14 x = (project.p3[1] + project.p4[1])/2   #?????//
    15 z, easting, northing  = redfearn(y, x)
     15lat = project.north
     16lon = (project.p3[1] + project.p4[1])/2
     17z, easting, northing  = redfearn(lat, lon)
    1618
    1719
    1820#From Neil
    19 gauges = [[easting, northing - 10],
     21gauges = [[easting, northing - 5],
     22          [easting, northing - 10],
     23          [easting, northing - 20],
     24          [easting, northing - 30],
     25          [easting, northing - 50],
     26          [easting, northing - 100],                                       
    2027          [470882, 7717952],
    2128          [469390, 7715426],
     
    2532          [470285, 7717470]]
    2633         
    27 
    28    
    29 #gauge_names = ['MOST',
    30 #               'Gauge_1', 'Gauge_2', 'Gauge_3',
    31 #               'Gauge_4', 'Gauge_5', 'Gauge_6']
    32  
    3334
    3435#Read model output
     
    4041
    4142f = cache(file_function, swwfile,
    42           {'quantities': 'stage',
     43          {'quantities': ['stage', 'elevation', 'xmomentum', 'ymomentum'],
    4344           'interpolation_points': gauges,
    4445           'verbose': True},
     
    5152
    5253
     54from math import sqrt
     55for k, g in enumerate(gauges):
    5356
    54 for k, g in enumerate(gauges):
    55     model_values = []
     57    stages = []
     58    elevations = []
     59    momenta = []
     60   
    5661    for i, t in enumerate(f.T):
    57         val = f(t, point_id = k)[0]
    58         model_values.append(val)
     62        w = f(t, point_id = k)[0]
     63        z = f(t, point_id = k)[1]
     64        uh = f(t, point_id = k)[2]
     65        vh = f(t, point_id = k)[3]       
     66       
     67        stages.append(w)
     68        elevations.append(z)  #Should be constant over time
     69        momenta.append( sqrt(uh*uh + vh*vh) )  #Absolute momentum
    5970
    6071    from pylab import *
     72   
    6173    ion()
    6274    hold(False)
    63     plot(model_time, model_values, '-b')
    64     name = 'Gauge_%d: %s' %(k, str(g))
     75
     76    if elevations[0] < -10:
     77        plot(model_time, stages, '-b')
     78    else:   
     79        plot(model_time, stages, '-b',
     80             model_time, elevations, '-k')
     81    name = 'Gauge_%d: (%.1f, %.1f)' %(k, g[0], g[1])
    6582    title(name)
    66    
    67     savefig('Gauge_%d' %k)
     83
     84    title('%s (stage)' %name)
     85    xlabel('time [s]')
     86    ylabel('elevation [m]')   
     87    legend(('Stage', 'Bed = %.1f' %elevations[0]),
     88           shadow=True,
     89           loc='upper right')
     90    savefig('Gauge_%d_stage' %k)
    6891
    6992    raw_input('Next')
     93
     94
     95    #Momentum plot
     96    ion()
     97    hold(False)
     98    plot(model_time, momenta, '-r')
     99    title(name)
     100
     101    title('%s (momentum)' %name)
     102    xlabel('time [s]')
     103    ylabel('sqrt( uh^2 + vh^2 ) [m^2/s]')   
     104    savefig('Gauge_%d_momentum' %k)
     105
     106    raw_input('Next')
     107   
     108
     109   
    70110show()
    71111
Note: See TracChangeset for help on using the changeset viewer.