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

Last change on this file since 4173 was 3563, checked in by ole, 18 years ago

Moved shallow water out from the old pyvolution directory.
All tests pass, most examples and okushiri works again.
File size: 17.1 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"""
6
7from math import sqrt
8from Scientific.Functions.Interpolation import InterpolatingFunction
9from Numeric import array, ravel, Float, zeros
10from random import choice
11
12try: 
13    import kinds 
14except ImportError: 
15    # Hand-built mockup of the things we need from the kinds package, since it
16    # was recently removed from the standard Numeric distro.  Some users may 
17    # not have it by default. 
18    class _bunch: 
19        pass 
20         
21    class _kinds(_bunch): 
22        default_float_kind = _bunch() 
23        default_float_kind.MIN = 2.2250738585072014e-308  #smallest +ve number
24        default_float_kind.MAX = 1.7976931348623157e+308 
25     
26    kinds = _kinds()
27   
28
29from anuga.utilities.numerical_tools import ensure_numeric
30from anuga.shallow_water.data_manager import Exposure_csv
31from anuga.abstract_2d_finite_volumes.util import file_function
32from anuga.geospatial_data.geospatial_data import ensure_absolute
33from anuga.utilities.numerical_tools import NAN
34from anuga_config import epsilon
35depth_epsilon = epsilon
36
37def inundation_damage(sww_file, exposure_file_in,
38                      exposure_file_out=None,
39                      ground_floor_height=0.3,
40                      overwrite=False, verbose=True,
41                                 use_cache = True):
42    """
43    This is the main function for calculating tsunami damage due to
44    inundation.  It gets the location of structures from the exposure
45    file and gets the inundation of these structures from the
46    sww file.
47
48    It then calculates the damage loss.
49
50    Note, structures outside of the sww file get the minimum inundation
51    (-ground_floor_height).
52    """
53
54    csv = Exposure_csv(exposure_file_in)
55    geospatial = csv.get_location()
56    geospatial = ensure_absolute(geospatial)
57    max_depths, max_momentums = calc_max_depth_and_momentum(sww_file,
58                                    geospatial,
59                                    ground_floor_height=ground_floor_height,
60                                    verbose=verbose,
61                                    use_cache=use_cache)
62    edm = EventDamageModel(max_depths,
63                           csv.get_column('SHORE_DIST'),
64                           csv.get_column('WALLS'),
65                           csv.get_column('STR_VALUE'),
66                           csv.get_column('C_VALUE')
67                           )
68    results_dic = edm.calc_damage_and_costs(verbose_csv=True, verbose=verbose)
69    for title, value in results_dic.iteritems():
70        csv.set_column(title, value, overwrite=overwrite)
71   
72    # Save info back to csv file
73    if exposure_file_out == None:
74        exposure_file_out = exposure_file_in
75    csv.save(exposure_file_out)
76   
77def add_depth_and_momentum2csv(sww_file, exposure_file_in,
78                      exposure_file_out=None,
79                      overwrite=False, verbose=True,
80                                 use_cache = True):
81    """
82    Calculate the maximum depth and momemtum in an sww file, for locations
83    specified in an csv exposure file.
84    """
85
86    csv = Exposure_csv(exposure_file_in)
87    geospatial = csv.get_location()
88    max_depths, max_momentums = calc_max_depth_and_momentum(sww_file,
89                                                          geospatial,
90                                                          verbose=verbose,
91                                                          use_cache=use_cache)
92    csv.set_column("MAX INUNDATION DEPTH (m)",max_depths, overwrite=overwrite)
93    csv.set_column("MOMENTUM (m^2/s) ",max_momentums, overwrite=overwrite)
94    csv.save(exposure_file_out)
95   
96def calc_max_depth_and_momentum(sww_file, points,
97                                ground_floor_height=0.0,
98                                verbose=True,
99                                 use_cache = True):
100    """
101    Calculate the maximum inundation height above ground floor for a list
102    of locations.
103
104    The inundation value is in the range -ground_floor_height to
105    overflow errors.
106    """
107    quantities =  ['stage', 'elevation', 'xmomentum', 'ymomentum']
108    #print "points",points
109    points = ensure_absolute(points)
110    point_count = len(points)
111    callable_sww = file_function(sww_file,
112                                 quantities=quantities,
113                                 interpolation_points=points,
114                                 verbose=verbose,
115                                 use_cache=use_cache)
116    # initialise the max lists
117    max_depths = [-ground_floor_height]*point_count
118    max_momentums = [-ground_floor_height]*point_count
119    for point_i, point in enumerate(points):
120        for time in callable_sww.get_time():
121            quantities = callable_sww(time,point_i)
122            #print "quantities", quantities
123           
124            w = quantities[0]
125            z = quantities[1]
126            uh = quantities[2]
127            vh = quantities[3]
128
129            #  -ground_floor_height is the minimum value.
130            depth = w - z - ground_floor_height
131             
132            if depth > max_depths[point_i]:
133                max_depths[point_i] = depth
134            if w == NAN or z == NAN or uh == NAN or vh == NAN:
135                continue
136            momentum = sqrt(uh*uh + vh*vh)
137            if momentum > max_momentums[point_i]:
138                max_momentums[point_i] = momentum
139    return max_depths, max_momentums
140
141class EventDamageModel:
142    """
143    Object for working out the damage and cost
144
145    """
146    STRUCT_LOSS_TITLE = "Structure Loss ($)"
147    CONTENTS_LOSS_TITLE = "Contents Loss ($)"
148    CONTENTS_DAMAGE_TITLE = "Contents damaged (fraction)"
149    STRUCT_DAMAGE_TITLE = "Structure damaged (fraction)"
150    COLLAPSE_CSV_INFO_TITLE = "Calculation notes"
151    MAX_DEPTH_TITLE = "Inundation height above ground floor (m)"
152    STRUCT_COLLAPSED_TITLE = "collapsed structure if 1"
153    STRUCT_INUNDATED_TITLE = "inundated structure if 1"
154    double_brick_damage_array = array([[-kinds.default_float_kind.MAX, 0.0],
155                                       [0.0-depth_epsilon, 0.0],
156                                       [0.0,0.016],
157                                       [0.1,0.150],
158                                       [0.3,0.425],
159                                       [0.5,0.449],
160                                       [1.0,0.572],
161                                       [1.5,0.582],
162                                       [2.0,0.587],
163                                       [2.5,0.647],
164                                       [kinds.default_float_kind.MAX,64.7]])
165    double_brick_damage_curve = InterpolatingFunction( \
166        (ravel(double_brick_damage_array[:,0:1]),),
167        ravel(double_brick_damage_array[:,1:]))
168   
169    brick_veeer_damage_array = array([[-kinds.default_float_kind.MAX, 0.0],
170                                       [0.0-depth_epsilon, 0.0],
171                                       [0.0,0.016],
172                                       [0.1,0.169],
173                                       [0.3,0.445],
174                                       [0.5,0.472],
175                                       [1.0,0.618],
176                                       [1.5,0.629],
177                                       [2.0,0.633],
178                                       [2.5,0.694],
179                                       [kinds.default_float_kind.MAX,69.4]])
180    brick_veeer_damage_curve = InterpolatingFunction( \
181        (ravel(brick_veeer_damage_array[:,0:1]),),
182        ravel(brick_veeer_damage_array[:,1:]))
183    struct_damage_curve = {'Double Brick':double_brick_damage_curve,
184                           'Brick Veneer':brick_veeer_damage_curve}
185    default_struct_damage_curve = brick_veeer_damage_curve
186
187    contents_damage_array = array([[-kinds.default_float_kind.MAX, 0.0],
188                                       [0.0-depth_epsilon, 0.0],
189                                       [0.0,0.013],
190                                       [0.1,0.102],
191                                       [0.3,0.381],
192                                       [0.5,0.500],
193                                       [1.0,0.970],
194                                       [1.5,0.976],
195                                       [2.0,0.986],
196                                       [kinds.default_float_kind.MAX,98.6]])
197    contents_damage_curve = InterpolatingFunction( \
198        (ravel(contents_damage_array[:,0:1]),),
199        ravel(contents_damage_array[:,1:]))
200
201    #building collapse probability
202    # inundation depth above ground floor, m
203    depth_upper_limits = [depth_epsilon, 1.0, 2.0, 3.0, 5.0,
204                          kinds.default_float_kind.MAX]
205    # shore mistance, m
206    shore_upper_limits = [125,200,250, kinds.default_float_kind.MAX]
207    # Building collapse probability
208    collapse_probability = [[0.0, 0.0, 0.0, 0.0], #Code below assumes 0.0
209                            [0.05, 0.02, 0.01, 0.0],
210                            [0.6, 0.3, 0.1, 0.05],
211                            [0.8, 0.4, 0.25, 0.15],
212                            [0.95, 0.7, 0.5, 0.3],
213                            [0.99, 0.9, 0.65, 0.45]]
214   
215    def __init__(self,max_depths, shore_distances, walls,
216                 struct_costs, content_costs):
217        """
218        max depth is Inundation height above ground floor (m), so
219                  the ground floor has been taken into account.
220        """
221        self.max_depths = [float(x) for x in max_depths]
222        self.shore_distances = [float(x) for x in shore_distances]
223        self.walls = walls
224        self.struct_costs = [float(x) for x in struct_costs]
225        self.content_costs = [float(x) for x in content_costs]
226
227        self.structure_count = len(self.max_depths)
228        #Fixme expand
229        assert self.structure_count == len(self.shore_distances)
230        assert  self.structure_count == len(self.walls)
231        assert  self.structure_count == len(self.struct_costs)
232        assert  self.structure_count == len(self.content_costs)
233        #assert  self.structure_count == len(self.)
234       
235    def calc_damage_and_costs(self, verbose_csv=False, verbose=False):
236        """
237        This is an overall method to calculate the % damage and collapsed
238        structures and then the $ loss.
239        """
240        self.calc_damage_percentages()
241        collapse_probability = self.calc_collapse_probability()
242        self._calc_collapse_structures(collapse_probability,
243                                  verbose_csv=verbose_csv)
244        self.calc_cost()
245        results_dict = {self.STRUCT_LOSS_TITLE:self.struct_loss
246                        ,self.STRUCT_DAMAGE_TITLE:self.struct_damage
247                        ,self.CONTENTS_LOSS_TITLE:self.contents_loss
248                        ,self.CONTENTS_DAMAGE_TITLE:self.contents_damage
249                        ,self.MAX_DEPTH_TITLE:self.max_depths
250                        ,self.STRUCT_COLLAPSED_TITLE:self.struct_collapsed
251                        ,self.STRUCT_INUNDATED_TITLE:self.struct_inundated
252                        }
253        if verbose_csv:
254           results_dict[self.COLLAPSE_CSV_INFO_TITLE] = self.collapse_csv_info
255        return results_dict
256           
257    def calc_damage_percentages(self):
258        """
259        Using stage curves calc the damage to structures and contents
260        """
261
262        # the data being created
263        struct_damage = zeros(self.structure_count,Float)
264        contents_damage = zeros(self.structure_count,Float)
265        self.struct_inundated = ['']* self.structure_count
266
267        for i,max_depth,shore_distance,wall in map(None,
268                                                   range(self.structure_count),
269                                                   self.max_depths,
270                                                   self.shore_distances,
271                                                   self.walls):
272            #print "i",i
273            #print "max_depth",max_depth
274            #print "shore_distance",shore_distance
275            #print "wall",wall
276            ## WARNING SKIP IF DEPTH < 0.0
277            if 0.0 > max_depth:
278                continue
279           
280            # The definition of inundated is if the max_depth is > 0.0
281            self.struct_inundated[i] = 1.0
282           
283            #calc structural damage %
284            damage_curve = self.struct_damage_curve.get(wall,
285                                              self.default_struct_damage_curve)
286            struct_damage[i] = damage_curve(max_depth)
287            contents_damage[i] = self.contents_damage_curve(max_depth)
288           
289        self.struct_damage = struct_damage
290        self.contents_damage = contents_damage
291           
292       
293    def calc_cost(self):
294        """
295        Once the damage has been calculated, determine the $ cost.
296        """
297        # ensure_numeric does not cut it.
298        self.struct_loss = self.struct_damage * \
299                           ensure_numeric(self.struct_costs)
300        self.contents_loss = self.contents_damage * \
301                           ensure_numeric(self.content_costs)
302       
303    def calc_collapse_probability(self):
304        """
305        return a dict of which structures have x probability of collapse.
306             key is collapse probability
307             value is list of struct indexes with key probability of collapse
308        """
309        # I could've done this is the calc_damage_percentages and
310        # Just had one loop.
311        # But for ease of testing and bug finding I'm seperating the loops.
312        # I'm make the outer loop for both of them the same though,
313        # so this loop can easily be folded into the other loop.
314       
315        # dict of which structures have x probability of collapse.
316        # key of collapse probability
317        # value of list of struct indexes
318        struct_coll_prob = {}
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 ASSUMING THE FIRST BIN OF DEPTHS GIVE A ZERO PROBABILITY
330            depth_upper_limits = self.depth_upper_limits
331            shore_upper_limits = self.shore_upper_limits
332            collapse_probability = self.collapse_probability
333            if max_depth <= depth_upper_limits[0]:
334                continue
335            start = 1
336            for i_depth, depth_limit in enumerate(depth_upper_limits[start:]):
337                #Have to change i_depth so it indexes into the lists correctly
338                i_depth += start
339                if max_depth <= depth_limit:
340                    for i_shore, shore_limit in enumerate(shore_upper_limits):
341                        if shore_distance <= shore_limit:
342                            coll_prob = collapse_probability[i_depth][i_shore]
343                            if 0.0 == collapse_probability[i_depth][i_shore]:
344                                break
345                            struct_coll_prob.setdefault(coll_prob,[]).append(i)
346                            break
347                    break
348                           
349        return struct_coll_prob
350   
351    def _calc_collapse_structures(self, collapse_probability,
352                                  verbose_csv=False):
353        """
354        Given the collapse probabilities, throw the dice
355        and collapse some houses
356        """
357       
358        self.struct_collapsed = ['']* self.structure_count
359        if verbose_csv:
360            self.collapse_csv_info = ['']* self.structure_count
361        #for a given 'bin', work out how many houses will collapse
362        for probability, house_indexes in collapse_probability.iteritems():
363            collapse_count = round(len(house_indexes) *probability)
364           
365            if verbose_csv:
366                for i in house_indexes:
367                    # This could be sped up I think
368                    self.collapse_csv_info[i] = str(probability) + ' prob.( ' \
369                           + str(int(collapse_count)) + ' collapsed out of ' \
370                           + str(len(house_indexes)) + ')'
371            for _ in range(int(collapse_count)):
372                house_index = choice(house_indexes)
373                self.struct_damage[house_index] = 1.0
374                self.contents_damage[house_index] = 1.0
375                house_indexes.remove(house_index)
376                self.struct_collapsed[house_index] = 1
377               
378            # Warning, the collapse_probability list now lists
379            # houses that did not collapse, (though not all of them)
380            #print "",self.collapse_csv_info
381           
382#############################################################################
383if __name__ == "__main__":
384    pass 
Note: See TracBrowser for help on using the repository browser.