source_modelling.scripts.plot_mw_contributions

  1from pathlib import Path
  2from typing import Annotated
  3
  4import numpy as np
  5import typer
  6from matplotlib import pyplot as plt
  7
  8from source_modelling import moment, rupture_propagation, srf
  9from workflow.realisations import (
 10    RealisationMetadata,
 11    RupturePropagationConfig,
 12    SourceConfig,
 13)
 14
 15app = typer.Typer()
 16
 17
 18@app.command(help="Plot segment magnitudes against the Leonard scaling relation.")
 19def plot_mw_contributions(
 20    srf_ffp: Annotated[
 21        Path, typer.Argument(help="Path to SRF file", exists=True, dir_okay=False)
 22    ],
 23    realisation_ffp: Annotated[
 24        Path,
 25        typer.Argument(help="Realisation filepath", dir_okay=False, exists=True),
 26    ],
 27    output_ffp: Annotated[
 28        Path, typer.Argument(help="Output plot path.", writable=True, dir_okay=False)
 29    ],
 30    dpi: Annotated[
 31        float, typer.Option(help="Output plot DPI (higher is better).")
 32    ] = 300,
 33) -> None:
 34    """Plot segment magnitudes against the Leonard scaling relation.
 35
 36    Parameters
 37    ----------
 38    srf_ffp : Path
 39        Path to SRF file.
 40    realisation_ffp : Path
 41        Realisation filepath.
 42    output_ffp : Path
 43        Output plot path.
 44    dpi : float, default 300
 45        Output plot DPI (higher is better).
 46    """
 47    source_config = SourceConfig.read_from_realisation(realisation_ffp)
 48    rupture_propogation_config = RupturePropagationConfig.read_from_realisation(
 49        realisation_ffp
 50    )
 51    realisation_metadata = RealisationMetadata.read_from_realisation(realisation_ffp)
 52    total_area = sum(fault.area() for fault in source_config.source_geometries.values())
 53    smallest_area = min(
 54        fault.area() for fault in source_config.source_geometries.values()
 55    )
 56    area = np.linspace(smallest_area, total_area)
 57    fig, ax = plt.subplots()
 58    # Mw = log(area) + 3.995 is the Leonard2014 magnitude scaling relation
 59    # for average rake.
 60    ax.plot(
 61        area, np.log10(area) + 3.995, label="Leonard 2014 Interplate (Average Rake)"
 62    )
 63
 64    srf_data = srf.read_srf(srf_ffp)
 65    total_magnitude = moment.moment_to_magnitude(
 66        moment.MU * (srf_data.points["area"] * srf_data.points["slip"] / (100**3)).sum()
 67    )
 68    ax.scatter(total_area, total_magnitude, label="Total Magnitude")
 69
 70    segment_counter = 0
 71    point_counter = 0
 72    for fault_name in rupture_propagation.tree_nodes_in_order(
 73        rupture_propogation_config.rupture_causality_tree
 74    ):
 75        plane_count = len(source_config.source_geometries[fault_name].planes)
 76        segments = srf_data.header.iloc[segment_counter : segment_counter + plane_count]
 77
 78        num_points = (segments["nstk"] * segments["ndip"]).sum()
 79        individual_area = source_config.source_geometries[fault_name].area()
 80
 81        # get all points associated with all segments in the current fault
 82        segment_points = srf_data.points.iloc[
 83            point_counter : point_counter + num_points
 84        ]
 85        individual_magnitude = moment.moment_to_magnitude(
 86            (
 87                moment.MU * segment_points["area"] * segment_points["slip"] / (100**3)
 88            ).sum()
 89        )
 90        ax.scatter(individual_area, individual_magnitude, label=fault_name)
 91
 92        # advance segment counter and point counter to skip all points from the current point
 93        segment_counter += plane_count
 94        point_counter += num_points
 95
 96    ax.set_xlabel("Area (m^2)")
 97    ax.set_ylabel("Mw")
 98    ax.set_xscale("log")
 99    ax.legend()
100    ax.set_title(f"Log Area vs Magnitude ({realisation_metadata.name})")
101    fig.savefig(output_ffp, dpi=dpi)
102
103
104if __name__ == "__main__":
105    app()