Source code for shakelib.sites

#!/usr/bin/env python

# stdlib imports

# third party imports
from mapio.reader import read, get_file_geodict
from mapio.grid2d import Grid2D
from mapio.geodict import GeoDict
from openquake.hazardlib.gsim.base import SitesContext
import numpy as np

# local imports
from shakelib.utils.exception import ShakeLibException


[docs]class Sites(object): """ An object to encapsulate information used to generate a GEM `SitesContext <https://github.com/gem/oq-hazardlib/blob/master/openquake/hazardlib/gsim/base.py>`__. """ # noqa def __init__( self, vs30grid, vs30measured_grid=None, backarc=None, defaultVs30=686.0 ): """ Construct a Sites object. Args: vs30grid: MapIO Grid2D object containing Vs30 values. vs30measured_grid: Boolean array indicating whether Vs30 values were measured or derived (i.e., from topographic slope). backarc: Boolean array indicating whether site is in the subduction `backarc <http://earthquake.usgs.gov/learn/glossary/?term=backarc>`__. defaultVs30: Default Vs30 value to use in locations where Vs30Grid is not specified. """ # noqa self._Vs30 = vs30grid if backarc is None: self._backarc = np.zeros_like(vs30grid.getData(), dtype=bool) else: self._backarc = backarc # Could add a check here that if backarc is provided that it's type # is bool and dimensions match vs30grid self._defaultVs30 = defaultVs30 self._vs30measured_grid = vs30measured_grid self._GeoDict = vs30grid.getGeoDict().copy() if self._GeoDict.xmin > self._GeoDict.xmax: xmin = self._GeoDict.xmin - 360 xmax = self._GeoDict.xmax else: xmin = self._GeoDict.xmin xmax = self._GeoDict.xmax self._lons = np.linspace(xmin, xmax, self._GeoDict.nx) self._lats = np.linspace( self._GeoDict.ymax, self._GeoDict.ymin, self._GeoDict.ny ) @classmethod def _create(cls, geodict, defaultVs30, vs30File, padding, resample): if vs30File is not None: fgeodict = get_file_geodict(vs30File) if not resample: if not padding: # we want something that is within and aligned geodict = fgeodict.getBoundsWithin(geodict) else: # we want something that is just aligned, since we're # padding edges geodict = fgeodict.getAligned(geodict) vs30grid = read( vs30File, samplegeodict=geodict, resample=resample, method="linear", doPadding=padding, padValue=defaultVs30, ) return vs30grid
[docs] @classmethod def fromBounds( cls, xmin, xmax, ymin, ymax, dx, dy, defaultVs30=686.0, vs30File=None, vs30measured_grid=None, backarc=None, padding=False, resample=False, ): """ Create a Sites object by defining a center point, resolution, extent, and Vs30 values. Args: xmin: X coordinate of left edge of bounds. xmax: X coordinate of right edge of bounds. ymin: Y coordinate of bottom edge of bounds. ymax: Y coordinate of top edge of bounds. dx: Resolution of desired grid in X direction. dy: Resolution of desired grid in Y direction. defaultVs30: Default Vs30 value to use if vs30File not specified. vs30File: Name of GMT or GDAL format grid file containing Vs30 values. vs30measured_grid: Boolean grid indicating whether Vs30 values were measured or derived (i.e., from slope). backarc: Boolean array indicating whether site is in the subduction `backarc <http://earthquake.usgs.gov/learn/glossary/?term=backarc>`__. padding: Boolean indicating whether or not to pad resulting Vs30 grid out to edges of input bounds. If False, grid will be clipped to the extent of the input file. resample: Boolean indicating whether or not the grid should be resampled. """ # noqa geodict = GeoDict.createDictFromBox(xmin, xmax, ymin, ymax, dx, dy) if vs30File is not None: vs30grid = cls._create(geodict, defaultVs30, vs30File, padding, resample) else: griddata = np.ones((geodict.ny, geodict.nx), dtype=np.float64) * defaultVs30 vs30grid = Grid2D(griddata, geodict) return cls( vs30grid, vs30measured_grid=vs30measured_grid, backarc=backarc, defaultVs30=defaultVs30, )
[docs] @classmethod def fromCenter( cls, cx, cy, xspan, yspan, dx, dy, defaultVs30=686.0, vs30File=None, vs30measured_grid=None, backarc=None, padding=False, resample=False, ): """ Create a Sites object by defining a center point, resolution, extent, and Vs30 values. Args: cx: X coordinate of desired center point. cy: Y coordinate of desired center point. xspan: Width of desired grid. yspan: Height of desired grid. dx: Resolution of desired grid in X direction. dy: Resolution of desired grid in Y direction. defaultVs30: Default Vs30 value to use if vs30File not specified. vs30File: Name of GMT or GDAL format grid file containing Vs30 values. vs30measured_grid: Boolean grid indicating whether Vs30 values were measured or derived (i.e., from slope). backarc: Boolean array indicating whether site is in the subduction `backarc <http://earthquake.usgs.gov/learn/glossary/?term=backarc>`__. padding: Boolean indicating whether or not to pad resulting Vs30 grid out to edges of input bounds. If False, grid will be clipped to the extent of the input file. resample: Boolean indicating whether or not the grid should be resampled. """ # noqa geodict = GeoDict.createDictFromCenter(cx, cy, dx, dy, xspan, yspan) if vs30File is not None: vs30grid = cls._create(geodict, defaultVs30, vs30File, padding, resample) else: griddata = np.ones((geodict.ny, geodict.nx), dtype=np.float64) * defaultVs30 vs30grid = Grid2D(griddata, geodict) return cls( vs30grid, vs30measured_grid=vs30measured_grid, backarc=backarc, defaultVs30=defaultVs30, )
[docs] def getSitesContext(self, lldict=None, rock_vs30=None): """ Create a SitesContext object by sampling the current Sites object. Args: lldict: Either - None, in which case the SitesContext for the complete Sites grid is returned, or - A location dictionary (elements are 'lats' and 'lons' and each is a numpy array). Each element must have the same shape. In this case the SitesContext for these locaitons is returned. rock_vs30: Either - None, in which case the SitesContext will reflect the Vs30 grid in the Sites instance, or - A float for the rock Vs30 value, in which case the SitesContext will be constructed for this constant Vs30 value. Returns: SitesContext object. Raises: ShakeLibException: When lat/lon input sequences do not share dimensionality. """ # noqa sctx = SitesContext() if lldict is not None: lats = lldict["lats"] lons = lldict["lons"] latshape = lats.shape lonshape = lons.shape if latshape != lonshape: msg = "Input lat/lon arrays must have the same dimensions" raise ShakeLibException(msg) if rock_vs30 is not None: tmp = self._Vs30.getValue(lats, lons, default=self._defaultVs30) sctx.vs30 = np.ones_like(tmp) * rock_vs30 else: sctx.vs30 = self._Vs30.getValue(lats, lons, default=self._defaultVs30) sctx.lats = lats sctx.lons = lons sctx.sids = np.array(range(np.size(lons))).reshape(lons.shape) else: sctx.lons = self._lons.copy() sctx.lats = self._lats.copy() if rock_vs30 is not None: sctx.vs30 = np.full_like(self._Vs30.getData(), rock_vs30) else: sctx.vs30 = self._Vs30.getData().copy() sctx.sids = np.array( range(np.size(sctx.lons) * np.size(sctx.lats)) ).reshape(sctx.vs30.shape) Sites._addDepthParameters(sctx) # For ShakeMap purposes, vs30 measured is always Fales sctx.vs30measured = np.zeros_like(sctx.vs30, dtype=bool) # Backarc should be a numpy array if lldict is not None: tgd = self._Vs30.getGeoDict() backarcgrid = Grid2D(self._backarc, tgd) sctx.backarc = backarcgrid.getValue(lats, lons, default=False) else: sctx.backarc = self._backarc.copy() return sctx
[docs] def getVs30Grid(self): """ Returns: Grid2D object containing Vs30 values for this Sites object. """ return self._Vs30
[docs] def getNxNy(self): """ Returns: The number of grid points in x and y. """ return self._GeoDict.nx, self._GeoDict.ny
@staticmethod def _load( vs30File, samplegeodict=None, resample=False, method="linear", doPadding=False, padValue=np.nan, ): try: vs30grid = read( vs30File, samplegeodict=samplegeodict, resample=resample, method=method, doPadding=doPadding, padValue=padValue, ) except Exception as msg1: msg = f'Load failure of {vs30File} - error message: "{str(msg1)}"' raise ShakeLibException(msg) if vs30grid.getData().dtype != np.float64: vs30grid.setData(vs30grid.getData().astype(np.float64)) return vs30grid @staticmethod def _getFileGeoDict(fname): geodict = None try: geodict = get_file_geodict(fname) except Exception as msg1: msg = f'File geodict failure with {fname} - error messages: "{str(msg1)}"' raise ShakeLibException(msg) return geodict @staticmethod def _addDepthParameters(sctx): """ Add the different depth parameters to a sites context from Vs30 values. Args: sctx: A sites context. Returns: A sites context with the depth parameters set. """ sctx.z1pt0_cy14_cal = Sites._z1pt0_from_vs30_cy14_cal(sctx.vs30) sctx.z1pt0_cy14_jpn = Sites._z1pt0_from_vs30_cy14_jpn(sctx.vs30) sctx.z1pt0_ask14_cal = Sites._z1pt0_from_vs30_ask14_cal(sctx.vs30) sctx.z1pt0_ask14_jpn = Sites._z1pt0_from_vs30_ask14_jpn(sctx.vs30) sctx.z2pt5_cb14_cal = Sites._z2pt5_from_vs30_cb14_cal(sctx.vs30) sctx.z2pt5_cb14_jpn = Sites._z2pt5_from_vs30_cb14_jpn(sctx.vs30) sctx.z1pt0_cy08 = Sites._z1pt0_from_vs30_cy08(sctx.vs30) sctx.z2pt5_cb07 = Sites._z2pt5_from_z1pt0_cb07(sctx.z1pt0_cy08) @staticmethod def _z1pt0_from_vs30_cy14_cal(vs30): """ Compute z1.0 using CY14 relationship for California. Args: vs30: Numpy array of Vs30 values in m/s. Returns: Numpy array of z1.0 in m. """ z1 = np.exp( -(7.15 / 4.0) * np.log((vs30 ** 4.0 + 571.0 ** 4) / (1360 ** 4.0 + 571.0 ** 4)) ) return z1 @staticmethod def _z1pt0_from_vs30_cy14_jpn(vs30): """ Compute z1.0 using CY14 relationship for Japan. Args: vs30: Numpy array of Vs30 values in m/s. Returns: Numpy array of z1.0 in m. """ z1 = np.exp( -(5.23 / 2.0) * np.log((vs30 ** 2.0 + 412.0 ** 2.0) / (1360 ** 2.0 + 412.0 ** 2.0)) ) return z1 @staticmethod def _z1pt0_from_vs30_ask14_cal(vs30): """ Calculate z1.0 using ASK14 relationship for California. Args: vs30: Numpy array of Vs30 values in m/s. Returns: Numpy array of z1.0 in m. """ # ASK14 define units as km, but implemented as m in OQ z1 = np.exp( -(7.67 / 4.0) * np.log((vs30 ** 4.0 + 610.0 ** 4) / (1360 ** 4.0 + 610.0 ** 4)) ) return z1 @staticmethod def _z1pt0_from_vs30_ask14_jpn(vs30): """ Calculate z1.0 using ASK14 relationship for Japan. Args: vs30: Numpy array of Vs30 values in m/s. Returns: Numpy array of z1.0 in m. """ # ASK14 define units as km, but implemented as m in OQ z1 = np.exp( -(5.32 / 2.0) * np.log((vs30 ** 2.0 + 412.0 ** 2) / (1360 ** 2.0 + 412.0 ** 2)) ) return z1 @staticmethod def _z2pt5_from_vs30_cb14_cal(vs30): """ Calculate z2.5 using CB14 relationship for California. Args: vs30: Numpy array of Vs30 values in m/s. Returns: Numpy array of z2.5 in m. *NOTE*: OQ's CampbellBozorgnia2014 class expects z2.5 to be in km! """ z2p5 = np.exp(7.089 - 1.144 * np.log(vs30)) return z2p5 @staticmethod def _z2pt5_from_vs30_cb14_jpn(vs30): """ Calculate z2.5 using CB14 relationship for Japan. Args: vs30: Numpy array of Vs30 values in m/s. Returns: Numpy array of z2.5 in m. *NOTE*: OQ's CampbellBozorgnia2014 class expects z2.5 to be in km! """ z2p5 = np.exp(5.359 - 1.102 * np.log(vs30)) return z2p5 @staticmethod def _z1pt0_from_vs30_cy08(vs30): """ Chiou and Youngs (2008) z1.0 equation. Args: vs30: Numpy array of Vs30 values in m/s. Returns: Numpy array of z1.0 in m. """ z1pt0 = np.exp(28.5 - (3.82 / 8.0) * np.log(vs30 ** 8 + 378.7 ** 8)) return z1pt0 @staticmethod def _z2pt5_from_z1pt0_cb07(z1pt0): """ Equations are from 2007 PEER report by Campbell and Bozorgnia. Args: z1pt0: Numpy array of z1.0 in m. Returns: Numpy array of z2.5 in m. """ z2pt5 = 519.0 + z1pt0 * 3.595 return z2pt5