Merge pull request #456 from GEMScienceTools/add_fault_geojsons

Add ability to create geojsons of the fault sources within a given source model + unit test

openquake/man/checks/plotting.py

@@ -1,6 +1,21 @@
+import os
 import numpy as np
+import pandas as pd
+import geopandas as gpd
 import matplotlib.pyplot as plt
+from glob import glob
+from shapely import Polygon as ShapelyPolygon, LineString
 
+from openquake.commonlib import readinput
+from openquake.hazardlib.nrml import GeometryModel
+from openquake.hazardlib.geo.mesh import RectangularMesh
+from openquake.hazardlib.geo.surface.simple_fault import SimpleFaultSurface
+from openquake.hazardlib.source.multi_fault import MultiFaultSource
+from openquake.hazardlib.source.complex_fault import ComplexFaultSource
+from openquake.hazardlib.source.simple_fault import SimpleFaultSource
+from openquake.hazardlib.source.kite_fault import KiteFaultSource
+from openquake.hazardlib.source.characteristic import CharacteristicFaultSource
+from openquake.hazardlib.source.complex_fault import ComplexFaultSurface
 from openquake.hazardlib.source.point import PointSource
 
 
@@ -58,3 +73,256 @@ def plot_polygon(poly):
 
 def plot_end():
     plt.show()
+
+
+def get_ssm_files(model_dir):
+    """
+    For a given source model get the xml ssm files
+    """
+    # Get seismic source model XMLs
+    base_path = os.path.join(model_dir, 'in', 'ssm')
+    try:
+        assert os.path.exists(base_path)
+    except:
+        raise ValueError("Admitted model does not abide to required directory "
+                         "hierarchy: model_dir/in/ssm/")
+    files = []
+
+    # Loop through each subfolder and get all XML files
+    for depth in range(5):  # Search up to 5 levels deep
+        fipath = os.path.join(base_path, *['**']*depth)
+        files.extend(glob(fipath + "/*.xml", recursive=True))
+
+    # Remove duplicates from recursive search
+    files = pd.Series(files).drop_duplicates().values
+
+    return files
+
+
+def get_sources(model_dir):
+    """
+    Load the sources in the given model and return them as a list.
+    """            
+    # Get source XMLs
+    files = get_ssm_files(model_dir)
+
+    # Fault classess
+    faults = [MultiFaultSource,
+              SimpleFaultSource,
+              ComplexFaultSource,
+              CharacteristicFaultSource,
+              KiteFaultSource]
+
+    # Read the XMLs for all srcs in the given model
+    ssm = readinput.read_source_models(files, 'tmp.hdf5', investigation_time=1,
+                                       rupture_mesh_spacing=5,
+                                       area_source_discretization=5,
+                                       width_of_mfd_bin=0.1)
+    # Get a list of the source objects
+    srcs = []
+    geom_model = []
+    for idx, ssm_obj in enumerate(ssm):
+        if isinstance(ssm_obj, GeometryModel):
+            geom_model.append(ssm_obj)
+            ssm = pd.Series(ssm).drop(idx).values
+        else:
+            for idx_sg, sg in enumerate(ssm_obj.src_groups):
+                for src in ssm_obj[idx_sg]:
+                    if type(src) in faults: # Only retain the fault sources
+                        srcs.append(src)
+    assert len(geom_model) < 2 # Should be one or none
+    if len(geom_model) == 1:
+        geom_model = geom_model[0]
+    else:
+        geom_model = None
+
+    return srcs, geom_model
+
+
+def get_boundary_2d(smsh):
+    """
+    Returns a polygon for the given surface.
+    """
+    coo = []
+
+    # Upper boundary + trace
+    idx = np.where(np.isfinite(smsh.lons[0, :]))[0]
+    tmp = [(smsh.lons[0, i], smsh.lats[0, i]) for i in idx]
+    trace = LineString(tmp)
+    coo.extend(tmp)
+
+    # Right boundary
+    idx = np.where(np.isfinite(smsh.lons[:, -1]))[0]
+    tmp = [(smsh.lons[i, -1], smsh.lats[i, -1]) for i in idx]
+    coo.extend(tmp)
+
+    # Lower boundary
+    idx = np.where(np.isfinite(smsh.lons[-1, :]))[0]
+    tmp = [(smsh.lons[-1, i], smsh.lats[-1, i]) for i in np.flip(idx)]
+    coo.extend(tmp)
+
+    # Left boundary
+    idx = idx = np.where(np.isfinite(smsh.lons[:, 0]))[0]
+    tmp = [(smsh.lons[i, 0], smsh.lats[i, 0]) for i in np.flip(idx)]
+    coo.extend(tmp)
+
+    return trace, ShapelyPolygon(coo)
+
+
+def get_complex_mesh(src):
+    """
+    Get the mesh of a ComplexFaultSource
+    """
+    # Get the surface
+    sfc = ComplexFaultSurface.from_fault_data(src.edges, mesh_spacing=5)
+
+    # Get mesh
+    mesh = RectangularMesh(sfc.mesh.lons, sfc.mesh.lats, sfc.mesh.depths)
+
+    return mesh
+
+
+def get_characteristic_mesh(src):
+    """
+    Get the mesh of a CharacteristicFaultSource
+    """
+    # Get the surface
+    sfc = src.surface.mesh
+
+    # Get mesh
+    mesh = RectangularMesh(sfc.mesh.lons, sfc.mesh.lats, sfc.mesh.depths)
+
+    return mesh
+
+
+def get_simple_mesh(src):
+    """
+    Get the mesh of a SimpleFaultSource
+    """
+    # Get the surface
+    sfc = SimpleFaultSurface.from_fault_data(
+       src.fault_trace, src.upper_seismogenic_depth,
+       src.lower_seismogenic_depth, src.dip, mesh_spacing=5)
+
+    # Get mesh
+    mesh = RectangularMesh(sfc.mesh.lons, sfc.mesh.lats, sfc.mesh.depths)
+
+    return mesh
+
+
+def get_kite_mesh(src):
+    """
+    Get the mesh of a KiteFaultSource
+    """
+    # Get the sfc (idl set to false as not automated yet in engine)
+    sfc = src.surface.from_profiles(src.profiles, src.profiles_sampling,
+                                    src.rupture_mesh_spacing, idl=False,
+                                    align=False)
+
+    # Get mesh
+    mesh = RectangularMesh(sfc.mesh.lons, sfc.mesh.lats, sfc.mesh.depths)
+
+    return mesh
+
+
+def get_geoms(srcs, geom_model):
+    """
+    Extract the geometry of each fault source and write to a geoJSON
+    """
+    traces = []
+    polys = []
+    suids = []
+
+    # First non-MFS sources
+    for i, src in enumerate(srcs):
+        if isinstance(src, MultiFaultSource): # Skip MFS as use the geom model
+            continue                          # object to get these geometries
+        if isinstance(src, ComplexFaultSource):
+            surf = get_complex_mesh(src)
+        elif isinstance(src, CharacteristicFaultSource):
+            surf = get_characteristic_mesh(src)
+        elif isinstance(src, SimpleFaultSource):
+            surf = get_simple_mesh(src)
+        elif isinstance(src, KiteFaultSource):
+            surf = get_kite_mesh(src)
+        else:
+            raise ValueError(f"Unknown source typology admitted: ({type(src)})")
+        trace, poly = get_boundary_2d(surf)        
+        traces.append(trace)
+        polys.append(poly)
+        suids.append(i)
+
+    # Get geometries of the MFS sources too
+    if geom_model:
+        for i, key in enumerate(geom_model.sections):
+            surf = geom_model.sections[key]
+            trace, poly = get_boundary_2d(surf.mesh)
+            traces.append(trace)
+            polys.append(poly)
+            suids.append(i)
+
+    daf = pd.DataFrame({'suid': suids, 'geometry': polys})
+    gdaf_polys = gpd.GeoDataFrame(daf, geometry='geometry') 
+    daf = pd.DataFrame({'suid': suids, 'geometry': traces})
+    gdaf_traces = gpd.GeoDataFrame(daf, geometry='geometry') 
+
+    return gdaf_polys, gdaf_traces
+
+
+def export_faults(gdaf_polys, gdaf_traces, model_dir):
+    """
+    Export the GeoDataFrames of the fault traces and fault polygons
+    """
+    out_traces = os.path.join(model_dir, 'fault_traces.geojson')
+    out_polys = os.path.join(model_dir, 'fault_sections.geojson')
+    gdaf_traces.to_file(out_traces)
+    gdaf_polys.to_file(out_polys)
+    print(f"Fault traces exported to {out_traces}")
+    print(f"Fault sections exported to {out_polys}")
+
+    return out_polys, out_traces 
+
+
+def plot_faults(gdaf_polys, gdaf_traces, region, model_dir):
+    """
+    Plot the faults within the geoJSONs
+    """
+    import pygmt 
+    out = os.path.join(model_dir, "faults_plot.png")
+    fig = pygmt.Figure()
+    fig.basemap(region=region, projection="M15c", frame=True)
+    fig.coast(land="grey", water="skyblue")
+    fig.plot(gdaf_polys, pen="0.5p,red")
+    fig.plot(gdaf_traces, pen="1p,green")
+    fig.savefig(out)
+    print(f"Faults plotted and saved to {out}")
+    fig.show()
+
+
+def get_fault_geojsons(model_dir, plotting=False, plotting_region=None):
+    """
+    Write the fault sections and fault traces within the given hazard model
+    model to geojsons.
+
+    :param model_dir: directory containing the required hazard model
+
+    :param plotting: Boolean which if True creates a plot using the geoJSONs
+                     of the faults in the given hazard model
+
+    :param plotting_region: list of [min lon, max lon, min lat, max lat] used
+                            to define axis limits of the plotted geoJSONs
+    """
+    # Get the sources in the given model
+    srcs, geom_model = get_sources(model_dir)
+
+    # Now get the geometries
+    gdaf_polys, gdaf_traces = get_geoms(srcs, geom_model)
+
+    # Export into geoJSONs
+    out_polys, out_traces = export_faults(gdaf_polys, gdaf_traces, model_dir)
+
+    # Plot them if specified
+    if plotting:
+        plot_faults(gdaf_polys, gdaf_traces, plotting_region, model_dir)
+
+    return gdaf_polys, gdaf_traces, out_polys, out_traces

openquake/man/tests/data/checks/expected/fault_sections.geojson

@@ -0,0 +1,6 @@
+{
+"type": "FeatureCollection",
+"features": [
+{ "type": "Feature", "properties": { "suid": 0 }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 1.021395468520495, -0.274887161684799 ], [ 1.061615435442728, -0.254777787522597 ], [ 1.101835125551899, -0.234668373923735 ], [ 1.142054573532541, -0.21455891841192 ], [ 1.182273814068017, -0.194449418510723 ], [ 1.222492881840751, -0.174339871743573 ], [ 1.262711811532374, -0.154230275633758 ], [ 1.302930637823961, -0.134120627704419 ], [ 1.343149395396204, -0.114010925478548 ], [ 1.383368118929616, -0.093901166478988 ], [ 1.424116359046308, -0.074887187950001 ], [ 1.469082479230706, -0.074887358892875 ], [ 1.514048659583617, -0.074887468242194 ], [ 1.559014927801107, -0.074887515997653 ], [ 1.603981311579775, -0.074887502158835 ], [ 1.648947838616836, -0.074887426725208 ], [ 1.693914536610364, -0.074887289696126 ], [ 1.738881433259303, -0.074887091070829 ], [ 1.738881433259303, -0.074887091070829 ], [ 1.749579076381868, -0.112330675622099 ], [ 1.760276746916998, -0.149774256257488 ], [ 1.770974454002376, -0.187217831671684 ], [ 1.781672206775785, -0.224661400559353 ], [ 1.781672206775785, -0.224661400559353 ], [ 1.73670531012761, -0.224661599184538 ], [ 1.691738612134616, -0.224661736213544 ], [ 1.646772085097836, -0.224661811647129 ], [ 1.601805701319218, -0.22466182548594 ], [ 1.55683943310155, -0.224661777730507 ], [ 1.511873252748232, -0.224661668381249 ], [ 1.466907132563147, -0.224661497438471 ], [ 1.426158969044339, -0.243675465362055 ], [ 1.385940331651136, -0.263785213144987 ], [ 1.345721665490477, -0.283894904154221 ], [ 1.305502935881849, -0.304004540866913 ], [ 1.265284108144636, -0.324114125760065 ], [ 1.225065147597765, -0.344223661310532 ], [ 1.184846019559678, -0.364333149995023 ], [ 1.144626689347953, -0.384442594290103 ], [ 1.104407122279304, -0.404551996672199 ], [ 1.064187283669246, -0.424661359617599 ], [ 1.064187283669246, -0.424661359617599 ], [ 1.05348919727219, -0.387217829077592 ], [ 1.042791205374247, -0.349774285038472 ], [ 1.032093298836565, -0.312330728805726 ], [ 1.021395468520495, -0.274887161684799 ] ] ] } }
+]
+}

openquake/man/tests/data/checks/in/ssm/source_model.xml

@@ -0,0 +1,47 @@
+<?xml version="1.0" encoding="utf-8"?>
+<nrml
+xmlns="http://openquake.org/xmlns/nrml/0.5"
+xmlns:gml="http://www.opengis.net/gml"
+>
+    <sourceModel
+    name="Example Source Model Containing a Simple Fault Source"
+    >
+        <sourceGroup
+        name="group 1"
+        tectonicRegion="Active Shallow Crust"
+        >
+            <simpleFaultSource
+            id="3"
+            name="Simple Fault Source"
+            tectonicRegion="Active Shallow Crust"
+            >
+                <simpleFaultGeometry>
+                    <gml:LineString>
+                        <gml:posList>
+                            1.0000000E+00 -2.0000000E-01 1.4000000E+00 0.0000000E+00 1.7000000E+00 0.0000000E+00
+                        </gml:posList>
+                    </gml:LineString>
+                    <dip>
+                        3.0000000E+01
+                    </dip>
+                    <upperSeismoDepth>
+                        5.0000000E+00
+                    </upperSeismoDepth>
+                    <lowerSeismoDepth>
+                        1.5000000E+01
+                    </lowerSeismoDepth>
+                </simpleFaultGeometry>
+                <magScaleRel>
+                    WC1994
+                </magScaleRel>
+                <ruptAspectRatio>
+                    2.0000000E+00
+                </ruptAspectRatio>
+                <truncGutenbergRichterMFD aValue="4.2000000E+00" bValue="9.0000000E-01" maxMag="7.0000000E+00" minMag="6.5000000E+00"/>
+                <rake>
+                    9.0000000E+01
+                </rake>
+            </simpleFaultSource>
+        </sourceGroup>
+    </sourceModel>
+</nrml>