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

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

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