Source code for openquake.hazardlib.site

# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2012-2017 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/>.

"""
Module :mod:`openquake.hazardlib.site` defines :class:`Site`.
"""
import numpy

from openquake.baselib.python3compat import range
from openquake.baselib.slots import with_slots
from openquake.baselib.general import split_in_blocks
from openquake.hazardlib.geo.mesh import Mesh


@with_slots
[docs]class Site(object): """ Site object represents a geographical location defined by its position as well as its soil characteristics. :param location: Instance of :class:`~openquake.hazardlib.geo.point.Point` representing where the site is located. :param vs30: Average shear wave velocity in the top 30 m, in m/s. :param vs30measured: Boolean value, ``True`` if ``vs30`` was measured on that location and ``False`` if it was inferred. :param z1pt0: Vertical distance from earth surface to the layer where seismic waves start to propagate with a speed above 1.0 km/sec, in meters. :param z2pt5: Vertical distance from earth surface to the layer where seismic waves start to propagate with a speed above 2.5 km/sec, in km. :param backarc": Boolean value, ``True`` if the site is in the subduction backarc and ``False`` if it is in the subduction forearc or is unknown :raises ValueError: If any of ``vs30``, ``z1pt0`` or ``z2pt5`` is zero or negative. .. note:: :class:`Sites <Site>` are pickleable """ _slots_ = 'location vs30 vs30measured z1pt0 z2pt5 backarc'.split() def __init__(self, location, vs30, vs30measured, z1pt0, z2pt5, backarc=False): if not vs30 > 0: raise ValueError('vs30 must be positive') if not z1pt0 > 0: raise ValueError('z1pt0 must be positive') if not z2pt5 > 0: raise ValueError('z2pt5 must be positive') self.location = location self.vs30 = vs30 self.vs30measured = vs30measured self.z1pt0 = z1pt0 self.z2pt5 = z2pt5 self.backarc = backarc def __str__(self): """ >>> import openquake.hazardlib >>> loc = openquake.hazardlib.geo.point.Point(1, 2, 3) >>> str(Site(loc, 760.0, True, 100.0, 5.0)) '<Location=<Latitude=2.000000, Longitude=1.000000, Depth=3.0000>, \ Vs30=760.0000, Vs30Measured=True, Depth1.0km=100.0000, Depth2.5km=5.0000, \ Backarc=False>' """ return ( "<Location=%s, Vs30=%.4f, Vs30Measured=%r, Depth1.0km=%.4f, " "Depth2.5km=%.4f, Backarc=%r>") % ( self.location, self.vs30, self.vs30measured, self.z1pt0, self.z2pt5, self.backarc) def __hash__(self): return hash((self.location.x, self.location.y)) def __eq__(self, other): return (self.location.x, self.location.y) == ( other.location.x, other.location.y) def __repr__(self): """ >>> import openquake.hazardlib >>> loc = openquake.hazardlib.geo.point.Point(1, 2, 3) >>> site = Site(loc, 760.0, True, 100.0, 5.0) >>> str(site) == repr(site) True """ return self.__str__()
def _extract(array_or_float, indices): try: # if array return array_or_float[indices] except TypeError: # if float return array_or_float @with_slots
[docs]class SiteCollection(object): """ A collection of :class:`sites <Site>`. Instances of this class are intended to represent a large collection of sites in a most efficient way in terms of memory usage. .. note:: If a :class:`SiteCollection` is created from sites containing only lon and lat, iterating over the collection will yield :class:`Sites <Site>` with a reference depth of 0.0 (the sea level). Otherwise, it is possible to model the sites on a realistic topographic surface by specifying the `depth` of each site. :param sites: A list of instances of :class:`Site` class. """ dtype = numpy.dtype([ ('sids', numpy.uint32), ('lons', numpy.float64), ('lats', numpy.float64), ('depths', numpy.float64), ('_vs30', numpy.float64), ('_vs30measured', numpy.bool), ('_z1pt0', numpy.float64), ('_z2pt5', numpy.float64), ('_backarc', numpy.bool), ]) _slots_ = dtype.names @classmethod
[docs] def from_points(cls, lons, lats, depths, sitemodel): """ Build the site collection from :param lons: a sequence of longitudes :param lats: a sequence of latitudes :param depths: a sequence of depths :param sitemodel: an object containing the attributes reference_vs30_value, reference_vs30_type, reference_depth_to_1pt0km_per_sec, reference_depth_to_2pt5km_per_sec, reference_backarc """ if depths is None: depths = numpy.zeros(len(lons)) assert len(lons) == len(lats) == len(depths), (len(lons), len(lats), len(depths)) self = cls.__new__(cls) self.complete = self self.total_sites = len(lons) self.sids = numpy.arange(len(lons), dtype=numpy.uint32) self.lons = numpy.array(lons) self.lats = numpy.array(lats) self.depths = numpy.array(depths) self._vs30 = sitemodel.reference_vs30_value self._vs30measured = sitemodel.reference_vs30_type == 'measured' self._z1pt0 = sitemodel.reference_depth_to_1pt0km_per_sec self._z2pt5 = sitemodel.reference_depth_to_2pt5km_per_sec self._backarc = sitemodel.reference_backarc return self
def __init__(self, sites): self.complete = self self.total_sites = n = len(sites) self.sids = numpy.zeros(n, dtype=int) self.lons = numpy.zeros(n, dtype=float) self.lats = numpy.zeros(n, dtype=float) self.depths = numpy.zeros(n, dtype=float) self._vs30 = numpy.zeros(n, dtype=float) self._vs30measured = numpy.zeros(n, dtype=bool) self._z1pt0 = numpy.zeros(n, dtype=float) self._z2pt5 = numpy.zeros(n, dtype=float) self._backarc = numpy.zeros(n, dtype=bool) for i in range(n): self.sids[i] = i self.lons[i] = sites[i].location.longitude self.lats[i] = sites[i].location.latitude self.depths[i] = sites[i].location.depth self._vs30[i] = sites[i].vs30 self._vs30measured[i] = sites[i].vs30measured self._z1pt0[i] = sites[i].z1pt0 self._z2pt5[i] = sites[i].z2pt5 self._backarc[i] = sites[i].backarc # protect arrays from being accidentally changed. it is useful # because we pass these arrays directly to a GMPE through # a SiteContext object and if a GMPE is implemented poorly it could # modify the site values, thereby corrupting site and all the # subsequent calculation. note that this doesn't protect arrays from # being changed by calling itemset() for arr in (self._vs30, self._vs30measured, self._z1pt0, self._z2pt5, self.lons, self.lats, self.depths, self._backarc, self.sids): arr.flags.writeable = False def __toh5__(self): array = numpy.zeros(self.total_sites, self.dtype) for slot in self._slots_: array[slot] = getattr(self, slot) attrs = dict(total_sites=self.total_sites) return array, attrs def __fromh5__(self, array, attrs): for slot in self._slots_: setattr(self, slot, array[slot]) vars(self).update(attrs) self.complete = self @property def mesh(self): """Return a mesh with the given lons, lats, and depths""" return Mesh(self.lons, self.lats, self.depths) @property def indices(self): """The full set of indices from 0 to total_sites - 1""" return numpy.arange(0, self.total_sites)
[docs] def at_sea_level(self): """True if all depths are zero""" return (self.depths == 0).all()
[docs] def split_in_tiles(self, hint): """ Split a SiteCollection into a set of tiles (SiteCollection instances). :param hint: hint for how many tiles to generate """ tiles = [] for seq in split_in_blocks(range(len(self)), hint or 1): indices = numpy.array(seq, int) sc = SiteCollection.__new__(SiteCollection) sc.complete = sc sc.total_sites = len(indices) sc.sids = self.sids[indices] sc.lons = self.lons[indices] sc.lats = self.lats[indices] sc.depths = self.depths[indices] sc._vs30 = _extract(self._vs30, indices) sc._vs30measured = _extract(self._vs30measured, indices) sc._z1pt0 = _extract(self._z1pt0, indices) sc._z2pt5 = _extract(self._z2pt5, indices) sc._backarc = _extract(self._backarc, indices) tiles.append(sc) return tiles
def __iter__(self): """ Iterate through all :class:`sites <Site>` in the collection, yielding one at a time. """ if isinstance(self.vs30, float): # from points for i, location in enumerate(self.mesh): yield Site(location, self._vs30, self._vs30measured, self._z1pt0, self._z2pt5, self._backarc) else: # from sites for i, location in enumerate(self.mesh): yield Site(location, self.vs30[i], self.vs30measured[i], self.z1pt0[i], self.z2pt5[i], self.backarc[i])
[docs] def filter(self, mask): """ Create a FilteredSiteCollection with only a subset of sites from this one. :param mask: Numpy array of boolean values of the same length as this sites collection. ``True`` values should indicate that site with that index should be included into the filtered collection. :returns: A new :class:`FilteredSiteCollection` instance, unless all the values in ``mask`` are ``True``, in which case this site collection is returned, or if all the values in ``mask`` are ``False``, in which case method returns ``None``. New collection has data of only those sites that were marked for inclusion in mask. """ assert len(mask) == len(self), (len(mask), len(self)) if mask.all(): # all sites satisfy the filter, return # this collection unchanged return self if not mask.any(): # no sites pass the filter, return None return None # extract indices of Trues from the mask [indices] = mask.nonzero() return FilteredSiteCollection(indices, self)
def __len__(self): """ Return the number of sites in the collection. """ return self.total_sites def __repr__(self): return '<SiteCollection with %d sites>' % self.total_sites
# adding a number of properties for the site model data for name in 'vs30 vs30measured z1pt0 z2pt5 backarc'.split():
[docs] def getarray(sc, name=name): # sc is a SiteCollection value = getattr(sc, '_' + name) if isinstance(value, (float, bool)): arr = numpy.array([value] * len(sc), dtype=type(value)) arr.flags.writeable = False return arr else: return value
setattr(SiteCollection, name, property(getarray, doc='%s array' % name)) @with_slots
[docs]class FilteredSiteCollection(object): """ A class meant to store proper subsets of a complete collection of sites in a memory-efficient way. :param indices: an array of indices referring to the complete site collection :param complete: the complete site collection the filtered collection was derived from Notice that if you filter a FilteredSiteCollection `fsc`, you will get a different FilteredSiteCollection referring to the complete SiteCollection `fsc.complete`, not to the filtered collection `fsc`. """ _slots_ = 'indices complete'.split() def __init__(self, indices, complete): if complete is not complete.complete: raise ValueError( 'You should pass a full site collection, not %s' % complete) self.indices = indices self.complete = complete @property def total_sites(self): """The total number of the original sites, without filtering""" return self.complete.total_sites @property def mesh(self): """Return a mesh with the given lons, lats, and depths""" return Mesh(self.lons, self.lats, self.depths)
[docs] def filter(self, mask): """ Create a FilteredSiteCollection with only a subset of sites from this one. :param mask: Numpy array of boolean values of the same length as this filtered sites collection. ``True`` values should indicate that site with that index should be included into the filtered collection. :returns: A new :class:`FilteredSiteCollection` instance, unless all the values in ``mask`` are ``True``, in which case this site collection is returned, or if all the values in ``mask`` are ``False``, in which case method returns ``None``. New collection has data of only those sites that were marked for inclusion in mask. """ assert len(mask) == len(self), (len(mask), len(self)) if mask.all(): return self elif not mask.any(): return None indices = self.indices.take(mask.nonzero()[0]) return FilteredSiteCollection(indices, self.complete)
def __iter__(self): """ Iterate through all :class:`sites <Site>` in the collection, yielding one at a time. """ for i, location in enumerate(self.mesh): yield Site(location, self.vs30[i], self.vs30measured[i], self.z1pt0[i], self.z2pt5[i], self.backarc[i]) def __len__(self): """Return the number of filtered sites""" return len(self.indices) def __toh5__(self): n = len(self.complete) array = numpy.zeros(n, self.complete.dtype) for slot in self.complete._slots_: array[slot] = getattr(self.complete, slot) attrs = dict(total_sites=n, indices=self.indices) return array, attrs def __fromh5__(self, array, attrs): complete = object.__new__(SiteCollection) complete.complete = complete complete.total_sites = attrs['total_sites'] for slot in complete._slots_: setattr(complete, slot, array[slot]) self.indices = attrs['indices'] self.complete = complete def __repr__(self): return '<FilteredSiteCollection with %d of %d sites>' % ( len(self.indices), self.total_sites)
def _extract_site_param(fsc, name): # extract the site parameter 'name' from the filtered site collection return getattr(fsc.complete, name)[fsc.indices] # attach a number of properties filtering the arrays for name in 'vs30 vs30measured z1pt0 z2pt5 backarc lons lats depths \ sids'.split(): prop = property( lambda fsc, name=name: _extract_site_param(fsc, name), doc='Extract %s array from FilteredSiteCollection' % name) setattr(FilteredSiteCollection, name, prop)