curves.py

import numpy

from miniquake import *


SITES = sites(boundary=[(-111, 29), (-127, 29), (-127, 43.5), (-111, 43.5)],
              discretization=10, vs30=760, vs30measured=False, z1pt0=100,
              z2pt5=5)

SOURCES = sources('/tmp/NRML04/CaliforniaGriddedSeismicity.xml',
                  mesh_spacing=2, bin_width=1, area_src_disc=30)

IMTS = {PGA(): list(10 ** (1 + numpy.arange(30.) / 300) - 10 + 0.000001)}

TIME_SPAN = 100

GSIMS = dict((getattr(TRT, trt), ChiouYoungs2008())
             for trt in vars(TRT) if not trt.startswith('_'))
GSIMS['Stable Continental Crust'] = ChiouYoungs2008()

TRUNCATION_LEVEL = 3

INTEGRATION_DISTANCE = 200

CURVES = calc(SOURCES, SITES, IMTS, TIME_SPAN, GSIMS, TRUNCATION_LEVEL,
              INTEGRATION_DISTANCE)


for imt in CURVES:
    numpy.save('%s.npy' % type(imt).__name__, CURVES[imt])

miniquake.py

import multiprocessing
import time

import numpy
from shapely import wkt

from nrml import parsers as nrml_parsers
from nrml import models as nrml_models

from nhlib import geo
from nhlib import mfd
from nhlib import pmf
from nhlib import scalerel
from nhlib import source
from nhlib.geo import Point, Polygon
from nhlib.site import SiteCollection
from nhlib.calc import hazard_curves_poissonian
from nhlib.calc.filters import rupture_site_distance_filter, \
                               source_site_distance_filter

from nhlib.gsim.sadigh_1997 import SadighEtAl1997
from nhlib.gsim.chiou_youngs_2008 import ChiouYoungs2008
from nhlib.imt import *
from nhlib.const import TRT


__all__ = ('SadighEtAl1997', 'ChiouYoungs2008', 'sites', 'calc', 'TRT',
           'sources', 'PGA', 'PGD', 'PGV', 'SA', 'IA', 'RSD', 'MMI')


def sites(boundary, discretization, vs30, vs30measured, z1pt0, z2pt5):
    boundary = [Point(*coords) for coords in boundary]
    sites_mesh = Polygon(boundary).discretize(discretization)
    sitecol = object.__new__(SiteCollection)
    sitecol.indices = None
    sitecol.mesh = sites_mesh
    sitecol.vs30 = numpy.repeat(float(vs30), len(sites_mesh))
    sitecol.vs30measured = numpy.repeat(bool(vs30measured), len(sites_mesh))
    sitecol.z1pt0 = numpy.repeat(float(z1pt0), len(sites_mesh))
    sitecol.z2pt5 = numpy.repeat(float(z2pt5), len(sites_mesh))
    return sitecol


def sources(model, mesh_spacing, bin_width, area_src_disc):
    nrml_sources = nrml_parsers.SourceModelParser(model).parse()
    for nrml_source in nrml_sources:
        nhlib_source = nrml_to_nhlib(nrml_source, mesh_spacing, bin_width,
                                     area_src_disc)
        yield nhlib_source


def calc(sources, sites, imts, time_span, gsims, truncation_level,
         integration_distance=None):
    kwargs = {'sites': sites, 'imts': imts, 'time_span': time_span,
              'gsims': gsims, 'truncation_level': truncation_level}
    if integration_distance is not None:
        kwargs['source_site_filter'] = source_site_distance_filter(
            integration_distance
        )
        kwargs['rupture_site_filter'] = rupture_site_distance_filter(
            integration_distance
        )

    print 'cntrl: considering %s sites' % len(sites.vs30)

    num_procs = multiprocessing.cpu_count()
    source_queue = multiprocessing.Queue(maxsize=num_procs * 4)
    result_queue = multiprocessing.Queue()
    loglock = multiprocessing.Lock()
    processes = []
    for i in xrange(num_procs):
        proc = _CalcProcess(kwargs, source_queue, result_queue, loglock)
        proc.start()
        processes.append(proc)

    for i, nhlib_source in enumerate(sources):
        source_queue.put(nhlib_source)
        with loglock:
            print 'cntrl: sent %s sources to the queue' % i

    for i in xrange(num_procs):
        with loglock:
            print 'cntrl: no more sources, sending exit command'
        source_queue.put(None)

    res = result_queue.get()
    for i in xrange(1, num_procs):
        with loglock:
            print 'cntrl: got results from %d sources' % i
        res = _merge_poes(res, result_queue.get())

    return res


class _CalcProcess(multiprocessing.Process):
    def __init__(self, kwargs, source_queue, result_queue, loglock):
        super(_CalcProcess, self).__init__()
        self.kwargs = kwargs
        self.source_queue = source_queue
        self.result_queue = result_queue
        self.loglock = loglock

    def run(self):
        curves = dict((imt, numpy.zeros([len(self.kwargs['sites'].vs30),
                                        len(self.kwargs['imts'][imt])]))
                      for imt in self.kwargs['imts'])
        sources_processed = 0
        self._log('ready')

        while True:
            source = self.source_queue.get()
            if source is None:
                self._log('got exit command')
                break
            self._log('got source %r' % source.name)
            start = time.time()
            source_curves = hazard_curves_poissonian(sources=[source],
                                                     **self.kwargs)
            stop = time.time()
            curves = _merge_poes(curves, source_curves)
            sources_processed += 1
            self._log('processed source %r in %.3f seconds' % (source.name,
                                                               stop - start))
            self._log('processed %d sources so far' % sources_processed)

        self.result_queue.put(curves)

    def _log(self, msg):
        with self.loglock:
            print '%r: %s' % (self.ident, msg)


def _merge_poes(curves1, curves2):
    return dict(
        (imt, 1 - (1 - curves1[imt]) * (1 - curves2[imt]))
        for imt in curves1
    )


# Copyright (c) 2010-2012, 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 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/>.


_SCALE_REL_MAP = {
    'PeerMSR': scalerel.PeerMSR,
    'WC1994': scalerel.WC1994,
}


def nrml_to_nhlib(src, mesh_spacing, bin_width, area_src_disc):
    """Convert a seismic source object from the NRML representation to the
    NHLib representation. Inputs can be point, area, simple fault, or complex
    fault sources.

    See :mod:`nrml.models` and :mod:`nhlib.source`.

    :param src:
        :mod:`nrml.models` seismic source instance.
    :param float mesh_spacing:
        Rupture mesh spacing, in km.
    :param float bin_width:
        Truncated Gutenberg-Richter MFD (Magnitude Frequency Distribution) bin
        width.
    :param float area_src_disc:
        Area source discretization, in km. Applies only to area sources.
        If the input source is known to be a type other than an area source,
        you can specify `area_src_disc=None`.
    :returns:
        The NHLib representation of the input source.
    """
    # The ordering of the switch here matters because:
    #   - AreaSource inherits from PointSource
    #   - ComplexFaultSource inherits from SimpleFaultSource
    if isinstance(src, nrml_models.AreaSource):
        return _area_to_nhlib(src, mesh_spacing, bin_width, area_src_disc)
    elif isinstance(src, nrml_models.PointSource):
        return _point_to_nhlib(src, mesh_spacing, bin_width)
    elif isinstance(src, nrml_models.ComplexFaultSource):
        return _complex_to_nhlib(src, mesh_spacing, bin_width)
    elif isinstance(src, nrml_models.SimpleFaultSource):
        return _simple_to_nhlib(src, mesh_spacing, bin_width)


def _point_to_nhlib(src, mesh_spacing, bin_width):
    """Convert a NRML point source to the NHLib equivalent.

    See :mod:`nrml.models` and :mod:`nhlib.source`.

    :param src:
        :class:`nrml.models.PointSource` instance.
    :param float mesh_spacing:
        Rupture mesh spacing, in km.
    :param float bin_width:
        Truncated Gutenberg-Richter MFD (Magnitude Frequency Distribution) bin
        width.
    :returns:
        The NHLib representation of the input source.
    """
    shapely_pt = wkt.loads(src.geometry.wkt)

    mf_dist = _mfd_to_nhlib(src.mfd, bin_width)

    # nodal plane distribution:
    npd = pmf.PMF(
        [(x.probability,
          geo.NodalPlane(strike=x.strike, dip=x.dip, rake=x.rake))
         for x in src.nodal_plane_dist]
    )

    # hypocentral depth distribution:
    hd = pmf.PMF([(x.probability, x.depth) for x in src.hypo_depth_dist])

    point = source.PointSource(
        source_id=src.id,
        name=src.name,
        tectonic_region_type=src.trt,
        mfd=mf_dist,
        rupture_mesh_spacing=mesh_spacing,
        magnitude_scaling_relationship=_SCALE_REL_MAP[src.mag_scale_rel](),
        rupture_aspect_ratio=src.rupt_aspect_ratio,
        upper_seismogenic_depth=src.geometry.upper_seismo_depth,
        lower_seismogenic_depth=src.geometry.lower_seismo_depth,
        location=geo.Point(shapely_pt.x, shapely_pt.y),
        nodal_plane_distribution=npd,
        hypocenter_distribution=hd
    )

    return point


def _area_to_nhlib(src, mesh_spacing, bin_width, area_src_disc):
    """Convert a NRML area source to the NHLib equivalent.

    See :mod:`nrml.models` and :mod:`nhlib.source`.

    :param src:
        :class:`nrml.models.PointSource` instance.
    :param float mesh_spacing:
        Rupture mesh spacing, in km.
    :param float bin_width:
        Truncated Gutenberg-Richter MFD (Magnitude Frequency Distribution) bin
        width.
    :param float area_src_disc:
        Area source discretization, in km. Applies only to area sources.
    :returns:
        The NHLib representation of the input source.
    """
    shapely_polygon = wkt.loads(src.geometry.wkt)
    nhlib_polygon = geo.Polygon(
        # We ignore the last coordinate in the sequence here, since it is a
        # duplicate of the first. nhlib will close the loop for us.
        [geo.Point(*x) for x in list(shapely_polygon.exterior.coords)[:-1]]
    )

    mf_dist = _mfd_to_nhlib(src.mfd, bin_width)

    # nodal plane distribution:
    npd = pmf.PMF(
        [(x.probability,
          geo.NodalPlane(strike=x.strike, dip=x.dip, rake=x.rake))
         for x in src.nodal_plane_dist]
    )

    # hypocentral depth distribution:
    hd = pmf.PMF([(x.probability, x.depth) for x in src.hypo_depth_dist])

    area = source.AreaSource(
        source_id=src.id,
        name=src.name,
        tectonic_region_type=src.trt,
        mfd=mf_dist,
        rupture_mesh_spacing=mesh_spacing,
        magnitude_scaling_relationship=_SCALE_REL_MAP[src.mag_scale_rel](),
        rupture_aspect_ratio=src.rupt_aspect_ratio,
        upper_seismogenic_depth=src.geometry.upper_seismo_depth,
        lower_seismogenic_depth=src.geometry.lower_seismo_depth,
        nodal_plane_distribution=npd, hypocenter_distribution=hd,
        polygon=nhlib_polygon,
        area_discretization=area_src_disc
    )

    return area


def _simple_to_nhlib(src, mesh_spacing, bin_width):
    """Convert a NRML simple fault source to the NHLib equivalent.

    See :mod:`nrml.models` and :mod:`nhlib.source`.

    :param src:
        :class:`nrml.models.PointSource` instance.
    :param float mesh_spacing:
        Rupture mesh spacing, in km.
    :param float bin_width:
        Truncated Gutenberg-Richter MFD (Magnitude Frequency Distribution) bin
        width.
    :returns:
        The NHLib representation of the input source.
    """
    shapely_line = wkt.loads(src.geometry.wkt)
    fault_trace = geo.Line([geo.Point(*x) for x in shapely_line.coords])

    mf_dist = _mfd_to_nhlib(src.mfd, bin_width)

    simple = source.SimpleFaultSource(
        source_id=src.id,
        name=src.name,
        tectonic_region_type=src.trt,
        mfd=mf_dist,
        rupture_mesh_spacing=mesh_spacing,
        magnitude_scaling_relationship=_SCALE_REL_MAP[src.mag_scale_rel](),
        rupture_aspect_ratio=src.rupt_aspect_ratio,
        upper_seismogenic_depth=src.geometry.upper_seismo_depth,
        lower_seismogenic_depth=src.geometry.lower_seismo_depth,
        fault_trace=fault_trace,
        dip=src.geometry.dip,
        rake=src.rake
    )

    return simple


def _complex_to_nhlib(src, mesh_spacing, bin_width):
    """Convert a NRML complex fault source to the NHLib equivalent.

    See :mod:`nrml.models` and :mod:`nhlib.source`.

    :param src:
        :class:`nrml.models.PointSource` instance.
    :param float mesh_spacing:
        Rupture mesh spacing, in km.
    :param float bin_width:
        Truncated Gutenberg-Richter MFD (Magnitude Frequency Distribution) bin
        width.
    :returns:
        The NHLib representation of the input source.
    """
    edges_wkt = []
    edges_wkt.append(src.geometry.top_edge_wkt)
    edges_wkt.extend(src.geometry.int_edges)
    edges_wkt.append(src.geometry.bottom_edge_wkt)

    edges = []

    for edge in edges_wkt:
        shapely_line = wkt.loads(edge)
        line = geo.Line([geo.Point(*x) for x in shapely_line.coords])
        edges.append(line)

    mf_dist = _mfd_to_nhlib(src.mfd, bin_width)

    cmplx = source.ComplexFaultSource(
        source_id=src.id,
        name=src.name,
        tectonic_region_type=src.trt,
        mfd=mf_dist,
        rupture_mesh_spacing=mesh_spacing,
        magnitude_scaling_relationship=_SCALE_REL_MAP[src.mag_scale_rel](),
        rupture_aspect_ratio=src.rupt_aspect_ratio,
        edges=edges,
        rake=src.rake,
    )

    return cmplx


def _mfd_to_nhlib(src_mfd, bin_width):
    """Convert a NRML MFD to an NHLib MFD.

    :param src_mfd:
        :class:`nrml.models.IncrementalMFD` or :class:`nrml.models.TGRMFD`
        instance.
    :param float bin_width:
        Optional. Required only for Truncated Gutenberg-Richter MFDs.
    :returns:
        The NHLib representation of the MFD. See :mod:`nhlib.mfd`.
    """
    if isinstance(src_mfd, nrml_models.TGRMFD):
        assert bin_width is not None
        return mfd.TruncatedGRMFD(
            a_val=src_mfd.a_val, b_val=src_mfd.b_val, min_mag=src_mfd.min_mag,
            max_mag=src_mfd.max_mag, bin_width=bin_width
        )
    elif isinstance(src_mfd, nrml_models.IncrementalMFD):
        return mfd.EvenlyDiscretizedMFD(
            min_mag=src_mfd.min_mag, bin_width=src_mfd.bin_width,
            occurrence_rates=src_mfd.occur_rates
        )