Source code for openquake.hmtk.parsers.source_model.nrml04_parser

# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4

#
# LICENSE
#
# Copyright (C) 2015-2023 GEM Foundation
#
# The Hazard Modeller's Toolkit 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.
#
# You should have received a copy of the GNU Affero General Public License
# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>
#
# DISCLAIMER
#
# The software Hazard Modeller's Toolkit (openquake.hmtk) provided herein
# is released as a prototype implementation on behalf of
# scientists and engineers working within the GEM Foundation (Global
# Earthquake Model).
#
# It is distributed for the purpose of open collaboration and in the
# hope that it will be useful to the scientific, engineering, disaster
# risk and software design communities.
#
# The software is NOT distributed as part of GEM's OpenQuake suite
# (https://www.globalquakemodel.org/tools-products) and must be considered as a
# separate entity. The software provided herein is designed and implemented
# by scientific staff. It is not developed to the design standards, nor
# subject to same level of critical review by professional software
# developers, as GEM's OpenQuake software suite.
#
# Feedback and contribution to the software is welcome, and can be
# directed to the hazard scientific staff of the GEM Model Facility
# (hazard@globalquakemodel.org).
#
# The Hazard Modeller's Toolkit (openquake.hmtk) is therefore distributed
# 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.
#
# The GEM Foundation, and the authors of the software, assume no
# liability for use of the software.

"""
Parser for input in a NRML format, with partial validation
"""
import re
from copy import copy
from openquake.baselib.node import node_from_xml
from openquake.baselib.general import deprecated
from openquake.hazardlib.geo.point import Point
from openquake.hazardlib.geo.line import Line
from openquake.hazardlib.geo.polygon import Polygon
from openquake.hazardlib import mfd, valid
from openquake.hazardlib.pmf import PMF
from openquake.hazardlib.geo.nodalplane import NodalPlane
from openquake.hmtk.sources.source_model import mtkSourceModel
from openquake.hmtk.sources.point_source import mtkPointSource
from openquake.hmtk.sources.area_source import mtkAreaSource
from openquake.hmtk.sources.simple_fault_source import mtkSimpleFaultSource
from openquake.hmtk.sources.complex_fault_source import mtkComplexFaultSource
from openquake.hmtk.parsers.source_model.base import BaseSourceModelParser


[docs]def string_(string):
    """
    Returns string or None
    """
    if string:
        return string
    else:
        return None


[docs]def float_(value):
    """
    Returns float of a value, or None
    """
    if value:
        return float(value)
    else:
        return None


[docs]def int_(value):
    """
    Returns int or None
    """
    if value:
        return float(value)
    else:
        return None


[docs]def get_taglist(node):
    """
    Return a list of tags (with NRML namespace removed) representing the
    order of the nodes within a node
    """
    return [
        re.sub(r"\{[^}]*\}", "", copy(subnode.tag)) for subnode in node.nodes
    ]


[docs]def linestring_node_to_line(node, with_depth=False):
    """
    Returns an instance of a Linestring node to :class:
    openquake.hazardlib.geo.line.Line
    """
    assert "LineString" in node.tag
    crds = [float(x) for x in node.nodes[0].text.split()]
    if with_depth:
        return Line(
            [
                Point(crds[iloc], crds[iloc + 1], crds[iloc + 2])
                for iloc in range(0, len(crds), 3)
            ]
        )
    else:
        return Line(
            [
                Point(crds[iloc], crds[iloc + 1])
                for iloc in range(0, len(crds), 2)
            ]
        )


[docs]def node_to_point_geometry(node):
    """
    Reads the node and returns the point geometry, upper depth and lower depth
    """
    assert "pointGeometry" in node.tag
    for subnode in node.nodes:
        if "Point" in subnode.tag:
            # Position
            lon, lat = map(float, subnode.nodes[0].text.split())
            point = Point(lon, lat)
        elif "upperSeismoDepth" in subnode.tag:
            upper_depth = float_(subnode.text)
        elif "lowerSeismoDepth" in subnode.tag:
            lower_depth = float_(subnode.text)
        else:
            # Redundent
            pass
    assert lower_depth > upper_depth
    return point, upper_depth, lower_depth


[docs]def node_to_area_geometry(node):
    """
    Reads an area geometry node and returns the polygon, upper depth and lower
    depth
    """
    assert "areaGeometry" in node.tag
    for subnode in node.nodes:
        if "Polygon" in subnode.tag:
            crds = [
                float(x)
                for x in subnode.nodes[0].nodes[0].nodes[0].text.split()
            ]
            polygon = Polygon(
                [
                    Point(crds[iloc], crds[iloc + 1])
                    for iloc in range(0, len(crds), 2)
                ]
            )
        elif "upperSeismoDepth" in subnode.tag:
            upper_depth = float_(subnode.text)
        elif "lowerSeismoDepth" in subnode.tag:
            lower_depth = float_(subnode.text)
        else:
            # Redundent
            pass
    assert lower_depth > upper_depth
    return polygon, upper_depth, lower_depth


[docs]def node_to_simple_fault_geometry(node):
    """
    Reads a simple fault geometry node and returns an OpenQuake representation

    :returns:
        trace - Trace of fault as instance
    """
    assert "simpleFaultGeometry" in node.tag
    for subnode in node.nodes:
        if "LineString" in subnode.tag:
            trace = linestring_node_to_line(subnode, with_depth=False)
        elif "dip" in subnode.tag:
            dip = float(subnode.text)
        elif "upperSeismoDepth" in subnode.tag:
            upper_depth = float(subnode.text)
        elif "lowerSeismoDepth" in subnode.tag:
            lower_depth = float(subnode.text)
        else:
            # Redundent
            pass
    assert lower_depth > upper_depth
    return trace, dip, upper_depth, lower_depth


[docs]def node_to_complex_fault_geometry(node):
    """
    Reads a complex fault geometry node and returns an
    """
    assert "complexFaultGeometry" in node.tag
    intermediate_edges = []
    for subnode in node.nodes:
        if "faultTopEdge" in subnode.tag:
            top_edge = linestring_node_to_line(
                subnode.nodes[0], with_depth=True
            )
        elif "intermediateEdge" in subnode.tag:
            int_edge = linestring_node_to_line(
                subnode.nodes[0], with_depth=True
            )
            intermediate_edges.append(int_edge)
        elif "faultBottomEdge" in subnode.tag:
            bottom_edge = linestring_node_to_line(
                subnode.nodes[0], with_depth=True
            )
        else:
            # Redundent
            pass
    return [top_edge] + intermediate_edges + [bottom_edge]


[docs]def node_to_scalerel(node):
    """
    Parses a node to an instance of a supported scaling relation class
    """
    if not node.text:
        return None
    return valid.mag_scale_rel(node.text)


[docs]def node_to_truncated_gr(node, bin_width=0.1):
    """
    Parses truncated GR node to an instance of the
    :class: openquake.hazardlib.mfd.truncated_gr.TruncatedGRMFD
    """
    # Parse to float dictionary
    if not all(
        [node.attrib[key] for key in ["minMag", "maxMag", "aValue", "bValue"]]
    ):
        return None
    tgr = dict((key, float_(node.attrib[key])) for key in node.attrib)
    return mfd.truncated_gr.TruncatedGRMFD(
        min_mag=tgr["minMag"],
        max_mag=tgr["maxMag"],
        bin_width=bin_width,
        a_val=tgr["aValue"],
        b_val=tgr["bValue"],
    )


[docs]def node_to_evenly_discretized(node):
    """
    Parses the evenly discretized mfd node to an instance of the
    :class: openquake.hazardlib.mfd.evenly_discretized.EvenlyDiscretizedMFD,
    or to None if not all parameters are available
    """
    if not all(
        [node.attrib["minMag"], node.attrib["binWidth"], node.nodes[0].text]
    ):
        return None
    # Text to float
    rates = [float(x) for x in node.nodes[0].text.split()]
    return mfd.evenly_discretized.EvenlyDiscretizedMFD(
        float(node.attrib["minMag"]), float(node.attrib["binWidth"]), rates
    )


[docs]def node_to_mfd(node, taglist):
    """
    Reads the node to return a magnitude frequency distribution
    """
    if "incrementalMFD" in taglist:
        mfd = node_to_evenly_discretized(
            node.nodes[taglist.index("incrementalMFD")]
        )
    elif "truncGutenbergRichterMFD" in taglist:
        mfd = node_to_truncated_gr(
            node.nodes[taglist.index("truncGutenbergRichterMFD")]
        )
    else:
        mfd = None
    return mfd


[docs]def node_to_nodal_planes(node):
    """
    Parses the nodal plane distribution to a PMF
    """
    if not len(node):
        return None
    npd_pmf = []
    for plane in node.nodes:
        if not all(plane.attrib[key] for key in plane.attrib):
            # One plane fails - return None
            return None
        npd = NodalPlane(
            float(plane.attrib["strike"]),
            float(plane.attrib["dip"]),
            float(plane.attrib["rake"]),
        )
        npd_pmf.append((float(plane.attrib["probability"]), npd))
    return PMF(npd_pmf)


[docs]def node_to_hdd(node):
    """
    Parses the node to a hpyocentral depth distribution PMF
    """
    if not len(node):
        return None
    hdds = []
    for subnode in node.nodes:
        if not all([subnode.attrib[key] for key in ["depth", "probability"]]):
            return None
        hdds.append(
            (
                float(subnode.attrib["probability"]),
                float(subnode.attrib["depth"]),
            )
        )
    return PMF(hdds)


[docs]def parse_point_source_node(node, mfd_spacing=0.1):
    """
    Returns an "areaSource" node into an instance of the :class:
    openquake.hmtk.sources.area.mtkAreaSource
    """
    assert "pointSource" in node.tag
    pnt_taglist = get_taglist(node)
    # Get metadata
    point_id, name, trt = (
        node.attrib["id"],
        node.attrib["name"],
        node.attrib["tectonicRegion"],
    )
    assert point_id  # Defensive validation!
    # Process geometry
    location, upper_depth, lower_depth = node_to_point_geometry(
        node.nodes[pnt_taglist.index("pointGeometry")]
    )
    # Process scaling relation
    msr = node_to_scalerel(node.nodes[pnt_taglist.index("magScaleRel")])
    # Process aspect ratio
    aspect = float_(node.nodes[pnt_taglist.index("ruptAspectRatio")].text)
    # Process MFD
    mfd = node_to_mfd(node, pnt_taglist)
    # Process nodal planes
    npds = node_to_nodal_planes(
        node.nodes[pnt_taglist.index("nodalPlaneDist")]
    )
    # Process hypocentral depths
    hdds = node_to_hdd(node.nodes[pnt_taglist.index("hypoDepthDist")])
    return mtkPointSource(
        point_id,
        name,
        trt,
        geometry=location,
        upper_depth=upper_depth,
        lower_depth=lower_depth,
        mag_scale_rel=msr,
        rupt_aspect_ratio=aspect,
        mfd=mfd,
        nodal_plane_dist=npds,
        hypo_depth_dist=hdds,
    )


[docs]def parse_area_source_node(node, mfd_spacing=0.1):
    """
    Returns an "areaSource" node into an instance of the :class:
    openquake.hmtk.sources.area.mtkAreaSource
    """
    assert "areaSource" in node.tag
    area_taglist = get_taglist(node)
    # Get metadata
    area_id, name, trt = (
        node.attrib["id"],
        node.attrib["name"],
        node.attrib["tectonicRegion"],
    )
    assert area_id  # Defensive validation!
    # Process geometry
    polygon, upper_depth, lower_depth = node_to_area_geometry(
        node.nodes[area_taglist.index("areaGeometry")]
    )
    # Process scaling relation
    msr = node_to_scalerel(node.nodes[area_taglist.index("magScaleRel")])
    # Process aspect ratio
    aspect = float_(node.nodes[area_taglist.index("ruptAspectRatio")].text)
    # Process MFD
    mfd = node_to_mfd(node, area_taglist)
    # Process nodal planes
    npds = node_to_nodal_planes(
        node.nodes[area_taglist.index("nodalPlaneDist")]
    )
    # Process hypocentral depths
    hdds = node_to_hdd(node.nodes[area_taglist.index("hypoDepthDist")])
    return mtkAreaSource(
        area_id,
        name,
        trt,
        geometry=polygon,
        upper_depth=upper_depth,
        lower_depth=lower_depth,
        mag_scale_rel=msr,
        rupt_aspect_ratio=aspect,
        mfd=mfd,
        nodal_plane_dist=npds,
        hypo_depth_dist=hdds,
    )


[docs]def parse_simple_fault_node(node, mfd_spacing=0.1, mesh_spacing=1.0):
    """
    Parses a "simpleFaultSource" node and returns an instance of the :class:
    openquake.hmtk.sources.simple_fault.mtkSimpleFaultSource
    """
    assert "simpleFaultSource" in node.tag
    sf_taglist = get_taglist(node)
    # Get metadata
    sf_id, name, trt = (
        node.attrib["id"],
        node.attrib["name"],
        node.attrib["tectonicRegion"],
    )
    # Process geometry
    trace, dip, upper_depth, lower_depth = node_to_simple_fault_geometry(
        node.nodes[sf_taglist.index("simpleFaultGeometry")]
    )
    # Process scaling relation
    msr = node_to_scalerel(node.nodes[sf_taglist.index("magScaleRel")])
    # Process aspect ratio
    aspect = float_(node.nodes[sf_taglist.index("ruptAspectRatio")].text)
    # Process MFD
    mfd = node_to_mfd(node, sf_taglist)
    # Process rake
    rake = float_(node.nodes[sf_taglist.index("rake")].text)
    simple_fault = mtkSimpleFaultSource(
        sf_id,
        name,
        trt,
        geometry=None,
        dip=dip,
        upper_depth=upper_depth,
        lower_depth=lower_depth,
        mag_scale_rel=msr,
        rupt_aspect_ratio=aspect,
        mfd=mfd,
        rake=rake,
    )
    simple_fault.create_geometry(
        trace, dip, upper_depth, lower_depth, mesh_spacing
    )
    return simple_fault


[docs]def parse_complex_fault_node(node, mfd_spacing=0.1, mesh_spacing=4.0):
    """
    Parses a "complexFaultSource" node and returns an instance of the :class:
    openquake.hmtk.sources.complex_fault.mtkComplexFaultSource
    """
    assert "complexFaultSource" in node.tag
    sf_taglist = get_taglist(node)
    # Get metadata
    sf_id, name, trt = (
        node.attrib["id"],
        node.attrib["name"],
        node.attrib["tectonicRegion"],
    )
    # Process geometry
    edges = node_to_complex_fault_geometry(
        node.nodes[sf_taglist.index("complexFaultGeometry")]
    )
    # Process scaling relation
    msr = node_to_scalerel(node.nodes[sf_taglist.index("magScaleRel")])
    # Process aspect ratio
    aspect = float_(node.nodes[sf_taglist.index("ruptAspectRatio")].text)
    # Process MFD
    mfd = node_to_mfd(node, sf_taglist)
    # Process rake
    rake = float_(node.nodes[sf_taglist.index("rake")].text)
    complex_fault = mtkComplexFaultSource(
        sf_id,
        name,
        trt,
        geometry=None,
        mag_scale_rel=msr,
        rupt_aspect_ratio=aspect,
        mfd=mfd,
        rake=rake,
    )
    complex_fault.create_geometry(edges, mesh_spacing)
    return complex_fault


[docs]class nrmlSourceModelParser(BaseSourceModelParser):
    """
    Parser for a source model in NRML format, permitting partial validation
    such that not all fields need to be specified for the file to be parsed
    """

[docs]    @deprecated(msg="Use openquake.hazardlib.nrml.to_python instead")
    def read_file(
        self,
        identifier,
        mfd_spacing=0.1,
        simple_mesh_spacing=1.0,
        complex_mesh_spacing=4.0,
        area_discretization=10.0,
    ):
        """
        Reads in the source model in returns an instance of the :class:
        openquake.hmtk.sourcs.source_model.mtkSourceModel
        """
        sm_node = node_from_xml(self.input_file)[0]
        nrml04 = "xmlns/nrml/0.4" in sm_node[0].tag
        if nrml04:
            node_sets = [sm_node]
            sm_name = sm_node.get("name", "")
        else:  # format NRML 0.5+
            node_sets = sm_node
            sm_name = sm_node["name"]
        source_model = mtkSourceModel(identifier, name=sm_name)
        for node_set in node_sets:
            for node in node_set:
                if not nrml04:  # get the TRT from the sourceGroup
                    node.attrib["tectonicRegion"] = node_set["tectonicRegion"]
                if "pointSource" in node.tag:
                    source_model.sources.append(
                        parse_point_source_node(node, mfd_spacing)
                    )
                elif "areaSource" in node.tag:
                    source_model.sources.append(
                        parse_area_source_node(node, mfd_spacing)
                    )
                elif "simpleFaultSource" in node.tag:
                    source_model.sources.append(
                        parse_simple_fault_node(
                            node, mfd_spacing, simple_mesh_spacing
                        )
                    )
                elif "complexFaultSource" in node.tag:
                    source_model.sources.append(
                        parse_complex_fault_node(
                            node, mfd_spacing, complex_mesh_spacing
                        )
                    )
                # TODO: multiPointSource are not supported
                else:
                    print(
                        "Source typology %s not recognised - skipping!"
                        % node.tag
                    )
        return source_model