source_modelling.scripts.plot_srf_moment

  1"""Utility script to plot moment over time for an SRF."""
  2
  3from pathlib import Path
  4from typing import Annotated, Optional
  5
  6import typer
  7from matplotlib import pyplot as plt
  8
  9from source_modelling import moment, rupture_propagation, srf
 10
 11app = typer.Typer()
 12
 13
 14@app.command(help="Plot released moment for an SRF over time.")
 15def plot_srf_moment(
 16    srf_ffp: Annotated[
 17        Path,
 18        typer.Argument(
 19            help="SRF filepath to plot", exists=True, readable=True, dir_okay=False
 20        ),
 21    ],
 22    output_png_ffp: Annotated[
 23        Path, typer.Argument(help="Output plot path", writable=True, dir_okay=False)
 24    ],
 25    dpi: Annotated[
 26        int, typer.Option(help="Plot image pixel density (higher = better)", min=300)
 27    ] = 300,
 28    realisation_ffp: Annotated[
 29        Optional[Path],
 30        typer.Option(
 31            help="Path to realisation, used to plot individual fault contribution."
 32        ),
 33    ] = None,
 34) -> None:
 35    """Plot released moment for an SRF over time.
 36
 37    Parameters
 38    ----------
 39    srf_ffp : Path
 40        SRF filepath to plot.
 41    output_png_ffp : Path
 42        Output plot path.
 43    dpi : float, default 300
 44        Plot image pixel density (higher = better).
 45    realisation_ffp : Optional[Path], default None
 46        Path to realisation, used to plot individual fault contribution.
 47    """
 48    srf_data = srf.read_srf(srf_ffp)
 49
 50    magnitude = moment.moment_to_magnitude(
 51        moment.MU * (srf_data.points["area"] * srf_data.points["slip"] / (100**3)).sum()
 52    )
 53
 54    dt = srf_data.points["dt"].iloc[0]
 55
 56    overall_moment_rate = moment.moment_rate_over_time_from_slip(
 57        srf_data.points["area"], srf_data.slip, dt, srf_data.nt
 58    )
 59    fig, ax = plt.subplots()
 60    ax.plot(
 61        overall_moment_rate.index.values,
 62        overall_moment_rate["moment_rate"],
 63        label="Overall Moment Rate",
 64    )
 65
 66    if realisation_ffp: # pragma: no cover
 67        # NOTE: this import is here because the workflow is, as yet,
 68        # not ready to be installed along-side source modelling.
 69        from workflow.realisations import RupturePropagationConfig, SourceConfig
 70
 71        source_config = SourceConfig.read_from_realisation(realisation_ffp)
 72        rupture_propogation_config = RupturePropagationConfig.read_from_realisation(
 73            realisation_ffp
 74        )
 75        segment_counter = 0
 76        point_counter = 0
 77        for fault_name in rupture_propagation.tree_nodes_in_order(
 78            rupture_propogation_config.rupture_causality_tree
 79        ):
 80            plane_count = len(source_config.source_geometries[fault_name].planes)
 81            segments = srf_data.header.iloc[
 82                segment_counter : segment_counter + plane_count
 83            ]
 84            num_points = (segments["nstk"] * segments["ndip"]).sum()
 85            individual_moment_rate = moment.moment_rate_over_time_from_slip(
 86                srf_data.points["area"]
 87                .iloc[point_counter : point_counter + num_points]
 88                .to_numpy(),
 89                srf_data.slip[point_counter : point_counter + num_points],
 90                dt,
 91                srf_data.nt,
 92            )
 93            ax.plot(
 94                individual_moment_rate.index.values,
 95                individual_moment_rate["moment_rate"],
 96                label=fault_name,
 97            )
 98            segment_counter += plane_count
 99            point_counter += num_points
100
101    ax.set_ylabel("Moment Rate (Nm/s)")
102    ax.set_xlabel("Time (s)")
103    ax.legend()
104    ax.set_title(f"Moment over Time (Total Mw: {magnitude:.2f})")
105
106    fig.savefig(output_png_ffp, dpi=dpi)
107
108
109if __name__ == "__main__":
110    app()