Meshing fault surfaces in Python (without Leapfrog Geo)

This notebook demonstrates how to build cut, depth-trimmed triangular fault meshes entirely in Python, without using Leapfrog Geo for the meshing/cutting step. It picks up where the “defining connected fault systems (without smoothing)” tutorial leaves off and re-uses the same shapefile, fault systems CSV, and cutting hierarchy CSV.

The workflow is:

  1. Re-build the curated fault network from the shapefile, fault systems CSV, and cutting hierarchy CSV.

  2. (Optional) Suggest, edit, and apply trace extensions so that fault traces extend far enough to be properly cut by neighbouring faults.

  3. Generate depth contours and a triangulated surface for every curated fault, and save them as VTK/OBJ/GeoJSON.

  4. (Optional) Merge the uncut meshes into a single VTK file for inspection in ParaView.

  5. Cut each fault mesh by all higher-priority faults in the cutting hierarchy and trim the result against a depth-of-faulting surface.

  6. (Optional) Merge the cut meshes into a single VTK file for inspection in ParaView.

Note: Step 5 needs a depth-of-faulting raster (here with_hannu_mods.tif). Drop your own raster into the working directory and update the filename if you are following along with a different dataset.

Imports and parameters

import os
import numpy as np

from fault_mesh.faults.leapfrog import LeapfrogMultiFault
from fault_mesh.faults.mesh import FaultMesh
from fault_mesh.io.array_operations import read_raster

# Output coordinate system EPSG code (set to None to disable reprojection)
epsg = 2193

# Distance tolerance (m) used when finding connections between segment ends
dist_tolerance = 1000.

Re-building the curated fault network

We read the same shapefile that was used in the connections tutorial, find segment-to-segment connections, and apply the manually-edited fault systems CSV and cutting hierarchy CSV.

fault_fname = "tutorial_gis/central_nz_minimal_data.shp"

fault_data = LeapfrogMultiFault.from_nz_cfm_shp(
    fault_fname,
    remove_colons=True,
    epsg=epsg,
    smoothing_n=None,
    dip_multiplier=1.0,
    exclude_zero=False,
)
fault_data.segment_distance_tolerance = dist_tolerance
fault_data.find_connections(verbose=False)

# Apply edited fault-system definitions and the cutting hierarchy
fault_data.read_fault_systems("./define_connections_data/central_gt1_5_connected_edited.csv")
fault_data.generate_curated_faults()
fault_data.read_cutting_hierarchy("./define_connections_data/central_gt1_5_hierarchy_edited.csv")

Trace extensions

Some faults need their surface traces extended at one or both ends so that, when projected to depth, they reach far enough laterally to be cut cleanly by neighbouring faults. This is especially relevant for short faults that terminate against a longer, faster-slipping fault.

Calling suggest_trace_extensions writes a CSV listing every end of every fault that the algorithm thinks might benefit from extending, together with default extend_distance and fit_distance values. Passing geojson_out_file additionally writes a GeoJSON containing the suggested extended trace geometries — drop it into QGIS to sanity-check the suggestions before committing to them. Edit the CSV by hand to keep only the rows you want, adjust the distances, and save it under a new name (here trace_extensions_edited.csv).

fault_data.suggest_trace_extensions(
    out_file="suggested_trace_extensions.csv",
    geojson_out_file="suggested_trace_extensions.geojson",
    fit_distance=5.e3,
    extend_distance=40.e3,
    proximity_threshold=1.e3,
)

Once you are happy with the edited CSV, read it in and apply the extensions to the curated fault traces.

fault_data.read_trace_extensions("trace_extensions_edited.csv")
fault_data.apply_trace_extensions()

Generating uncut meshes for every fault surface

Next we loop over every curated fault, build its depth contours (here from the surface to 32 km depth in 500 m steps), and triangulate them into a surface mesh. Each mesh is written out as both VTK (for ParaView) and OBJ (which we re-load in the cutting step), and the contours are written as GeoJSON for inspection in GIS software.

If a particular fault fails the standard meshing routine (for example because of awkward geometry), we fall back to mesh_simple_contours, which is more forgiving but produces a less faithful surface, and write the failing contours out for debugging.

for dir_name in ["test_vtks", "test_contours", "test_objs"]:
    if not os.path.exists(dir_name):
        os.mkdir(dir_name)

# Loop over every curated (multi-segment) fault and build its mesh
for fault in fault_data.curated_faults:
    print(fault.name)
    try:
        # Project the surface trace down-dip and slice into horizontal contours from
        # 0 m to 32 km depth at 500 m intervals; smoothing=False because the input
        # traces already have closely-spaced vertices
        fault.generate_depth_contours(np.arange(0., 32000., 500.), smoothing=False)
        # Triangulate the contours into a single fault surface; resolution sets the
        # target along-strike triangle edge length (m)
        mesh = fault.mesh_fault_surface(check_mesh=False, resolution=500.)
    except Exception as e:
        # If the standard mesher fails (usually due to awkward contour geometry),
        # fall back to the simpler contour-stitching mesher and dump the failing
        # contours to a GeoJSON so we can inspect them in GIS
        print(f"Failed to mesh {fault.name}: {e}")
        fault.contours.to_file(f"failed_contours_{fault.name}.geojson", driver="GeoJSON")
        mesh = fault.mesh_simple_contours(np.arange(0., 32000., 500.))

    # Build per-fault output paths
    vtk_file = os.path.join("test_vtks", f"{fault.name}_depth_contours.vtk")          # for ParaView
    obj_file = os.path.join("test_objs", f"{fault.name}_depth_contours.obj")          # re-loaded in the cutting step
    geojson_file = os.path.join("test_contours", f"{fault.name}_depth_contours.geojson")  # for GIS inspection

    # Write contours and mesh to disk
    fault.contours.to_file(geojson_file, driver="GeoJSON")
    mesh.write(vtk_file)
    mesh.write(obj_file)

(Optional) Merge the uncut meshes for ParaView

If you want to look at every uncut fault in a single ParaView scene, the easiest thing to do is merge them into one VTK file.

import meshio
import pyvista as pv
from glob import glob

uncut_meshes = [pv.from_meshio(meshio.read(f)) for f in glob("test_vtks/*.vtk")]
if uncut_meshes:
    pv.merge(uncut_meshes).save("merged_mesh_uncut.vtk")

Cutting the meshes against each other and trimming by a depth surface

Now that every fault has an uncut mesh, we can cut each one against all higher-priority faults in the cutting hierarchy. The cutting step also trims the result against a depth-of-faulting surface so that the meshes don’t extend below the seismogenic zone.

Inputs:

  • with_hannu_mods.tif: a raster of depth of faulting (Z values are negative-down, in metres). Replace this with your own raster as needed.

  • additional_cuts.csv: optional header-less CSV of fault_to_cut,cutting_fault pairs that should always be cut. Use it when the slip-rate-based hierarchy puts two faults in the wrong order, or (more often) for edge cases the cutting logic doesn’t cover. A cutting fault listed here is tried even if it ranks below the fault it cuts.

  • excluded_cuts.csv: optional CSV in the same format, listing pairs that must never be cut. Use it when the cutting logic makes a cut that shouldn’t happen, typically one that leaves a jagged outline. This tutorial doesn’t need one; the code below reads it only if the file exists.

Tuning knobs:

  • min_distance (m): minimum trace-to-trace separation below which we skip a cut. This avoids spurious cuts driven by KDTree artefacts or tiny faults that happen to sit right next to a much larger one.

  • bottom_depth (m, negative-down): lower depth bound used when generating the cutting fragment.

  • fancy_cutting=True: enables the more careful cutting routine that takes the already-cut state of faults higher in the hierarchy into account.

The Mesh cutting logic page explains how the code decides whether to cut and which side of a cut to keep.

# Read depth-of-faulting raster as a PyVista surface
depth_pyvista = read_raster("with_hannu_mods.tif", use_z=True, out_crs=f"EPSG:{epsg}")

# Optional manual cut overrides (fault_to_cut,cutting_fault pairs)
fault_data.read_additional_cuts("additional_cuts.csv")
if os.path.exists("excluded_cuts.csv"):
    fault_data.read_excluded_cuts("excluded_cuts.csv")

# Re-attach the uncut meshes we just wrote to disk
for fault in fault_data.curated_faults:
    obj_file = os.path.join("test_objs", f"{fault.name}_depth_contours.obj")
    if os.path.exists(obj_file):
        fault.mesh = FaultMesh.from_file(obj_file)

cutting_dict = {fault.name: fault for fault in fault_data.curated_faults if fault.mesh is not None}

output_dir = "test_final_meshes"
if not os.path.exists(output_dir):
    os.mkdir(output_dir)

# Threshold is currently unused, kept for forward-compatibility
threshold = 0.7
# Skip cuts where the two faults are closer than this (m)
min_distance = 2.e3

for fault_name in fault_data.cutting_hierarchy:
    print(f"Processing {fault_name}...")
    fault = fault_data.name_dict[fault_name]
    if fault.mesh is None:
        continue
    # Faults higher than this one in the hierarchy are candidate cutters
    faults_to_cut = fault_data.cutting_hierarchy[:fault_data.cutting_hierarchy.index(fault_name)]
    # Cutters forced by additional_cuts are tried even if they rank lower
    forced_cutters = [cutter for (cut_fault, cutter) in fault_data.additional_cuts
                      if cut_fault == fault_name and cutter not in faults_to_cut]

    for cut_name in faults_to_cut + forced_cutters:
        if cut_name not in cutting_dict:
            continue
        cutting_mesh = cutting_dict[cut_name].mesh

        # Manual overrides first: True forces a cut, False skips it, None falls
        # through to the geometric test
        decide_to_cut = fault_data.should_cut(fault_name, cut_name)
        if decide_to_cut is None:
            higher_faults_to_cut = fault_data.cutting_hierarchy[:fault_data.cutting_hierarchy.index(cut_name)]
            higher_meshes = [cutting_dict[name].mesh for name in higher_faults_to_cut if name in cutting_dict]
            decide_to_cut = fault.mesh.decide_whether_to_cut(
                cutting_mesh,
                threshold=threshold,
                min_distance=min_distance,
                higher_meshes=higher_meshes,
                bottom_depth=-31500.,
                fancy_cutting=True,
            )
        if decide_to_cut:
            print(f"Cutting {fault.name} by {cut_name}")
            cutting_fragment = cutting_mesh.generate_cutting_mesh(fault.mesh, max_distance=5.e3)
            new_mesh = fault.mesh.cut_mesh(
                cutting_mesh,
                fault_trace=fault.original_nztm_trace_array,
                cutting_fragment=cutting_fragment,
                fancy_cutting=True,
            )
            fault_data.name_dict[fault.name].mesh = new_mesh
            cutting_dict[fault.name].mesh = new_mesh

    # Trim the cut mesh against the depth-of-faulting surface and save
    try:
        depth_trimmed = fault_data.name_dict[fault.name].mesh.cut_mesh_pv(
            depth_pyvista,
            fault_trace=fault.original_nztm_trace_array,
        )
        fault_data.name_dict[fault.name].mesh = depth_trimmed
    except Exception as e:
        print(f"Error trimming depth for {fault.name}: {e}")
        continue

    output_file = os.path.join(output_dir, f"{fault.name}_cut.obj")
    fault_data.name_dict[fault.name].mesh.mesh.write(output_file)
    print(f"Done: {fault.name}")

(Optional) Merge the cut meshes for ParaView

Same idea as above, but for the final cut and depth-trimmed meshes. The result is a single VTK file you can drag into ParaView.

cut_meshes = [pv.from_meshio(meshio.read(f)) for f in glob("test_final_meshes/*.obj")]
if cut_meshes:
    pv.merge(cut_meshes).save("merged_mesh_cut.vtk")

Next steps

For a full fault model, these steps are run as a script, and the cut meshes are then checked, reviewed in 3-D with PyVista, and remeshed. See Whole-model workflow without Leapfrog.