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

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

Getting inundation_damage to work with scipy

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