Changeset 2069


Ignore:
Timestamp:
Nov 24, 2005, 6:12:22 PM (18 years ago)
Author:
ole
Message:

Karratha work

Files:
3 edited

Legend:

Unmodified
Added
Removed
  • inundation/pyvolution/data_manager.py

    r2064 r2069  
    362362        from os import stat
    363363
    364         minimum_allowed_depth = 0.001
    365         #minimum_allowed_depth = 0.0  #FIXME pass in or read from domain
     364        #minimum_allowed_depth = 0.001
     365        minimum_allowed_depth = 0.0  #FIXME pass in or read from domain
    366366        from Numeric import choose
    367367
  • production/karratha_2005/get_timeseries.py

    r1987 r2069  
    1 """Read in sww file, interpolate at specified locations and plot
     1"""Read in sww file, interpolate at specified locations and plot time series
     2Additionally, store augmented building database with max inundation depths
    23"""
    34
     
    67from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
    78from pylab import *
    8    
    99
    1010
    11 #swwfile = project.outputdir + project.basename + '.sww'
    12 #swwfile = project.outputname + '_0.0tide' + '.sww'
     11plot = False
     12   
    1313swwfile = project.outputname + '.sww'
    14 #swwfile = project.outputname + '_notsunami.sww'
    15 #swwfile = project.outputdir + 'karratha_100m_12m.sww'
    1614
    1715#MOST Boundary directly north of Legendre island
     
    3230    for line in lines[1:]:
    3331        fields = line.split(',')
    34         lon = float(fields[5])
    35         lat = float(fields[6])
     32        lon = float(fields[4])
     33        lat = float(fields[5])
    3634
    3735        z, easting, northing  = redfearn(lat, lon)       
    3836        gauges.append([easting, northing])
    3937
    40     return gauges
     38    #Return gauges and raw data for subsequent storage   
     39    return gauges, lines
    4140
    42 gauges = get_gauges_from_file(project.gauge_filename)
    43 
    44 
    45 #gauges = gauges[-10:]
    46 #print gauges
    47 #import sys; sys.exit()
    48 
    49 #From Neil
    50 #gauges = [[easting, northing - 5],
    51 #          [easting, northing - 10],
    52 #          [easting, northing - 50],
    53 #          [easting, northing - 100],                                       
    54 #          [470882, 7717952],
    55 #          [469390, 7715426],
    56 #          [469214, 7714979],
    57 #          [468899, 7715177],
    58 #          [469038, 7715725],
    59 #          [470285, 7717470]]
    60          
     41gauges, buildings = get_gauges_from_file(project.gauge_filename)
    6142
    6243#Read model output
     
    6849                  use_cache = True)
    6950
    70 print f
    71 #model_time = f.T
    7251
    73 
     52buildings[0] = buildings[0].strip() +\
     53               ',MAX INUNDATION DEPTH (m)' +\
     54               ',MAX MOMENTUM (m^2/s)\n'           
     55               
    7456from math import sqrt
     57N = len(gauges)
    7558for k, g in enumerate(gauges):
     59    if k%((N+10)/10)==0:
     60        print 'Doing row %d of %d' %(k, N)
    7661
    7762    model_time = []
     
    8065    momenta = []
    8166
     67    max_depth = 0
     68    max_momentum = 0   
     69    for i, t in enumerate(f.T):
     70        w = f(t, point_id = k)[0]
     71        z = f(t, point_id = k)[1]
     72        uh = f(t, point_id = k)[2]
     73        vh = f(t, point_id = k)[3]       
    8274
    83     max_depth = 0
    84     for i, t in enumerate(f.T):
    85         if tmin < t < tmax:
    86             w = f(t, point_id = k)[0]
    87             z = f(t, point_id = k)[1]
    88             uh = f(t, point_id = k)[2]
    89             vh = f(t, point_id = k)[3]       
     75        m = sqrt(uh*uh + vh*vh)   #Absolute momentum
     76       
     77        model_time.append(t)       
     78        stages.append(w)
     79        elevations.append(z)  #Should be constant over time
     80        momenta.append(m)
    9081
    91             model_time.append(t)       
    92             stages.append(w)
    93             elevations.append(z)  #Should be constant over time
    94             momenta.append( sqrt(uh*uh + vh*vh) )  #Absolute momentum
     82        if w-z > max_depth:
     83            max_depth = w-z
     84        if m > max_momentum:           
     85            max_momentum = m   
    9586
    96             if w-z > max_depth:
    97                 max_depth = w-z
     87    #Augment building data
     88    buildings[k+1] = buildings[k+1].strip() +\
     89                     ',%f' %max_depth +\
     90                     ',%f\n' %max_momentum
     91   
     92    #print buildings[k]
     93   
     94    if plot is True:
     95        #Plot only those gauges that have been inundated by more than a threshold
     96        if max_depth < 0.2:
     97            print 'Skipping gauge %d' %k
     98            continue
    9899
    99                
     100        #FIXME Reduce to tmin and tmax
     101        #i0 = model_timefind
     102   
     103        ion()
     104        hold(False)
    100105
    101     #Plot only those gauges that have been inundated by more than a threshold
    102     if max_depth < 0.2:
    103         print 'Skipping gauge %d' %k
    104         continue
    105    
    106     ion()
    107     hold(False)
     106        if elevations[0] < -10:
     107            plot(model_time, stages, '-b')
     108        else:   
     109            plot(model_time, stages, '-b',
     110                 model_time, elevations, '-k')
     111        name = 'Gauge_%d: (%.1f, %.1f)' %(k, g[0], g[1])
     112        title(name)
    108113
    109     if elevations[0] < -10:
    110         plot(model_time, stages, '-b')
    111     else:   
    112         plot(model_time, stages, '-b',
    113              model_time, elevations, '-k')
    114     name = 'Gauge_%d: (%.1f, %.1f)' %(k, g[0], g[1])
    115     title(name)
     114        title('%s (stage)' %name)
     115        xlabel('time [s]')
     116        ylabel('elevation [m]')   
     117        legend(('Stage', 'Bed = %.1f' %elevations[0]),
     118               shadow=True,
     119               loc='upper right')
     120        savefig('Gauge_%d_stage' %k)
    116121
    117     title('%s (stage)' %name)
    118     xlabel('time [s]')
    119     ylabel('elevation [m]')   
    120     legend(('Stage', 'Bed = %.1f' %elevations[0]),
    121            shadow=True,
    122            loc='upper right')
    123     savefig('Gauge_%d_stage' %k)
    124 
    125     raw_input('Next')
     122        raw_input('Next')
    126123
    127124
    128     #Momentum plot
    129     ion()
    130     hold(False)
    131     plot(model_time, momenta, '-r')
    132     title(name)
     125        #Momentum plot
     126        ion()
     127        hold(False)
     128        plot(model_time, momenta, '-r')
     129        title(name)
    133130
    134     title('%s (momentum)' %name)
    135     xlabel('time [s]')
    136     ylabel('sqrt( uh^2 + vh^2 ) [m^2/s]')   
    137     savefig('Gauge_%d_momentum' %k)
     131        title('%s (momentum)' %name)
     132        xlabel('time [s]')
     133        ylabel('sqrt( uh^2 + vh^2 ) [m^2/s]')   
     134        savefig('Gauge_%d_momentum' %k)
    138135
    139     raw_input('Next')
     136        raw_input('Next')
    140137   
    141138
    142    
    143 show()
     139#Store new building file with mak depths added
    144140
     141FN = 'inundation_augmented_' + project.gauge_filename
     142fid = open(FN, 'w')
     143for line in buildings:
     144    fid.write(line)
     145fid.close()   
     146
     147
     148if plot is True:   
     149    show()
     150
  • production/karratha_2005/project.py

    r2017 r2069  
    3131outputname = outputdir + basename  #Used by post processing
    3232
    33 gauge_filename = 'all_bld_ind.csv'
     33#gauge_filename = 'all_bld_ind.csv'
     34gauge_filename = 'karratha_buildings_23112005.csv'
    3435
    3536
Note: See TracChangeset for help on using the changeset viewer.