Source code for shakelib.conversions.imc.boore_kishida_2017

"""
Module implements BooreKishida2017 class to convert between various
horizontal intensity measure components.
"""
# Standard imports
import glob
import logging
import os.path

# This is what we want for the future (Python 2.9 an beyond):
# import importlib.resources as importlib_resources
# but right now we have to support Py 3.8, so we have this:
import pkg_resources

# Third party imports
import numpy as np
from openquake.hazardlib.const import IMC
from openquake.hazardlib.imt import PGA, PGV
import pandas as pd

# Local imports
from shakelib.conversions.convert_imc import ComponentConverter


[docs]class BooreKishida2017(ComponentConverter): """ This class implements the Boore and Kishida (2017) conversions for horizontal intensity measure components. This class explicitly supports the following subset of the conversions provided by B&K (specified as OpenQuake IMCs): - RotD50 <=> GMRotI50 - RotD50 <=> GEOMETRIC_MEAN (i.e., Geometric Mean "as recorded") - RotD100 <=> RotD50 - RotD50 <=> GREATER_OF_TWO_HORIZONTAL - RotD100 <=> GREATER_OF_TWO_HORIZONTAL - GMRotI50 <=> GREATER_OF_TWO_HORIZONTAL - GEOMETRIC_MEAN <=> GREATER_OF_TWO_HORIZONTAL Chain conversions are supported when using `convertAmps`. Otherwise conversions must be done in two+ steps using `convertAmpsOnce`. For IMCs not explicitly supported by B&K, we assume the IMC is equivalent to the geometric mean (which B&K call GM_AR). Notes - Assumes ALL unknown IMC types are GEOMETRIC_MEAN. References Boore, D.M. and T. Kishida (2017). Relations between some horizontal-component ground-motion intensity measures used in practice. Bulletin of the Seismological Society of America, 107(1), 334-343, doi: 10.1785/0120160250. Beyer, K., & Bommer, J. J. (2006). Relationships between median values and between aleatory variabilities for different definitions of the horizontal component of motion. Bulletin of the Seismological Society of America, 96(4A), 1512-1522. """ def __init__(self, imc_in, imc_out): super().__init__() if type(imc_in) == IMC: self.imc_in = imc_in elif type(imc_in) == str: imc_from_str = self.imc_from_str(imc_in) if imc_from_str is None: logging.warning( "Unknown IMC string '%s' in input, using Geometric Mean", imc_in, ) self.imc_in = IMC.GEOMETRIC_MEAN else: self.imc_in = imc_from_str else: logging.warning( "Unknown object passed as input %s, using Geometric Mean", type(imc_in), ) self.imc_in = IMC.GEOMETRIC_MEAN if type(imc_out) == IMC: self.imc_out = imc_out elif type(imc_out) == str: imc_from_str = self.imc_from_str(imc_out) if imc_from_str is None: logging.warning( "Unknown IMC string '%s' in output, using Geometric Mean", imc_out, ) self.imc_out = IMC.GEOMETRIC_MEAN else: self.imc_out = imc_from_str else: logging.warning( "Unknown object passed as opuput %s, using Geometric Mean", type(imc_out), ) self.imc_out = IMC.GEOMETRIC_MEAN # Possible conversions self.conversion_graph = { IMC.RotD50: set( [ IMC.GMRotI50, IMC.GEOMETRIC_MEAN, IMC.RotD100, IMC.GREATER_OF_TWO_HORIZONTAL, IMC.RANDOM_HORIZONTAL, IMC.HORIZONTAL, IMC.MEDIAN_HORIZONTAL, ] ), IMC.GMRotI50: set([IMC.RotD50, IMC.GREATER_OF_TWO_HORIZONTAL]), IMC.GEOMETRIC_MEAN: set([IMC.RotD50, IMC.GREATER_OF_TWO_HORIZONTAL]), IMC.RotD100: set([IMC.RotD50, IMC.GREATER_OF_TWO_HORIZONTAL]), IMC.GREATER_OF_TWO_HORIZONTAL: set( [ IMC.RotD50, IMC.GMRotI50, IMC.GEOMETRIC_MEAN, IMC.RotD100, IMC.RANDOM_HORIZONTAL, IMC.HORIZONTAL, IMC.MEDIAN_HORIZONTAL, ] ), IMC.HORIZONTAL: set([IMC.GREATER_OF_TWO_HORIZONTAL, IMC.RotD50]), IMC.MEDIAN_HORIZONTAL: set([IMC.GREATER_OF_TWO_HORIZONTAL, IMC.RotD50]), IMC.RANDOM_HORIZONTAL: set([IMC.GREATER_OF_TWO_HORIZONTAL, IMC.RotD50]), } # Check if any imc values are unknown. If they are, convert # to GEOMETRIC_MEAN self.checkUnknown() # Get shortest conversion "path" between imc_in and imc_out self.path = self.getShortestPath( self.conversion_graph, self.imc_in, self.imc_out )
[docs] def convertAmpsOnce(self, imt, amps, rrups=None, mag=None): """ Return an array of amps converted from one IMC to another. Args: imt (OpenQuake IMT): The intensity measure type of the input ground motions. Valid IMTs are PGA, PGV, and SA. amps (array): A numpy array of the (logged) ground motions to be converted. rrups (array): A numpy array of the same shape as amps, containing the rupture distances of the ground motions. mag (float): The earthquake magnitude. Returns: array: A numpy array of converted ground motions (logged). Raises: ValueError: If mag and rrup are none or the IMT is not an allowed type. """ # Check if mag and rrups are real values if mag is None or rrups is None: raise ValueError("No magnitude or rupture distances specified.") # Verify that the conversion is possible self._verifyConversion(self.imc_in, self.imc_out) # Return original amps if imc_in and imc_out are the same if self.pars is None: return amps.copy() # Get coeffecients (sigma, c0, r1, m1, m2) = self._getParamsFromIMT(imt) # Limit magnitude and rupture distances rrups_clipped = np.clip(rrups, 1e-2, 400) if mag < 2: mag = 2.0 elif mag > 9: mag = 9.0 # Calculate conversion variable ln_ratio = ( c0 + r1 * np.log(rrups_clipped / 50) + m1 * (mag - 5.5) + m2 * (mag - 5.5) ** 2 ) # # The B&K file naming convention has things like D100D50, which # means the parameters give the (log) ratio of RotD100/RotD50, # but we use the convention that RotD100 would be the input IMC # and RotD50 would be the output IMC, so we reverse the sense # of the conversion here. # if self.forward: amps = amps - ln_ratio else: amps = amps + ln_ratio return amps
[docs] def convertSigmasOnce(self, imt, sigmas): """ Return an array of standard deviations converted from one IMC to another. Note that the action of this method is to always increase the input standard deviations. Thus, while converting from one IMC to another and then back again will yield the original ground motions via convertAmps(), the standard deviations will be inflated by both conversions via this method. Args: imt (OpenQuake IMT): The intensity measure type of the input ground motions. Valid IMTs are PGA, PGV, and SA. sigmas (array): A numpy array of the standard deviations of the logged ground motions. Returns: array: A numpy array of converted standard deviations. Raises: ValueError: If mag and rrup are none or the IMT is not an allowed type. """ # Verify that the conversion is possible self._verifyConversion(self.imc_in, self.imc_out) # Return original sigmas if imc_in and imc_out are the same if self.pars is None: return sigmas.copy() # Get coeffecients (sigma, c0, r1, m1, m2) = self._getParamsFromIMT(imt) # Calculate conversion sigmas = np.sqrt(sigmas ** 2 + sigma ** 2) return sigmas
def _getParamsFromIMT(self, imt): """ Helper function to return (possibly interpolated) conversion parameters for a given IMT. Args: imt (OpenQuake IMT): The intensity measure type of the input ground motions. Valid IMTs are PGA, PGV, and SA. Returns: (float, float, float, float, float): Coeffients for conversion. """ if imt == PGA(): sigma = self.pars["sigma"][0] c0 = self.pars["c0smooth"][0] r1 = self.pars["r1smooth"][0] m1 = self.pars["m1smooth"][0] m2 = self.pars["m2smooth"][0] elif imt == PGV(): sigma = self.pars["sigma"][1] c0 = self.pars["c0smooth"][1] r1 = self.pars["r1smooth"][1] m1 = self.pars["m1smooth"][1] m2 = self.pars["m2smooth"][1] elif "SA" in imt.string: imt_per = imt.period pa = self.pars["per"][2:] sigma = np.interp(imt_per, pa, self.pars["sigma"][2:]) c0 = np.interp(imt_per, pa, self.pars["c0smooth"][2:]) r1 = np.interp(imt_per, pa, self.pars["r1smooth"][2:]) m1 = np.interp(imt_per, pa, self.pars["m1smooth"][2:]) m2 = np.interp(imt_per, pa, self.pars["m2smooth"][2:]) else: raise ValueError(f"Unknown IMT: {str(imt)}") return (sigma, c0, r1, m1, m2) @staticmethod def _imcPairToFile(imc_in, imc_out): """ Helper function to find the name of the file representing the conversion. Returns: (str, bool): The filename and a boolean 'forward' indicating whether the conversion should be done in the forward (True) or inverse (False) direction. If filename is None, then no appropriate conversion file coule be found; if it is the string 'Null', then imc_in and imc_out evaluated to be the same. """ # This is what we want for the future (Python 2.9 an beyond): # datadir = importlib_resources.files("shakelib.conversions.imc") / "data" # but right now we have to support Py 3.8, so we have this: datadir = pkg_resources.resource_filename("shakelib.conversions.imc", "data") conv_files = glob.glob(os.path.join(datadir, "*.csv")) stub1 = BooreKishida2017._imcToFilestr(imc_in) stub2 = BooreKishida2017._imcToFilestr(imc_out) if stub1 == stub2: # No conversion necessary return ("Null", True) # # Look for the conversion from imc_in -> imc_out # stub = stub1 + stub2 filelist = [name for name in conv_files if stub in name] if len(filelist) == 1: return (filelist[0], True) # # Now try the conversion from imc_out -> imc_in # stub = stub2 + stub1 filelist = [name for name in conv_files if stub in name] if len(filelist) == 1: return (filelist[0], False) # # Can't find anything # return (None, None) @staticmethod def _imcToFilestr(oq_imc): """ Helper function to convert an OpenQuake IMC into part of the Boore & Kishida file name. """ if oq_imc == IMC.RotD50: return "D50" elif oq_imc == IMC.RotD100: return "D100" elif oq_imc == IMC.GMRotI50: return "GM50" elif ( oq_imc == IMC.GEOMETRIC_MEAN or oq_imc == IMC.HORIZONTAL or oq_imc == IMC.RANDOM_HORIZONTAL or oq_imc == IMC.MEDIAN_HORIZONTAL ): return "GMAR" elif oq_imc == IMC.GREATER_OF_TWO_HORIZONTAL: return "Larger" else: # # For less common IMCs, Beyer & Bommer (2006) found most # of them to be more or less equivalent to geometric mean # logging.warning(f"Can't handle IMC {oq_imc}, using GMAR") return "GMAR" def _verifyConversion(self, imc_in, imc_out=None): """ Helper method to ensure that the conversion is possible. Args: imc_in (IMC): OpenQuake IMC type of the input amp array. imc_out (IMC): Desired OpenQuake IMC type of the output amps. Default is None. Raises: ValueError if imc_in or imc_out are not valid.. """ filename, forward = self._imcPairToFile(imc_in, imc_out) if filename is None: raise ValueError(f"Can't find a conversion file for {imc_in} and {imc_out}") self.forward = forward if filename == "Null": # Null conversion -- imc_in and imc_out are either identical # or at least functionally equivalent self.pars = None else: self.pars = pd.read_csv(filename)