source: production/pt_hedland_2006/get_building_inundation.py @ 3290

""" Script for calculating maximum inundation at building locations

    Inputs:

    swwfile             - name of sww file
                        - assume that all conserved quantities have been stored
                    
    buildings_filename  - name of file containing building data, sourced
                          from NBED

    time_min            - beginning of user defined time range 
                          - default will be first available time found in swwfile
                        
    time_max            - end of user defined time range
                          - default will be last available time found in swwfile
                        
    Output:
    
    - building_filename (.csv) augmented to include maximum depth, momentum and
      velocity and stored in same location as swwfile
      Name = augmented_buildings.csv
      
    """

from pyvolution.util import file_function
import project
from os import sep
from math import sqrt

# Inputs
#timestampdir = '20060707_001859' # MSL
#timestampdir = '20060707_003301' # HAT
timestampdir = '20060707_003424' # LAT
file_loc = project.outputdir + timestampdir + sep 
swwfile = file_loc + project.basename + '.sww'
buildings_filename = project.buildings_filename
time_min = None
time_max = None

def get_buildings_from_file(filename):
    from coordinate_transforms.redfearn import redfearn
    fid = open(filename)
    lines = fid.readlines()
    fid.close()
    buildings = []
    line1 = lines[0]
    line11 = line1.split(',')
    for i in range(len(line11)):
        if line11[i].strip('\n').strip(' ') == 'LATITUDE': lat_index = i
        if line11[i].strip('\n').strip(' ') == 'LONGITUDE': lon_index = i
        if line11[i].strip('\n').strip(' ') == 'BUILDING_N': name_index = i
    for line in lines[1:]:
        fields = line.split(',')
        lat = float(fields[lat_index])
        lon = float(fields[lon_index])
        z, easting, northing = redfearn(lat,lon)
        buildings.append([easting, northing])
        loc = fields[name_index]
    
    return buildings, lines

print '\n Buildings obtained from: %s \n' %buildings_filename
buildings, lines = get_buildings_from_file(buildings_filename)

sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum']

f = file_function(swwfile,
                  quantities = sww_quantity,
                  interpolation_points = buildings,
                  verbose = True,
                  use_cache = True)

T = f.get_time()

if time_min is None:
    time_min = min(T)
else:
    if time_min < min(T):
        msg = 'Minimum time entered not correct - please try again'
        raise Exception, msg

if time_max is None:
    time_max = max(T)
else:
    if time_max > max(T):
        msg = 'Maximum time entered not correct - please try again'
        raise Exception, msg

lines[0] = lines[0].strip() +\
           ',MAX INUNDATION DEPTH (m)' +\
           ',MAX MOMENTUM (m^2/s)' +\
           ',MAX SPEED (m/s) \n'

for k, g in enumerate(buildings):
    print 'Building %d of %d' %(k, len(buildings))
    max_depth = 0
    max_momentum = 0
    max_velocity = 0      
    for i, t in enumerate(T): 
            w = f(t, point_id = k)[0]
            z = f(t, point_id = k)[1]
            uh = f(t, point_id = k)[2]
            vh = f(t, point_id = k)[3]
            depth = w-z
            m = sqrt(uh*uh + vh*vh)   
            vel = m / (depth + 1.e-30) 
            if depth > max_depth: max_depth = w-z
            if m > max_momentum: max_momentum = m
            if vel > max_velocity: max_velocity = vel

    lines[k+1] = lines[k+1].strip() +\
                 ',%f' %max_depth +\
                 ',%f' %max_momentum +\
                 ',%f\n' %max_velocity

augbuildingsfile = file_loc + 'augmented_buildings.csv'

fid = open(augbuildingsfile, 'w')
for line in lines:
    fid.write(line)
fid.close()

print '\n Augmented building file written to %s \n' %augbuildingsfile