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

Last change on this file since 8125 was 8125, checked in by wilsonr, 13 years ago

Changes to address ticket 360.

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