| 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 | |
|---|
| 25 | from pyvolution.util import file_function |
|---|
| 26 | import project |
|---|
| 27 | from os import sep |
|---|
| 28 | from math import sqrt |
|---|
| 29 | |
|---|
| 30 | # Inputs |
|---|
| 31 | #timestampdir = '20060704_063005' # HAT |
|---|
| 32 | #timestampdir = '20060706_235246' # LAT |
|---|
| 33 | timestampdir = '20060704_063234' # MSL |
|---|
| 34 | #timestampdir = '20060515_001733' # MSL DTED |
|---|
| 35 | file_loc = project.outputdir + timestampdir + sep |
|---|
| 36 | swwfile = file_loc + project.basename + '.sww' |
|---|
| 37 | buildings_filename = project.buildings_filename |
|---|
| 38 | time_min = None |
|---|
| 39 | time_max = None |
|---|
| 40 | |
|---|
| 41 | def get_buildings_from_file(filename,bounding_polygon): |
|---|
| 42 | from coordinate_transforms.redfearn import redfearn |
|---|
| 43 | fid = open(filename) |
|---|
| 44 | lines = fid.readlines() |
|---|
| 45 | fid.close() |
|---|
| 46 | buildings = [] |
|---|
| 47 | line1 = lines[0] |
|---|
| 48 | line11 = line1.split(',') |
|---|
| 49 | for i in range(len(line11)): |
|---|
| 50 | if line11[i].strip('\n').strip(' ') == 'LATITUDE': lat_index = i |
|---|
| 51 | if line11[i].strip('\n').strip(' ') == 'LONGITUDE': lon_index = i |
|---|
| 52 | if line11[i].strip('\n').strip(' ') == 'BUILDING_N': name_index = i |
|---|
| 53 | for line in lines[1:]: |
|---|
| 54 | fields = line.split(',') |
|---|
| 55 | lat = float(fields[lat_index]) |
|---|
| 56 | lon = float(fields[lon_index]) |
|---|
| 57 | z, easting, northing = redfearn(lat,lon) |
|---|
| 58 | utm_pt = [easting, northing] |
|---|
| 59 | buildings.append(utm_pt) |
|---|
| 60 | loc = fields[name_index] |
|---|
| 61 | |
|---|
| 62 | return buildings, lines |
|---|
| 63 | |
|---|
| 64 | print '\n Buildings obtained from: %s \n' %buildings_filename |
|---|
| 65 | buildings, lines = get_buildings_from_file(buildings_filename,project.polyAll) |
|---|
| 66 | |
|---|
| 67 | sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum'] |
|---|
| 68 | |
|---|
| 69 | f = file_function(swwfile, |
|---|
| 70 | quantities = sww_quantity, |
|---|
| 71 | interpolation_points = buildings, |
|---|
| 72 | verbose = True, |
|---|
| 73 | use_cache = True) |
|---|
| 74 | |
|---|
| 75 | T = f.get_time() |
|---|
| 76 | |
|---|
| 77 | if time_min is None: |
|---|
| 78 | time_min = min(T) |
|---|
| 79 | else: |
|---|
| 80 | if time_min < min(T): |
|---|
| 81 | msg = 'Minimum time entered not correct - please try again' |
|---|
| 82 | raise Exception, msg |
|---|
| 83 | |
|---|
| 84 | if time_max is None: |
|---|
| 85 | time_max = max(T) |
|---|
| 86 | else: |
|---|
| 87 | if time_max > max(T): |
|---|
| 88 | msg = 'Maximum time entered not correct - please try again' |
|---|
| 89 | raise Exception, msg |
|---|
| 90 | |
|---|
| 91 | lines[0] = lines[0].strip() +\ |
|---|
| 92 | ',MAX INUNDATION DEPTH (m)' +\ |
|---|
| 93 | ',MAX MOMENTUM (m^2/s) \n' |
|---|
| 94 | |
|---|
| 95 | from utilities.polygon import inside_polygon |
|---|
| 96 | N = len(buildings) |
|---|
| 97 | for k, g in enumerate(buildings): |
|---|
| 98 | if k%((N+10)/10)==0: print 'Building %d of %d' %(k, N) |
|---|
| 99 | max_depth = 0.0 |
|---|
| 100 | max_momentum = 0.0 |
|---|
| 101 | max_velocity = 0.0 |
|---|
| 102 | zero_depth = 0.0 |
|---|
| 103 | zero_momentum = 0.0 |
|---|
| 104 | #if inside_polygon(g,project.polyAll) == True: |
|---|
| 105 | for i, t in enumerate(T): |
|---|
| 106 | w = f(t, point_id = k)[0] |
|---|
| 107 | z = f(t, point_id = k)[1] |
|---|
| 108 | uh = f(t, point_id = k)[2] |
|---|
| 109 | vh = f(t, point_id = k)[3] |
|---|
| 110 | depth = w-z |
|---|
| 111 | m = sqrt(uh*uh + vh*vh) |
|---|
| 112 | #vel = m / (depth + 1.e-30) |
|---|
| 113 | if depth > max_depth: max_depth = depth |
|---|
| 114 | if m > max_momentum: max_momentum = m |
|---|
| 115 | #else: |
|---|
| 116 | # max_depth = 0.0 |
|---|
| 117 | # max_momentum = 0.0 |
|---|
| 118 | print 'gauge ', g |
|---|
| 119 | print 'max depth', max_depth |
|---|
| 120 | lines[k+1] = lines[k+1].strip() +\ |
|---|
| 121 | ',%f' %max_depth +\ |
|---|
| 122 | ',%f\n' %max_momentum |
|---|
| 123 | |
|---|
| 124 | |
|---|
| 125 | augbuildingsfile = file_loc + 'augmented_buildings.csv' |
|---|
| 126 | |
|---|
| 127 | fid = open(augbuildingsfile, 'w') |
|---|
| 128 | for line in lines: |
|---|
| 129 | fid.write(line) |
|---|
| 130 | fid.close() |
|---|
| 131 | |
|---|
| 132 | print '\n Augmented building file written to %s \n' %augbuildingsfile |
|---|
