source: branches/numpy/anuga/damage_modelling/inundation_damage.py @ 9053

Last change on this file since 9053 was 6304, checked in by rwilson, 16 years ago

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