source_modelling.scripts.plot_srf_cumulative_moment

  1"""Utility script to plot cumulative moment over time for an SRF."""
  2
  3from pathlib import Path
  4from typing import Annotated, Optional
  5
  6import numpy as np
  7import typer
  8from matplotlib import pyplot as plt
  9
 10from source_modelling import moment, rupture_propagation, srf
 11
 12app = typer.Typer()
 13
 14
 15@app.command(help="Plot cumulative moment for an SRF over time.")
 16def plot_srf_cumulative_moment(
 17    srf_ffp: Annotated[
 18        Path,
 19        typer.Argument(
 20            help="SRF filepath to plot", exists=True, readable=True, dir_okay=False
 21        ),
 22    ],
 23    output_png_ffp: Annotated[
 24        Path, typer.Argument(help="Output plot path", writable=True, dir_okay=False)
 25    ],
 26    dpi: Annotated[
 27        int, typer.Option(help="Plot image pixel density (higher = better)", min=300)
 28    ] = 300,
 29    realisation_ffp: Annotated[
 30        Optional[Path],
 31        typer.Option(
 32            help="Path to realisation, used to plot individual fault contribution."
 33        ),
 34    ] = None,
 35    min_shade_cutoff: Annotated[
 36        float, typer.Option(help="Minimum shading cutoff", min=0, max=1)
 37    ] = 0.05,
 38    max_shade_cutoff: Annotated[
 39        float, typer.Option(help="Maximum shading cutoff", min=0, max=1)
 40    ] = 0.95,
 41):
 42    """Plot cumulative moment for an SRF over time.
 43
 44    Parameters
 45    ----------
 46    srf_ffp : Annotated[ Path, typer.Argument( help
 47        SRF filepath to plot.
 48    output_png_ffp : Annotated[ Path, typer.Argument(help
 49        Output plot path.
 50    dpi : Annotated[ int, typer.Option(help
 51        Plot image pixel density (higher = better).
 52    realisation_ffp : Annotated[ Optional[Path], typer.Option( help
 53        Path to realisation, used to plot individual fault
 54        contribution.
 55    min_shade_cutoff : Annotated[ float, typer.Option(help
 56        Minimum shading cutoff.
 57    max_shade_cutoff : Annotated[ float, typer.Option(help
 58        Maximum shading cutoff.
 59    """
 60    srf_data = srf.read_srf(srf_ffp)
 61
 62    overall_moment = moment.moment_over_time_from_moment_rate(
 63        moment.moment_rate_over_time_from_slip(
 64            srf_data.points["area"],
 65            srf_data.slip,
 66            srf_data.dt,
 67            srf_data.nt,
 68        )
 69    )
 70    total_moment = overall_moment["moment"].iloc[-1]
 71
 72    shaded_moments = overall_moment[
 73        (overall_moment["moment"] >= total_moment * min_shade_cutoff)
 74        & (overall_moment["moment"] <= total_moment * max_shade_cutoff)
 75    ]
 76    fig, ax = plt.subplots()
 77    ax.fill_between(shaded_moments.index.values, shaded_moments["moment"], alpha=0.2)
 78    ax.plot(
 79        overall_moment.index.values, overall_moment["moment"], label="Overall Moment"
 80    )
 81
 82    if realisation_ffp: # pragma: no cover
 83        # NOTE: this import is here because the workflow is, as yet,
 84        # not ready to be installed along-side source modelling.
 85        from workflow.realisations import RupturePropagationConfig, SourceConfig
 86
 87        source_config = SourceConfig.read_from_realisation(realisation_ffp)
 88        rupture_propogation_config = RupturePropagationConfig.read_from_realisation(
 89            realisation_ffp
 90        )
 91        segment_counter = 0
 92        point_counter = 0
 93        for fault_name in rupture_propagation.tree_nodes_in_order(
 94            rupture_propogation_config.rupture_causality_tree
 95        ):
 96            plane_count = len(source_config.source_geometries[fault_name].planes)
 97            segments = srf_data.header.iloc[
 98                segment_counter : segment_counter + plane_count
 99            ]
100            num_points = (segments["nstk"] * segments["ndip"]).sum()
101            segment_points = srf_data.points.iloc[
102                point_counter : point_counter + num_points
103            ]
104
105            individual_moment = moment.moment_over_time_from_moment_rate(
106                moment.moment_rate_over_time_from_slip(
107                    segment_points["area"],
108                    srf_data.slip[point_counter : point_counter + num_points],
109                    srf_data.dt,
110                    srf_data.nt,
111                )
112            )
113
114            ax.plot(
115                individual_moment.index.values,
116                individual_moment["moment"],
117                label=fault_name,
118            )
119            total_moment = individual_moment["moment"].iloc[-1]
120            shaded_moments = individual_moment[
121                (individual_moment["moment"] >= total_moment * min_shade_cutoff)
122                & (individual_moment["moment"] <= total_moment * max_shade_cutoff)
123            ]
124            ax.fill_between(
125                shaded_moments.index.values, shaded_moments["moment"], alpha=0.2
126            )
127            segment_counter += plane_count
128            point_counter += num_points
129
130    ax.set_ylabel("Cumulative Moment (Nm)")
131    ax.set_xlabel("Time (s)")
132    ax.legend()
133    min_shade_percent = int(np.round(min_shade_cutoff * 100))
134    max_shade_percent = int(np.round(max_shade_cutoff * 100))
135    ax.set_title(
136        f"Cumulative Moment over Time (Shaded Area: {min_shade_percent}% - {max_shade_percent}%)"
137    )
138
139    fig.savefig(output_png_ffp, dpi=dpi)
140
141
142if __name__ == "__main__":
143    app()