Source code for openquake.calculators.scenario_damage

# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2014-2018 GEM Foundation
#
# OpenQuake is free software: you can redistribute it and/or modify it
# under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# OpenQuake is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>.

import itertools
import logging
import numpy

from openquake.risklib import scientific, riskmodels
from openquake.calculators import base

U16 = numpy.uint16
U64 = numpy.uint64
F32 = numpy.float32
F64 = numpy.float64


[docs]def dist_by_asset(data, multi_stat_dt, number): """ :param data: array of shape (N, R, L, 2, ...) :param multi_stat_dt: numpy dtype for statistical outputs :param number: expected number of units per asset :returns: array of shape (N, R) with records of type multi_stat_dt """ N, R, L = data.shape[:3] out = numpy.zeros((N, R), multi_stat_dt) for l, lt in enumerate(multi_stat_dt.names): out_lt = out[lt] for n, r in itertools.product(range(N), range(R)): mean, stddev = data[n, r, l] out_lt[n, r] = (mean, stddev) # sanity check on the sum over all damage states if abs(mean.sum() / number[n] - 1) > 1E-3: logging.warn( 'Asset #%d, rlz=%d, expected %s, got %s for %s damage', n, r, number[n], mean.sum(), lt) return out
[docs]def scenario_damage(riskinput, riskmodel, param, monitor): """ Core function for a damage computation. :param riskinput: a :class:`openquake.risklib.riskinput.RiskInput` object :param riskmodel: a :class:`openquake.risklib.riskinput.CompositeRiskModel` instance :param monitor: :class:`openquake.baselib.performance.Monitor` instance :param param: dictionary of extra parameters :returns: a dictionary {'d_asset': [(l, r, a, mean-stddev), ...], 'd_event': damage array of shape R, L, E, D, 'c_asset': [(l, r, a, mean-stddev), ...], 'c_event': damage array of shape R, L, E} `d_asset` and `d_tag` are related to the damage distributions whereas `c_asset` and `c_tag` are the consequence distributions. If there is no consequence model `c_asset` is an empty list and `c_tag` is a zero-valued array. """ c_models = param['consequence_models'] L = len(riskmodel.loss_types) R = riskinput.hazard_getter.num_rlzs D = len(riskmodel.damage_states) E = param['number_of_ground_motion_fields'] result = dict(d_asset=[], d_event=numpy.zeros((E, R, L, D), F64), c_asset=[], c_event=numpy.zeros((E, R, L), F64)) for outputs in riskmodel.gen_outputs(riskinput, monitor): r = outputs.rlzi for l, damages in enumerate(outputs): loss_type = riskmodel.loss_types[l] c_model = c_models.get(loss_type) for a, fraction in enumerate(damages): asset = outputs.assets[a] taxo = riskmodel.taxonomy[asset.taxonomy] damages = fraction * asset.number result['d_event'][:, r, l] += damages # shape (E, D) if c_model: # compute consequences means = [par[0] for par in c_model[taxo].params] # NB: we add a 0 in front for nodamage state c_ratio = numpy.dot(fraction, [0] + means) consequences = c_ratio * asset.value(loss_type) result['c_asset'].append( (l, r, asset.ordinal, scientific.mean_std(consequences))) result['c_event'][:, r, l] += consequences # TODO: consequences for the occupants result['d_asset'].append( (l, r, asset.ordinal, scientific.mean_std(damages))) result['gmdata'] = riskinput.gmdata return result
[docs]@base.calculators.add('scenario_damage') class ScenarioDamageCalculator(base.RiskCalculator): """ Scenario damage calculator """ pre_calculator = 'scenario' core_task = scenario_damage is_stochastic = True
[docs] def pre_execute(self): if 'gmfs' in self.oqparam.inputs: self.pre_calculator = None base.RiskCalculator.pre_execute(self) E = self.oqparam.number_of_ground_motion_fields self.param['number_of_ground_motion_fields'] = E self.param['consequence_models'] = riskmodels.get_risk_models( self.oqparam, 'consequence') self.riskinputs = self.build_riskinputs('gmf', num_events=E) self.param['tags'] = list(self.assetcol.tagcol)
[docs] def post_execute(self, result): """ Compute stats for the aggregated distributions and save the results on the datastore. """ dstates = self.riskmodel.damage_states ltypes = self.riskmodel.loss_types L = len(ltypes) R = len(self.rlzs_assoc.realizations) D = len(dstates) N = len(self.assetcol) E = self.oqparam.number_of_ground_motion_fields # damage distributions dt_list = [] for ltype in ltypes: dt_list.append((ltype, numpy.dtype([('mean', (F32, D)), ('stddev', (F32, D))]))) multi_stat_dt = numpy.dtype(dt_list) d_asset = numpy.zeros((N, R, L, 2, D), F32) for (l, r, a, stat) in result['d_asset']: d_asset[a, r, l] = stat self.datastore['dmg_by_asset'] = dist_by_asset( d_asset, multi_stat_dt, self.assetcol.array['number']) dmg_dt = [(ds, F32) for ds in self.riskmodel.damage_states] d_event = numpy.zeros((E, R, L), dmg_dt) for d, ds in enumerate(self.riskmodel.damage_states): d_event[ds] = result['d_event'][:, :, :, d] self.datastore['dmg_by_event'] = d_event # consequence distributions if result['c_asset']: dtlist = [('eid', U64), ('rlzi', U16), ('loss', (F32, L))] stat_dt = numpy.dtype([('mean', F32), ('stddev', F32)]) c_asset = numpy.zeros((N, R, L), stat_dt) for (l, r, a, stat) in result['c_asset']: c_asset[a, r, l] = stat multi_stat_dt = self.oqparam.loss_dt(stat_dt) self.datastore['losses_by_asset'] = c_asset self.datastore['losses_by_event'] = numpy.fromiter( ((eid, rlzi, F32(result['c_event'][eid, rlzi])) for rlzi in range(R) for eid in range(E)), dtlist) # save gmdata self.gmdata = result['gmdata'] for arr in self.gmdata.values(): arr[-2] = self.oqparam.number_of_ground_motion_fields # events base.save_gmdata(self, R)