Source code for safe.impact_functions.earthquake.pager_earthquake_fatality_model.impact_function

# coding=utf-8
"""InaSAFE Disaster risk tool by Australian Aid - PAGER Earthquake
Impact Function on Population.

Contact : [email protected]

.. note:: This program is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
     the Free Software Foundation; either version 2 of the License, or
     (at your option) any later version.

"""

__author__ = 'lucernae'
__date__ = '24/03/15'
__revision__ = '$Format:%H$'
__copyright__ = ('Copyright 2012, Australia Indonesia Facility for '
                 'Disaster Reduction')

import math
import numpy
import heapq

from safe.common.utilities import OrderedDict
from safe.impact_functions.earthquake.\
    itb_earthquake_fatality_model.impact_function import ITBFatalityFunction
from safe.impact_functions.earthquake\
    .pager_earthquake_fatality_model.metadata_definitions import \
    PAGFatalityMetadata


[docs]class PAGFatalityFunction(ITBFatalityFunction): # noinspection PyUnresolvedReferences """USGS Pager fatality estimation model. Fatality rate(MMI) = cum. standard normal dist(1/BETA * ln(MMI/THETA)). Reference: Jaiswal, K. S., Wald, D. J., and Hearne, M. (2009a). Estimating casualties for large worldwide earthquakes using an empirical approach. U.S. Geological Survey Open-File Report 2009-1136. v1.0: Theta: 14.05, Beta: 0.17, Zeta 2.15 Jaiswal, K, and Wald, D (2010) An Empirical Model for Global Earthquake Fatality Estimation Earthquake Spectra, Volume 26, No. 4, pages 1017–1037 v2.0: Theta: 13.249, Beta: 0.151, Zeta: 1.641) (http://pubs.usgs.gov/of/2009/1136/pdf/ PAGER%20Implementation%20of%20Empirical%20model.xls) """ _metadata = PAGFatalityMetadata() def __init__(self): super(PAGFatalityFunction, self).__init__() self.hardcoded_parameters = OrderedDict([ ('Theta', 13.249), ('Beta', 0.151), ('Zeta', 1.641), # Model coefficients # Rates of people displaced for each MMI level # FIXME: should be independent from fatality model - Hyeuk ('displacement_rate', { 2: 0.0, 3: 0.0, 4: 0.0, 5: 0.0, 6: 1.0, 7: 1.0, 8: 1.0, 9: 1.0, 10: 1.0}), ('mmi_range', range(2, 11)), ('step', 0.5), # Threshold below which layer should be transparent ('tolerance', 0.01), ('calculate_displaced_people', True), ('magnitude_bin', numpy.power(10, range(0, 6), dtype=float)) ]) # noinspection PyPep8Naming
[docs] def cdf_normal(self, x): """Cumulative distribution function of standard normal distribution. Logic based on http://en.wikipedia.org/wiki/Normal_distribution :param x: random variable x :type x: float :returns: phi of (x) :rtype: float """ return (1.0 + math.erf(x / math.sqrt(2.0))) / 2.0
[docs] def fatality_rate(self, mmi): """Pager method to compute fatality rate. :param mmi: MMI :type mmi: float :returns: Fatality rate calculated as: lognorm.cdf(mmi, shape=Beta, scale=Theta) :rtype: float """ theta = self.hardcoded_parameters['Theta'] beta = self.hardcoded_parameters['Beta'] x = math.log(mmi / theta) / beta return self.cdf_normal(x)
@staticmethod
[docs] def round_to_sum(l, r): """ Round a list of numbers while maintaining the sum. http://stackoverflow.com/questions/15769948/ round-a-python-list-of-numbers-and-maintain-the-sum :param l: array :type l: list(float) :param r: decimal place :type r: int :returns: A list of rounded numbers whose sum is equal to the sum of the list of input numbers. :rtype: list """ q = 10 ** (-r) d = (round(sum(l), r) - sum([round(x, r) for x in l])) * (10 ** r) d = int(d) if d == 0: return [round(x, r) for x in l] elif d in [-1, 1]: c, _ = max(enumerate(l), key=lambda x: math.copysign( 1, d) * math.fmod(x[1] - 0.5 * q, q)) return [round(x, r) + q * math.copysign(1, d) if i == c else round( x, r) for (i, x) in enumerate(l)] else: c = [i for i, _ in heapq.nlargest(abs(d), enumerate( l), key=lambda x: math.copysign(1, d) * math.fmod( x[1] - 0.5 * q, q))] return [round(x, r) + q * math.copysign( 1, d) if i in c else round(x, r) for (i, x) in enumerate(l)]
[docs] def compute_probability(self, total_fatalities): """Pager method compute probaility of fatality in each magnitude bin. (0,1), (1,10), (10,10^2), (10^2,10^3), (10^3, 10^4), (10^4, 10^5), (10^5, 10^6+) :param total_fatalities: List of total fatalities in each MMI class. :type total_fatalities: int, float :returns: Probability of fatality magnitude bin from lognorm.cdf(bin, shape=Zeta, scale=total_fatalities) :rtype: list(float) """ zeta = self.hardcoded_parameters['Zeta'] magnitude_bin = self.hardcoded_parameters['magnitude_bin'] cprob = numpy.ones_like(magnitude_bin, dtype=float) if total_fatalities > 1: for (i, mbin) in enumerate(magnitude_bin): x = math.log(mbin / total_fatalities) / zeta cprob[i] = self.cdf_normal(x) cprob = numpy.append(cprob, 1.0) # 1000+ # percentage prob = numpy.hstack((cprob[0], numpy.diff(cprob))) * 100.0 return self.round_to_sum(prob, 0) # rounding to decimal