source: production/pt_hedland_2006/get_building_inundation.py @ 3289

Last change on this file since 3289 was 3289, checked in by sexton, 18 years ago

update export coords for pthedland

File size: 4.0 KB
Line 
1""" Script for calculating maximum inundation at building locations
2
3    Inputs:
4
5    swwfile             - name of sww file
6                        - assume that all conserved quantities have been stored
7                   
8    buildings_filename  - name of file containing building data, sourced
9                          from NBED
10
11    time_min            - beginning of user defined time range
12                          - default will be first available time found in swwfile
13                       
14    time_max            - end of user defined time range
15                          - default will be last available time found in swwfile
16                       
17    Output:
18   
19    - building_filename (.csv) augmented to include maximum depth, momentum and
20      velocity and stored in same location as swwfile
21      Name = augmented_buildings.csv
22     
23    """
24
25from pyvolution.util import file_function
26import project
27from os import sep
28from math import sqrt
29
30# Inputs
31#timestampdir = '20060707_001859' # MSL
32#timestampdir = '20060707_003301' # HAT
33timestampdir = '20060707_003424' # LAT
34file_loc = project.outputdir + timestampdir + sep
35swwfile = file_loc + project.basename + '.sww'
36buildings_filename = project.buildings_filename
37time_min = None
38time_max = None
39
40def get_buildings_from_file(filename):
41    from coordinate_transforms.redfearn import redfearn
42    fid = open(filename)
43    lines = fid.readlines()
44    fid.close()
45    buildings = []
46    line1 = lines[0]
47    line11 = line1.split(',')
48    for i in range(len(line11)):
49        if line11[i].strip('\n').strip(' ') == 'LATITUDE': lat_index = i
50        if line11[i].strip('\n').strip(' ') == 'LONGITUDE': lon_index = i
51        if line11[i].strip('\n').strip(' ') == 'BUILDING_N': name_index = i
52    for line in lines[1:]:
53        fields = line.split(',')
54        lat = float(fields[lat_index])
55        lon = float(fields[lon_index])
56        z, easting, northing = redfearn(lat,lon)
57        buildings.append([easting, northing])
58        loc = fields[name_index]
59   
60    return buildings, lines
61
62print '\n Buildings obtained from: %s \n' %buildings_filename
63buildings, lines = get_buildings_from_file(buildings_filename)
64
65sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum']
66
67f = file_function(swwfile,
68                  quantities = sww_quantity,
69                  interpolation_points = buildings,
70                  verbose = True,
71                  use_cache = True)
72
73T = f.get_time()
74
75if time_min is None:
76    time_min = min(T)
77else:
78    if time_min < min(T):
79        msg = 'Minimum time entered not correct - please try again'
80        raise Exception, msg
81
82if time_max is None:
83    time_max = max(T)
84else:
85    if time_max > max(T):
86        msg = 'Maximum time entered not correct - please try again'
87        raise Exception, msg
88
89lines[0] = lines[0].strip() +\
90           ',MAX INUNDATION DEPTH (m)' +\
91           ',MAX MOMENTUM (m^2/s)' +\
92           ',MAX SPEED (m/s) \n'
93
94for k, g in enumerate(buildings):
95    print 'Building %d of %d' %(k, len(buildings))
96    max_depth = 0
97    max_momentum = 0
98    max_velocity = 0     
99    for i, t in enumerate(T): 
100            w = f(t, point_id = k)[0]
101            z = f(t, point_id = k)[1]
102            uh = f(t, point_id = k)[2]
103            vh = f(t, point_id = k)[3]
104            depth = w-z
105            m = sqrt(uh*uh + vh*vh)   
106            vel = m / (depth + 1.e-30) 
107            if depth > max_depth: max_depth = w-z
108            if m > max_momentum: max_momentum = m
109            if vel > max_velocity: max_velocity = vel
110
111    lines[k+1] = lines[k+1].strip() +\
112                 ',%f' %max_depth +\
113                 ',%f' %max_momentum +\
114                 ',%f\n' %max_velocity
115
116augbuildingsfile = file_loc + 'augmented_buildings.csv'
117
118fid = open(augbuildingsfile, 'w')
119for line in lines:
120    fid.write(line)
121fid.close()
122
123print '\n Augmented building file written to %s \n' %augbuildingsfile
Note: See TracBrowser for help on using the repository browser.