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

Last change on this file since 8249 was 8249, checked in by steve, 11 years ago

Got rid of those annoying double_scalar warnings in the windows code (just divide by zero warnings)
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
10
11import numpy as num
12
13
14try: 
15    import kinds 
16except ImportError: 
17    # Hand-built mockup of the things we need from the kinds package, since it
18    # was recently removed from the standard numeric distro.  Some users may 
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
32from exposure import Exposure
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
36import anuga.utilities.log as log
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    """
69    if isinstance(exposure_files_in, basestring):
70        exposure_files_in = [exposure_files_in]
71
72
73    for exposure_file_in in exposure_files_in:
74        csv = Exposure(exposure_file_in,
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)
104        if verbose: log.critical('Augmented building file written to %s'
105                                 % exposure_file_out)
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
119    csv = Exposure(exposure_file_in)
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
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
165    for this_sww_file in interate_over:
166        callable_sww = file_function(dir+sep+this_sww_file,
167                                     quantities=quantities,
168                                     interpolation_points=points,
169                                     verbose=verbose,
170                                     use_cache=use_cache)
171
172        for point_i, point in enumerate(points):
173            for time in callable_sww.get_time():
174                quantity_values = callable_sww(time,point_i)
175                w = quantity_values[0]
176                z = quantity_values[1]
177                uh = quantity_values[2] 
178                vh = quantity_values[3]
179
180                #print w,z,uh,vh
181                if w == NAN or z == NAN or uh == NAN or vh == NAN:
182                    continue
183                   
184                #  -ground_floor_height is the minimum value.
185                depth = w - z - ground_floor_height
186             
187                if depth > max_depths[point_i]:
188                    max_depths[point_i] = depth
189               
190                momentum = sqrt(uh*uh + vh*vh)
191                if momentum > max_momentums[point_i]:
192                    max_momentums[point_i] = momentum
193
194
195    return max_depths, max_momentums
196
197class EventDamageModel:
198    """
199    Object for working out the damage and cost
200
201    """
202    STRUCT_LOSS_TITLE = "STRUCT_LOSS_$"#"Structure Loss ($)"
203    CONTENTS_LOSS_TITLE = "CONTENTS_LOSS_$"#"Contents Loss ($)"
204    CONTENTS_DAMAGE_TITLE = "CONTENTS_DAMAGE_fraction"#"Contents damaged (fraction)"
205    STRUCT_DAMAGE_TITLE = "STRUCT_DAMAGE_fraction" #"Structure damaged (fraction)"
206    COLLAPSE_CSV_INFO_TITLE = "COLLAPSE_CSV_INFO"#"Calculation notes"
207    MAX_DEPTH_TITLE = "MAX_DEPTH_m" #"Inundation height above ground floor (m)"
208    STRUCT_COLLAPSED_TITLE = "STRUCT_COLLAPSED"#"collapsed structure if 1"
209    STRUCT_INUNDATED_TITLE = "STRUCT_INUNDATED"#"inundated structure if 1"
210    double_brick_damage_array = num.array([[-kinds.default_float_kind.MAX, 0.0],
211                                           [0.0-depth_epsilon, 0.0],
212                                           [0.0,0.016],
213                                           [0.1,0.150],
214                                           [0.3,0.425],
215                                           [0.5,0.449],
216                                           [1.0,0.572],
217                                           [1.5,0.582],
218                                           [2.0,0.587],
219                                           [2.5,0.647],
220                                           [kinds.default_float_kind.MAX,64.7]])
221    double_brick_damage_curve = InterpolatingFunction( \
222        (num.ravel(double_brick_damage_array[:,0:1]),),
223        num.ravel(double_brick_damage_array[:,1:]))
224   
225    brick_veeer_damage_array = num.array([[-kinds.default_float_kind.MAX, 0.0],
226                                          [0.0-depth_epsilon, 0.0],
227                                          [0.0,0.016],
228                                          [0.1,0.169],
229                                          [0.3,0.445],
230                                          [0.5,0.472],
231                                          [1.0,0.618],
232                                          [1.5,0.629],
233                                          [2.0,0.633],
234                                          [2.5,0.694],
235                                          [kinds.default_float_kind.MAX,69.4]])
236    brick_veeer_damage_curve = InterpolatingFunction( \
237        (num.ravel(brick_veeer_damage_array[:,0:1]),),
238        num.ravel(brick_veeer_damage_array[:,1:]))
239    struct_damage_curve = {'Double Brick':double_brick_damage_curve,
240                           'Brick Veneer':brick_veeer_damage_curve}
241    default_struct_damage_curve = brick_veeer_damage_curve
242
243    contents_damage_array = num.array([[-kinds.default_float_kind.MAX, 0.0],
244                                       [0.0-depth_epsilon, 0.0],
245                                       [0.0,0.013],
246                                       [0.1,0.102],
247                                       [0.3,0.381],
248                                       [0.5,0.500],
249                                       [1.0,0.970],
250                                       [1.5,0.976],
251                                       [2.0,0.986],
252                                       [kinds.default_float_kind.MAX,98.6]])
253    contents_damage_curve = InterpolatingFunction( \
254        (num.ravel(contents_damage_array[:,0:1]),),
255        num.ravel(contents_damage_array[:,1:]))
256
257    #building collapse probability
258    # inundation depth above ground floor, m
259    depth_upper_limits = [depth_epsilon, 1.0, 2.0, 3.0, 5.0,
260                          kinds.default_float_kind.MAX]
261    # shore mistance, m
262    shore_upper_limits = [125,200,250, kinds.default_float_kind.MAX]
263    # Building collapse probability
264    collapse_probability = [[0.0, 0.0, 0.0, 0.0], #Code below assumes 0.0
265                            [0.05, 0.02, 0.01, 0.0],
266                            [0.6, 0.3, 0.1, 0.05],
267                            [0.8, 0.4, 0.25, 0.15],
268                            [0.95, 0.7, 0.5, 0.3],
269                            [0.99, 0.9, 0.65, 0.45]]
270   
271    def __init__(self,max_depths, shore_distances, walls,
272                 struct_costs, content_costs):
273        """
274        max depth is Inundation height above ground floor (m), so
275                  the ground floor has been taken into account.
276        """
277        self.max_depths = [float(x) for x in max_depths]
278        self.shore_distances = [float(x) for x in shore_distances]
279        self.walls = walls
280        self.struct_costs = [float(x) for x in struct_costs]
281        self.content_costs = [float(x) for x in content_costs]
282
283        self.structure_count = len(self.max_depths)
284        #Fixme expand
285        assert self.structure_count == len(self.shore_distances)
286        assert  self.structure_count == len(self.walls)
287        assert  self.structure_count == len(self.struct_costs)
288        assert  self.structure_count == len(self.content_costs)
289        #assert  self.structure_count == len(self.)
290       
291    def calc_damage_and_costs(self, verbose_csv=False, verbose=False):
292        """
293        This is an overall method to calculate the % damage and collapsed
294        structures and then the $ loss.
295        """
296        self.calc_damage_percentages()
297        collapse_probability = self.calc_collapse_probability()
298        self._calc_collapse_structures(collapse_probability,
299                                  verbose_csv=verbose_csv)
300        self.calc_cost()
301        results_dict = {self.STRUCT_LOSS_TITLE:self.struct_loss
302                        ,self.STRUCT_DAMAGE_TITLE:self.struct_damage
303                        ,self.CONTENTS_LOSS_TITLE:self.contents_loss
304                        ,self.CONTENTS_DAMAGE_TITLE:self.contents_damage
305                        ,self.MAX_DEPTH_TITLE:self.max_depths
306                        ,self.STRUCT_COLLAPSED_TITLE:self.struct_collapsed
307                        ,self.STRUCT_INUNDATED_TITLE:self.struct_inundated
308                        }
309        if verbose_csv:
310           results_dict[self.COLLAPSE_CSV_INFO_TITLE] = self.collapse_csv_info
311        return results_dict
312           
313    def calc_damage_percentages(self):
314        """
315        Using stage curves calc the damage to structures and contents
316        """
317
318        # the data being created
319        struct_damage = num.zeros(self.structure_count, num.float)
320        contents_damage = num.zeros(self.structure_count, num.float)
321        self.struct_inundated = ['']* self.structure_count
322
323        for i,max_depth,shore_distance,wall in map(None,
324                                                   range(self.structure_count),
325                                                   self.max_depths,
326                                                   self.shore_distances,
327                                                   self.walls):
328            ## WARNING SKIP IF DEPTH < 0.0
329            if 0.0 > max_depth:
330                continue
331           
332            # The definition of inundated is if the max_depth is > 0.0
333            self.struct_inundated[i] = 1.0
334           
335            #calc structural damage %
336            damage_curve = self.struct_damage_curve.get(wall,
337                                              self.default_struct_damage_curve)
338            struct_damage[i] = damage_curve(max_depth)
339            contents_damage[i] = self.contents_damage_curve(max_depth)
340           
341        self.struct_damage = struct_damage
342        self.contents_damage = contents_damage
343           
344       
345    def calc_cost(self):
346        """
347        Once the damage has been calculated, determine the $ cost.
348        """
349        # ensure_numeric does not cut it.
350        self.struct_loss = self.struct_damage * \
351                           ensure_numeric(self.struct_costs)
352        self.contents_loss = self.contents_damage * \
353                           ensure_numeric(self.content_costs)
354       
355    def calc_collapse_probability(self):
356        """
357        return a dict of which structures have x probability of collapse.
358             key is collapse probability
359             value is list of struct indexes with key probability of collapse
360        """
361        # I could've done this is the calc_damage_percentages and
362        # Just had one loop.
363        # But for ease of testing and bug finding I'm seperating the loops.
364        # I'm make the outer loop for both of them the same though,
365        # so this loop can easily be folded into the other loop.
366       
367        # dict of which structures have x probability of collapse.
368        # key of collapse probability
369        # value of list of struct indexes
370        struct_coll_prob = {}
371       
372        for i,max_depth,shore_distance,wall in map(None,
373                                                   range(self.structure_count),
374                                                   self.max_depths,
375                                                   self.shore_distances,
376                                                   self.walls):
377            # WARNING ASSUMING THE FIRST BIN OF DEPTHS GIVE A ZERO PROBABILITY
378            depth_upper_limits = self.depth_upper_limits
379            shore_upper_limits = self.shore_upper_limits
380            collapse_probability = self.collapse_probability
381            if max_depth <= depth_upper_limits[0]:
382                continue
383            start = 1
384            for i_depth, depth_limit in enumerate(depth_upper_limits[start:]):
385                #Have to change i_depth so it indexes into the lists correctly
386                i_depth += start
387                if max_depth <= depth_limit:
388                    for i_shore, shore_limit in enumerate(shore_upper_limits):
389                        if shore_distance <= shore_limit:
390                            coll_prob = collapse_probability[i_depth][i_shore]
391                            if 0.0 == collapse_probability[i_depth][i_shore]:
392                                break
393                            struct_coll_prob.setdefault(coll_prob,[]).append(i)
394                            break
395                    break
396                           
397        return struct_coll_prob
398   
399    def _calc_collapse_structures(self, collapse_probability,
400                                  verbose_csv=False):
401        """
402        Given the collapse probabilities, throw the dice
403        and collapse some houses
404        """
405       
406        self.struct_collapsed = ['']* self.structure_count
407        if verbose_csv:
408            self.collapse_csv_info = ['']* self.structure_count
409        #for a given 'bin', work out how many houses will collapse
410        for probability, house_indexes in collapse_probability.iteritems():
411            collapse_count = round(len(house_indexes) *probability)
412           
413            if verbose_csv:
414                for i in house_indexes:
415                    # This could be sped up I think
416                    self.collapse_csv_info[i] = str(probability) + ' prob.( ' \
417                           + str(int(collapse_count)) + ' collapsed out of ' \
418                           + str(len(house_indexes)) + ')'
419            for _ in range(int(collapse_count)):
420                house_index = choice(house_indexes)
421                self.struct_damage[house_index] = 1.0
422                self.contents_damage[house_index] = 1.0
423                house_indexes.remove(house_index)
424                self.struct_collapsed[house_index] = 1
425               
426            # Warning, the collapse_probability list now lists
427            # houses that did not collapse, (though not all of them)
428           
429#############################################################################
430if __name__ == "__main__":
431    pass 
Note: See TracBrowser for help on using the repository browser.