source: trunk/anuga_core/source/anuga/damage_modelling/inundation_damage.py @ 7753

Last change on this file since 7753 was 7743, checked in by hudson, 14 years ago

Removed more modules from data_handler: code to do with building destruction.
File size: 18.8 KB
Line 
1"""Classes for implementing damage curves and calculating financial damage
2
3   Duncan Gray, Ole Nielsen, Jane Sexton, Nick Bartzis
4   Geoscience Australia, 2006
5"""
6import os
7from math import sqrt
8from Scientific.Functions.Interpolation import InterpolatingFunction
9from random import choice
10from types import StringType
11
12import numpy as num
13
14
15try: 
16    import kinds 
17except ImportError: 
18    # Hand-built mockup of the things we need from the kinds package, since it
19    # was recently removed from the standard numeric distro.  Some users may 
20    # not have it by default. 
21    class _bunch: 
22        pass 
23         
24    class _kinds(_bunch): 
25        default_float_kind = _bunch() 
26        default_float_kind.MIN = 2.2250738585072014e-308  #smallest +ve number
27        default_float_kind.MAX = 1.7976931348623157e+308 
28     
29    kinds = _kinds()
30   
31
32from anuga.utilities.numerical_tools import ensure_numeric
33from exposure import Exposure
34from anuga.abstract_2d_finite_volumes.util import file_function
35from anuga.geospatial_data.geospatial_data import ensure_absolute
36from anuga.utilities.numerical_tools import NAN
37import anuga.utilities.log as log
38from config import epsilon
39depth_epsilon = epsilon
40
41# Change these if the ouput from nexix changes
42SHORE_DIST_LABEL = 'SHORE_DIST'
43WALL_TYPE_LABEL = 'WALL_TYPE'
44STR_VALUE_LABEL = 'STR_VALUE'
45CONT_VALUE_LABEL = 'CONT_VALUE'
46
47def inundation_damage(sww_base_name, exposure_files_in,
48                      exposure_file_out_marker=None,
49                      ground_floor_height=0.3,
50                      overwrite=False, verbose=True,
51                                 use_cache = True):
52    """
53    This is the main function for calculating tsunami damage due to
54    inundation.  It gets the location of structures from the exposure
55    file and gets the inundation of these structures from the
56    sww file.
57
58    It then calculates the damage loss.
59
60    Note, structures outside of the sww file get the minimum inundation
61    (-ground_floor_height).
62   
63    These calculations are done over all the sww files with the sww_base_name
64    in the specified directory.
65
66    exposure_files_in - a file or a list of files to input from
67    exposure_file_out_marker -  this string will be added to the input file
68                                name to get the output file name
69    """
70    if type(exposure_files_in) == StringType:
71        exposure_files_in = [exposure_files_in]
72
73
74    for exposure_file_in in exposure_files_in:
75        csv = Exposure(exposure_file_in,
76                           title_check_list=[SHORE_DIST_LABEL,WALL_TYPE_LABEL,
77                                             STR_VALUE_LABEL,CONT_VALUE_LABEL])
78        geospatial = csv.get_location()
79        geospatial = ensure_absolute(geospatial)
80        max_depths, max_momentums = calc_max_depth_and_momentum(sww_base_name,
81                        geospatial,
82                        ground_floor_height=ground_floor_height,
83                        verbose=verbose,
84                        use_cache=use_cache)
85        edm = EventDamageModel(max_depths,
86                               csv.get_column(SHORE_DIST_LABEL),
87                               csv.get_column(WALL_TYPE_LABEL),
88                               csv.get_column(STR_VALUE_LABEL),
89                               csv.get_column(CONT_VALUE_LABEL)
90                               )
91        results_dic = edm.calc_damage_and_costs(verbose_csv=True,
92                                                verbose=verbose)
93        for title, value in results_dic.iteritems():
94            csv.set_column(title, value, overwrite=overwrite)
95   
96        # Save info back to csv file
97        if exposure_file_out_marker == None:
98            exposure_file_out = exposure_file_in
99        else:
100            # split off extension, in such a way to deal with more than one '.' in the name of file
101            split_name = exposure_file_in.split('.')
102            exposure_file_out =  '.'.join(split_name[:-1]) + exposure_file_out_marker + \
103                                '.' + split_name[-1]
104        csv.save(exposure_file_out)
105        if verbose: log.critical('Augmented building file written to %s'
106                                 % exposure_file_out)
107   
108def add_depth_and_momentum2csv(sww_base_name, exposure_file_in,
109                      exposure_file_out=None,
110                      overwrite=False, verbose=True,
111                                 use_cache = True):
112    """
113    Calculate the maximum depth and momemtum in an sww file, for locations
114    specified in an csv exposure file.
115   
116    These calculations are done over all the sww files with the sww_base_name
117    in the specified directory.
118    """
119
120    csv = Exposure(exposure_file_in)
121    geospatial = csv.get_location()
122    max_depths, max_momentums = calc_max_depth_and_momentum(sww_base_name,
123                                                          geospatial,
124                                                          verbose=verbose,
125                                                          use_cache=use_cache)
126    csv.set_column("MAX INUNDATION DEPTH (m)",max_depths, overwrite=overwrite)
127    csv.set_column("MOMENTUM (m^2/s) ",max_momentums, overwrite=overwrite)
128    csv.save(exposure_file_out)
129   
130def calc_max_depth_and_momentum(sww_base_name, points,
131                                ground_floor_height=0.0,
132                                verbose=True,
133                                 use_cache = True):
134    """
135    Calculate the maximum inundation height above ground floor for a list
136    of locations.
137
138    The inundation value is in the range -ground_floor_height to
139    overflow errors.
140
141    These calculations are done over all the sww files with the sww_base_name
142    in the specified directory.
143    """
144    quantities =  ['stage', 'elevation', 'xmomentum', 'ymomentum']
145    points = ensure_absolute(points)
146    point_count = len(points)
147
148    # initialise the max lists
149    max_depths = [-ground_floor_height]*point_count
150    max_momentums = [-ground_floor_height]*point_count
151   
152    # How many sww files are there?
153    dir, base = os.path.split(sww_base_name)
154    if base[-4:] == '.sww':
155        base = base[:-4]
156    if dir == "": dir = "." # Unix compatibility
157    dir_ls = os.listdir(dir)
158    interate_over = [x for x in dir_ls if base in x and x[-4:] == '.sww']
159    if len(interate_over) == 0:
160        msg = 'No files of the base name %s.'\
161              %(sww_base_name)
162        raise IOError, msg
163    from os import sep
164    for this_sww_file in interate_over:
165        callable_sww = file_function(dir+sep+this_sww_file,
166                                     quantities=quantities,
167                                     interpolation_points=points,
168                                     verbose=verbose,
169                                     use_cache=use_cache)
170   
171        for point_i, point in enumerate(points):
172            for time in callable_sww.get_time():
173                quantity_values = callable_sww(time,point_i)
174           
175                w = quantity_values[0]
176                z = quantity_values[1]
177                uh = quantity_values[2] 
178                vh = quantity_values[3]
179               
180                #  -ground_floor_height is the minimum value.
181                depth = w - z - ground_floor_height
182             
183                if depth > max_depths[point_i]:
184                    max_depths[point_i] = depth
185                if w == NAN or z == NAN or uh == NAN or vh == NAN:
186                    continue
187                momentum = sqrt(uh*uh + vh*vh)
188                if momentum > max_momentums[point_i]:
189                    max_momentums[point_i] = momentum
190    return max_depths, max_momentums
191
192class EventDamageModel:
193    """
194    Object for working out the damage and cost
195
196    """
197    STRUCT_LOSS_TITLE = "STRUCT_LOSS_$"#"Structure Loss ($)"
198    CONTENTS_LOSS_TITLE = "CONTENTS_LOSS_$"#"Contents Loss ($)"
199    CONTENTS_DAMAGE_TITLE = "CONTENTS_DAMAGE_fraction"#"Contents damaged (fraction)"
200    STRUCT_DAMAGE_TITLE = "STRUCT_DAMAGE_fraction" #"Structure damaged (fraction)"
201    COLLAPSE_CSV_INFO_TITLE = "COLLAPSE_CSV_INFO"#"Calculation notes"
202    MAX_DEPTH_TITLE = "MAX_DEPTH_m" #"Inundation height above ground floor (m)"
203    STRUCT_COLLAPSED_TITLE = "STRUCT_COLLAPSED"#"collapsed structure if 1"
204    STRUCT_INUNDATED_TITLE = "STRUCT_INUNDATED"#"inundated structure if 1"
205    double_brick_damage_array = num.array([[-kinds.default_float_kind.MAX, 0.0],
206                                           [0.0-depth_epsilon, 0.0],
207                                           [0.0,0.016],
208                                           [0.1,0.150],
209                                           [0.3,0.425],
210                                           [0.5,0.449],
211                                           [1.0,0.572],
212                                           [1.5,0.582],
213                                           [2.0,0.587],
214                                           [2.5,0.647],
215                                           [kinds.default_float_kind.MAX,64.7]])
216    double_brick_damage_curve = InterpolatingFunction( \
217        (num.ravel(double_brick_damage_array[:,0:1]),),
218        num.ravel(double_brick_damage_array[:,1:]))
219   
220    brick_veeer_damage_array = num.array([[-kinds.default_float_kind.MAX, 0.0],
221                                          [0.0-depth_epsilon, 0.0],
222                                          [0.0,0.016],
223                                          [0.1,0.169],
224                                          [0.3,0.445],
225                                          [0.5,0.472],
226                                          [1.0,0.618],
227                                          [1.5,0.629],
228                                          [2.0,0.633],
229                                          [2.5,0.694],
230                                          [kinds.default_float_kind.MAX,69.4]])
231    brick_veeer_damage_curve = InterpolatingFunction( \
232        (num.ravel(brick_veeer_damage_array[:,0:1]),),
233        num.ravel(brick_veeer_damage_array[:,1:]))
234    struct_damage_curve = {'Double Brick':double_brick_damage_curve,
235                           'Brick Veneer':brick_veeer_damage_curve}
236    default_struct_damage_curve = brick_veeer_damage_curve
237
238    contents_damage_array = num.array([[-kinds.default_float_kind.MAX, 0.0],
239                                       [0.0-depth_epsilon, 0.0],
240                                       [0.0,0.013],
241                                       [0.1,0.102],
242                                       [0.3,0.381],
243                                       [0.5,0.500],
244                                       [1.0,0.970],
245                                       [1.5,0.976],
246                                       [2.0,0.986],
247                                       [kinds.default_float_kind.MAX,98.6]])
248    contents_damage_curve = InterpolatingFunction( \
249        (num.ravel(contents_damage_array[:,0:1]),),
250        num.ravel(contents_damage_array[:,1:]))
251
252    #building collapse probability
253    # inundation depth above ground floor, m
254    depth_upper_limits = [depth_epsilon, 1.0, 2.0, 3.0, 5.0,
255                          kinds.default_float_kind.MAX]
256    # shore mistance, m
257    shore_upper_limits = [125,200,250, kinds.default_float_kind.MAX]
258    # Building collapse probability
259    collapse_probability = [[0.0, 0.0, 0.0, 0.0], #Code below assumes 0.0
260                            [0.05, 0.02, 0.01, 0.0],
261                            [0.6, 0.3, 0.1, 0.05],
262                            [0.8, 0.4, 0.25, 0.15],
263                            [0.95, 0.7, 0.5, 0.3],
264                            [0.99, 0.9, 0.65, 0.45]]
265   
266    def __init__(self,max_depths, shore_distances, walls,
267                 struct_costs, content_costs):
268        """
269        max depth is Inundation height above ground floor (m), so
270                  the ground floor has been taken into account.
271        """
272        self.max_depths = [float(x) for x in max_depths]
273        self.shore_distances = [float(x) for x in shore_distances]
274        self.walls = walls
275        self.struct_costs = [float(x) for x in struct_costs]
276        self.content_costs = [float(x) for x in content_costs]
277
278        self.structure_count = len(self.max_depths)
279        #Fixme expand
280        assert self.structure_count == len(self.shore_distances)
281        assert  self.structure_count == len(self.walls)
282        assert  self.structure_count == len(self.struct_costs)
283        assert  self.structure_count == len(self.content_costs)
284        #assert  self.structure_count == len(self.)
285       
286    def calc_damage_and_costs(self, verbose_csv=False, verbose=False):
287        """
288        This is an overall method to calculate the % damage and collapsed
289        structures and then the $ loss.
290        """
291        self.calc_damage_percentages()
292        collapse_probability = self.calc_collapse_probability()
293        self._calc_collapse_structures(collapse_probability,
294                                  verbose_csv=verbose_csv)
295        self.calc_cost()
296        results_dict = {self.STRUCT_LOSS_TITLE:self.struct_loss
297                        ,self.STRUCT_DAMAGE_TITLE:self.struct_damage
298                        ,self.CONTENTS_LOSS_TITLE:self.contents_loss
299                        ,self.CONTENTS_DAMAGE_TITLE:self.contents_damage
300                        ,self.MAX_DEPTH_TITLE:self.max_depths
301                        ,self.STRUCT_COLLAPSED_TITLE:self.struct_collapsed
302                        ,self.STRUCT_INUNDATED_TITLE:self.struct_inundated
303                        }
304        if verbose_csv:
305           results_dict[self.COLLAPSE_CSV_INFO_TITLE] = self.collapse_csv_info
306        return results_dict
307           
308    def calc_damage_percentages(self):
309        """
310        Using stage curves calc the damage to structures and contents
311        """
312
313        # the data being created
314        struct_damage = num.zeros(self.structure_count, num.float)
315        contents_damage = num.zeros(self.structure_count, num.float)
316        self.struct_inundated = ['']* self.structure_count
317
318        for i,max_depth,shore_distance,wall in map(None,
319                                                   range(self.structure_count),
320                                                   self.max_depths,
321                                                   self.shore_distances,
322                                                   self.walls):
323            ## WARNING SKIP IF DEPTH < 0.0
324            if 0.0 > max_depth:
325                continue
326           
327            # The definition of inundated is if the max_depth is > 0.0
328            self.struct_inundated[i] = 1.0
329           
330            #calc structural damage %
331            damage_curve = self.struct_damage_curve.get(wall,
332                                              self.default_struct_damage_curve)
333            struct_damage[i] = damage_curve(max_depth)
334            contents_damage[i] = self.contents_damage_curve(max_depth)
335           
336        self.struct_damage = struct_damage
337        self.contents_damage = contents_damage
338           
339       
340    def calc_cost(self):
341        """
342        Once the damage has been calculated, determine the $ cost.
343        """
344        # ensure_numeric does not cut it.
345        self.struct_loss = self.struct_damage * \
346                           ensure_numeric(self.struct_costs)
347        self.contents_loss = self.contents_damage * \
348                           ensure_numeric(self.content_costs)
349       
350    def calc_collapse_probability(self):
351        """
352        return a dict of which structures have x probability of collapse.
353             key is collapse probability
354             value is list of struct indexes with key probability of collapse
355        """
356        # I could've done this is the calc_damage_percentages and
357        # Just had one loop.
358        # But for ease of testing and bug finding I'm seperating the loops.
359        # I'm make the outer loop for both of them the same though,
360        # so this loop can easily be folded into the other loop.
361       
362        # dict of which structures have x probability of collapse.
363        # key of collapse probability
364        # value of list of struct indexes
365        struct_coll_prob = {}
366       
367        for i,max_depth,shore_distance,wall in map(None,
368                                                   range(self.structure_count),
369                                                   self.max_depths,
370                                                   self.shore_distances,
371                                                   self.walls):
372            # WARNING ASSUMING THE FIRST BIN OF DEPTHS GIVE A ZERO PROBABILITY
373            depth_upper_limits = self.depth_upper_limits
374            shore_upper_limits = self.shore_upper_limits
375            collapse_probability = self.collapse_probability
376            if max_depth <= depth_upper_limits[0]:
377                continue
378            start = 1
379            for i_depth, depth_limit in enumerate(depth_upper_limits[start:]):
380                #Have to change i_depth so it indexes into the lists correctly
381                i_depth += start
382                if max_depth <= depth_limit:
383                    for i_shore, shore_limit in enumerate(shore_upper_limits):
384                        if shore_distance <= shore_limit:
385                            coll_prob = collapse_probability[i_depth][i_shore]
386                            if 0.0 == collapse_probability[i_depth][i_shore]:
387                                break
388                            struct_coll_prob.setdefault(coll_prob,[]).append(i)
389                            break
390                    break
391                           
392        return struct_coll_prob
393   
394    def _calc_collapse_structures(self, collapse_probability,
395                                  verbose_csv=False):
396        """
397        Given the collapse probabilities, throw the dice
398        and collapse some houses
399        """
400       
401        self.struct_collapsed = ['']* self.structure_count
402        if verbose_csv:
403            self.collapse_csv_info = ['']* self.structure_count
404        #for a given 'bin', work out how many houses will collapse
405        for probability, house_indexes in collapse_probability.iteritems():
406            collapse_count = round(len(house_indexes) *probability)
407           
408            if verbose_csv:
409                for i in house_indexes:
410                    # This could be sped up I think
411                    self.collapse_csv_info[i] = str(probability) + ' prob.( ' \
412                           + str(int(collapse_count)) + ' collapsed out of ' \
413                           + str(len(house_indexes)) + ')'
414            for _ in range(int(collapse_count)):
415                house_index = choice(house_indexes)
416                self.struct_damage[house_index] = 1.0
417                self.contents_damage[house_index] = 1.0
418                house_indexes.remove(house_index)
419                self.struct_collapsed[house_index] = 1
420               
421            # Warning, the collapse_probability list now lists
422            # houses that did not collapse, (though not all of them)
423           
424#############################################################################
425if __name__ == "__main__":
426    pass 
Note: See TracBrowser for help on using the repository browser.