Source code for safe.impact_functions.earthquake.padang_building_impact_model

"""Impact function based on Padang 2009 post earthquake survey

This impact function estimates percentual damage to buildings as a
function of ground shaking measured in MMI.
Buildings are assumed to fall the 9 classes below as described in
the Geoscience Australia/ITB 2009 Padang earthquake
survey (http://trove.nla.gov.au/work/38470066).

Class Building Type                              Median (MMI)  Beta (MMI)
-------------------------------------------------------------------------
1     URM with river rock walls                        7.5     0.11
2     URM with Metal Roof                              8.3     0.1
3     Timber frame with masonry in-fill                8.8     0.11
4     RC medium rise Frame with Masonry in-fill walls  8.4     0.05
5     Timber frame with stucco in-fill                 9.2     0.11
6     Concrete Shear wall  high rise* Hazus C2H        9.7     0.15
7     RC low rise Frame with Masonry in-fill walls     9       0.08
8     Confined Masonry                                 8.9     0.07
9     Timber frame residential                        10.5     0.15
"""

from safe.impact_functions.core import (FunctionProvider,
                                        get_hazard_layer,
                                        get_exposure_layer,
                                        get_question,
                                        format_int)
from safe.storage.vector import Vector
from safe.common.utilities import ugettext as tr
from safe.common.numerics import lognormal_cdf
from safe.common.tables import Table, TableRow
from safe.impact_functions.mappings import osm2padang, sigab2padang
from safe.engine.interpolation import assign_hazard_values_to_exposure_data


# Damage curves for each of the nine classes derived from the Padang survey
damage_curves = {1: dict(median=7.5, beta=0.11),
                 2: dict(median=8.3, beta=0.1),
                 3: dict(median=8.8, beta=0.11),
                 4: dict(median=8.4, beta=0.05),
                 5: dict(median=9.2, beta=0.11),
                 6: dict(median=9.7, beta=0.15),
                 7: dict(median=9.0, beta=0.08),
                 8: dict(median=8.9, beta=0.07),
                 9: dict(median=10.5, beta=0.15)}


class PadangEarthquakeBuildingDamageFunction(FunctionProvider):
[docs] """Risk plugin for Padang earthquake damage to buildings :param requires category=='hazard' and \ subcategory=='earthquake' and \ layertype=='raster' and \ unit=='MMI' :param requires category=='exposure' and \ subcategory=='structure' and \ layertype=='vector' and \ datatype in ['osm', 'itb', 'sigab'] """ title = tr('Be damaged depending on building type') def run(self, layers):
[docs] """Risk plugin for Padang building survey """ # Extract data H = get_hazard_layer(layers) # Ground shaking E = get_exposure_layer(layers) # Building locations question = get_question(H.get_name(), E.get_name(), self) # Map from different kinds of datasets to Padang vulnerability classes datatype = E.get_keywords()['datatype'] vclass_tag = 'VCLASS' if datatype.lower() == 'osm': # Map from OSM attributes Emap = osm2padang(E) elif datatype.lower() == 'sigab': # Map from SIGAB attributes Emap = sigab2padang(E) else: Emap = E # Interpolate hazard level to building locations I = assign_hazard_values_to_exposure_data(H, Emap, attribute_name='MMI') # Extract relevant numerical data attributes = I.get_data() N = len(I) # Calculate building damage count_high = count_medium = count_low = count_none = 0 for i in range(N): mmi = float(attributes[i]['MMI']) building_type = Emap.get_data(vclass_tag, i) damage_params = damage_curves[building_type] beta = damage_params['beta'] median = damage_params['median'] percent_damage = lognormal_cdf(mmi, median=median, sigma=beta) * 100 # Add calculated impact to existing attributes attributes[i][self.target_field] = percent_damage # Calculate statistics if percent_damage < 10: count_none += 1 if 10 <= percent_damage < 33: count_low += 1 if 33 <= percent_damage < 66: count_medium += 1 if 66 <= percent_damage: count_high += 1 # Generate impact report table_body = [question, TableRow([tr('Buildings'), tr('Total')], header=True), TableRow([tr('All'), N]), TableRow([tr('No damage'), format_int(count_none)]), TableRow([tr('Low damage'), format_int(count_low)]), TableRow([tr('Medium damage'), format_int(count_medium)]), TableRow([tr('High damage'), format_int(count_high)])] table_body.append(TableRow(tr('Notes'), header=True)) table_body.append(tr('Levels of impact are defined by post 2009 ' 'Padang earthquake survey conducted by Geoscience ' 'Australia and Institute of Teknologi Bandung.')) table_body.append(tr('Unreinforced masonry is assumed where no ' 'structural information is available.')) impact_summary = Table(table_body).toNewlineFreeString() impact_table = impact_summary map_title = tr('Earthquake damage to buildings') # Create style style_classes = [dict(label=tr('No damage'), min=0, max=10, colour='#00ff00', transparency=0), dict(label=tr('Low damage'), min=10, max=33, colour='#ffff00', transparency=0), dict(label=tr('Medium damage'), min=33, max=66, colour='#ffaa00', transparency=0), dict(label=tr('High damage'), min=66, max=100, colour='#ff0000', transparency=0)] style_info = dict(target_field=self.target_field, style_classes=style_classes) # Create vector layer and return V = Vector(data=attributes, projection=E.get_projection(), geometry=E.get_geometry(), name='Estimated pct damage', keywords={'impact_summary': impact_summary, 'impact_table': impact_table, 'map_title': map_title, 'target_field': self.target_field}, style_info=style_info) return V