Source code for openquake.hazardlib.site

# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2012-2016 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
from openquake.hazardlib.geo.utils import cross_idl


@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:: Because calculations assume that :class:`Sites <Site>` are on the Earth's surface, all `depth` information in a :class:`SiteCollection` is discarded. The collection `mesh` will only contain lon and lat. So even if a :class:`SiteCollection` is created from sites containing `depth` in their geometry, iterating over the collection will yield :class:`Sites <Site>` with a reference depth of 0.0. :param sites: A list of instances of :class:`Site` class. """ dtype = numpy.dtype([ ('sids', numpy.uint32), ('lons', numpy.float64), ('lats', 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, sitemodel): """ Build the site collection from :param lons: a sequence of longitudes :param lats: a sequence of latitudes :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 """ assert len(lons) == len(lats), (len(lons), len(lats)) 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._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._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._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._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 and lats""" return Mesh(self.lons, self.lats, depths=None) @property def indices(self): """The full set of indices from 0 to total_sites - 1""" return numpy.arange(0, self.total_sites)
[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._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. See also :meth:`expand`. """ 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)
[docs] def expand(self, data, placeholder): """ For non-filtered site collections just checks that data has the right number of elements and returns it. It is here just for API compatibility with filtered site collections. """ assert len(data) == len(self), (len(data), len(self)) return data
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(): 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 and lats""" return Mesh(self.lons, self.lats, depths=None)
[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. See also :meth:`expand`. """ 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)
[docs] def expand(self, data, placeholder): """ Expand a short array `data` over a filtered site collection of the same length and return a long array of size `total_sites` filled with the placeholder. The typical workflow is the following: there is a whole site collection, the one that has an information about all the sites. Then it gets filtered for performing some calculation on a limited set of sites (like for instance filtering sites by their proximity to a rupture). That filtering process can be repeated arbitrary number of times, i.e. a collection that is already filtered can be filtered for further limiting the set of sites to compute on. Then the (supposedly expensive) computation is done on a limited set of sites which still appears as just a :class:`SiteCollection` instance, so that computation code doesn't need to worry about filtering, it just needs to handle site collection objects. The calculation result comes in a form of 1d or 2d numpy array (that is, either one value per site or one 1d array per site) with length equal to number of sites in a filtered collection. That result needs to be expanded to an array of similar structure but the one that holds values for all the sites in the original (unfiltered) collection. This is what :meth:`expand` is for. It creates a result array of ``total_sites`` length and puts values from ``data`` into appropriate places in it remembering indices of sites that were chosen for actual calculation and leaving ``placeholder`` value everywhere else. :param data: 1d or 2d numpy array with first dimension representing values computed for site from this collection. :param placeholder: A scalar value to be put in result array for those sites that were filtered out and no real calculation was performed for them. :returns: Array of length ``total_sites`` with values from ``data`` distributed in the appropriate places. """ len_data = data.shape[0] assert len_data == len(self), (len_data, len(self)) assert len_data <= self.total_sites assert self.indices[-1] < self.total_sites, ( self.indices[-1], self.total_sites) if data.ndim == 1: # single-dimensional array result = numpy.empty(self.total_sites) result.fill(placeholder) result.put(self.indices, data) return result assert data.ndim == 2 # two-dimensional array num_values = data.shape[1] result = numpy.empty((self.total_sites, num_values)) result.fill(placeholder) for i in range(num_values): result[:, i].put(self.indices, data[:, i]) return result
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 __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 sids'.split(): prop = property( lambda fsc, name=name: _extract_site_param(fsc, name), doc='Extract %s array from FilteredSiteCollection' % name) setattr(FilteredSiteCollection, name, prop)