source: anuga_core/source/anuga/damage_modelling/inundation_damage.py @ 6034

Last change on this file since 6034 was 6034, checked in by kristy, 16 years ago

Change titles so that the csv can go straight into Arc
File size: 19.5 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 Numeric import array, ravel, Float, zeros
10from random import choice
11from types import StringType
12
13try: 
14    import kinds 
15except ImportError: 
16    # Hand-built mockup of the things we need from the kinds package, since it
17    # was recently removed from the standard Numeric distro.  Some users may 
18    # not have it by default. 
19    class _bunch: 
20        pass 
21         
22    class _kinds(_bunch): 
23        default_float_kind = _bunch() 
24        default_float_kind.MIN = 2.2250738585072014e-308  #smallest +ve number
25        default_float_kind.MAX = 1.7976931348623157e+308 
26     
27    kinds = _kinds()
28   
29
30from anuga.utilities.numerical_tools import ensure_numeric
31from anuga.shallow_water.data_manager import Exposure_csv
32from anuga.abstract_2d_finite_volumes.util import file_function
33from anuga.geospatial_data.geospatial_data import ensure_absolute
34from anuga.utilities.numerical_tools import NAN
35from config import epsilon
36depth_epsilon = epsilon
37
38# Change these if the ouput from nexix changes
39SHORE_DIST_LABEL = 'SHORE_DIST'
40WALL_TYPE_LABEL = 'WALL_TYPE'
41STR_VALUE_LABEL = 'STR_VALUE'
42CONT_VALUE_LABEL = 'CONT_VALUE'
43
44def inundation_damage(sww_base_name, exposure_files_in,
45                      exposure_file_out_marker=None,
46                      ground_floor_height=0.3,
47                      overwrite=False, verbose=True,
48                                 use_cache = True):
49    """
50    This is the main function for calculating tsunami damage due to
51    inundation.  It gets the location of structures from the exposure
52    file and gets the inundation of these structures from the
53    sww file.
54
55    It then calculates the damage loss.
56
57    Note, structures outside of the sww file get the minimum inundation
58    (-ground_floor_height).
59   
60    These calculations are done over all the sww files with the sww_base_name
61    in the specified directory.
62
63    exposure_files_in - a file or a list of files to input from
64    exposure_file_out_marker -  this string will be added to the input file
65                                name to get the output file name
66    """
67    if type(exposure_files_in) == StringType:
68        exposure_files_in = [exposure_files_in]
69
70
71    for exposure_file_in in exposure_files_in:
72        csv = Exposure_csv(exposure_file_in,
73                           title_check_list=[SHORE_DIST_LABEL,WALL_TYPE_LABEL,
74                                             STR_VALUE_LABEL,CONT_VALUE_LABEL])
75        geospatial = csv.get_location()
76        geospatial = ensure_absolute(geospatial)
77        max_depths, max_momentums = calc_max_depth_and_momentum(sww_base_name,
78                        geospatial,
79                        ground_floor_height=ground_floor_height,
80                        verbose=verbose,
81                        use_cache=use_cache)
82        edm = EventDamageModel(max_depths,
83                               csv.get_column(SHORE_DIST_LABEL),
84                               csv.get_column(WALL_TYPE_LABEL),
85                               csv.get_column(STR_VALUE_LABEL),
86                               csv.get_column(CONT_VALUE_LABEL)
87                               )
88        results_dic = edm.calc_damage_and_costs(verbose_csv=True,
89                                                verbose=verbose)
90        for title, value in results_dic.iteritems():
91            csv.set_column(title, value, overwrite=overwrite)
92   
93        # Save info back to csv file
94        if exposure_file_out_marker == None:
95            exposure_file_out = exposure_file_in
96        else:
97            # split off extension, in such a way to deal with more than one '.' in the name of file
98            split_name = exposure_file_in.split('.')
99            exposure_file_out =  '.'.join(split_name[:-1]) + exposure_file_out_marker + \
100                                '.' + split_name[-1]
101        csv.save(exposure_file_out)
102        if verbose: print '\n Augmented building file written to %s \n' %exposure_file_out
103   
104def add_depth_and_momentum2csv(sww_base_name, exposure_file_in,
105                      exposure_file_out=None,
106                      overwrite=False, verbose=True,
107                                 use_cache = True):
108    """
109    Calculate the maximum depth and momemtum in an sww file, for locations
110    specified in an csv exposure file.
111   
112    These calculations are done over all the sww files with the sww_base_name
113    in the specified directory.
114    """
115
116    csv = Exposure_csv(exposure_file_in)
117    geospatial = csv.get_location()
118    max_depths, max_momentums = calc_max_depth_and_momentum(sww_base_name,
119                                                          geospatial,
120                                                          verbose=verbose,
121                                                          use_cache=use_cache)
122    csv.set_column("MAX INUNDATION DEPTH (m)",max_depths, overwrite=overwrite)
123    csv.set_column("MOMENTUM (m^2/s) ",max_momentums, overwrite=overwrite)
124    csv.save(exposure_file_out)
125   
126def calc_max_depth_and_momentum(sww_base_name, points,
127                                ground_floor_height=0.0,
128                                verbose=True,
129                                 use_cache = True):
130    """
131    Calculate the maximum inundation height above ground floor for a list
132    of locations.
133
134    The inundation value is in the range -ground_floor_height to
135    overflow errors.
136
137    These calculations are done over all the sww files with the sww_base_name
138    in the specified directory.
139    """
140    quantities =  ['stage', 'elevation', 'xmomentum', 'ymomentum']
141    #print "points",points
142    points = ensure_absolute(points)
143    point_count = len(points)
144
145    # initialise the max lists
146    max_depths = [-ground_floor_height]*point_count
147    max_momentums = [-ground_floor_height]*point_count
148   
149    # How many sww files are there?
150    dir, base = os.path.split(sww_base_name)
151    #print "basename_in",basename_in
152    if base[-4:] == '.sww':
153        base = base[:-4]
154    #print "base",base
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    #print "interate_over", interate_over
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 = 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        (ravel(double_brick_damage_array[:,0:1]),),
218        ravel(double_brick_damage_array[:,1:]))
219   
220    brick_veeer_damage_array = 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        (ravel(brick_veeer_damage_array[:,0:1]),),
233        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 = 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        (ravel(contents_damage_array[:,0:1]),),
250        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 = zeros(self.structure_count,Float)
315        contents_damage = zeros(self.structure_count,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            #print "i",i
324            #print "max_depth",max_depth
325            #print "shore_distance",shore_distance
326            #print "wall",wall
327            ## WARNING SKIP IF DEPTH < 0.0
328            if 0.0 > max_depth:
329                continue
330           
331            # The definition of inundated is if the max_depth is > 0.0
332            self.struct_inundated[i] = 1.0
333           
334            #calc structural damage %
335            damage_curve = self.struct_damage_curve.get(wall,
336                                              self.default_struct_damage_curve)
337            struct_damage[i] = damage_curve(max_depth)
338            contents_damage[i] = self.contents_damage_curve(max_depth)
339           
340        self.struct_damage = struct_damage
341        self.contents_damage = contents_damage
342           
343       
344    def calc_cost(self):
345        """
346        Once the damage has been calculated, determine the $ cost.
347        """
348        # ensure_numeric does not cut it.
349        self.struct_loss = self.struct_damage * \
350                           ensure_numeric(self.struct_costs)
351        self.contents_loss = self.contents_damage * \
352                           ensure_numeric(self.content_costs)
353       
354    def calc_collapse_probability(self):
355        """
356        return a dict of which structures have x probability of collapse.
357             key is collapse probability
358             value is list of struct indexes with key probability of collapse
359        """
360        # I could've done this is the calc_damage_percentages and
361        # Just had one loop.
362        # But for ease of testing and bug finding I'm seperating the loops.
363        # I'm make the outer loop for both of them the same though,
364        # so this loop can easily be folded into the other loop.
365       
366        # dict of which structures have x probability of collapse.
367        # key of collapse probability
368        # value of list of struct indexes
369        struct_coll_prob = {}
370       
371        for i,max_depth,shore_distance,wall in map(None,
372                                                   range(self.structure_count),
373                                                   self.max_depths,
374                                                   self.shore_distances,
375                                                   self.walls):
376            #print "i",i
377            #print "max_depth",max_depth
378            #print "shore_distance",shore_distance
379            #print "wall",wall
380            # WARNING ASSUMING THE FIRST BIN OF DEPTHS GIVE A ZERO PROBABILITY
381            depth_upper_limits = self.depth_upper_limits
382            shore_upper_limits = self.shore_upper_limits
383            collapse_probability = self.collapse_probability
384            if max_depth <= depth_upper_limits[0]:
385                continue
386            start = 1
387            for i_depth, depth_limit in enumerate(depth_upper_limits[start:]):
388                #Have to change i_depth so it indexes into the lists correctly
389                i_depth += start
390                if max_depth <= depth_limit:
391                    for i_shore, shore_limit in enumerate(shore_upper_limits):
392                        if shore_distance <= shore_limit:
393                            coll_prob = collapse_probability[i_depth][i_shore]
394                            if 0.0 == collapse_probability[i_depth][i_shore]:
395                                break
396                            struct_coll_prob.setdefault(coll_prob,[]).append(i)
397                            break
398                    break
399                           
400        return struct_coll_prob
401   
402    def _calc_collapse_structures(self, collapse_probability,
403                                  verbose_csv=False):
404        """
405        Given the collapse probabilities, throw the dice
406        and collapse some houses
407        """
408       
409        self.struct_collapsed = ['']* self.structure_count
410        if verbose_csv:
411            self.collapse_csv_info = ['']* self.structure_count
412        #for a given 'bin', work out how many houses will collapse
413        for probability, house_indexes in collapse_probability.iteritems():
414            collapse_count = round(len(house_indexes) *probability)
415           
416            if verbose_csv:
417                for i in house_indexes:
418                    # This could be sped up I think
419                    self.collapse_csv_info[i] = str(probability) + ' prob.( ' \
420                           + str(int(collapse_count)) + ' collapsed out of ' \
421                           + str(len(house_indexes)) + ')'
422            for _ in range(int(collapse_count)):
423                house_index = choice(house_indexes)
424                self.struct_damage[house_index] = 1.0
425                self.contents_damage[house_index] = 1.0
426                house_indexes.remove(house_index)
427                self.struct_collapsed[house_index] = 1
428               
429            # Warning, the collapse_probability list now lists
430            # houses that did not collapse, (though not all of them)
431            #print "",self.collapse_csv_info
432           
433#############################################################################
434if __name__ == "__main__":
435    pass 
Note: See TracBrowser for help on using the repository browser.