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

Last change on this file since 4855 was 4855, checked in by duncan, 17 years ago

bug fix

File size: 19.3 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            name, extension = exposure_file_in.split('.')
98            exposure_file_out = name + exposure_file_out_marker + \
99                                '.' + extension
100        csv.save(exposure_file_out)
101        if verbose: print '\n Augmented building file written to %s \n' %exposure_file_out
102   
103def add_depth_and_momentum2csv(sww_base_name, exposure_file_in,
104                      exposure_file_out=None,
105                      overwrite=False, verbose=True,
106                                 use_cache = True):
107    """
108    Calculate the maximum depth and momemtum in an sww file, for locations
109    specified in an csv exposure file.
110   
111    These calculations are done over all the sww files with the sww_base_name
112    in the specified directory.
113    """
114
115    csv = Exposure_csv(exposure_file_in)
116    geospatial = csv.get_location()
117    max_depths, max_momentums = calc_max_depth_and_momentum(sww_base_name,
118                                                          geospatial,
119                                                          verbose=verbose,
120                                                          use_cache=use_cache)
121    csv.set_column("MAX INUNDATION DEPTH (m)",max_depths, overwrite=overwrite)
122    csv.set_column("MOMENTUM (m^2/s) ",max_momentums, overwrite=overwrite)
123    csv.save(exposure_file_out)
124   
125def calc_max_depth_and_momentum(sww_base_name, points,
126                                ground_floor_height=0.0,
127                                verbose=True,
128                                 use_cache = True):
129    """
130    Calculate the maximum inundation height above ground floor for a list
131    of locations.
132
133    The inundation value is in the range -ground_floor_height to
134    overflow errors.
135
136    These calculations are done over all the sww files with the sww_base_name
137    in the specified directory.
138    """
139    quantities =  ['stage', 'elevation', 'xmomentum', 'ymomentum']
140    #print "points",points
141    points = ensure_absolute(points)
142    point_count = len(points)
143
144    # initialise the max lists
145    max_depths = [-ground_floor_height]*point_count
146    max_momentums = [-ground_floor_height]*point_count
147   
148    # How many sww files are there?
149    dir, base = os.path.split(sww_base_name)
150    #print "basename_in",basename_in
151    if base[-4:] == '.sww':
152        base = base[:-4]
153    #print "base",base
154    if dir == "": dir = "." # Unix compatibility
155    dir_ls = os.listdir(dir)
156    interate_over = [x for x in dir_ls if base in x and x[-4:] == '.sww']
157    if len(interate_over) == 0:
158        msg = 'No files of the base name %s.'\
159              %(sww_base_name)
160        raise IOError, msg
161    #print "interate_over", interate_over
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 = "Structure Loss ($)"
197    CONTENTS_LOSS_TITLE = "Contents Loss ($)"
198    CONTENTS_DAMAGE_TITLE = "Contents damaged (fraction)"
199    STRUCT_DAMAGE_TITLE = "Structure damaged (fraction)"
200    COLLAPSE_CSV_INFO_TITLE = "Calculation notes"
201    MAX_DEPTH_TITLE = "Inundation height above ground floor (m)"
202    STRUCT_COLLAPSED_TITLE = "collapsed structure if 1"
203    STRUCT_INUNDATED_TITLE = "inundated structure if 1"
204    double_brick_damage_array = 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        (ravel(double_brick_damage_array[:,0:1]),),
217        ravel(double_brick_damage_array[:,1:]))
218   
219    brick_veeer_damage_array = 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        (ravel(brick_veeer_damage_array[:,0:1]),),
232        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 = 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        (ravel(contents_damage_array[:,0:1]),),
249        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
313        struct_damage = zeros(self.structure_count,Float)
314        contents_damage = zeros(self.structure_count,Float)
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            #print "i",i
323            #print "max_depth",max_depth
324            #print "shore_distance",shore_distance
325            #print "wall",wall
326            ## WARNING SKIP IF DEPTH < 0.0
327            if 0.0 > max_depth:
328                continue
329           
330            # The definition of inundated is if the max_depth is > 0.0
331            self.struct_inundated[i] = 1.0
332           
333            #calc structural damage %
334            damage_curve = self.struct_damage_curve.get(wall,
335                                              self.default_struct_damage_curve)
336            struct_damage[i] = damage_curve(max_depth)
337            contents_damage[i] = self.contents_damage_curve(max_depth)
338           
339        self.struct_damage = struct_damage
340        self.contents_damage = contents_damage
341           
342       
343    def calc_cost(self):
344        """
345        Once the damage has been calculated, determine the $ cost.
346        """
347        # ensure_numeric does not cut it.
348        self.struct_loss = self.struct_damage * \
349                           ensure_numeric(self.struct_costs)
350        self.contents_loss = self.contents_damage * \
351                           ensure_numeric(self.content_costs)
352       
353    def calc_collapse_probability(self):
354        """
355        return a dict of which structures have x probability of collapse.
356             key is collapse probability
357             value is list of struct indexes with key probability of collapse
358        """
359        # I could've done this is the calc_damage_percentages and
360        # Just had one loop.
361        # But for ease of testing and bug finding I'm seperating the loops.
362        # I'm make the outer loop for both of them the same though,
363        # so this loop can easily be folded into the other loop.
364       
365        # dict of which structures have x probability of collapse.
366        # key of collapse probability
367        # value of list of struct indexes
368        struct_coll_prob = {}
369       
370        for i,max_depth,shore_distance,wall in map(None,
371                                                   range(self.structure_count),
372                                                   self.max_depths,
373                                                   self.shore_distances,
374                                                   self.walls):
375            #print "i",i
376            #print "max_depth",max_depth
377            #print "shore_distance",shore_distance
378            #print "wall",wall
379            # WARNING ASSUMING THE FIRST BIN OF DEPTHS GIVE A ZERO PROBABILITY
380            depth_upper_limits = self.depth_upper_limits
381            shore_upper_limits = self.shore_upper_limits
382            collapse_probability = self.collapse_probability
383            if max_depth <= depth_upper_limits[0]:
384                continue
385            start = 1
386            for i_depth, depth_limit in enumerate(depth_upper_limits[start:]):
387                #Have to change i_depth so it indexes into the lists correctly
388                i_depth += start
389                if max_depth <= depth_limit:
390                    for i_shore, shore_limit in enumerate(shore_upper_limits):
391                        if shore_distance <= shore_limit:
392                            coll_prob = collapse_probability[i_depth][i_shore]
393                            if 0.0 == collapse_probability[i_depth][i_shore]:
394                                break
395                            struct_coll_prob.setdefault(coll_prob,[]).append(i)
396                            break
397                    break
398                           
399        return struct_coll_prob
400   
401    def _calc_collapse_structures(self, collapse_probability,
402                                  verbose_csv=False):
403        """
404        Given the collapse probabilities, throw the dice
405        and collapse some houses
406        """
407       
408        self.struct_collapsed = ['']* self.structure_count
409        if verbose_csv:
410            self.collapse_csv_info = ['']* self.structure_count
411        #for a given 'bin', work out how many houses will collapse
412        for probability, house_indexes in collapse_probability.iteritems():
413            collapse_count = round(len(house_indexes) *probability)
414           
415            if verbose_csv:
416                for i in house_indexes:
417                    # This could be sped up I think
418                    self.collapse_csv_info[i] = str(probability) + ' prob.( ' \
419                           + str(int(collapse_count)) + ' collapsed out of ' \
420                           + str(len(house_indexes)) + ')'
421            for _ in range(int(collapse_count)):
422                house_index = choice(house_indexes)
423                self.struct_damage[house_index] = 1.0
424                self.contents_damage[house_index] = 1.0
425                house_indexes.remove(house_index)
426                self.struct_collapsed[house_index] = 1
427               
428            # Warning, the collapse_probability list now lists
429            # houses that did not collapse, (though not all of them)
430            #print "",self.collapse_csv_info
431           
432#############################################################################
433if __name__ == "__main__":
434    pass 
Note: See TracBrowser for help on using the repository browser.