Refactoring of the code that creates in-slab ruptures

openquake/mbi/sub/create_ruptures.py

@@ -9,6 +9,9 @@
 
 
 def main(config_fname: str, only_plt: bool = False):
+    """
+    Creates inslab ruptures
+    """
 
     # Parsing config
     model = toml.load(config_fname)
@@ -24,7 +27,7 @@ def main(config_fname: str, only_plt: bool = False):
         calculate_ruptures(ini_fname, ref_fdr=ref_fdr, agr=agr, bgr=bgr,
                            mmin=mmin, mmax=mmax)
 
-descr = 'The path to the .ini file containing info to build the ruptures'
+descr = 'The path to the .toml file containing info to build the ruptures'
 main.config_fname = descr
 
 if __name__ == '__main__':

openquake/mbi/sub/get_profiles_from_slab2pt0_geojson.py

@@ -0,0 +1,88 @@
+#!/usr/bin/env python
+# ------------------- The OpenQuake Model Building Toolkit --------------------
+# Copyright (C) 2022 GEM Foundation
+#           _______  _______        __   __  _______  _______  ___   _
+#          |       ||       |      |  |_|  ||  _    ||       ||   | | |
+#          |   _   ||   _   | ____ |       || |_|   ||_     _||   |_| |
+#          |  | |  ||  | |  ||____||       ||       |  |   |  |      _|
+#          |  |_|  ||  |_|  |      |       ||  _   |   |   |  |     |_
+#          |       ||      |       | ||_|| || |_|   |  |   |  |    _  |
+#          |_______||____||_|      |_|   |_||_______|  |___|  |___| |_|
+#
+# This program 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.
+#
+# This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+# -----------------------------------------------------------------------------
+# vim: tabstop=4 shiftwidth=4 softtabstop=4
+# coding: utf-8
+
+import toml
+from openquake.baselib import sap
+from openquake.wkf.utils import create_folder
+from openquake.sub.get_profiles_from_slab2pt0 import get_profiles_geojson
+
+
+def main(fname_conf: str):
+    """
+    Given a geojson file and the Slab2.0 depth .grd file, this code creates 
+    a set of profiles describing the surface of the slab.
+    The geojson file should contain cross-sections as line segments, with
+    dipdir and length (in m) columns for each cross-section
+    (see sub/tests/data/izu_slab2_css.geojson for example)
+
+    :param fname_conf:
+        Name of the .toml formatted file. If the configuration file defines
+        the `output_folder` variable, the profiles are saved in that folder
+        for further use. If the configuration file defines the `fname_fig`
+        variable, a figure with the traces is saved.
+
+        Example of .toml file
+        ```
+        fname_geojson ='izu_slab2_css.geojson'
+        fname_dep ='izu_slab2_dep_02.24.18.grd'
+        spacing = 200.0
+        folder_out = './cs'
+        fname_fig = './mar.pdf'
+        ```
+    """
+
+    # Read configuration file
+    conf = toml.load(fname_conf)
+
+    # Name of the .grd file with the depth values
+    fname_dep = conf.get('fname_dep', None)
+
+    # Name of the geojson file
+    fname_geojson = conf.get('fname_geojson', None)
+
+    # set spacing from configuration file
+    spacing = conf.get('spacing', 100.)
+
+    # Name of the folder where to save the profiles
+    folder_out = conf.get('folder_out', None)
+
+    # Name of the figure
+    fname_fig = conf.get('fname_fig', '')
+
+    # Get profiles
+    slb = get_profiles_geojson(fname_geojson, fname_dep, spacing, fname_fig)
+
+    # Save profiles
+    if folder_out is not None:
+        create_folder(folder_out)
+        slb.write_profiles(folder_out)
+
+
+main.fname_conf = 'Name of .toml file with configuration parameters'
+
+if __name__ == '__main__':
+    sap.run(main)

openquake/mbi/wkf/create_nrml_sources.py

@@ -10,12 +10,12 @@
 from glob import glob
 from openquake.wkf.utils import create_folder
 
-import importlib
 from openquake.baselib import sap
 from openquake.hazardlib.sourcewriter import write_source_model
 from openquake.hazardlib.source import PointSource, MultiPointSource
 from openquake.hazardlib.mfd import TruncatedGRMFD
-from openquake.hazardlib.mfd. multi_mfd import MultiMFD
+from openquake.hazardlib.mfd.multi_mfd import MultiMFD
+from openquake.hazardlib.scalerel import get_available_magnitude_scalerel
 
 from shapely.geometry import Point as PointShapely
 from openquake.hazardlib.geo.point import Point
@@ -40,14 +40,42 @@ def _get_hypocenter_distribution(data):
     return PMF(out)
 
 
-def write_as_multipoint_sources(df, model, src_id, module, subzones,
+def write_as_multipoint_sources(df, model, src_id, msr_dict, subzones,
                                 model_subz, mmin, bwid, rms, tom, folder_out):
+    """
+    Write a set of point sources to NRML as a multi-point
+
+    :param df:
+        A dataframe where each row is a point source
+    :param model:
+        A dictionary with the model representation
+    :param src_id:
+        A string with the ID of the source
+    :param msr_dict:
+        A dictionary created with  the `get_available_magnitude_scalerel`
+        function available in OQ Engine
+    :param subzones:
+        Must be false since we do not support this feature
+    :param model_subz:
+        ditto
+    :param mmin:
+        A float defining the minimum magnitude of the newly created source
+    :param bwid:
+        A float defining the width of the magnitude bins for the MFD of the
+        newly created source
+    :param rms:
+        A float specifying the rupture mesh spacing
+    :param tom:
+        An instance of :class:`openquake.hazardlib.tom.BaseTOM` subclasses
+    :param folder_out:
+        The output folder
+    """
 
-    # We do not support subzones in this case
+    # We do not support subzones in this case hence 'subzones' must be False
     assert subzones is False
     srcd = model['sources'][src_id]
 
-    # Prefix
+    # Get the prefix
     pfx = model.get("source_prefix", "")
     pfx += "_" if len(pfx) else pfx
 
@@ -170,10 +198,19 @@ def create_nrml_sources(fname_input_pattern: str, fname_config: str,
                         folder_out: str, as_multipoint: bool,
                         fname_subzone_shp: str = "",
                         fname_subzone_config: str = "",):
+    """
+    :param fname_input_pattern:
+    :param fname_config:
+    :param folder_out:
+    :param as_multipoint:
+    :param fname_subzone_shp:
+    :param fname_subzone_config:
+    """
 
+    # Create the output folder
     create_folder(folder_out)
 
-    # If true we take some of the information from subzones
+    # If `subzones` is true we take some of the information from subzones
     subzones = (len(fname_subzone_shp) > 0 and len(fname_subzone_config) > 0)
     model_subz = None
     if subzones:
@@ -182,6 +219,7 @@ def create_nrml_sources(fname_input_pattern: str, fname_config: str,
 
     # This is used to instantiate the MSR
     module = importlib.import_module('openquake.hazardlib.scalerel')
+    msr_dict = get_available_magnitude_scalerel
 
     # Parsing config
     model = toml.load(fname_config)
@@ -211,11 +249,11 @@ def create_nrml_sources(fname_input_pattern: str, fname_config: str,
             df = gpd.sjoin(gdf, tdf, op='within')
 
         if as_multipoint:
-            write_as_multipoint_sources(df, model, src_id, module, subzones,
+            write_as_multipoint_sources(df, model, src_id, msr_dict, subzones,
                                         model_subz, mmin, bwid, rms, tom,
                                         folder_out)
         else:
-            write_as_set_point_sources(df, model, src_id, module, subzones,
+            write_as_set_point_sources(df, model, src_id, msr_dict, subzones,
                                        model_subz, mmin, bwid, rms, tom,
                                        folder_out)
 

openquake/sub/cross_sections.py

@@ -46,7 +46,7 @@
 
 from openquake.hazardlib.geo.utils import OrthographicProjection
 from scipy.interpolate import LinearNDInterpolator
-
+from scipy.spatial import Delaunay
 from openquake.hmtk.seismicity.selector import CatalogueSelector
 from openquake.hmtk.parsers.catalogue.csv_catalogue_parser import CsvCatalogueParser
 from openquake.hmtk.parsers.catalogue.gcmt_ndk_parser import ParseNDKtoGCMT
@@ -113,11 +113,12 @@ def compute_profiles(self, bffer):
 
             # Get min and max longitude and latitude values
             minlo, maxlo, minla, maxla, qual = cs.get_mm(2.0)
-
+
+            #Sbreakpoint()
             # Find the nodes of the grid within a certain distance from the
             # plane of the cross-section
             if qual == 0:
-                minlo, maxlo, minla, maxla, _ = cs.get_mm(2.0)
+                minlo, maxlo, minla, maxla, _ = cs.get_mm(5.0)
                 idxslb, dsts = cs.get_grd_nodes_within_buffer(
                     pnts[:, 0], pnts[:, 1], bffer, minlo, maxlo, minla, maxla)
             if qual == 1:
@@ -137,11 +138,16 @@ def compute_profiles(self, bffer):
                                    cs.length[0], 0., num)
             p = pnts[idxslb, :]
 
-            try:
+            try: 
                 interp = LinearNDInterpolator(p[:, 0:2], p[:, 2])
                 z = interp(psec[0], psec[1])
+
             except:
-                breakpoint()
+                print("trying altered qhull for interpolation")
+                tri = Delaunay(numpy.c_[(p[:, 0], p[:,1])], qhull_options = "QJ")
+                ip = LinearNDInterpolator(tri, p[:,2])
+                z = ip(psec[0], psec[1])
+
 
             iii = numpy.isfinite(z)
             pro = numpy.concatenate((numpy.expand_dims(psec[0][iii], axis=1),