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

Last change on this file since 4474 was 4326, checked in by duncan, 18 years ago

damage modelling change. Get a message quickly if the csv label is not present.

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