Merge pull request #395 from GEMScienceTools/slab_geojson

Slab geojson

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,11 @@ def compute_profiles(self, bffer):
 
             # Get min and max longitude and latitude values
             minlo, maxlo, minla, maxla, qual = cs.get_mm(2.0)
-
+            
             # 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 +137,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),

openquake/sub/get_profiles_from_slab2pt0.py

@@ -25,11 +25,13 @@
 # coding: utf-8
 
 
+import pyproj
 import netCDF4
 import numpy as np
+import geopandas as gpd
 from numba import njit
-from openquake.hazardlib.geo.geodetic import (point_at, npoints_towards,
-        geodetic_distance, azimuth)
+from openquake.hazardlib.geo.geodetic import (
+    point_at, npoints_towards, geodetic_distance, azimuth)
 from openquake.sub.cross_sections import CrossSection, Slab2pt0
 
 pygmt_available = False
@@ -281,6 +283,134 @@ def tmp_get_initial_traces(bb, dip_dir, spacing):
     return profiles, distance
 
 
+def get_profiles_geojson(geojson: str, fname_dep: str, spacing: float,
+                         fname_fig: str = ''):
+    """
+    :param fname_str:
+        The name of the Slab2.0 .grd file with the values of strike
+    :param fname_dep:
+        The name of the Slab2.0 .grd file with the values of depth
+    :param spacing:
+        The separation distance between traces
+    :param fname_fig:
+        String specifiying location in which to save output figure
+    """
+
+    # Reading file with depth values
+    f_dep = netCDF4.Dataset(fname_dep)
+    depths = np.array(f_dep.variables['z'])
+    mask = np.where(np.isfinite(depths))
+
+    # Mesh
+    x = np.array(f_dep.variables['x'])
+    y = np.array(f_dep.variables['y'])
+    xx, yy = np.meshgrid(x, y)
+
+    css = []
+    gdf = gpd.read_file(geojson)
+    gdf['coords'] = gdf.geometry.apply(lambda geom: list(geom.coords))
+
+    # Create cross-sections
+    min_lo = 180.0
+    min_la = 90.
+    max_lo = -180.0
+    max_la = -90.0
+    for index, row in gdf.iterrows():
+        coo = np.array(row.coords)
+        min_lo = np.min([min_lo, np.min(coo[:, 0])])
+        min_la = np.min([min_la, np.min(coo[:, 1])])
+        max_lo = np.max([max_lo, np.max(coo[:, 0])])
+        max_la = np.max([max_la, np.max(coo[:, 1])])
+    lon_c = min_lo + (max_lo - min_lo) / 2
+    lat_c = min_la + (max_la - min_la) / 2            
+
+    # Define the forward projection
+    aeqd = pyproj.Proj(proj='aeqd', ellps='WGS84', datum='WGS84',
+                       lat_0=lat_c, lon_0=lon_c).srs
+    gdf_pro = gdf.to_crs(crs=aeqd)
+
+    # Create cross-sections
+    for index, row in gdf.iterrows():
+        # calculate dipdir and lengths from profiles directly
+        azim = azimuth(gdf.coords[index][0][0], gdf.coords[index][0][1], gdf.coords[index][-1][0], gdf.coords[index][-1][1])
+        length = geodetic_distance(gdf.coords[index][0][0], gdf.coords[index][0][1], gdf.coords[index][-1][0], gdf.coords[index][-1][1])
+        cs = CrossSection(gdf.coords[index][0][0], gdf.coords[index][0][1],
+                          length, azim)
+        css.append(cs)
+
+    # Filter
+    depths = depths[mask]
+    xx = xx[mask]
+    yy = yy[mask]
+
+    # Coords
+    tmp = zip(xx.flatten(), yy.flatten(), depths.flatten())
+    depths = [[x, y, z] for x, y, z in tmp]
+    depths = np.array(depths)
+    mask = depths[:, 0] > 180
+    depths[mask, 0] = depths[mask, 0] - 360
+    milo = np.min(depths[:, 0])
+    mila = np.min(depths[:, 1])
+    print(f'Min lon {milo:.2f} Max lon {np.max(depths[:, 0]):.2f}')
+    print(f'Min lat {mila:.2f} Max lat {np.max(depths[:, 1]):.2f}')
+
+    # Slab 2.0
+    slb = Slab2pt0(depths, css)
+    slb.compute_profiles(spacing)
+
+
+    if len(str(fname_fig)) > 0:
+        bb = np.array([min_lo, min_la, max_lo, max_la])
+        dlt = 5.0
+        reg = [bb[0] - dlt, bb[2] + dlt, bb[1] - dlt, bb[3] + dlt]
+        clo = np.mean([bb[0], bb[1]])
+        cla = np.mean([bb[2], bb[3]])
+
+        if pygmt_available:
+            fig = pygmt.Figure()
+            pygmt.makecpt(cmap="jet", series=[0.0, 800])
+            fig.basemap(region=reg, projection=f"T{clo}/{cla}/12c", frame=True)
+            fig.coast(land="gray", water="skyblue")
+
+            fig.plot(x=depths[:, 0], y=depths[:, 1],
+                     color=-depths[:, 2],
+                     style='c0.025c',
+                     cmap=True)
+            # Profiles
+            for i, key in enumerate(slb.profiles):
+                pro = slb.profiles[key]
+                if pro.shape[0] > 0:
+                    fig.plot(x=pro[:, 0],
+                             y=pro[:, 1],
+                             color=pro[:, 2],
+                             cmap=True,
+                             style="h0.025c",
+                             pen='black')
+                    fig.text(x=(pro[0, 0] + 0.3), y=pro[0, 1], 
+                             text=f'{i}', font="4p") 
+            fig.savefig(fname_fig)
+            fig.show()
+        else:
+            from matplotlib import pyplot as plt
+            plt.scatter(depths[:, 0], depths[:, 1], c=-depths[:, 2])
+            for i, pro in enumerate(traces):
+                plt.plot(pro[:, 0], pro[:, 1], 'k')
+                plt.text(pro[0, 0], pro[0, 1], f'{i}')
+            for key in slb.profiles:
+                pro = slb.profiles[key]
+                if pro.shape[0] > 0:
+                    plt.plot(pro[:, 0], pro[:, 1], c='r')
+            if max(reg[0], reg[1]) > 180:
+                xmin = reg[0]-360; xmax = reg[1]-360
+            else:
+                xmin = reg[0]; xmax = reg[1]
+            plt.xlim([xmin, xmax])
+            plt.colorbar(label='depth to slab (km)')
+            plt.savefig(fname_fig)
+
+    return slb
+
+
 def get_profiles(fname_str: str, fname_dep: str, spacing: float, fname_fig:
                  str = ''):
     """
@@ -422,7 +552,7 @@ def get_profiles(fname_str: str, fname_dep: str, spacing: float, fname_fig:
             plt.colorbar(label='depth to slab (km)')
             plt.savefig(fname_fig)
 
-            
+
 
     return slb