"""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