Source code for shakelib.gmice.ak07

# third party imports
import numpy as np

# stdlib imports
from openquake.hazardlib.imt import PGA, PGV, SA
from shakelib.gmice.gmice import GMICE


[docs]class AK07(GMICE): """ Implements the ground motion intensity conversion equations (GMICE) of Atkinson and Kaka (2007). This module implements a simplified version. References: Atkinson, Gail & Kaka, SanLinn. (2007). Relationships between Felt Intensity and Instrumental Ground Motion in the Central United States and California. Bulletin of The Seismological Society of America - BULL SEISMOL SOC AMER. 97. 497-510. 10.1785/0120060154. """ # ----------------------------------------------------------------------- # # MMI = C1 + C2 * log(Y) for c2->T1 < log(Y) <= T1 # MMI = C3 + C4 * log(Y) for log(Y) > T1 # # or # # MMI = c2->C1 + c2->C2 * log(Y) + C5 + C6 * M + C7 * log(D) # for log(Y) <= c2->T1 # MMI = C1 + C2 * log(Y) + C5 + C6 * log(D) + C7 * M # for c2->T1 < log(Y) <= T1 # MMI = C3 + C4 * log(Y) + C5 + C6 * M + C7 * log(D) # for log(Y) > T1 # # Limit the distance residuals to between 10 and 300 km. # Limit the magnitude residuals to between M3.0 and M7.3. # # Constants taken from Table 4 and 5 in the reference material # # ----------------------------------------------------------------------- def __init__(self): super().__init__() self.min_max = (1.0, 10.0) self.name = "Atkinson and Kaka (2007)" self.scale = "scale_ak07.ps" self._constants = { self._pga: { "C1": 2.65, "C2": 1.39, "C3": -1.91, "C4": 4.09, "C5": -1.96, "C6": 0.02, "C7": 0.98, "T1": 1.69, "T2": 5, "SMMI": 0.89, }, self._pgv: { "C1": 4.37, "C2": 1.32, "C3": 3.54, "C4": 3.03, "C5": 0.47, "C6": -0.19, "C7": 0.26, "T1": 0.48, "T2": 5, "SMMI": 0.76, }, self._sa03: { "C1": 2.4, "C2": 1.36, "C3": -1.83, "C4": 3.56, "C5": -0.11, "C6": -0.2, "C7": 0.64, "T1": 1.92, "T2": 5, "SMMI": 0.79, }, self._sa10: { "C1": 3.23, "C2": 1.18, "C3": 0.57, "C4": 2.95, "C5": 1.92, "C6": -0.39, "C7": 0.04, "T1": 1.5, "T2": 5, "SMMI": 0.73, }, self._sa30: { "C1": 3.72, "C2": 1.29, "C3": 1.99, "C4": 3.0, "C5": 2.24, "C6": -0.33, "C7": -0.31, "T1": 1, "T2": 5, "SMMI": 0.72, }, } self.DEFINED_FOR_INTENSITY_MEASURE_TYPES = set([PGA, PGV, SA]) self.DEFINED_FOR_SA_PERIODS = set([0.3, 1.0, 3.0])
[docs] def getMIfromGM(self, amps, imt, dists=None, mag=None): """ Function to compute macroseismic intensity from ground-motion intensity. Supported ground-motion IMTs are PGA, PGV and PSA at 0.3, 1.0, and 2.0 sec periods. Args: amps (ndarray): Ground motion amplitude; natural log units; g for PGA and PSA, cm/s for PGV. imt (OpenQuake IMT): Type the input amps (must be one of PGA, PGV, or SA). Supported SA periods are 0.3, 1.0, and 3.0 sec. `[link] <http://docs.openquake.org/oq-hazardlib/master/imt.html>` dists (ndarray): Numpy array of distances from rupture (km). mag (float): Earthquake magnitude. Returns: ndarray of Modified Mercalli Intensity and ndarray of dMMI / dln(amp) (i.e., the slope of the relationship at the point in question). """ lfact = np.log10(np.e) c = self._getConsts(imt) if dists is not None and mag is not None: doresid = True # Limit distances and take the log10 ldd = np.log10(np.clip(dists, 10, 300)) # Llimit magnitudes lmm = np.clip(mag, 3.0, 7.3) else: doresid = False # Account for ln form of amplitudes # Convert (for accelerations) from ln(g) to cm/s^2 # then take the log10 # if imt != self._pgv: units = 981.0 else: units = 1.0 # # Math: log10(981 * exp(amps)) = log10(981) + log10(exp(amps)) # = log10(981) + amps * log10(e) # For PGV, just convert ln(amp) to log10(amp) by multiplying # by log10(e) # lamps = np.log10(units) + amps * lfact mmi = np.zeros_like(amps) dmmi_damp = np.zeros_like(amps) # # This is the lower segment of the bi-linear fit # where log(Y) is less than T1 # idx = lamps < c["T1"] mmi[idx] = c["C1"] + c["C2"] * lamps[idx] dmmi_damp[idx] = c["C2"] * lfact # # This is the upper segment of the bi-linear fit # where log(Y) is greater than or equal to T1 # idx = lamps >= c["T1"] mmi[idx] = c["C3"] + c["C4"] * lamps[idx] dmmi_damp[idx] = c["C4"] * lfact # Inclusion of residuals if magnitude and # distance information is available if doresid: mmi += c["C5"] + c["C6"] * lmm + c["C7"] * ldd # Limit mmi values mmi = np.clip(mmi, 1.0, 10.0) return mmi, dmmi_damp
[docs] def getGMfromMI(self, mmi, imt, dists=None, mag=None): """ Function to tcompute ground-motion intensity from macroseismic intensity. Supported IMTs are PGA, PGV and PSA for 0.3, 1.0, and 3.0 sec periods. Args: mmi (ndarray): Macroseismic intensity. imt (OpenQuake IMT): IMT of the requested ground-motions intensities (must be one of PGA, PGV, or SA). `[link] <http://docs.openquake.org/oq-hazardlib/master/imt.html>` dists (ndarray): Rupture distances (km) to the corresponding MMIs. mag (float): Earthquake magnitude. Returns: Ndarray of ground motion intensity in natural log of g for PGA and PSA, and natural log cm/s for PGV; ndarray of dln(amp) / dMMI (i.e., the slope of the relationship at the point in question). """ lfact = np.log10(np.e) c = self._getConsts(imt) mmi = mmi.copy() # Set nan values to 1 ix_nan = np.isnan(mmi) mmi[ix_nan] = 1.0 if dists is not None and mag is not None: doresid = True # Limit distances and take the log10 ldd = np.log10(np.clip(dists, 10, 300)) # Limit magnitudes lmm = np.clip(mag, 3.0, 7.3) else: doresid = False # Inclusion of residuals if magnitude and # distance information is available if doresid: mmi -= c["C5"] + c["C6"] * lmm + c["C7"] * ldd pgm = np.zeros_like(mmi) dpgm_dmmi = np.zeros_like(mmi) # # This is the lower segment of the bi-linear fit # where MMI is less than I5 # idx = mmi < c["T2"] pgm[idx] = np.power(10, (mmi[idx] - c["C1"]) / c["C2"]) dpgm_dmmi[idx] = 1.0 / (c["C2"] * lfact) # # This is the upper segment of the bi-linear fit # where MMI is greater than or equal to I5 # idx = mmi >= c["T2"] pgm[idx] = np.power(10, (mmi[idx] - c["C3"]) / c["C4"]) dpgm_dmmi[idx] = 1.0 / (c["C4"] * lfact) if imt != self._pgv: units = 981.0 else: units = 1.0 # Return a ln(amp) value. Convert PGA to from cm/s^2 to g pgm /= units pgm = np.log(pgm) # Set nan values back from 1 to nan pgm[ix_nan] = np.nan dpgm_dmmi[ix_nan] = np.nan return pgm, dpgm_dmmi
[docs] def getGM2MIsd(self): """ Return a dictionary of standard deviations for the ground-motion to MMI conversion. The keys are the ground motion types. Returns: Dictionary of GM to MI sigmas (in MMI units). """ return { self._pga: self._constants[self._pga]["SMMI"], self._pgv: self._constants[self._pgv]["SMMI"], self._sa03: self._constants[self._sa03]["SMMI"], self._sa10: self._constants[self._sa10]["SMMI"], self._sa30: self._constants[self._sa30]["SMMI"], }
[docs] def getMI2GMsd(self): """ Return a dictionary of standard deviations for the MMI to ground-motion conversion. The keys are the ground motion types. Returns: Dictionary of MI to GM sigmas (ln(PGM) units). """ # # Need to convert log10 to ln units # lfact = np.log(10.0) return { self._pga: lfact * 0.57, self._pgv: lfact * 0.52, self._sa03: lfact * 0.63, self._sa10: lfact * 0.57, self._sa30: lfact * 0.81, }
def _getConsts(self, imt): """ Helper function to get the constants. """ if ( imt != self._pga and imt != self._pgv and imt != self._sa03 and imt != self._sa10 and imt != self._sa30 ): raise ValueError("Invalid IMT " + str(imt)) c = self._constants[imt] return c