#!/usr/bin/env python3 # -*- coding: utf-8 -*- __version__ = "1.0.0" """ff_slip_dist_latlon.py Read combined reg_inv output with columns: lon lat depth slip moment Key conventions --------------- * The input moment column is scaled. param_M0 converts it to dyne-cm: M0_dyne_cm = moment * param_M0 Default: param_M0 = 1.0e20. * Subfault area is computed from grid spacing: subfault_area_km2 = (x_length_km/xgrid_num) * (z_length_km/zgrid_num) subfault_area_cm2 = subfault_area_km2 * (100000^2) * --use_moment_for_slip affects only the active slip used for slip_mag and plotting. The original input slip is always preserved as slip_input. * mu_gpa_rake1/2 is always computed from original input slip. * For two-rake output, moment_mag at each subfault is sqrt(moment_rake1^2 + moment_rake2^2). * Plot x/y scale is equal by default: 1 km along strike has the same length as 1 km along dip. """ import argparse import json import math import re import sys import subprocess import traceback from pathlib import Path import numpy as np import pandas as pd import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from ff_catalog_common import ( add_catalog_model_depth, apply_event_filters, find_column, load_polygon_points, read_and_filter_repeating, read_and_filter_seismicity, read_nca_eqddrt_catalog, read_repeating_catalog, resolve_repeating_locations_from_seis_catalog, ) from ff_map_common import ( map_tick_values, read_map_line_segments, map_midpoint_cells, ) from matplotlib.colors import Normalize from matplotlib.patches import Circle from matplotlib.ticker import MultipleLocator from matplotlib.collections import PolyCollection from matplotlib.path import Path as MplPath COLS = ["lon", "lat", "depth", "slip", "moment"] CM_PER_KM = 100000.0 ONE_KM2_CM2 = CM_PER_KM ** 2 PARAM_M0_DEFAULT = 1.0e20 DYNE_CM_TO_N_M = 1.0e-7 CATALOG_COLUMNS = [ "year", "month", "day", "hour", "minute", "second", "lat", "lon", "dep", "ex", "ey", "ez", "mag", "id", "ver", "base", "meth" ] REPEATING_COLUMNS = [ "year_decimal", "rep_lat", "rep_lon", "rep_dep", "cumD_cm", "evid", "csid", "rep_mag", "D_cm" ] LAT_CANDIDATES = ["lat", "latitude", "Lat", "Latitude", "LAT", "y", "Y"] LON_CANDIDATES = ["lon", "longitude", "lng", "Lon", "Longitude", "LON", "x", "X"] def fail(msg): raise RuntimeError("\nERROR: " + msg + "\n") def clamp(v, lo, hi): return lo if v < lo else hi if v > hi else v def mu_pa_from_depth(depth): return 10.0 ** (0.2607 * np.log(depth) + 9.8387) def dyne_cm_to_N_m(M0_dyne_cm): """Convert scalar moment from dyne-cm to N-m. 1 dyne = 1e-5 N and 1 cm = 1e-2 m, so 1 dyne-cm = 1e-7 N-m. """ return M0_dyne_cm * DYNE_CM_TO_N_M def moment_magnitude_from_N_m(M0_N_m): """Return Mw from scalar moment in N-m using log10(M0_N_m)=1.5*Mw+9.05.""" if M0_N_m <= 0.0 or not np.isfinite(M0_N_m): return float("nan") return (math.log10(M0_N_m) - 9.05) / 1.5 def mu_gpa_from_moment_and_slip(moment, slip_cm, subfault_area_cm2, param_M0): """Infer mu in GPa from moment, slip_cm, subfault area, and param_M0. M0_dyne_cm = moment * param_M0 mu_gpa = M0_dyne_cm / slip_cm / subfault_area_cm2 / 10 / 1e9 Non-positive slip or area gives NaN. """ moment_arr = np.asarray(moment, dtype=float) slip_arr = np.asarray(slip_cm, dtype=float) area_arr = np.asarray(subfault_area_cm2, dtype=float) out = np.full_like(moment_arr, np.nan, dtype=float) mask = (slip_arr > 0.0) & (area_arr > 0.0) out[mask] = moment_arr[mask] * param_M0 / slip_arr[mask] / area_arr[mask] / 10.0 / 1.0e9 return out def read_combined_file(path): path = Path(path) if not path.is_file(): fail(f"combined file not found: {path}") rows = [] with path.open("r", encoding="utf-8") as fin: for line_number, line in enumerate(fin, start=1): stripped = line.strip() if not stripped: continue parts = stripped.split() if len(parts) != 5: fail(f"bad format in {path} at line {line_number}: expected 5 columns, got {len(parts)}\nline: {line.rstrip()}") try: rows.append([float(v) for v in parts]) except ValueError as exc: fail(f"non-numeric value in {path} at line {line_number}\nline: {line.rstrip()}\n{exc}") if not rows: fail(f"combined file has no data rows: {path}") return pd.DataFrame(rows, columns=COLS) def lon_to_0360(lon): lon_arr = np.asarray(lon, dtype=float) return np.where(lon_arr < 0.0, lon_arr + 360.0, lon_arr) def find_slab_ids(slabdir): return sorted({f.name.split("_")[0] for f in sorted(Path(slabdir).glob("*_slab2_dep*.grd"))}) def find_grid_file(slabdir, slab_id, quantity): matches = sorted(Path(slabdir).glob(f"{slab_id}_slab2_{quantity}*.grd")) return matches[0] if matches else None def run_grdtrack_batch(lons, lats, grdfile, gmt_cmd="gmt"): cmd = [gmt_cmd, "grdtrack", f"-G{grdfile}", "-Z"] input_text = "".join(f"{float(lon):.8f} {float(lat):.8f}\n" for lon, lat in zip(lons, lats)) try: result = subprocess.run(cmd, input=input_text, text=True, capture_output=True, check=False) except FileNotFoundError: raise RuntimeError(f"Could not find GMT command '{gmt_cmd}'. Make sure GMT is installed and available in your PATH.") if result.returncode != 0: fail( "GMT grdtrack failed.\n" f"command: {' '.join(cmd)}\n" f"return code: {result.returncode}\n" f"stderr: {result.stderr.strip() or ''}\n" f"stdout: {result.stdout.strip() or ''}" ) lines = result.stdout.splitlines() if len(lines) != len(lons): fail( "GMT grdtrack returned an unexpected number of rows.\n" f"grid: {grdfile}\nrequested points: {len(lons)}\n" f"returned rows: {len(lines)}" ) vals = [] for row_number, line in enumerate(lines, start=1): fields = line.split() if not fields: fail(f"GMT grdtrack returned an empty row at output row {row_number}") try: vals.append(float(fields[0])) except ValueError as exc: fail( "GMT grdtrack returned a non-numeric value.\n" f"output row: {row_number}\nline: {line}\nerror: {exc}" ) return np.asarray(vals, dtype=float) def _xarray_dataarray_from_grid(grdfile): import xarray as xr try: da = xr.open_dataarray(grdfile) except Exception: ds = xr.open_dataset(grdfile) data_vars = list(ds.data_vars) if not data_vars: raise RuntimeError(f"No data variables found in grid file: {grdfile}") da = ds[data_vars[0]] if da.ndim > 2: da = da.squeeze(drop=True) return da def _detect_xy_coord_names(da): lower = {str(c).lower(): c for c in da.coords} x_name = next((lower[c] for c in ["lon", "longitude", "x"] if c in lower), None) y_name = next((lower[c] for c in ["lat", "latitude", "y"] if c in lower), None) dims = list(da.dims) if x_name is None and len(dims) >= 1: x_name = dims[-1] if y_name is None and len(dims) >= 2: y_name = dims[-2] if x_name is None or y_name is None: raise RuntimeError(f"Could not detect x/y coordinates for grid with coords={list(da.coords)} dims={list(da.dims)}") return x_name, y_name def sample_grid_xarray(lons, lats, grdfile): import xarray as xr da = _xarray_dataarray_from_grid(grdfile) x_name, y_name = _detect_xy_coord_names(da) sampled = da.interp({x_name: xr.DataArray(np.asarray(lons, dtype=float), dims="points"), y_name: xr.DataArray(np.asarray(lats, dtype=float), dims="points")}) return np.asarray(sampled.to_numpy(), dtype=float) def sample_grid_pygmt(lons, lats, grdfile): import pygmt pts = pd.DataFrame({"lon": np.asarray(lons, dtype=float), "lat": np.asarray(lats, dtype=float)}) out = pygmt.grdtrack(points=pts, grid=str(grdfile), newcolname="value") return pd.to_numeric(out["value"], errors="coerce").to_numpy(dtype=float) if "value" in out.columns else np.full(len(pts), np.nan) def sample_slab_grid(lons, lats, grdfile, args): backend = args.slab_sample_backend if backend == "auto": fail( "--slab_sample_backend auto is disabled because backend fallback can " "hide dependency or grid-reading errors. Choose xarray, gmt, or pygmt explicitly." ) try: if backend == "xarray": return sample_grid_xarray(lons, lats, grdfile), "xarray" if backend == "gmt": return run_grdtrack_batch(lons, lats, grdfile, gmt_cmd=args.gmt_cmd), "gmt" if backend == "pygmt": return sample_grid_pygmt(lons, lats, grdfile), "pygmt" except Exception as exc: fail( "Slab-grid sampling failed without fallback.\n" f"backend: {backend}\ngrid: {grdfile}\n" f"exception: {type(exc).__name__}: {exc}" ) fail(f"unsupported slab sampling backend: {backend}") def slab_ids_for_query(args): return [s.strip() for s in args.slab_id.split(",") if s.strip()] if args.slab_id else find_slab_ids(args.slabdir) def top_trace_samples_for_segment(seg, args): top = seg[seg["z_index"].astype(int) == int(seg["z_index"].min())].copy() if top.empty: top = seg.copy() top = top.sort_values("x_km").drop_duplicates(subset=["x_km"]) xmin, xmax = float(top["x_km"].min()), float(top["x_km"].max()) dx = float(args.slab_sample_dx_km) if dx <= 0: fail("--slab_sample_dx_km must be positive") xs = np.arange(xmin, xmax + 0.5 * dx, dx) if len(xs) == 0 or xs[-1] < xmax: xs = np.append(xs, xmax) lon = np.interp(xs, top["x_km"].to_numpy(float), top["lon"].to_numpy(float)) lat = np.interp(xs, top["x_km"].to_numpy(float), top["lat"].to_numpy(float)) return xs, lon, lat def _print_slab_range(prefix, values, unit=""): arr = np.asarray(values, dtype=float) finite = arr[np.isfinite(arr)] if finite.size: print(f"# {prefix:<27s} = {finite.min():.3f} to {finite.max():.3f}{unit}") else: print(f"# {prefix:<27s} = no finite values") def query_slab_profile_for_segment(seg, args): seg_id = int(seg["segment_id"].iloc[0]) xs, lons, lats = top_trace_samples_for_segment(seg, args) lons_grid = lon_to_0360(lons) slab_offset = float(args.slab_depth_offset_km) projection_mode = "depth_only_collapse" if args.ignore_across_strike_distance_depth_projection else "full_across_strike_projection" print(f"# slab query segment {seg_id}") print(f"# requested samples = {len(xs)}") print(f"# sample dx = {args.slab_sample_dx_km:g} km") print(f"# x trace range = {xs.min():.3f} to {xs.max():.3f} km") print(f"# lon range = {np.nanmin(lons):.5f} to {np.nanmax(lons):.5f}") print(f"# lat range = {np.nanmin(lats):.5f} to {np.nanmax(lats):.5f}") print(f"# lon grid range = {np.nanmin(lons_grid):.5f} to {np.nanmax(lons_grid):.5f}") print(f"# depth offset = {slab_offset:g} km") print(f"# projection mode = {projection_mode}") rows = [] for sid in slab_ids_for_query(args): dep_file = find_grid_file(args.slabdir, sid, "dep") str_file = find_grid_file(args.slabdir, sid, "str") dip_file = find_grid_file(args.slabdir, sid, "dip") thk_file = find_grid_file(args.slabdir, sid, "thk") unc_file = find_grid_file(args.slabdir, sid, "unc") print(f"# slab query segment {seg_id} slab_id={sid}") print(f"# dep file = {dep_file if dep_file else 'missing'}") print(f"# str/dip/thk/unc files = {bool(str_file)}/{bool(dip_file)}/{bool(thk_file)}/{bool(unc_file)}") if dep_file is None: fail( f"required Slab2 depth grid not found for slab_id={sid} in " f"{args.slabdir}; expected a file matching {sid}_slab2_dep*.grd" ) raw_dep, dep_backend = sample_slab_grid(lons_grid, lats, dep_file, args) depth_km = -np.asarray(raw_dep, dtype=float) finite_dep = np.isfinite(depth_km) print(f"# sample backend = {dep_backend}") print(f"# finite depth samples = {finite_dep.sum()} / {len(xs)}") print(f"# NaN depth samples = {len(xs) - int(finite_dep.sum())}") if finite_dep.any(): print(f"# valid x range = {xs[finite_dep].min():.3f} to {xs[finite_dep].max():.3f} km") _print_slab_range("raw depth grid", raw_dep[finite_dep], "") _print_slab_range("positive depth", depth_km[finite_dep], " km") _print_slab_range("model depth", depth_km[finite_dep] + slab_offset, " km") else: message = ( f"all Slab2 depth samples are NaN/invalid for segment {seg_id}, " f"slab_id={sid}, grid={dep_file}; skipping this segment/model" ) if args.require_slab_profile: fail(message) print(f"# WARNING: {message}", file=sys.stderr) continue slab_str = np.full(len(xs), np.nan) slab_dip = np.full(len(xs), np.nan) slab_thk = np.full(len(xs), np.nan) slab_unc = np.full(len(xs), np.nan) if str_file is not None: slab_str, _ = sample_slab_grid(lons_grid, lats, str_file, args) print(f"# finite slab strike samples = {np.isfinite(slab_str).sum()} / {len(xs)}") if dip_file is not None: slab_dip, _ = sample_slab_grid(lons_grid, lats, dip_file, args) print(f"# finite slab dip samples = {np.isfinite(slab_dip).sum()} / {len(xs)}") if thk_file is not None: slab_thk, _ = sample_slab_grid(lons_grid, lats, thk_file, args) print(f"# finite slab thk samples = {np.isfinite(slab_thk).sum()} / {len(xs)}") if unc_file is not None: slab_unc, _ = sample_slab_grid(lons_grid, lats, unc_file, args) print(f"# finite slab unc samples = {np.isfinite(slab_unc).sum()} / {len(xs)}") for i in range(len(xs)): if not np.isfinite(depth_km[i]): continue rows.append({ "segment_id": seg_id, "slab_id": sid, "x_km_trace": xs[i], "lon": lons[i], "lat": lats[i], "grid_lon": lons_grid[i], "dep": depth_km[i], "dep_model": depth_km[i] + slab_offset, "depth_offset_km": slab_offset, "raw_depth_grid_value": raw_dep[i], "slab_strike_deg": slab_str[i], "slab_dip_deg": slab_dip[i], "slab_thickness_km": slab_thk[i], "slab_depth_unc_km": slab_unc[i], "sample_backend": dep_backend, "dep_file": str(dep_file), }) if not rows: print(f"# slab query segment {seg_id}: no valid rows after sampling") return pd.DataFrame() prof = pd.DataFrame(rows) projected = [] for sid, sub in prof.groupby("slab_id", sort=True): proj = project_events_to_segment(sub, seg, args) proj["x_km"] = proj["x_km_trace"] print(f"# slab projection segment {seg_id} slab_id={sid}") print(f"# projected rows = {len(proj)}") _print_slab_range("projected x", proj["x_km"], " km") _print_slab_range("projected z", proj["z_km"], " km") _print_slab_range("projection down-dip", proj["projection_down_dip_distance_km"], " km") projected.append(proj) return pd.concat(projected, ignore_index=True) if projected else pd.DataFrame() def build_slab_profiles(outputs, args): combined = outputs["combined_slip_dist_latlon"] segment_ids = sorted(combined["segment_id"].dropna().astype(int).unique()) profiles = [] print("# slab profile setup") print(f"# slabdir = {args.slabdir}") print(f"# slab_id option = {args.slab_id if args.slab_id else 'auto/all'}") print(f"# sample backend option = {args.slab_sample_backend}") print(f"# gmt command = {args.gmt_cmd}") print(f"# available slab IDs = {', '.join(find_slab_ids(args.slabdir)) if Path(args.slabdir).is_dir() else 'slabdir not found'}") for seg_id in segment_ids: seg = combined[combined["segment_id"].astype(int) == seg_id].copy() prof = query_slab_profile_for_segment(seg, args) print(f"# slab profile segment {seg_id}: rows={len(prof)}") if not prof.empty: profiles.append(prof) if not profiles: message = ( "Slab plotting was requested, but no valid Slab2 profile was produced " "for any segment/model. The Slab2 curve will be omitted; other figures " "and outputs will continue. Check --slabdir, --slab_id, grid coverage, " "and the explicitly selected backend." ) if args.require_slab_profile: fail(message) print(f"# WARNING: {message}", file=sys.stderr) return pd.DataFrame() out = pd.concat(profiles, ignore_index=True) print(f"# slab profile total rows = {len(out)}") print("# slab profile valid x range by segment/slab") for (seg_id, sid), sub in out.groupby(["segment_id", "slab_id"], sort=True): print(f"# seg {int(seg_id)} slab {sid}: count={len(sub)} x={sub['x_km'].min():.3f}-{sub['x_km'].max():.3f} km z={sub['z_km'].min():.3f}-{sub['z_km'].max():.3f} km") return out def add_bool_flag(parser, name, default=True, help_text=""): group = parser.add_mutually_exclusive_group() group.add_argument(f"--{name}", dest=name, action="store_true", help=help_text) group.add_argument(f"--no_{name}", dest=name, action="store_false", help=f"Disable {name}") parser.set_defaults(**{name: default}) def parse_multime_count(rupture_config): """Return the multimeN count encoded in rupture_config, or 1.""" matches = re.findall(r"(?:^|_)multime(\d+)(?:_|$)", str(rupture_config)) if len(matches) > 1: fail(f"rupture_config contains more than one multimeN field: {rupture_config}") count = int(matches[0]) if matches else 1 if count < 1: fail(f"time-window count must be at least 1: rupture_config={rupture_config}") return count def resolve_num_time_windows(args): """Resolve explicit/configured window count and reject contradictions.""" encoded = parse_multime_count(args.rupture_config) explicit = args.num_time_windows if explicit is None: return encoded if explicit < 1: fail(f"--num_time_windows must be at least 1, got {explicit}") has_encoded_multime = bool(re.search(r"(?:^|_)multime\d+(?:_|$)", args.rupture_config)) if has_encoded_multime and explicit != encoded: fail( "inconsistent time-window configuration\n" f"--rupture_config = {args.rupture_config}\n" f"encoded window count = {encoded}\n" f"--num_time_windows = {explicit}" ) return explicit def build_frw_layout(args, num_time_windows): """Build the confirmed segment -> rake -> time-window FRW layout.""" layout = [] for block_name, segment_id, rake_id in rupture_block_plan(args.rupture_config): xnum, znum, _, _ = segment_grid_params(segment_id, args) layout.append({ "block_name": block_name, "segment_id": segment_id, "rake_id": rake_id, "xgrid_num": xnum, "zgrid_num": znum, "num_time_windows": num_time_windows, }) return layout def read_frw_file(path, layout): """Read FRW grids in segment -> rake -> window -> z -> x order.""" path = Path(path) if not path.is_file(): fail(f"FRW file not found: {path}") physical_rows = [] with path.open("r", encoding="utf-8") as fin: for line_number, line in enumerate(fin, start=1): stripped = line.strip() if not stripped: continue values = [] for token in stripped.split(): try: values.append(float(token.replace("D", "E").replace("d", "e"))) except ValueError as exc: fail(f"non-numeric FRW value in {path} at line {line_number}: {token}\n{exc}") physical_rows.append((line_number, values)) expected_rows = sum(item["zgrid_num"] * item["num_time_windows"] for item in layout) if len(physical_rows) != expected_rows: fail( f"FRW row-count mismatch: {path}\n" f"expected nonblank rows = {expected_rows}\n" f"actual nonblank rows = {len(physical_rows)}" ) windows = {} cursor = 0 for item in layout: block_name = item["block_name"] xnum = item["xgrid_num"] znum = item["zgrid_num"] windows[block_name] = [] for window_id in range(1, item["num_time_windows"] + 1): rows = [] for z_index in range(znum): line_number, values = physical_rows[cursor] cursor += 1 if len(values) != xnum: fail( f"FRW column-count mismatch: {path}\n" f"block={block_name} window={window_id} z_index={z_index}\n" f"physical line={line_number}, expected columns={xnum}, actual columns={len(values)}" ) rows.append(values) grid = np.asarray(rows, dtype=float) if not np.isfinite(grid).all(): bad = np.argwhere(~np.isfinite(grid))[0] fail(f"non-finite FRW value: block={block_name} window={window_id} z_index={bad[0]} x_index={bad[1]}") if (grid < 0.0).any(): bad = np.argwhere(grid < 0.0)[0] fail(f"negative FRW value: block={block_name} window={window_id} z_index={bad[0]} x_index={bad[1]} value={grid[tuple(bad)]}") windows[block_name].append(grid) return windows def sum_frw_time_windows(frw_windows): return {name: np.sum(np.stack(grids, axis=0), axis=0) for name, grids in frw_windows.items()} def validate_frw_against_combined(frw_totals, combined_blocks, args): """Validate cumulative FRW moments against each combined-file block.""" summaries = {} for block_name, total_grid in frw_totals.items(): if block_name not in combined_blocks: fail(f"FRW block is absent from combined-file blocks: {block_name}") expected = combined_blocks[block_name]["moment"].to_numpy(dtype=float) actual = total_grid.ravel(order="C") if actual.size != expected.size: fail(f"FRW/combined size mismatch for {block_name}: FRW={actual.size}, combined={expected.size}") zero_mismatch = (actual == 0.0) != (expected == 0.0) close = np.isclose(actual, expected, rtol=args.frw_rtol, atol=args.frw_atol) bad = zero_mismatch | ~close abs_diff = np.abs(actual - expected) rel_diff = np.zeros_like(abs_diff) nonzero = expected != 0.0 rel_diff[nonzero] = abs_diff[nonzero] / np.abs(expected[nonzero]) summaries[block_name] = { "max_abs_diff": float(abs_diff.max(initial=0.0)), "max_rel_diff": float(rel_diff[nonzero].max(initial=0.0)), "zero_mask_mismatches": int(zero_mismatch.sum()), "window_count": len(frw_windows_global[block_name]), } if bad.any(): idx = int(np.flatnonzero(bad)[0]) xnum, _, _, _ = segment_grid_params(int(combined_blocks[block_name]["segment_id"].iloc[0]), args) z_index, x_index = divmod(idx, xnum) fail( f"FRW validation failed for {block_name}\n" f"z_index={z_index}, x_index={x_index}\n" f"summed FRW moment={actual[idx]:.12g}\n" f"combined moment={expected[idx]:.12g}\n" f"absolute difference={abs_diff[idx]:.12g}\n" f"relative difference={rel_diff[idx]:.12g}\n" f"rtol={args.frw_rtol}, atol={args.frw_atol}\n" f"zero-mask mismatch count={int(zero_mismatch.sum())}" ) return summaries def _parse_numbered_frw_patches(args, layout): """Normalize --frw_patch1..3 and enforce contiguous numbering.""" raw = [args.frw_patch1, args.frw_patch2, args.frw_patch3] if raw[1] is not None and raw[0] is None: fail("--frw_patch2 requires --frw_patch1; patch definitions must be consecutive") if raw[2] is not None and raw[1] is None: fail("--frw_patch3 requires --frw_patch2; patch definitions must be consecutive") if any(item is not None for item in raw) and not args.frw_file: fail("FRW patch generation requires --frw_file") valid_segments = {int(item["segment_id"]) for item in layout} outputs = [args.frw_patch1_output, args.frw_patch2_output, args.frw_patch3_output] patches = [] for i, values in enumerate(raw, start=1): if values is None: continue try: segment_f = float(values[0]) segment_id = int(segment_f) if segment_f != segment_id: raise ValueError("segment ID is not an integer") xmin, xmax, zmin, zmax = (float(v) for v in values[1:]) except ValueError as exc: fail(f"invalid --frw_patch{i} values {values}: {exc}") bounds = np.asarray([xmin, xmax, zmin, zmax], dtype=float) if not np.isfinite(bounds).all(): fail(f"--frw_patch{i} bounds must be finite: {values}") if segment_id not in valid_segments: fail(f"--frw_patch{i} requests segment {segment_id}; available segments are {sorted(valid_segments)}") if xmin > xmax or zmin > zmax: fail(f"--frw_patch{i} requires XMIN<=XMAX and ZMIN<=ZMAX; got {values}") output = Path(outputs[i - 1]) if output.resolve() == Path(args.frw_file).resolve(): fail(f"--frw_patch{i}_output must not overwrite the source FRW file: {output}") patches.append({ "patch_id": i, "segment_id": segment_id, "xmin_km": xmin, "xmax_km": xmax, "zmin_km": zmin, "zmax_km": zmax, "output": output, }) resolved_outputs = [p["output"].resolve() for p in patches] if len(set(resolved_outputs)) != len(resolved_outputs): fail("FRW patch output paths must be unique") return patches def _frw_patch_mask(patch, args): xnum, znum, xlen, zlen = segment_grid_params(patch["segment_id"], args) dx, dz = xlen / xnum, zlen / znum xs = np.arange(xnum, dtype=float) * dx zs = np.arange(znum, dtype=float) * dz mask = ((zs[:, None] >= patch["zmin_km"]) & (zs[:, None] <= patch["zmax_km"]) & (xs[None, :] >= patch["xmin_km"]) & (xs[None, :] <= patch["xmax_km"])) if not mask.any(): fail( f"--frw_patch{patch['patch_id']} selects no grid coordinates on segment {patch['segment_id']}; " f"requested x={patch['xmin_km']}..{patch['xmax_km']} km, " f"z={patch['zmin_km']}..{patch['zmax_km']} km, dx={dx}, dz={dz}" ) zi, xi = np.where(mask) meta = { "xnum": xnum, "znum": znum, "dx_km": dx, "dz_km": dz, "x_index_min": int(xi.min()), "x_index_max": int(xi.max()), "z_index_min": int(zi.min()), "z_index_max": int(zi.max()), "x_realized_min_km": float(xs[xi].min()), "x_realized_max_km": float(xs[xi].max()), "z_realized_min_km": float(zs[zi].min()), "z_realized_max_km": float(zs[zi].max()), "x_edge_min_km": float(xs[xi].min() - 0.5 * dx), "x_edge_max_km": float(xs[xi].max() + 0.5 * dx), "z_edge_min_km": float(zs[zi].min() - 0.5 * dz), "z_edge_max_km": float(zs[zi].max() + 0.5 * dz), "selected_x_count": int(np.unique(xi).size), "selected_z_count": int(np.unique(zi).size), "spatial_cell_count": int(mask.sum()), } return mask, meta def print_frw_patch_coordinate_guide(layout, args): """Print available local FRW grid-center and cell-edge coordinates.""" segment_ids = sorted({int(item["segment_id"]) for item in layout}) print("# FRW patch coordinate guide") print("# patch bounds are inclusive local subfault-center coordinates") for segment_id in segment_ids: xnum, znum, xlen, zlen = segment_grid_params(segment_id, args) dx, dz = xlen / xnum, zlen / znum x_center_max = (xnum - 1) * dx z_center_max = (znum - 1) * dz print(f"# segment {segment_id}") print(f"# grid shape (x,z) = {xnum} x {znum}") print(f"# along-strike spacing dx = {dx:.3f} km") print(f"# along-dip spacing dz = {dz:.3f} km") print(f"# along-strike centers = 0.000 to {x_center_max:.3f} km") print(f"# along-dip centers = 0.000 to {z_center_max:.3f} km") print(f"# along-strike cell edges = {-0.5 * dx:.3f} to {x_center_max + 0.5 * dx:.3f} km") print(f"# along-dip cell edges = {-0.5 * dz:.3f} to {z_center_max + 0.5 * dz:.3f} km") print(f"# valid center sequence = 0, {dx:g}, {2.0 * dx:g}, ... (x); 0, {dz:g}, {2.0 * dz:g}, ... (z)") def _warn_for_overlapping_frw_patches(prepared): for i, left in enumerate(prepared): for right in prepared[i + 1:]: if left["segment_id"] == right["segment_id"] and np.any(left["mask"] & right["mask"]): overlap = int(np.sum(left["mask"] & right["mask"])) print( f"# WARNING: FRW patches {left['patch_id']} and {right['patch_id']} overlap " f"on segment {left['segment_id']} in {overlap} spatial cells; overlapping moment " "is retained independently in both files.", file=sys.stderr ) def _make_frw_patch_windows(source_windows, layout, patch, mask): segment_by_block = {item["block_name"]: int(item["segment_id"]) for item in layout} patched = {} for block_name, grids in source_windows.items(): patched[block_name] = [] for source in grids: output = np.zeros_like(source) if segment_by_block[block_name] == patch["segment_id"]: output[mask] = source[mask] patched[block_name].append(output) return patched def write_frw_file(path, windows, layout): """Write FRW with no header in segment->rake->window->z->x order.""" path = Path(path) if not path.parent.is_dir(): fail(f"FRW patch output directory does not exist: {path.parent}") with path.open("w", encoding="utf-8") as fout: for item in layout: block_name = item["block_name"] grids = windows[block_name] if len(grids) != item["num_time_windows"]: fail(f"internal FRW window-count mismatch for {block_name}") for grid in grids: if grid.shape != (item["zgrid_num"], item["xgrid_num"]): fail(f"internal FRW grid-shape mismatch for {block_name}: {grid.shape}") for row in grid: fout.write(" ".join(f"{float(value):.3e}" for value in row) + "\n") def _rounded_frw_grid(grid): return np.asarray([[float(f"{float(value):.3e}") for value in row] for row in grid], dtype=float) def _validate_written_frw_patch(path, written_windows, layout, patch, mask): reread = read_frw_file(path, layout) segment_by_block = {item["block_name"]: int(item["segment_id"]) for item in layout} for block_name, grids in written_windows.items(): for window_index, expected_grid in enumerate(grids): expected = _rounded_frw_grid(expected_grid) actual = reread[block_name][window_index] if not np.array_equal(actual, expected): fail(f"FRW patch read-back mismatch: {path}, block={block_name}, window={window_index + 1}") if segment_by_block[block_name] != patch["segment_id"]: if np.count_nonzero(actual): fail(f"FRW patch has nonzero values on unselected segment: {path}, block={block_name}") elif np.count_nonzero(actual[~mask]): fail(f"FRW patch has nonzero values outside selected rectangle: {path}, block={block_name}") return reread def _normalize_angle_deg(angle): """Normalize an angle to [-180, 180).""" return (float(angle) + 180.0) % 360.0 - 180.0 def _circular_mean_deg(angles_deg, weights=None): angles = np.asarray(angles_deg, dtype=float) finite = np.isfinite(angles) if weights is None: weights_arr = np.ones_like(angles) else: weights_arr = np.asarray(weights, dtype=float) finite &= np.isfinite(weights_arr) & (weights_arr >= 0.0) angles = angles[finite] weights_arr = weights_arr[finite] if angles.size == 0 or weights_arr.sum() <= 0.0: return float("nan") radians = np.deg2rad(angles) result = math.degrees(math.atan2( float(np.sum(weights_arr * np.sin(radians))), float(np.sum(weights_arr * np.cos(radians))), )) return _normalize_angle_deg(result) def _circular_median_deg(angles_deg): """Return a circular L1 median chosen from observed finite angles.""" angles = np.asarray(angles_deg, dtype=float) angles = angles[np.isfinite(angles)] if angles.size == 0: return float("nan") candidates = np.unique([_normalize_angle_deg(v) for v in angles]) def objective(candidate): delta = np.abs((angles - candidate + 180.0) % 360.0 - 180.0) return float(delta.sum()) return float(min(candidates, key=lambda candidate: (objective(candidate), candidate))) def _segment_strike_dip(segment_id, args): if int(segment_id) == 1: return float(args.strike1_deg), float(args.dip1_deg) if int(segment_id) == 2: return float(args.strike2_deg), float(args.dip2_deg) fail(f"unsupported segment for focal mechanism: {segment_id}") def _double_couple_tensor_ned(strike_deg, dip_deg, rake_deg): """Return a normalized double-couple moment tensor in NED coordinates.""" strike, dip, rake = np.deg2rad([strike_deg, dip_deg, rake_deg]) mnn = -(np.sin(dip) * np.cos(rake) * np.sin(2.0 * strike) + np.sin(2.0 * dip) * np.sin(rake) * np.sin(strike) ** 2) mee = (np.sin(dip) * np.cos(rake) * np.sin(2.0 * strike) - np.sin(2.0 * dip) * np.sin(rake) * np.cos(strike) ** 2) mdd = np.sin(2.0 * dip) * np.sin(rake) mne = (np.sin(dip) * np.cos(rake) * np.cos(2.0 * strike) + 0.5 * np.sin(2.0 * dip) * np.sin(rake) * np.sin(2.0 * strike)) mnd = -(np.cos(dip) * np.cos(rake) * np.cos(strike) + np.cos(2.0 * dip) * np.sin(rake) * np.sin(strike)) med = -(np.cos(dip) * np.cos(rake) * np.sin(strike) - np.cos(2.0 * dip) * np.sin(rake) * np.cos(strike)) return np.asarray([[mnn, mne, mnd], [mne, mee, med], [mnd, med, mdd]], dtype=float) def _draw_beachball(ax, strike_deg, dip_deg, rake_deg, radius): """Draw a lower-hemisphere equal-area double-couple beach ball.""" samples = 420 axis_values = np.linspace(-radius, radius, samples) x, y = np.meshgrid(axis_values, axis_values) rr = np.sqrt(x * x + y * y) inside = rr <= radius rho = np.clip(rr / radius, 0.0, 1.0) theta = 2.0 * np.arcsin(rho / math.sqrt(2.0)) azimuth = np.arctan2(x, y) north = np.sin(theta) * np.cos(azimuth) east = np.sin(theta) * np.sin(azimuth) down = np.cos(theta) tensor = _double_couple_tensor_ned(strike_deg, dip_deg, rake_deg) amplitude = ( tensor[0, 0] * north * north + tensor[1, 1] * east * east + tensor[2, 2] * down * down + 2.0 * tensor[0, 1] * north * east + 2.0 * tensor[0, 2] * north * down + 2.0 * tensor[1, 2] * east * down ) amplitude = np.ma.array(amplitude, mask=~inside) ax.contourf(x, y, amplitude, levels=[-1.0e9, 0.0, 1.0e9], colors=["white", "black"], antialiased=True) ax.contour(x, y, amplitude, levels=[0.0], colors="black", linewidths=1.1) ax.add_patch(Circle((0.0, 0.0), radius, fill=False, edgecolor="black", linewidth=1.4)) ax.set_xlim(-1.12, 1.12) ax.set_ylim(-1.12, 1.12) ax.set_aspect("equal") ax.axis("off") def plot_frw_patch_focal_mechanisms(summaries, args): if not summaries: return output = Path(args.frw_patch_beachball_file) if not output.parent.is_dir(): fail(f"FRW patch focal-mechanism output directory does not exist: {output.parent}") patches = sorted(summaries, key=lambda row: int(row["patch_id"])) max_percent = max(float(row["seismic_moment_percent_of_total"]) for row in patches) fig, axes = plt.subplots(1, len(patches), figsize=(3.45 * len(patches), 3.65), squeeze=False) axes = axes[0] for ax, row in zip(axes, patches): percent = float(row["seismic_moment_percent_of_total"]) radius = 0.38 + 0.57 * math.sqrt(percent / max_percent) if max_percent > 0.0 else 0.65 _draw_beachball(ax, row["strike_deg"], row["dip_deg"], row["mean_rake_deg"], radius) ax.text(0.0, 1.08, f"{percent:.2f}% of total M0\nMw {float(row['moment_magnitude_Mw']):.2f}", ha="center", va="bottom", fontsize=11, fontweight="bold", linespacing=1.15) ax.text(0.0, -1.08, f"Patch {int(row['patch_id'])} - Segment {int(row['segment_id'])}\n" f"{float(row['strike_deg']):.0f} deg / {float(row['dip_deg']):.0f} deg / " f"{float(row['mean_rake_deg']):.1f} deg", ha="center", va="top", fontsize=8.5, linespacing=1.15) fig.suptitle("FRW patch focal mechanisms", fontsize=14, y=0.985) fig.text(0.5, 0.018, "Beach-ball area scales with M0 percentage; labels below show strike / dip / moment-weighted mean rake.", ha="center", fontsize=8.5) fig.tight_layout(rect=[0.0, 0.075, 1.0, 0.93]) fig.savefig(output, format="pdf", bbox_inches="tight", pad_inches=0.08) plt.close(fig) print(f"# wrote FRW patch focal mechanisms: {output}") def _frw_patch_physical_metrics(patch, mask, outputs, args): """Return scalar moment and slip metrics, consistent with global summaries.""" combined = outputs["combined_slip_dist_latlon"] seg = combined[combined["segment_id"].astype(int) == int(patch["segment_id"])].copy() if seg.empty: fail(f"no combined output rows for FRW patch segment {patch['segment_id']}") x_index = pd.to_numeric(seg["x_index"], errors="coerce").to_numpy() z_index = pd.to_numeric(seg["z_index"], errors="coerce").to_numpy() valid_index = ( np.isfinite(x_index) & np.isfinite(z_index) & (x_index >= 0) & (x_index < mask.shape[1]) & (z_index >= 0) & (z_index < mask.shape[0]) ) selected = np.zeros(len(seg), dtype=bool) selected[valid_index] = mask[ z_index[valid_index].astype(int), x_index[valid_index].astype(int) ] patch_rows = seg.loc[selected].copy() if patch_rows.empty: fail(f"no combined output rows selected for FRW patch {patch['patch_id']}") moment = ( pd.to_numeric(patch_rows["moment_mag"], errors="coerce") .replace([np.inf, -np.inf], np.nan).fillna(0.0).abs() ) moment_sum_input_units = float(moment.sum()) seismic_moment_dyne_cm = moment_sum_input_units * args.param_M0 seismic_moment_N_m = dyne_cm_to_N_m(seismic_moment_dyne_cm) moment_magnitude_Mw = moment_magnitude_from_N_m(seismic_moment_N_m) total_moment = ( pd.to_numeric(combined["moment_mag"], errors="coerce") .replace([np.inf, -np.inf], np.nan).fillna(0.0).abs() ) total_moment_sum_input_units = float(total_moment.sum()) total_seismic_moment_dyne_cm = total_moment_sum_input_units * args.param_M0 total_seismic_moment_N_m = dyne_cm_to_N_m(total_seismic_moment_dyne_cm) total_moment_magnitude_Mw = moment_magnitude_from_N_m(total_seismic_moment_N_m) seismic_moment_percent_of_total = ( 100.0 * seismic_moment_dyne_cm / total_seismic_moment_dyne_cm if total_seismic_moment_dyne_cm > 0.0 else float("nan") ) mw_percent_of_total = ( 100.0 * moment_magnitude_Mw / total_moment_magnitude_Mw if np.isfinite(total_moment_magnitude_Mw) and total_moment_magnitude_Mw != 0.0 else float("nan") ) slip = ( pd.to_numeric(patch_rows["slip_mag"], errors="coerce") .replace([np.inf, -np.inf], np.nan).dropna() ) max_slip_cm = float(slip.max()) if not slip.empty else float("nan") rake = pd.to_numeric(patch_rows["rake_angle_deg"], errors="coerce").to_numpy(dtype=float) rake_weights = moment.to_numpy(dtype=float) rake_valid = np.isfinite(rake) & np.isfinite(rake_weights) & (rake_weights > 0.0) mean_rake_deg = _circular_mean_deg(rake[rake_valid], rake_weights[rake_valid]) median_rake_deg = _circular_median_deg(rake[rake_valid]) strike_deg, dip_deg = _segment_strike_dip(patch["segment_id"], args) return { "moment_sum_input_units": moment_sum_input_units, "seismic_moment_dyne_cm": seismic_moment_dyne_cm, "seismic_moment_N_m": seismic_moment_N_m, "moment_magnitude_Mw": moment_magnitude_Mw, "total_moment_sum_input_units": total_moment_sum_input_units, "total_seismic_moment_dyne_cm": total_seismic_moment_dyne_cm, "total_seismic_moment_N_m": total_seismic_moment_N_m, "total_moment_magnitude_Mw": total_moment_magnitude_Mw, "seismic_moment_percent_of_total": seismic_moment_percent_of_total, "mw_percent_of_total": mw_percent_of_total, "max_slip_cm": max_slip_cm, "mean_rake_deg": mean_rake_deg, "median_rake_deg": median_rake_deg, "rake_sample_count": int(rake_valid.sum()), "strike_deg": strike_deg, "dip_deg": dip_deg, "selected_subfault_count": int(len(patch_rows)), } def _write_frw_patch_info(patch, meta, source_windows, patch_windows, physical, args): output = patch["output"] info_path = Path(str(output) + ".info.txt") source_sum = sum(float(np.sum(g)) for grids in source_windows.values() for g in grids) retained_sum = sum(float(np.sum(g)) for grids in patch_windows.values() for g in grids) fraction = retained_sum / source_sum if source_sum > 0.0 else float("nan") lines = [ f"source_frw = {args.frw_file}", f"output_frw = {output}", f"patch_id = {patch['patch_id']}", f"segment_id = {patch['segment_id']}", f"requested_x_km = {patch['xmin_km']} {patch['xmax_km']}", f"requested_z_km = {patch['zmin_km']} {patch['zmax_km']}", f"grid_shape = {meta['znum']} {meta['xnum']}", f"grid_spacing_km = {meta['dx_km']} {meta['dz_km']}", f"selected_x_indices = {meta['x_index_min']} {meta['x_index_max']}", f"selected_z_indices = {meta['z_index_min']} {meta['z_index_max']}", f"realized_x_km = {meta['x_realized_min_km']} {meta['x_realized_max_km']}", f"realized_z_km = {meta['z_realized_min_km']} {meta['z_realized_max_km']}", f"spatial_cell_count = {meta['spatial_cell_count']}", f"num_time_windows = {args.num_time_windows_resolved}", f"source_frw_value_sum = {source_sum:.12g}", f"retained_frw_value_sum = {retained_sum:.12g}", f"retained_fraction = {fraction:.12g}", f"moment_sum_input_units = {physical['moment_sum_input_units']:.12g}", f"seismic_moment_dyne_cm = {physical['seismic_moment_dyne_cm']:.12g}", f"seismic_moment_N_m = {physical['seismic_moment_N_m']:.12g}", f"moment_magnitude_Mw = {physical['moment_magnitude_Mw']:.6g}", f"total_moment_magnitude_Mw = {physical['total_moment_magnitude_Mw']:.6g}", f"mw_percent_of_total = {physical['mw_percent_of_total']:.12g}", f"seismic_moment_percent_of_total = {physical['seismic_moment_percent_of_total']:.12g}", f"max_slip_cm = {physical['max_slip_cm']:.12g}", f"mean_rake_deg = {physical['mean_rake_deg']:.12g}", f"median_rake_deg = {physical['median_rake_deg']:.12g}", f"rake_sample_count = {physical['rake_sample_count']}", f"strike_deg = {physical['strike_deg']:.12g}", f"dip_deg = {physical['dip_deg']:.12g}", f"selected_subfault_count = {physical['selected_subfault_count']}", "readback_validation = passed", ] info_path.write_text("\n".join(lines) + "\n", encoding="utf-8") return { "patch_id": patch["patch_id"], "segment_id": patch["segment_id"], "output_frw": str(output), "info_file": str(info_path), "requested_xmin_km": patch["xmin_km"], "requested_xmax_km": patch["xmax_km"], "requested_zmin_km": patch["zmin_km"], "requested_zmax_km": patch["zmax_km"], "realized_xmin_km": meta["x_realized_min_km"], "realized_xmax_km": meta["x_realized_max_km"], "realized_zmin_km": meta["z_realized_min_km"], "realized_zmax_km": meta["z_realized_max_km"], "x_edge_min_km": meta["x_edge_min_km"], "x_edge_max_km": meta["x_edge_max_km"], "z_edge_min_km": meta["z_edge_min_km"], "z_edge_max_km": meta["z_edge_max_km"], "selected_x_count": meta["selected_x_count"], "selected_z_count": meta["selected_z_count"], "spatial_cell_count": meta["spatial_cell_count"], "dx_km": meta["dx_km"], "dz_km": meta["dz_km"], "source_frw_value_sum": source_sum, "retained_frw_value_sum": retained_sum, "retained_fraction": fraction, **physical, } def generate_numbered_frw_patches(args, layout, source_windows, outputs): patches = _parse_numbered_frw_patches(args, layout) if not patches: return [] print_frw_patch_coordinate_guide(layout, args) prepared = [] for patch in patches: mask, meta = _frw_patch_mask(patch, args) prepared.append({**patch, "mask": mask, "meta": meta}) _warn_for_overlapping_frw_patches(prepared) summaries = [] for patch in prepared: patched = _make_frw_patch_windows(source_windows, layout, patch, patch["mask"]) write_frw_file(patch["output"], patched, layout) _validate_written_frw_patch(patch["output"], patched, layout, patch, patch["mask"]) physical = _frw_patch_physical_metrics(patch, patch["mask"], outputs, args) summary = _write_frw_patch_info( patch, patch["meta"], source_windows, patched, physical, args ) summaries.append(summary) print( f"# wrote FRW patch {patch['patch_id']}: {patch['output']} " f"segment={patch['segment_id']} cells={patch['meta']['spatial_cell_count']} " f"realized_x={patch['meta']['x_realized_min_km']}..{patch['meta']['x_realized_max_km']} km " f"realized_z={patch['meta']['z_realized_min_km']}..{patch['meta']['z_realized_max_km']} km" ) print(f"# FRW patch {patch['patch_id']} coordinate selection") print(f"# requested center bounds = x {patch['xmin_km']:.3f}..{patch['xmax_km']:.3f}, z {patch['zmin_km']:.3f}..{patch['zmax_km']:.3f} km") print(f"# realized center bounds = x {patch['meta']['x_realized_min_km']:.3f}..{patch['meta']['x_realized_max_km']:.3f}, z {patch['meta']['z_realized_min_km']:.3f}..{patch['meta']['z_realized_max_km']:.3f} km") print(f"# realized cell-edge bounds = x {patch['meta']['x_edge_min_km']:.3f}..{patch['meta']['x_edge_max_km']:.3f}, z {patch['meta']['z_edge_min_km']:.3f}..{patch['meta']['z_edge_max_km']:.3f} km") print(f"# selected x indices = {patch['meta']['x_index_min']}..{patch['meta']['x_index_max']} ({patch['meta']['selected_x_count']} centers)") print(f"# selected z indices = {patch['meta']['z_index_min']}..{patch['meta']['z_index_max']} ({patch['meta']['selected_z_count']} centers)") print(f"# selected grid dimensions = {patch['meta']['selected_x_count']} x {patch['meta']['selected_z_count']}") print(f"# FRW patch {patch['patch_id']} source metrics") print(f"# moment sum input units = {physical['moment_sum_input_units']:.12g}") print(f"# seismic moment = {physical['seismic_moment_dyne_cm']:.6e} dyne-cm") print(f"# seismic moment = {physical['seismic_moment_N_m']:.6e} N-m") print(f"# moment magnitude = {physical['moment_magnitude_Mw']:.4f} Mw") print(f"# total moment magnitude = {physical['total_moment_magnitude_Mw']:.4f} Mw") print(f"# Mw percentage of total = {physical['mw_percent_of_total']:.4f} %") print(f"# M0 percentage of total = {physical['seismic_moment_percent_of_total']:.4f} %") print(f"# maximum slip = {physical['max_slip_cm']:.6f} cm") print(f"# mean rake (M0-weighted) = {physical['mean_rake_deg']:.4f} deg") print(f"# median rake (circular) = {physical['median_rake_deg']:.4f} deg") print(f"# rake samples = {physical['rake_sample_count']}") print(f"# selected subfaults = {physical['selected_subfault_count']}") pd.DataFrame(summaries).to_csv(args.frw_patch_summary, index=False) print(f"# wrote FRW patch summary: {args.frw_patch_summary}") plot_frw_patch_focal_mechanisms(summaries, args) return summaries def rupture_block_plan(rupture_config): rc = rupture_config if rc == "sing": return [("seg1_rake1", 1, 1)] if rc == "seg2": return [("seg1_rake1", 1, 1), ("seg2_rake1", 2, 1)] if "rakeval" in rc and not rc.startswith("seg2_rakeval"): return [("seg1_rake1", 1, 1), ("seg1_rake2", 1, 2)] if rc.startswith("seg2_rakeval"): return [ ("seg1_rake1", 1, 1), ("seg1_rake2", 1, 2), ("seg2_rake1", 2, 1), ("seg2_rake2", 2, 2), ] fail("unknown rupture_config: " + rc + "\nAllowed examples: sing, seg2, rakeval, rakeval_multime3, seg2_rakeval*") def segment_grid_params(segment_id, args): if segment_id == 1: xnum, znum = args.xgrid_num, args.zgrid_num xlen, zlen = args.x_length_km, args.z_length_km else: missing = [ name for name, value in [ ("--xgrid_num2", args.xgrid_num2), ("--zgrid_num2", args.zgrid_num2), ("--x_length_km2", args.x_length_km2), ("--z_length_km2", args.z_length_km2), ] if value is None ] if missing: fail( "segment 2 geometry must be specified explicitly; missing: " + ", ".join(missing) ) xnum, znum = args.xgrid_num2, args.zgrid_num2 xlen, zlen = args.x_length_km2, args.z_length_km2 if xnum <= 0 or znum <= 0: fail(f"grid numbers must be positive for segment {segment_id}: x={xnum}, z={znum}") if xlen <= 0.0 or zlen <= 0.0: fail(f"grid lengths must be positive for segment {segment_id}: x={xlen}, z={zlen}") return xnum, znum, xlen, zlen def add_local_grid_columns(block, segment_id, args): xnum, znum, xlen, zlen = segment_grid_params(segment_id, args) expected = xnum * znum if len(block) != expected: fail(f"cannot assign grid columns for segment {segment_id}: expected {expected} rows, got {len(block)} rows") x_grid_space_km = xlen / xnum z_grid_space_km = zlen / znum subfault_area_km2 = x_grid_space_km * z_grid_space_km subfault_area_cm2 = subfault_area_km2 * ONE_KM2_CM2 local_index = np.arange(len(block)) block["local_index"] = local_index block["x_index"] = local_index % xnum block["z_index"] = local_index // xnum block["x_grid_space_km"] = x_grid_space_km block["z_grid_space_km"] = z_grid_space_km block["subfault_area_km2"] = subfault_area_km2 block["subfault_area_cm2"] = subfault_area_cm2 block["x_km"] = block["x_index"] * x_grid_space_km # Offset grid z by the shallowest subfault depth so z_km is on an absolute-depth-like scale. z_origin_depth_km = float(block["depth"].min()) block["z_origin_depth_km"] = z_origin_depth_km block["z_km"] = z_origin_depth_km + block["z_index"] * z_grid_space_km return block def compute_mu_and_slip_from_moment(block, args): if (block["depth"] <= 0.0).any(): bad = block[block["depth"] <= 0.0].head(5) fail("depth must be positive to compute mu from log(depth). " f"Example bad rows:\n{bad.to_string(index=False)}") block["mu_pa_from_depth"] = mu_pa_from_depth(block["depth"].to_numpy(dtype=float)) block["mu_gpa_from_depth"] = block["mu_pa_from_depth"] / 1.0e9 block["mu_cgs_from_depth"] = block["mu_pa_from_depth"] * 10.0 block["slip_cm_from_moment"] = ( block["moment"] * args.param_M0 / block["mu_cgs_from_depth"] / block["subfault_area_cm2"] ) block["mu_gpa_from_input_slip"] = mu_gpa_from_moment_and_slip( block["moment"], block["slip"], block["subfault_area_cm2"], args.param_M0 ) return block def apply_slip_source(block, args): block["slip_input"] = block["slip"] block = compute_mu_and_slip_from_moment(block, args) if args.use_moment_for_slip: block["slip"] = block["slip_cm_from_moment"] block["slip_source"] = "moment" else: block["slip_source"] = "input" block["mu_gpa_from_active_slip"] = mu_gpa_from_moment_and_slip( block["moment"], block["slip"], block["subfault_area_cm2"], args.param_M0 ) return block def split_blocks(df, args): blocks = {} cursor = 0 for block_name, segment_id, rake_id in rupture_block_plan(args.rupture_config): xnum, znum, _, _ = segment_grid_params(segment_id, args) n = xnum * znum start, stop = cursor, cursor + n if stop > len(df): fail(f"not enough rows for block {block_name}: need rows {start}:{stop}, but file has {len(df)} rows") block = df.iloc[start:stop].copy().reset_index(drop=True) block["block_name"] = block_name block["segment_id"] = segment_id block["rake_id"] = rake_id block = add_local_grid_columns(block, segment_id, args) block = apply_slip_source(block, args) blocks[block_name] = block cursor = stop if cursor != len(df): if args.allow_extra_rows: print(f"# WARNING: file has {len(df) - cursor} extra rows after expected blocks. They are ignored.") else: fail(f"file has extra rows after expected blocks: expected {cursor}, got {len(df)}. Use --allow_extra_rows to ignore extras.") return blocks def assert_same_geometry(block1, block2, label, tol): if len(block1) != len(block2): fail(f"geometry mismatch for {label}: row counts differ ({len(block1)} vs {len(block2)})") for col in ["lon", "lat", "depth", "x_index", "z_index", "x_km", "z_km"]: max_abs = (block1[col] - block2[col]).abs().max() if max_abs > tol: fail(f"geometry mismatch for {label}: max |delta {col}| = {max_abs}, tolerance={tol}") def compute_slip_mag_az_from_rake_blocks(block1, block2, args, label): assert_same_geometry(block1, block2, label=label, tol=args.geometry_tol) theta1 = math.radians(args.theta1_deg) theta2 = math.radians(args.theta2_deg) rows = [] for row1, row2 in zip(block1.itertuples(index=False), block2.itertuples(index=False)): active_slip1, active_slip2 = float(row1.slip), float(row2.slip) input_slip1, input_slip2 = float(row1.slip_input), float(row2.slip_input) m1, m2 = float(row1.moment), float(row2.moment) area1, area2 = float(row1.subfault_area_cm2), float(row2.subfault_area_cm2) slip_mag = math.sqrt(active_slip1 * active_slip1 + active_slip2 * active_slip2) moment_mag = math.sqrt(m1 * m1 + m2 * m2) if slip_mag == 0.0: rake_angle_deg = float("nan") else: c = clamp((active_slip1 * math.cos(theta1) + active_slip2 * math.cos(theta2)) / slip_mag, -1.0, 1.0) rake_angle_deg = math.degrees(math.atan2(math.sqrt(max(0.0, 1.0 - c * c)), c)) * args.theta_sign rows.append({ "lon": row1.lon, "lat": row1.lat, "depth": row1.depth, "x_index": int(row1.x_index), "z_index": int(row1.z_index), "x_km": row1.x_km, "z_km": row1.z_km, "x_grid_space_km": row1.x_grid_space_km, "z_grid_space_km": row1.z_grid_space_km, "subfault_area_km2": row1.subfault_area_km2, "subfault_area_cm2": row1.subfault_area_cm2, "mu_gpa_rake1": mu_gpa_from_moment_and_slip([m1], [input_slip1], [area1], args.param_M0)[0], "mu_gpa_rake2": mu_gpa_from_moment_and_slip([m2], [input_slip2], [area2], args.param_M0)[0], "mu_gpa_from_depth_rake1": row1.mu_gpa_from_depth, "mu_gpa_from_depth_rake2": row2.mu_gpa_from_depth, "mu_gpa_from_active_slip_rake1": mu_gpa_from_moment_and_slip([m1], [active_slip1], [area1], args.param_M0)[0], "mu_gpa_from_active_slip_rake2": mu_gpa_from_moment_and_slip([m2], [active_slip2], [area2], args.param_M0)[0], "slip_source": row1.slip_source, "slip_mag": slip_mag, "rake_angle_deg": rake_angle_deg, "moment_mag": moment_mag, "slip_rake1": active_slip1, "slip_rake2": active_slip2, "slip_input_rake1": input_slip1, "slip_input_rake2": input_slip2, "slip_cm_from_moment_rake1": row1.slip_cm_from_moment, "slip_cm_from_moment_rake2": row2.slip_cm_from_moment, "moment_rake1": m1, "moment_rake2": m2, "segment_id": int(row1.segment_id), "local_index": int(row1.local_index), }) return pd.DataFrame(rows) def single_rake_block_to_slip(block): return pd.DataFrame({ "lon": block["lon"], "lat": block["lat"], "depth": block["depth"], "x_index": block["x_index"].astype(int), "z_index": block["z_index"].astype(int), "x_km": block["x_km"], "z_km": block["z_km"], "x_grid_space_km": block["x_grid_space_km"], "z_grid_space_km": block["z_grid_space_km"], "subfault_area_km2": block["subfault_area_km2"], "subfault_area_cm2": block["subfault_area_cm2"], "mu_gpa_rake1": block["mu_gpa_from_input_slip"], "mu_gpa_rake2": float("nan"), "mu_gpa_from_depth_rake1": block["mu_gpa_from_depth"], "mu_gpa_from_depth_rake2": float("nan"), "mu_gpa_from_active_slip_rake1": block["mu_gpa_from_active_slip"], "mu_gpa_from_active_slip_rake2": float("nan"), "slip_source": block["slip_source"], "slip_mag": block["slip"], "rake_angle_deg": float("nan"), "moment_mag": block["moment"].abs(), "slip_rake1": block["slip"], "slip_rake2": float("nan"), "slip_input_rake1": block["slip_input"], "slip_input_rake2": float("nan"), "slip_cm_from_moment_rake1": block["slip_cm_from_moment"], "slip_cm_from_moment_rake2": float("nan"), "moment_rake1": block["moment"], "moment_rake2": float("nan"), "segment_id": block["segment_id"].astype(int), "local_index": block["local_index"].astype(int), }) def compute_outputs(blocks, args): rc = args.rupture_config outputs = {} if rc == "sing": outputs["seg1_slip"] = single_rake_block_to_slip(blocks["seg1_rake1"]) elif rc == "seg2": outputs["seg1_slip"] = single_rake_block_to_slip(blocks["seg1_rake1"]) outputs["seg2_slip"] = single_rake_block_to_slip(blocks["seg2_rake1"]) elif "rakeval" in rc and not rc.startswith("seg2_rakeval"): outputs["seg1_slip_mag_az"] = compute_slip_mag_az_from_rake_blocks(blocks["seg1_rake1"], blocks["seg1_rake2"], args, "seg1 rake pair") elif rc.startswith("seg2_rakeval"): outputs["seg1_slip_mag_az"] = compute_slip_mag_az_from_rake_blocks(blocks["seg1_rake1"], blocks["seg1_rake2"], args, "seg1 rake pair") outputs["seg2_slip_mag_az"] = compute_slip_mag_az_from_rake_blocks(blocks["seg2_rake1"], blocks["seg2_rake2"], args, "seg2 rake pair") else: fail(f"unsupported rupture_config in compute_outputs: {rc}") outputs["combined_slip_dist_latlon"] = pd.concat(outputs.values(), ignore_index=True) return outputs def compute_summary_metrics(outputs, args): combined = outputs["combined_slip_dist_latlon"] moment_mag = pd.to_numeric(combined["moment_mag"], errors="coerce").replace([np.inf, -np.inf], np.nan).fillna(0.0) total_moment_input_units = float(moment_mag.abs().sum()) total_M0_dyne_cm = total_moment_input_units * args.param_M0 total_M0_N_m = dyne_cm_to_N_m(total_M0_dyne_cm) Mw = moment_magnitude_from_N_m(total_M0_N_m) slip = pd.to_numeric(combined["slip_mag"], errors="coerce").replace([np.inf, -np.inf], np.nan).dropna() if slip.empty: max_slip_cm = min_slip_for_stats_cm = mean_slip_cm = median_slip_cm = float("nan") n_slip_total = n_slip_used = 0 else: max_slip_cm = float(slip.max()) min_slip_for_stats_cm = max_slip_cm * args.min_slip_thresh slip_for_stats = slip[slip >= min_slip_for_stats_cm] n_slip_total = int(slip.size) n_slip_used = int(slip_for_stats.size) mean_slip_cm = float(slip_for_stats.mean()) if n_slip_used else float("nan") median_slip_cm = float(slip_for_stats.median()) if n_slip_used else float("nan") return { "param_M0": args.param_M0, "total_moment_input_units": total_moment_input_units, "total_M0_dyne_cm": total_M0_dyne_cm, "total_M0_N_m": total_M0_N_m, "Mw": Mw, "max_slip_cm": max_slip_cm, "min_slip_thresh": args.min_slip_thresh, "min_slip_for_stats_cm": min_slip_for_stats_cm, "mean_slip_cm": mean_slip_cm, "median_slip_cm": median_slip_cm, "n_slip_total": n_slip_total, "n_slip_used": n_slip_used, } def summary_text(metrics): return ( f"M0 = {metrics['total_M0_dyne_cm']:.3e} dyne-cm\n" f"M0 = {metrics['total_M0_N_m']:.3e} N-m, Mw = {metrics['Mw']:.2f}\n" f"Max slip = {metrics['max_slip_cm']:.1f} cm\n" f"Mean/Median slip = {metrics['mean_slip_cm']:.1f}/{metrics['median_slip_cm']:.1f} cm\n" f"Slip stats threshold = {metrics['min_slip_for_stats_cm']:.1f} cm ({metrics['min_slip_thresh']:.2g}*max)" ) def add_summary_box(fig, metrics, x=0.5, y=0.985): """Place summary metrics at the top center without a surrounding box.""" fig.text(x, y, summary_text(metrics), ha="center", va="top", fontsize=8) def strike_for_segment(seg_id, args): if int(seg_id) == 1: return float(args.strike1_deg) if int(seg_id) == 2: return float(args.strike2_deg) return float(args.strike1_deg) def x_axis_flip_for_segment(seg_id, args): if int(seg_id) == 1: return bool(args.x_axis_flip) if int(seg_id) == 2: return bool(args.x_axis_flip2) return bool(args.x_axis_flip) def data_span_km(seg): xpad = max(float(seg["x_grid_space_km"].median()) * 0.75, 0.5) zpad = max(float(seg["z_grid_space_km"].median()) * 0.75, 0.5) xspan = float(seg["x_km"].max() - seg["x_km"].min()) + 2.0 * xpad zspan = float(seg["z_km"].max() - seg["z_km"].min()) + 2.0 * zpad return max(xspan, 1.0), max(zspan, 1.0), xpad, zpad def effective_repeating_value(args, repeating_name, seis_name): val = getattr(args, repeating_name) return getattr(args, seis_name) if val is None else val def lonlat_to_local_en_km(lon, lat, lon0, lat0): lon = np.asarray(lon, dtype=float) lat = np.asarray(lat, dtype=float) lat0_rad = math.radians(float(lat0)) east_km = (lon - float(lon0)) * 111.32 * math.cos(lat0_rad) north_km = (lat - float(lat0)) * 111.32 return east_km, north_km def seismicity_strike_for_segment(seg_id, args): if int(seg_id) == 1: return float(args.seis_strike1_deg) if int(seg_id) == 2: return float(args.seis_strike2_deg) return float(args.seis_strike1_deg) def dip_for_segment(seg_id, args): if int(seg_id) == 1: return float(args.dip1_deg) if int(seg_id) == 2: return float(args.dip2_deg) return float(args.dip1_deg) def segment_trace_origin_lonlat(seg): top = seg[seg["z_index"].astype(int) == int(seg["z_index"].min())].copy() if top.empty: top = seg.copy() return float(top["lon"].mean()), float(top["lat"].mean()) def segment_top_depth_km(seg): if "z_origin_depth_km" in seg.columns: return float(seg["z_origin_depth_km"].iloc[0]) return float(seg["depth"].min()) def add_along_across_for_segment(eq_df, seg, seg_id, args): lon0, lat0 = segment_trace_origin_lonlat(seg) east_km, north_km = lonlat_to_local_en_km(eq_df["lon"], eq_df["lat"], lon0, lat0) az = math.radians(seismicity_strike_for_segment(seg_id, args)) x_along_km = east_km * math.sin(az) + north_km * math.cos(az) y_across_km = east_km * math.cos(az) - north_km * math.sin(az) out = eq_df.copy() out["seis_segment_id"] = int(seg_id) out["seis_x_along_km"] = x_along_km out["seis_y_across_km"] = y_across_km out["seis_abs_across_km"] = np.abs(y_across_km) out["seis_strike_deg"] = seismicity_strike_for_segment(seg_id, args) out["seis_origin_lon"] = lon0 out["seis_origin_lat"] = lat0 return out def event_identity_columns(eq_df): """Return stable columns used to identify one catalog event.""" if "id" in eq_df.columns and eq_df["id"].notna().any(): return ["id"] if "repeating_evid" in eq_df.columns: return ["csid", "repeating_evid"] return ["lat", "lon", "dep", "mag"] def build_event_segment_diagnostics(eq_df, outputs, args, label="seismicity"): """Preserve and diagnose every event x segment candidate. The table records the corridor test, the current nearest-line assignment, the independent projection onto every segment, panel-bound tests, and the final condition required for the event to be drawn on that segment. """ if eq_df is None or eq_df.empty: return pd.DataFrame() combined = outputs["combined_slip_dist_latlon"] key_cols = event_identity_columns(eq_df) candidates = [] for seg_id in sorted(combined["segment_id"].dropna().astype(int).unique()): seg = combined[combined["segment_id"].astype(int) == seg_id].copy() cand = add_along_across_for_segment(eq_df, seg, seg_id, args) projected = project_events_to_segment(cand, seg, args) _, _, xpad, zpad = data_span_km(seg) panel_xmin = float(seg["x_km"].min()) - xpad - args.seis_plot_buffer_km panel_xmax = float(seg["x_km"].max()) + xpad + args.seis_plot_buffer_km panel_zmin = float(seg["z_km"].min()) - zpad - args.seis_plot_buffer_km panel_zmax = float(seg["z_km"].max()) + zpad + args.seis_plot_buffer_km projected["diagnostic_catalog"] = label projected["candidate_segment_id"] = int(seg_id) projected["across_filter_enabled"] = bool(args.seis_across_filter) projected["passed_corridor"] = ( True if not args.seis_across_filter else projected["seis_abs_across_km"] <= args.seis_max_across_km ) projected["panel_xmin_km"] = panel_xmin projected["panel_xmax_km"] = panel_xmax projected["panel_zmin_km"] = panel_zmin projected["panel_zmax_km"] = panel_zmax projected["projected_inside_x"] = ( (projected["x_km"] >= panel_xmin) & (projected["x_km"] <= panel_xmax) ) projected["projected_inside_z"] = ( (projected["z_km"] >= panel_zmin) & (projected["z_km"] <= panel_zmax) ) projected["remained_inside_projected_panel"] = ( projected["projected_inside_x"] & projected["projected_inside_z"] ) candidates.append(projected) diagnostics = pd.concat(candidates, ignore_index=True) diagnostics["segment_assignment_mode"] = args.seis_segment_assignment_mode diagnostics["survived_assignment"] = False passed = diagnostics[diagnostics["passed_corridor"]].copy() if not passed.empty: if args.seis_segment_assignment_mode == "all_candidates": diagnostics.loc[passed.index, "survived_assignment"] = True else: winner_index = ( passed.sort_values( key_cols + ["seis_abs_across_km", "candidate_segment_id"], kind="mergesort", ) .drop_duplicates(subset=key_cols, keep="first") .index ) diagnostics.loc[winner_index, "survived_assignment"] = True diagnostics["would_plot_on_segment"] = ( diagnostics["passed_corridor"] & diagnostics["survived_assignment"] & diagnostics["remained_inside_projected_panel"] ) diagnostics["diagnostic_reason"] = np.select( [ ~diagnostics["passed_corridor"], diagnostics["passed_corridor"] & ~diagnostics["survived_assignment"], diagnostics["survived_assignment"] & ~diagnostics["projected_inside_x"], diagnostics["survived_assignment"] & diagnostics["projected_inside_x"] & ~diagnostics["projected_inside_z"], diagnostics["would_plot_on_segment"], ], [ "failed_across_corridor", "lost_nearest_segment_assignment", "assigned_but_outside_projected_x", "assigned_but_outside_projected_z", "plotted", ], default="not_plotted_other", ) ordered = key_cols + [ "candidate_segment_id", "passed_corridor", "survived_assignment", "projected_inside_x", "projected_inside_z", "remained_inside_projected_panel", "would_plot_on_segment", "diagnostic_reason", "seis_abs_across_km", "seis_x_along_km", "seis_y_across_km", "x_km", "z_km", "panel_xmin_km", "panel_xmax_km", "panel_zmin_km", "panel_zmax_km", "seis_origin_lon", "seis_origin_lat", "seis_strike_deg", "projection_mode", "projection_depth_col", "projection_dip_deg", "projection_y_right_km", "projection_dz_from_top_km", "projection_down_dip_distance_km", ] ordered += [col for col in diagnostics.columns if col not in ordered] diagnostics = diagnostics[ordered] print(f"# {label} event x segment diagnostics") print(f"# source events = {len(eq_df)}") print(f"# candidate rows = {len(diagnostics)}") for (seg_id, reason), count in diagnostics.groupby( ["candidate_segment_id", "diagnostic_reason"], sort=True ).size().items(): print(f"# segment {int(seg_id)} {reason:<38s} = {int(count)}") return diagnostics def write_event_segment_diagnostics(eq_df, outputs, args, path, label): if not path: return None diagnostics = build_event_segment_diagnostics(eq_df, outputs, args, label=label) diagnostics.to_csv(path, index=False) print(f"# wrote {path}") return diagnostics def filter_events_by_across_distance(eq_df, outputs, args, label="seismicity"): if eq_df is None or eq_df.empty: return eq_df combined = outputs["combined_slip_dist_latlon"] candidates = [] for seg_id in sorted(combined["segment_id"].dropna().astype(int).unique()): seg = combined[combined["segment_id"].astype(int) == seg_id].copy() tmp = add_along_across_for_segment(eq_df, seg, seg_id, args) if args.seis_across_filter: tmp = tmp[tmp["seis_abs_across_km"] <= args.seis_max_across_km].copy() print(f"# {label} across filter segment {seg_id}: strike={seismicity_strike_for_segment(seg_id, args):g} deg, dip={dip_for_segment(seg_id, args):g} deg, within +/-{args.seis_max_across_km:g} km = {len(tmp)}") else: print(f"# {label} across filter segment {seg_id}: DISABLED; all {len(tmp)} events are assignment candidates") if not tmp.empty: candidates.append(tmp) if not candidates: if args.seis_across_filter: fail( f"no {label} events passed the +/-{args.seis_max_across_km:g} km " "across-fault corridor for any segment" ) fail(f"no {label} events are available for segment assignment") near = pd.concat(candidates, ignore_index=True) key_cols = event_identity_columns(near) if args.seis_segment_assignment_mode == "nearest": near = ( near.sort_values( key_cols + ["seis_abs_across_km", "seis_segment_id"], kind="mergesort", ) .drop_duplicates(subset=key_cols, keep="first") .reset_index(drop=True) ) else: near = near.sort_values( key_cols + ["seis_segment_id"], kind="mergesort" ).reset_index(drop=True) print(f"# {label} segment assignment") print(f"# assignment_mode = {args.seis_segment_assignment_mode}") print(f"# across_filter_enabled = {args.seis_across_filter}") print(f"# max_across_km = {args.seis_max_across_km:g}" if args.seis_across_filter else "# max_across_km = disabled") print(f"# retained event-segment rows = {len(near)}") return near def filter_seismicity_by_across_distance(eq_df, outputs, args): return filter_events_by_across_distance(eq_df, outputs, args, label="seismicity") def filter_repeating_by_across_distance(eq_df, outputs, args): return filter_events_by_across_distance(eq_df, outputs, args, label="repeating earthquakes") def fit_fault_plane_projection(seg): use = seg[["lon", "lat", "depth", "x_km"]].replace([np.inf, -np.inf], np.nan).dropna() if len(use) < 4: fail("not enough finite subfault points to fit event projection") A = np.column_stack([use["lon"].to_numpy(float), use["lat"].to_numpy(float), use["depth"].to_numpy(float), np.ones(len(use))]) coef_x, _, _, _ = np.linalg.lstsq(A, use["x_km"].to_numpy(float), rcond=None) return coef_x def project_events_to_segment(eq_df, seg, args): seg_id = int(seg["segment_id"].iloc[0]) depth_col = "dep_model" if "dep_model" in eq_df.columns else "dep" coef_x = fit_fault_plane_projection(seg) Aeq = np.column_stack([eq_df["lon"].to_numpy(float), eq_df["lat"].to_numpy(float), eq_df[depth_col].to_numpy(float), np.ones(len(eq_df))]) lon0, lat0 = segment_trace_origin_lonlat(seg) east_km, north_km = lonlat_to_local_en_km(eq_df["lon"], eq_df["lat"], lon0, lat0) az = math.radians(seismicity_strike_for_segment(seg_id, args)) y_right_km = east_km * math.cos(az) - north_km * math.sin(az) dip_deg = dip_for_segment(seg_id, args) dip_rad = math.radians(dip_deg) sin_dip = math.sin(dip_rad) if abs(sin_dip) < 1.0e-8: fail(f"dip angle is too close to 0 deg for depth-to-dip projection: segment {seg_id}, dip={dip_deg}") top_depth = segment_top_depth_km(seg) dz_km = eq_df[depth_col].to_numpy(float) - top_depth if args.ignore_across_strike_distance_depth_projection: s_dip_km = dz_km / sin_dip projection_mode = "depth_only_collapse" else: s_dip_km = y_right_km * math.cos(dip_rad) + dz_km * sin_dip projection_mode = "full_across_strike_projection" out = eq_df.copy() out["x_km"] = Aeq @ coef_x out["z_km"] = top_depth + s_dip_km out["projection_depth_col"] = depth_col out["projection_dip_deg"] = dip_deg out["projection_mode"] = projection_mode out["projection_y_right_km"] = y_right_km out["projection_dz_from_top_km"] = dz_km out["projection_down_dip_distance_km"] = s_dip_km return out def target_event_dataframe(args): """Return the CLI target hypocenter as a one-row model-coordinate table.""" values = [ args.target_event_lat, args.target_event_lon, args.target_event_model_depth_km, ] if not all(np.isfinite(value) for value in values): fail("target-event latitude, longitude, and model depth must be finite") if args.target_event_model_depth_km < 0.0: fail("--target_event_model_depth_km must be non-negative") return pd.DataFrame({ "lat": [float(args.target_event_lat)], "lon": [float(args.target_event_lon)], "dep": [float(args.target_event_model_depth_km)], "dep_model": [float(args.target_event_model_depth_km)], "catalog_type": ["target_event"], }) def _closest_point_on_polyline_xy(point_x, point_y, line_x, line_y, line_value): """Return closest polyline point, interpolated value, and distance.""" line_x = np.asarray(line_x, dtype=float) line_y = np.asarray(line_y, dtype=float) line_value = np.asarray(line_value, dtype=float) if len(line_x) == 0: fail("cannot project target onto an empty model-depth line") if len(line_x) == 1: distance = math.hypot(point_x - line_x[0], point_y - line_y[0]) return line_x[0], line_y[0], line_value[0], distance best = None for index in range(len(line_x) - 1): x0, y0 = line_x[index], line_y[index] x1, y1 = line_x[index + 1], line_y[index + 1] dx, dy = x1 - x0, y1 - y0 length2 = dx * dx + dy * dy if length2 <= 0.0: fraction = 0.0 else: fraction = ((point_x - x0) * dx + (point_y - y0) * dy) / length2 fraction = clamp(fraction, 0.0, 1.0) closest_x = x0 + fraction * dx closest_y = y0 + fraction * dy distance = math.hypot(point_x - closest_x, point_y - closest_y) value = line_value[index] + fraction * ( line_value[index + 1] - line_value[index] ) candidate = (distance, closest_x, closest_y, value) if best is None or candidate[0] < best[0]: best = candidate return best[1], best[2], best[3], best[0] def target_event_surface_projection_for_segment(target, seg, args): """Project target lon/lat at a model-z or vertical-depth surface.""" target_value = float(target["dep_model"].iloc[0]) if args.target_event_depth_mode == "vertical_depth": coordinate = "depth" tolerance = max(float(seg["depth"].max() - seg["depth"].min()) * 1.0e-9, 1.0e-8) else: coordinate = "z_km" tolerance = max(float(seg["z_grid_space_km"].median()) * 1.0e-6, 1.0e-8) coord_min = float(seg[coordinate].min()) coord_max = float(seg[coordinate].max()) inside_z = coord_min - tolerance <= target_value <= coord_max + tolerance line_rows = [] for x_index, column in seg.groupby("x_index", sort=True): column = column.sort_values(coordinate).drop_duplicates(subset=[coordinate]) values = column[coordinate].to_numpy(float) if target_value < values.min() - tolerance or target_value > values.max() + tolerance: continue line_rows.append({ "x_index": int(x_index), "x_km": float(np.interp(target_value, values, column["x_km"].to_numpy(float))), "z_km": float(np.interp(target_value, values, column["z_km"].to_numpy(float))), "lon": float(np.interp(target_value, values, column["lon"].to_numpy(float))), "lat": float(np.interp(target_value, values, column["lat"].to_numpy(float))), }) if not line_rows: return { "x_km": float("nan"), "z_km": target_value, "surface_location_residual_km": float("inf"), "inside_z_range": False, "inside_x_range": False, "nearest_surface_lon": float("nan"), "nearest_surface_lat": float("nan"), "constant_depth_line_points": 0, } line = pd.DataFrame(line_rows).sort_values("x_km") reference_lon = float(line["lon"].mean()) reference_lat = float(line["lat"].mean()) line_east, line_north = lonlat_to_local_en_km( line["lon"], line["lat"], reference_lon, reference_lat ) target_east, target_north = lonlat_to_local_en_km( [float(target["lon"].iloc[0])], [float(target["lat"].iloc[0])], reference_lon, reference_lat, ) closest_east, closest_north, target_x, residual = _closest_point_on_polyline_xy( float(target_east[0]), float(target_north[0]), line_east, line_north, line["x_km"].to_numpy(float), ) nearest_lon = reference_lon + closest_east / ( 111.32 * math.cos(math.radians(reference_lat)) ) nearest_lat = reference_lat + closest_north / 111.32 xmin = float(seg["x_km"].min()) xmax = float(seg["x_km"].max()) x_tolerance = max(float(seg["x_grid_space_km"].median()) * 1.0e-6, 1.0e-8) inside_x = target_x >= xmin - x_tolerance and target_x <= xmax + x_tolerance target_z = float(np.interp(target_x, line["x_km"], line["z_km"])) return { "x_km": float(target_x), "z_km": target_z, "surface_location_residual_km": float(residual), "inside_z_range": bool(inside_z), "inside_x_range": bool(inside_x), "nearest_surface_lon": float(nearest_lon), "nearest_surface_lat": float(nearest_lat), "constant_depth_line_points": int(len(line)), } def target_event_projection_by_segment(outputs, args): """Project target using lon/lat plus an existing fault-plane z coordinate.""" target = target_event_dataframe(args) combined = outputs["combined_slip_dist_latlon"] rows = [] for seg_id in sorted(combined["segment_id"].dropna().astype(int).unique()): seg = combined[combined["segment_id"].astype(int) == seg_id].copy() projection = target_event_surface_projection_for_segment(target, seg, args) projection["segment_id"] = int(seg_id) rows.append(projection) result = pd.DataFrame(rows) if result.empty: fail("no fault segments available for target-event projection") requested_mode = args.target_event_segment_mode mode = "nearest_surface" if requested_mode == "nearest_plane" else requested_mode available = set(result["segment_id"].astype(int)) valid = result["inside_z_range"] & result["inside_x_range"] & np.isfinite( result["surface_location_residual_km"] ) if mode == "nearest_surface": candidates = result[valid].copy() if candidates.empty: fail( "target model depth is outside all segment ranges or no valid " "constant-model-depth surface line could be constructed" ) selected_id = int( candidates.sort_values( ["surface_location_residual_km", "segment_id"] ).iloc[0]["segment_id"] ) result["selected"] = result["segment_id"].astype(int) == selected_id elif mode == "all": result["selected"] = valid elif mode == "segment1": if 1 not in available: fail("--target_event_segment_mode segment1 requested but unavailable") result["selected"] = (result["segment_id"].astype(int) == 1) & valid else: if 2 not in available: fail("--target_event_segment_mode segment2 requested but unavailable") result["selected"] = (result["segment_id"].astype(int) == 2) & valid segment_count = int(result["segment_id"].nunique()) print("# target-event cross-section projection") print(f"# fault segment count = {segment_count}") print( f"# input lon/lat/model_depth = {args.target_event_lon:g} / " f"{args.target_event_lat:g} / {args.target_event_model_depth_km:g} km" ) depth_description = ( "vertical hypocentral depth; interpolated onto fault geometry" if args.target_event_depth_mode == "vertical_depth" else "fault-plane z; no dip conversion" ) print(f"# target depth interpretation= {depth_description}") print(f"# segment mode = {requested_mode}") for row in result.itertuples(index=False): residual_text = ( f"{row.surface_location_residual_km:.6g}" if np.isfinite(row.surface_location_residual_km) else "invalid" ) print( f"# segment {int(row.segment_id)}: " f"x={row.x_km:.6g} km z={row.z_km:.6g} km " f"surface_location_residual={residual_text} km " f"inside_x={bool(row.inside_x_range)} " f"inside_z={bool(row.inside_z_range)} " f"line_points={int(row.constant_depth_line_points)} " f"selected={bool(row.selected)}" ) return result def overlay_target_event_on_axis(ax, seg, target_projection, args): """Draw the selected projected target hypocenter on one segment panel.""" if target_projection is None or target_projection.empty: return None, 0 seg_id = int(seg["segment_id"].iloc[0]) selected = target_projection[ (target_projection["segment_id"].astype(int) == seg_id) & target_projection["selected"].astype(bool) ].copy() if selected.empty: return None, 0 if args.target_event_cross_clip: _, _, xpad, zpad = data_span_km(seg) xmin = float(seg["x_km"].min()) - xpad - args.seis_plot_buffer_km xmax = float(seg["x_km"].max()) + xpad + args.seis_plot_buffer_km zmin = float(seg["z_km"].min()) - zpad - args.seis_plot_buffer_km zmax = float(seg["z_km"].max()) + zpad + args.seis_plot_buffer_km selected = selected[ (selected["x_km"] >= xmin) & (selected["x_km"] <= xmax) & (selected["z_km"] >= zmin) & (selected["z_km"] <= zmax) ] if selected.empty: print(f"# target hypocenter selected for segment {seg_id} but outside panel") return None, 0 row = selected.iloc[0] handle = ax.scatter( [row["x_km"]], [row["z_km"]], s=args.target_event_cross_size, marker=args.target_event_cross_marker, facecolors=args.target_event_cross_facecolor, edgecolors=args.target_event_cross_edgecolor, linewidths=args.target_event_cross_linewidth, alpha=args.target_event_cross_alpha, zorder=args.target_event_cross_zorder, label=args.target_event_cross_label, ) return handle, 1 def project_seismicity_to_segment(eq_df, seg, args): return project_events_to_segment(eq_df, seg, args) def event_marker_sizes(eq, base_size, scale_by_mag, min_size, max_size, mag_ref, mag_scale): if not scale_by_mag: return np.full(len(eq), base_size, dtype=float) mag = pd.to_numeric(eq["mag"], errors="coerce") if "mag" in eq.columns else pd.Series([], dtype=float) if mag.dropna().empty: return np.full(len(eq), base_size, dtype=float) sizes = base_size * np.power(mag_scale, mag.fillna(mag.median()) - mag_ref) return np.clip(sizes, min_size, max_size).to_numpy(dtype=float) def seismicity_marker_sizes(eq, args): return event_marker_sizes(eq, args.seis_marker_size, args.seis_scale_by_mag, args.seis_marker_size_min, args.seis_marker_size_max, args.seis_marker_mag_ref, args.seis_marker_mag_scale) def repeating_marker_sizes(eq, args): return event_marker_sizes(eq, effective_repeating_value(args, "repeating_marker_size", "seis_marker_size"), effective_repeating_value(args, "repeating_scale_by_mag", "seis_scale_by_mag"), effective_repeating_value(args, "repeating_marker_size_min", "seis_marker_size_min"), effective_repeating_value(args, "repeating_marker_size_max", "seis_marker_size_max"), effective_repeating_value(args, "repeating_marker_mag_ref", "seis_marker_mag_ref"), effective_repeating_value(args, "repeating_marker_mag_scale", "seis_marker_mag_scale")) def seismicity_type_key(method): method = "" if pd.isna(method) else str(method).strip() return 1 if method == "HINV" else 2 if method == "DDRT-3.0" else 0 def seismicity_type_label(method): method = "" if pd.isna(method) else str(method).strip() return method if method in {"HINV", "DDRT-3.0"} else "other" def seismicity_marker_for_method(method, args): method = "" if pd.isna(method) else str(method).strip() return args.seis_marker_ddrt30 if method == "DDRT-3.0" else args.seis_marker_hinv if method == "HINV" else args.seis_marker_other def clip_projected_events_to_panel(eq, seg, args): _, _, xpad, zpad = data_span_km(seg) xmin = float(seg["x_km"].min()) - xpad - args.seis_plot_buffer_km xmax = float(seg["x_km"].max()) + xpad + args.seis_plot_buffer_km zmin = float(seg["z_km"].min()) - zpad - args.seis_plot_buffer_km zmax = float(seg["z_km"].max()) + zpad + args.seis_plot_buffer_km return eq[(eq["x_km"] >= xmin) & (eq["x_km"] <= xmax) & (eq["z_km"] >= zmin) & (eq["z_km"] <= zmax)].copy() def draw_method_sorted_events(ax, eq, args, color_mode="seismicity", label_prefix=""): if eq.empty: return None, 0 if "meth" not in eq.columns: eq["meth"] = "other" eq["_plot_type_key"] = eq["meth"].map(seismicity_type_key) eq["_plot_type_label"] = eq["meth"].map(seismicity_type_label) eq["_plot_mag"] = pd.to_numeric(eq["mag"], errors="coerce").fillna(-np.inf) eq = eq.sort_values(["_plot_type_key", "_plot_mag"], ascending=[True, True]).reset_index(drop=True) handles, used_labels = [], set() for _, row in eq.iterrows(): base_label = row["_plot_type_label"] if args.seis_label_by_method else args.seis_label label = f"{label_prefix}{base_label}" if label_prefix else base_label legend_label = label if label not in used_labels else "_nolegend_" used_labels.add(label) single = pd.DataFrame([row]) if color_mode == "repeating": facecolor, edgecolor = args.repeating_facecolor, args.repeating_edgecolor linewidth, alpha, zorder = args.repeating_linewidth, args.repeating_alpha, args.repeating_zorder sizes = repeating_marker_sizes(single, args) else: facecolor, edgecolor = args.seis_facecolor, args.seis_edgecolor linewidth, alpha, zorder = args.seis_linewidth, args.seis_alpha, args.seis_zorder sizes = seismicity_marker_sizes(single, args) handles.append(ax.scatter([row["x_km"]], [row["z_km"]], s=sizes, marker=seismicity_marker_for_method(row["meth"], args), facecolors=facecolor, edgecolors=edgecolor, linewidths=linewidth, alpha=alpha, zorder=zorder, label=legend_label)) return handles[0] if handles else None, len(eq) def overlay_seismicity_on_axis(ax, seg, eq_df, args): if eq_df is None or eq_df.empty: return None, 0 seg_id = int(seg["segment_id"].iloc[0]) if "seis_segment_id" in eq_df.columns: eq_df = eq_df[eq_df["seis_segment_id"].astype(int) == seg_id].copy() if eq_df.empty: return None, 0 eq = clip_projected_events_to_panel(project_seismicity_to_segment(eq_df, seg, args), seg, args) return draw_method_sorted_events(ax, eq, args, color_mode="seismicity", label_prefix="") def overlay_repeating_on_axis(ax, seg, rep_df, args): if rep_df is None or rep_df.empty: return None, 0 seg_id = int(seg["segment_id"].iloc[0]) if "seis_segment_id" in rep_df.columns: rep_df = rep_df[rep_df["seis_segment_id"].astype(int) == seg_id].copy() if rep_df.empty: return None, 0 rep = clip_projected_events_to_panel(project_events_to_segment(rep_df, seg, args), seg, args) return draw_method_sorted_events(ax, rep, args, color_mode="repeating", label_prefix=args.repeating_label_prefix) def overlay_slab_on_axis(ax, seg, slab_df, args): if slab_df is None or slab_df.empty: return None, 0 seg_id = int(seg["segment_id"].iloc[0]) sub = slab_df[slab_df["segment_id"].astype(int) == seg_id].copy() if sub.empty: return None, 0 handles = [] for slab_id, prof in sub.groupby("slab_id", sort=True): prof = prof.sort_values("x_km") label = args.slab_label if args.slab_label else f"Slab2 {slab_id}" if len(sub["slab_id"].dropna().unique()) > 1: label = f"{label} ({slab_id})" h, = ax.plot(prof["x_km"], prof["z_km"], color=args.slab_color, linewidth=args.slab_linewidth, linestyle=args.slab_linestyle, alpha=args.slab_alpha, zorder=args.slab_zorder, label=label) handles.append(h) return handles[0] if handles else None, len(sub) def map_segment_outline_linewidth(seg_id, args): """Return segment-specific map outline width with common fallback.""" if int(seg_id) == 1 and args.map_segment1_outline_linewidth is not None: return float(args.map_segment1_outline_linewidth) if int(seg_id) == 2 and args.map_segment2_outline_linewidth is not None: return float(args.map_segment2_outline_linewidth) return float(args.map_segment_outline_linewidth) def map_segment_outline_linestyle(seg_id, args): """Return segment-specific map outline style with common fallback.""" if int(seg_id) == 1 and args.map_segment1_outline_linestyle is not None: return args.map_segment1_outline_linestyle if int(seg_id) == 2 and args.map_segment2_outline_linestyle is not None: return args.map_segment2_outline_linestyle return args.map_segment_outline_linestyle def map_legend_display_label(raw_label, args): """Rewrite one map legend label from semantic/internal text.""" text = "" if raw_label is None else str(raw_label) target_internal = str(args.map_target_event_label) repeating_prefix = str(args.repeating_label_prefix) if text == "segment 1 outline": return args.map_legend_label_segment1 if text == "segment 2 outline": return args.map_legend_label_segment2 if text == "segment 1 event-selection corridor": return args.map_legend_label_selection_segment1 if text == "segment 2 event-selection corridor": return args.map_legend_label_selection_segment2 if text == "HINV": return args.map_legend_label_hinv if text == "DDRT-3.0": return args.map_legend_label_ddrt30 if text in {"repeating HINV", repeating_prefix + "HINV"}: return args.map_legend_label_repeating_hinv if text in {"repeating DDRT-3.0", repeating_prefix + "DDRT-3.0"}: return args.map_legend_label_repeating_ddrt30 if text == target_internal or text.startswith(target_internal + " ("): return args.map_legend_label_target if text.startswith("Slab2 sampled trace") or text.startswith("Slab2 interface"): return args.map_legend_label_slab return text def map_legend_order(raw_label): """Return stable map-legend order for recognized categories.""" text = "" if raw_label is None else str(raw_label) if text == "segment 1 outline": return 10 if text == "segment 2 outline": return 20 if text == "segment 1 event-selection corridor": return 21 if text == "segment 2 event-selection corridor": return 22 if text == "HINV": return 30 if text == "DDRT-3.0": return 40 if text.endswith("HINV") and text.startswith("repeating"): return 50 if text.endswith("DDRT-3.0") and text.startswith("repeating"): return 60 if text.startswith("target earthquake"): return 70 if text.startswith("Slab2"): return 80 return 100 def map_selection_corridor_polygon(seg, seg_id, args, map_bounds): """Return a map-spanning polygon for the current +/- across filter. The event filter uses perpendicular distance to an infinite strike line through the mean top-row segment center. This polygon depicts that exact criterion over the displayed map extent; it is not a finite fault buffer. """ min_lon, max_lon, min_lat, max_lat = map_bounds lon0, lat0 = segment_trace_origin_lonlat(seg) corner_lon = np.array([min_lon, min_lon, max_lon, max_lon], dtype=float) corner_lat = np.array([min_lat, max_lat, min_lat, max_lat], dtype=float) east, north = lonlat_to_local_en_km(corner_lon, corner_lat, lon0, lat0) az = math.radians(seismicity_strike_for_segment(seg_id, args)) along = east * math.sin(az) + north * math.cos(az) half_width = float(args.seis_max_across_km) margin = max(half_width, 1.0) along_min = float(np.min(along)) - margin along_max = float(np.max(along)) + margin polygon = [] for along_km, across_km in [ (along_min, -half_width), (along_max, -half_width), (along_max, half_width), (along_min, half_width), ]: east_km = along_km * math.sin(az) + across_km * math.cos(az) north_km = along_km * math.cos(az) - across_km * math.sin(az) polygon.append(_offset_lonlat_km(lon0, lat0, east_km, north_km)) polygon.append(polygon[0]) return np.asarray(polygon, dtype=float), lon0, lat0 def draw_map_selection_corridors(ax, combined, args, map_bounds, map_crs=None): """Draw segment-specific polygons showing the actual event filter strips.""" from matplotlib.patches import Polygon as MplPolygon colors = [args.map_selection_segment1_color, args.map_selection_segment2_color] count = 0 for seg_id, seg in combined.groupby("segment_id", sort=True): polygon, lon0, lat0 = map_selection_corridor_polygon( seg, int(seg_id), args, map_bounds ) color = colors[int(seg_id) - 1] if 1 <= int(seg_id) <= 2 else args.map_selection_other_color patch_kwargs = {} if map_crs is None else {"transform": map_crs} patch = MplPolygon( polygon, closed=True, facecolor=color, edgecolor=color, linewidth=args.map_selection_polygon_linewidth, linestyle=args.map_selection_polygon_linestyle, alpha=args.map_selection_polygon_alpha, zorder=args.map_selection_polygon_zorder, label=f"segment {int(seg_id)} event-selection corridor", **patch_kwargs, ) ax.add_patch(patch) if args.map_selection_draw_center: center_kwargs = {} if map_crs is None else {"transform": map_crs} ax.scatter( [lon0], [lat0], marker=args.map_selection_center_marker, s=args.map_selection_center_size, color=color, linewidths=1.2, zorder=args.map_selection_polygon_zorder + 0.1, label="_nolegend_", **center_kwargs, ) print( f"# map event-selection corridor segment {int(seg_id)}: " f"origin={lon0:.6f}/{lat0:.6f}, " f"strike={seismicity_strike_for_segment(seg_id, args):g} deg, " f"half_width={args.seis_max_across_km:g} km" ) count += 1 return count def map_segment_outline_color(seg_id, args): if int(seg_id) == 1: return args.map_segment1_outline_color if int(seg_id) == 2: return args.map_segment2_outline_color return args.map_segment_outline_color def map_segment_outline_points(segment): """Return a closed segment outline from the actual lon/lat grid centers.""" required = {"lon", "lat", "x_index", "z_index"} missing = required.difference(segment.columns) if missing: fail("cannot construct segment outline; missing columns: " + ", ".join(sorted(missing))) seg = segment.copy() xmin = int(seg["x_index"].min()) xmax = int(seg["x_index"].max()) zmin = int(seg["z_index"].min()) zmax = int(seg["z_index"].max()) top = seg[seg["z_index"].astype(int) == zmin].sort_values("x_index") right = seg[seg["x_index"].astype(int) == xmax].sort_values("z_index") bottom = seg[seg["z_index"].astype(int) == zmax].sort_values( "x_index", ascending=False ) left = seg[seg["x_index"].astype(int) == xmin].sort_values( "z_index", ascending=False ) if top.empty or right.empty or bottom.empty or left.empty: fail("cannot construct a complete segment outline from grid edges") pieces = [top, right.iloc[1:], bottom.iloc[1:], left.iloc[1:]] outline = pd.concat(pieces, ignore_index=True) lon = outline["lon"].to_numpy(dtype=float) lat = outline["lat"].to_numpy(dtype=float) lon = np.append(lon, lon[0]) lat = np.append(lat, lat[0]) return lon, lat def map_dip_direction_for_segment(seg_id, args): """Return map dip-direction azimuth; default is right-hand strike + 90.""" if int(seg_id) == 1: value = args.map_dip_direction1_deg elif int(seg_id) == 2: value = args.map_dip_direction2_deg else: value = None if value is None: value = strike_for_segment(seg_id, args) + 90.0 return float(value) % 360.0 def _offset_lonlat_km(lon, lat, east_km, north_km): """Apply a small local east/north offset to lon/lat.""" cos_lat = math.cos(math.radians(float(lat))) if abs(cos_lat) < 1.0e-8: fail("cannot construct map polygon too close to a geographic pole") return ( float(lon) + float(east_km) / (111.32 * cos_lat), float(lat) + float(north_km) / 111.32, ) def map_dip_rectangle_outlines(combined, args): """Build closed segment outlines matching dip_rectangles geometry. The top and bottom boundaries follow the top-row subfault centers. The bottom boundary is shifted down dip by the full reconstructed grid width. End boundaries include half-cell extensions, matching rectangle edges. """ outlines = [] for seg_id, seg in combined.groupby("segment_id", sort=True): seg_id = int(seg_id) strike = math.radians(strike_for_segment(seg_id, args)) dip = math.radians(dip_for_segment(seg_id, args)) dip_direction = math.radians( map_dip_direction_for_segment(seg_id, args) ) zmin = int(seg["z_index"].min()) zmax = int(seg["z_index"].max()) top = seg[seg["z_index"].astype(int) == zmin].copy() top = top.sort_values("x_index").drop_duplicates("x_index") if top.empty: fail(f"segment {seg_id} has no top row for dip-rectangle outline") dx = float(seg["x_grid_space_km"].iloc[0]) dz_horizontal = ( float(seg["z_grid_space_km"].iloc[0]) * abs(math.cos(dip)) ) number_z = zmax - zmin + 1 full_horizontal_width = number_z * dz_horizontal strike_e = math.sin(strike) strike_n = math.cos(strike) dip_e = math.sin(dip_direction) dip_n = math.cos(dip_direction) top_left = top.iloc[0] top_right = top.iloc[-1] # Extend from center locations to the outside rectangle edges. left_e = -0.5 * dx * strike_e left_n = -0.5 * dx * strike_n right_e = 0.5 * dx * strike_e right_n = 0.5 * dx * strike_n top_e = -0.5 * dz_horizontal * dip_e top_n = -0.5 * dz_horizontal * dip_n bottom_e = (full_horizontal_width - 0.5 * dz_horizontal) * dip_e bottom_n = (full_horizontal_width - 0.5 * dz_horizontal) * dip_n top_boundary = [] for row in top.itertuples(index=False): top_boundary.append( _offset_lonlat_km(row.lon, row.lat, top_e, top_n) ) bottom_boundary = [] for row in reversed(list(top.itertuples(index=False))): bottom_boundary.append( _offset_lonlat_km(row.lon, row.lat, bottom_e, bottom_n) ) top_left_corner = _offset_lonlat_km( top_left.lon, top_left.lat, left_e + top_e, left_n + top_n, ) top_right_corner = _offset_lonlat_km( top_right.lon, top_right.lat, right_e + top_e, right_n + top_n, ) bottom_right_corner = _offset_lonlat_km( top_right.lon, top_right.lat, right_e + bottom_e, right_n + bottom_n, ) bottom_left_corner = _offset_lonlat_km( top_left.lon, top_left.lat, left_e + bottom_e, left_n + bottom_n, ) # Use top/bottom center-derived boundaries plus explicit end corners. points = [top_left_corner] points.extend(top_boundary) points.append(top_right_corner) points.append(bottom_right_corner) points.extend(bottom_boundary) points.append(bottom_left_corner) points.append(top_left_corner) lon = np.asarray([point[0] for point in points], dtype=float) lat = np.asarray([point[1] for point in points], dtype=float) outlines.append((seg_id, lon, lat, full_horizontal_width)) return outlines def map_dip_rectangles(combined, args): """Build slip rectangles from top-row geometry, strike, and dip. Rectangle centers for z_index > zmin are reconstructed from the matching top-row x_index center. This avoids using potentially inconsistent deeper lon/lat coordinates when the requested map geometry must follow dip. """ rectangles = [] slips = [] diagnostics = [] for seg_id, seg in combined.groupby("segment_id", sort=True): seg_id = int(seg_id) strike_deg = strike_for_segment(seg_id, args) dip_deg = dip_for_segment(seg_id, args) dip_direction_deg = map_dip_direction_for_segment(seg_id, args) strike = math.radians(strike_deg) dip = math.radians(dip_deg) dip_direction = math.radians(dip_direction_deg) horizontal_dz_factor = abs(math.cos(dip)) zmin = int(seg["z_index"].min()) top = seg[seg["z_index"].astype(int) == zmin].copy() top = top.sort_values("x_index").drop_duplicates("x_index") top_by_x = {int(row.x_index): row for row in top.itertuples(index=False)} if not top_by_x: fail(f"segment {seg_id} has no top-row points for dip rectangles") strike_e = math.sin(strike) strike_n = math.cos(strike) dip_e = math.sin(dip_direction) dip_n = math.cos(dip_direction) first_dx = float(seg["x_grid_space_km"].iloc[0]) first_dz = float(seg["z_grid_space_km"].iloc[0]) diagnostics.append({ "segment_id": seg_id, "strike_deg": strike_deg, "dip_deg": dip_deg, "dip_direction_deg": dip_direction_deg, "cell_length_km": first_dx, "cell_horizontal_width_km": first_dz * horizontal_dz_factor, }) for row in seg.itertuples(index=False): x_index = int(row.x_index) z_index = int(row.z_index) if x_index not in top_by_x: fail( f"segment {seg_id} has no top-row anchor for x_index={x_index}" ) anchor_row = top_by_x[x_index] dx = float(row.x_grid_space_km) dz_horizontal = float(row.z_grid_space_km) * horizontal_dz_factor # The top-row input points are treated as subfault centers. Move the # center down dip by one horizontal cell width per z-index increment. down_dip_center_km = (z_index - zmin) * dz_horizontal center_e = down_dip_center_km * dip_e center_n = down_dip_center_km * dip_n center_lon, center_lat = _offset_lonlat_km( anchor_row.lon, anchor_row.lat, center_e, center_n ) half_strike = 0.5 * dx half_dip = 0.5 * dz_horizontal polygon = [] for strike_sign, dip_sign in [(-1, -1), (1, -1), (1, 1), (-1, 1)]: east = strike_sign * half_strike * strike_e + dip_sign * half_dip * dip_e north = strike_sign * half_strike * strike_n + dip_sign * half_dip * dip_n polygon.append( _offset_lonlat_km(center_lon, center_lat, east, north) ) rectangles.append(polygon) slips.append(float(row.slip_mag)) return rectangles, np.asarray(slips, dtype=float), diagnostics def map_subfault_polygons(combined, args): """Build strike-aligned map polygons using subfault physical dimensions.""" polygons = [] slips = [] metadata = [] for row in combined.itertuples(index=False): seg_id = int(row.segment_id) strike = math.radians(strike_for_segment(seg_id, args)) dip = math.radians(dip_for_segment(seg_id, args)) dip_direction_deg = map_dip_direction_for_segment(seg_id, args) dip_direction = math.radians(dip_direction_deg) half_strike = 0.5 * float(row.x_grid_space_km) # Only the horizontal projection of the down-dip cell is visible in map view. half_dip_horizontal = 0.5 * float(row.z_grid_space_km) * abs(math.cos(dip)) strike_e = math.sin(strike) strike_n = math.cos(strike) dip_e = math.sin(dip_direction) dip_n = math.cos(dip_direction) polygon = [] for strike_sign, dip_sign in [(-1, -1), (1, -1), (1, 1), (-1, 1)]: east = strike_sign * half_strike * strike_e + dip_sign * half_dip_horizontal * dip_e north = strike_sign * half_strike * strike_n + dip_sign * half_dip_horizontal * dip_n polygon.append(_offset_lonlat_km(row.lon, row.lat, east, north)) polygons.append(polygon) slips.append(float(row.slip_mag)) metadata.append((seg_id, dip_direction_deg, 2.0 * half_dip_horizontal)) return polygons, np.asarray(slips, dtype=float), metadata def draw_map_scale_bar(ax, min_lon, max_lon, min_lat, max_lat, args): """Draw a geographic scale bar at a named or custom map position.""" length_km = float(args.map_scale_bar_km) if length_km <= 0.0: return if not (0.0 <= args.map_scale_bar_x <= 1.0): fail("--map_scale_bar_x must be between 0 and 1") if not (0.0 <= args.map_scale_bar_y <= 1.0): fail("--map_scale_bar_y must be between 0 and 1") lon_span = max_lon - min_lon lat_span = max_lat - min_lat location = args.map_scale_bar_loc # Determine the latitude first because longitude length depends on latitude. if location in {"upper_left", "upper_right"}: y = max_lat - args.map_scale_bar_y * lat_span label_sign = -1.0 label_va = "top" else: y = min_lat + args.map_scale_bar_y * lat_span label_sign = 1.0 label_va = "bottom" cos_lat = math.cos(math.radians(y)) if abs(cos_lat) < 1.0e-8: fail("cannot draw map scale bar too close to a geographic pole") dlon = length_km / (111.32 * cos_lat) if location in {"lower_right", "upper_right"}: x1 = max_lon - args.map_scale_bar_x * lon_span x0 = x1 - dlon else: # custom and left-side locations use x as the left-edge fraction. x0 = min_lon + args.map_scale_bar_x * lon_span x1 = x0 + dlon if x0 < min_lon or x1 > max_lon: fail( "map scale bar does not fit at the selected location; reduce " "--map_scale_bar_km or adjust --map_scale_bar_x" ) ax.plot( [x0, x1], [y, y], color=args.map_scale_bar_color, linewidth=args.map_scale_bar_linewidth, solid_capstyle="butt", zorder=20, ) text_y = y + label_sign * 0.018 * lat_span ax.text( 0.5 * (x0 + x1), text_y, f"{length_km:g} km", ha="center", va=label_va, fontsize=args.map_scale_bar_fontsize, color=args.map_scale_bar_color, zorder=20, ) print( f"# map scale bar: length={length_km:g} km " f"location={location} x_fraction={args.map_scale_bar_x:g} " f"y_fraction={args.map_scale_bar_y:g}" ) def _map_scatter_catalog(ax, catalog, args, repeating=False): """Plot catalog rows directly in geographic coordinates.""" if catalog is None or catalog.empty: return 0 data = catalog.copy() if "meth" not in data.columns: data["meth"] = "other" data["_plot_type_key"] = data["meth"].map(seismicity_type_key) data["_plot_mag"] = pd.to_numeric( data["mag"], errors="coerce" ).fillna(-np.inf) if "mag" in data.columns else -np.inf data = data.sort_values(["_plot_type_key", "_plot_mag"]) used_labels = set() for method, group in data.groupby("meth", dropna=False, sort=False): method_label = seismicity_type_label(method) if repeating: label = f"{args.repeating_label_prefix}{method_label}" facecolor = args.repeating_facecolor edgecolor = args.repeating_edgecolor linewidth = args.repeating_linewidth alpha = args.repeating_alpha zorder = args.repeating_zorder sizes = repeating_marker_sizes(group, args) else: label = method_label if args.seis_label_by_method else args.seis_label facecolor = args.seis_facecolor edgecolor = args.seis_edgecolor linewidth = args.seis_linewidth alpha = args.seis_alpha zorder = args.seis_zorder sizes = seismicity_marker_sizes(group, args) if label in used_labels: label = "_nolegend_" else: used_labels.add(label) ax.scatter( group["lon"], group["lat"], s=np.asarray(sizes, dtype=float) * args.map_event_marker_scale, marker=seismicity_marker_for_method(method, args), facecolors=facecolor, edgecolors=edgecolor, linewidths=linewidth, alpha=alpha, zorder=zorder, label=label, ) return len(data) def _map_bounds(combined, seismicity_map, repeating_map, slab_profile, args): """Return explicit or automatically derived lon/lat map bounds.""" lon_values = [pd.to_numeric(combined["lon"], errors="coerce").to_numpy()] lat_values = [pd.to_numeric(combined["lat"], errors="coerce").to_numpy()] optional = [] if args.map_plot_seismicity and seismicity_map is not None: optional.append(seismicity_map) if args.map_plot_repeating and repeating_map is not None: optional.append(repeating_map) if args.map_plot_slab_trace and slab_profile is not None: optional.append(slab_profile) for data in optional: if not data.empty and "lon" in data.columns and "lat" in data.columns: lon_values.append(pd.to_numeric(data["lon"], errors="coerce").to_numpy()) lat_values.append(pd.to_numeric(data["lat"], errors="coerce").to_numpy()) lon = np.concatenate(lon_values) lat = np.concatenate(lat_values) lon = lon[np.isfinite(lon)] lat = lat[np.isfinite(lat)] if lon.size == 0 or lat.size == 0: fail("no finite longitude/latitude values available for map bounds") if args.map_padding_deg < 0.0: fail("--map_padding_deg must be non-negative") min_lon = float(lon.min()) - args.map_padding_deg if args.map_min_lon is None else args.map_min_lon max_lon = float(lon.max()) + args.map_padding_deg if args.map_max_lon is None else args.map_max_lon min_lat = float(lat.min()) - args.map_padding_deg if args.map_min_lat is None else args.map_min_lat max_lat = float(lat.max()) + args.map_padding_deg if args.map_max_lat is None else args.map_max_lat if not (min_lon < max_lon and min_lat < max_lat): fail("invalid map bounds: minimum values must be smaller than maximum values") return min_lon, max_lon, min_lat, max_lat def plot_latlon_slip_map(outputs, args, metrics=None, seismicity_map=None, repeating_map=None, slab_profile=None): """Plot slip and optional catalogs in a longitude/latitude framework.""" combined = outputs["combined_slip_dist_latlon"].copy() map_render_mode = ( args.map_fault_render if args.map_plot_slip else "outline" ) if not args.map_plot_slip: print("# map slip layer disabled; plotting fault geometry and overlays only") slip = pd.to_numeric(combined["slip_mag"], errors="coerce") finite_slip = slip.replace([np.inf, -np.inf], np.nan).dropna() plotted_slip_max = float(finite_slip.max()) if not finite_slip.empty else float("nan") colorbar_extend = "max" if plotted_slip_max > args.slip_cmax else "neither" min_lon, max_lon, min_lat, max_lat = _map_bounds( combined, seismicity_map, repeating_map, slab_profile, args ) # Include target epicenter in automatic map bounds. if args.map_plot_target_event: if not ( np.isfinite(args.target_event_lat) and np.isfinite(args.target_event_lon) and np.isfinite(args.target_event_model_depth_km) ): fail("target-event latitude, longitude, and model depth must be finite") if args.target_event_model_depth_km < 0.0: fail("--target_event_model_depth_km must be non-negative") if args.map_min_lon is None: min_lon = min(min_lon, args.target_event_lon - args.map_padding_deg) if args.map_max_lon is None: max_lon = max(max_lon, args.target_event_lon + args.map_padding_deg) if args.map_min_lat is None: min_lat = min(min_lat, args.target_event_lat - args.map_padding_deg) if args.map_max_lat is None: max_lat = max(max_lat, args.target_event_lat + args.map_padding_deg) if args.map_fig_width <= 0.0 or args.map_fig_height <= 0.0: fail("--map_fig_width and --map_fig_height must be positive") use_cartopy = args.map_coastline_source == "cartopy" map_crs = None if use_cartopy: try: import cartopy.crs as ccrs except ImportError: fail( "--map_coastline_source cartopy requires Cartopy. " "Use file or none if Cartopy is unavailable." ) map_crs = ccrs.PlateCarree() fig, ax = plt.subplots( figsize=(args.map_fig_width, args.map_fig_height), subplot_kw={"projection": map_crs}, ) else: fig, ax = plt.subplots(figsize=(args.map_fig_width, args.map_fig_height)) norm = Normalize(vmin=args.slip_cmin, vmax=args.slip_cmax) # Render fault geometry using the selected map mode. if map_render_mode == "midpoint_cells": midpoint_polygons, midpoint_slip, midpoint_outlines, midpoint_info = \ map_midpoint_cells(combined) map_mappable = PolyCollection( midpoint_polygons, array=midpoint_slip, cmap=args.map_cmap, norm=norm, edgecolors=args.map_midpoint_cell_edgecolor, linewidths=args.map_midpoint_cell_linewidth, alpha=args.map_midpoint_cell_alpha, zorder=args.map_midpoint_cell_zorder, label="_nolegend_", transform=map_crs, ) ax.add_collection(map_mappable) if args.map_midpoint_draw_outlines: for outline_seg_id, outline_lon, outline_lat in midpoint_outlines: ax.plot( outline_lon, outline_lat, color=map_segment_outline_color(outline_seg_id, args), linewidth=map_segment_outline_linewidth(outline_seg_id, args), linestyle=map_segment_outline_linestyle(outline_seg_id, args), zorder=args.map_segment_outline_zorder, label=f"segment {int(outline_seg_id)} outline", transform=map_crs, ) for info_seg_id, info_nx, info_nz in midpoint_info: print( f"# map segment {info_seg_id}: midpoint cells = " f"{info_nx} x {info_nz}; outer boundaries pinned to grid endpoints" ) elif map_render_mode == "dip_rectangles": rectangles, rectangle_slip, rectangle_diagnostics = map_dip_rectangles( combined, args ) map_mappable = PolyCollection( rectangles, array=rectangle_slip, cmap=args.map_cmap, norm=norm, edgecolors=args.map_dip_rectangle_edgecolor, linewidths=args.map_dip_rectangle_linewidth, alpha=args.map_dip_rectangle_alpha, zorder=args.map_dip_rectangle_zorder, label="_nolegend_", transform=map_crs, ) ax.add_collection(map_mappable) for item in rectangle_diagnostics: print( f"# map segment {item['segment_id']}: " f"strike={item['strike_deg']:g} deg " f"dip={item['dip_deg']:g} deg " f"dip_direction={item['dip_direction_deg']:g} deg " f"rectangle={item['cell_length_km']:g} x " f"{item['cell_horizontal_width_km']:g} km" ) # Draw outlines matching the reconstructed dip rectangles. for outline_seg_id, outline_lon, outline_lat, outline_width_km in \ map_dip_rectangle_outlines(combined, args): ax.plot( outline_lon, outline_lat, color=map_segment_outline_color(outline_seg_id, args), linewidth=map_segment_outline_linewidth(outline_seg_id, args), linestyle=map_segment_outline_linestyle(outline_seg_id, args), zorder=args.map_segment_outline_zorder, label=f"segment {int(outline_seg_id)} outline", transform=map_crs, ) print( f"# map segment {int(outline_seg_id)} dip-rectangle " f"outline width = {outline_width_km:g} km" ) elif map_render_mode == "polygons": polygons, polygon_slip, polygon_metadata = map_subfault_polygons(combined, args) map_mappable = PolyCollection( polygons, array=polygon_slip, cmap=args.map_cmap, norm=norm, edgecolors=args.map_fault_edgecolor, linewidths=args.map_fault_linewidth, alpha=args.map_polygon_alpha, zorder=3, label="slip subfault polygons", transform=map_crs, ) ax.add_collection(map_mappable) for seg_id in sorted(set(item[0] for item in polygon_metadata)): sub = [item for item in polygon_metadata if item[0] == seg_id] print( f"# map segment {seg_id}: dip direction = {sub[0][1]:g} deg; " f"horizontal cell width = {sub[0][2]:g} km" ) elif map_render_mode in {"centers", "outline_centers"}: marker_size = ( args.map_outline_center_size if map_render_mode == "outline_centers" else args.map_marker_size ) marker_alpha = ( args.map_outline_center_alpha if map_render_mode == "outline_centers" else 1.0 ) marker_symbol = ( args.map_outline_center_marker if map_render_mode == "outline_centers" else "s" ) map_mappable = ax.scatter( combined["lon"], combined["lat"], c=slip, cmap=args.map_cmap, norm=norm, marker=marker_symbol, s=marker_size, edgecolors=args.map_fault_edgecolor, linewidths=args.map_fault_linewidth, alpha=marker_alpha, zorder=4, label="slip subfault centers", transform=map_crs, ) else: # Outline-only mode still needs a ScalarMappable for the common slip # colorbar, even though no slip-colored artist is drawn on the map. map_mappable = matplotlib.cm.ScalarMappable(norm=norm, cmap=args.map_cmap) map_mappable.set_array(slip.to_numpy(dtype=float)) if map_render_mode in {"outline", "outline_centers"}: for seg_id, seg in combined.groupby("segment_id", sort=True): outline_lon, outline_lat = map_segment_outline_points(seg) ax.plot( outline_lon, outline_lat, color=map_segment_outline_color(seg_id, args), linewidth=map_segment_outline_linewidth(seg_id, args), linestyle=map_segment_outline_linestyle(seg_id, args), zorder=args.map_segment_outline_zorder, label=f"segment {int(seg_id)} outline", transform=map_crs, ) if args.map_draw_top_trace: for seg_id, seg in combined.groupby("segment_id", sort=True): top = seg[seg["z_index"].astype(int) == int(seg["z_index"].min())].copy() top = top.sort_values("x_km") ax.plot( top["lon"], top["lat"], color=args.map_top_trace_color, linewidth=args.map_top_trace_linewidth, linestyle=args.map_top_trace_linestyle, zorder=args.map_top_trace_zorder, label=f"segment {int(seg_id)} top trace", ) if args.map_plot_selection_polygons: draw_map_selection_corridors( ax, combined, args, (min_lon, max_lon, min_lat, max_lat), map_crs=map_crs, ) n_seis = 0 n_rep = 0 if args.map_plot_seismicity: n_seis = _map_scatter_catalog(ax, seismicity_map, args, repeating=False) if args.map_plot_repeating: n_rep = _map_scatter_catalog(ax, repeating_map, args, repeating=True) # Target-earthquake epicenter. Model depth is metadata in map view. if args.map_plot_target_event: target_label = args.map_target_event_label if args.map_target_event_label_depth: target_label = ( f"{target_label} ({args.target_event_model_depth_km:g} km depth)" ) ax.scatter( [args.target_event_lon], [args.target_event_lat], s=args.map_target_event_size, marker=args.map_target_event_marker, facecolors=args.map_target_event_facecolor, edgecolors=args.map_target_event_edgecolor, linewidths=args.map_target_event_linewidth, alpha=args.map_target_event_alpha, zorder=args.map_target_event_zorder, label=target_label, transform=map_crs, ) print( f"# target earthquake epicenter: " f"lat={args.target_event_lat:g} lon={args.target_event_lon:g} " f"model_depth={args.target_event_model_depth_km:g} km" ) n_slab = 0 if args.map_plot_slab_trace and slab_profile is not None and not slab_profile.empty: for slab_id, profile in slab_profile.groupby("slab_id", sort=True): profile = profile.sort_values(["segment_id", "x_km_trace"]) ax.plot( profile["lon"], profile["lat"], color=args.slab_color, linewidth=args.slab_linewidth, linestyle=args.slab_linestyle, alpha=args.slab_alpha, zorder=args.slab_zorder, label=f"Slab2 sampled trace ({slab_id})", ) n_slab += len(profile) if args.map_coastline_source == "cartopy": try: ax.coastlines( resolution=args.map_coastline_resolution, linewidth=args.map_coastline_linewidth, color=args.map_coastline_color, zorder=args.map_coastline_zorder, ) except Exception as exc: fail( "Cartopy coastline rendering failed. Natural Earth data may " "not be locally available. Use --map_coastline_source file " "or none. Original error: %s" % exc ) print( f"# map coastline source = cartopy; " f"resolution = {args.map_coastline_resolution}" ) elif args.map_coastline_source == "file": if not args.map_coastline_file: fail("--map_coastline_source file requires --map_coastline_file") coastline_segments = read_map_line_segments(args.map_coastline_file) for segment in coastline_segments: coast_lon, coast_lat = zip(*segment) ax.plot( coast_lon, coast_lat, color=args.map_coastline_color, linewidth=args.map_coastline_linewidth, zorder=args.map_coastline_zorder, ) print( f"# map coastline source = file; segments = " f"{len(coastline_segments)}; file = {args.map_coastline_file}" ) elif args.map_coastline_file: print( "# WARNING: --map_coastline_file ignored because " "--map_coastline_source is none" ) if use_cartopy: ax.set_extent( [min_lon, max_lon, min_lat, max_lat], crs=map_crs, ) else: ax.set_xlim(min_lon, max_lon) ax.set_ylim(min_lat, max_lat) mean_lat = 0.5 * (min_lat + max_lat) cos_lat = math.cos(math.radians(mean_lat)) if abs(cos_lat) < 1.0e-8: fail("map mean latitude is too close to a pole for aspect correction") ax.set_aspect(1.0 / cos_lat, adjustable="box") draw_map_scale_bar(ax, min_lon, max_lon, min_lat, max_lat, args) # Configure geographic tick labels and guide lines. lon_ticks = map_tick_values(min_lon, max_lon, args.map_lon_tick_interval) lat_ticks = map_tick_values(min_lat, max_lat, args.map_lat_tick_interval) if use_cartopy: import matplotlib.ticker as mticker from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER gridliner = ax.gridlines( crs=map_crs, draw_labels=True, linewidth=(args.map_gridline_linewidth if args.map_gridlines else 0.0), color=args.map_gridline_color, alpha=(args.map_gridline_alpha if args.map_gridlines else 0.0), linestyle=args.map_gridline_linestyle, x_inline=False, y_inline=False, ) gridliner.top_labels = False gridliner.right_labels = False gridliner.xformatter = LONGITUDE_FORMATTER gridliner.yformatter = LATITUDE_FORMATTER gridliner.xlabel_style = {"size": args.map_tick_label_size} gridliner.ylabel_style = {"size": args.map_tick_label_size} if lon_ticks is not None and len(lon_ticks): gridliner.xlocator = mticker.FixedLocator(lon_ticks) if lat_ticks is not None and len(lat_ticks): gridliner.ylocator = mticker.FixedLocator(lat_ticks) else: if lon_ticks is not None and len(lon_ticks): ax.set_xticks(lon_ticks) if lat_ticks is not None and len(lat_ticks): ax.set_yticks(lat_ticks) ax.tick_params( axis="both", which="both", direction="out", bottom=True, left=True, top=True, right=True, labelbottom=True, labelleft=True, labeltop=False, labelright=False, labelsize=args.map_tick_label_size, ) ax.grid( args.map_gridlines, color=args.map_gridline_color, linewidth=args.map_gridline_linewidth, alpha=args.map_gridline_alpha, linestyle=args.map_gridline_linestyle, ) ax.set_xlabel("Longitude (deg)") ax.set_ylabel("Latitude (deg)") if use_cartopy: ax.tick_params(axis="both", which="both", direction="out") title_lines = [] if args.map_title: title_lines.append(args.map_title) elif args.title: title_lines.append(args.title) if metrics is not None and args.map_plot_slip: title_lines.append(summary_text(metrics)) if title_lines: ax.set_title("\n".join(title_lines), fontsize=9, pad=8, linespacing=1.05) if args.map_plot_slip: cbar = fig.colorbar( map_mappable, ax=ax, orientation="horizontal", pad=0.12, fraction=0.05, extend=colorbar_extend, ) cbar.set_label("Slip (cm)") cbar.set_ticks(range(int(args.slip_cmin), int(args.slip_cmax) + 1, 100)) if args.map_legend: handles, raw_labels = ax.get_legend_handles_labels() if handles: entries = [] seen = set() for appearance_index, (handle, raw_label) in enumerate( zip(handles, raw_labels) ): display_label = map_legend_display_label(raw_label, args) if ( not display_label or display_label == "_nolegend_" or display_label in seen ): continue seen.add(display_label) proxy = legend_scatter_proxy( handle, raw_label, args.map_legend_marker_size, args, ) entries.append( (map_legend_order(raw_label), appearance_index, proxy, display_label) ) entries.sort(key=lambda item: (item[0], item[1])) map_handles = [item[2] for item in entries] map_labels = [item[3] for item in entries] ax.legend( map_handles, map_labels, loc=args.map_legend_loc, fontsize=args.map_legend_fontsize, markerscale=args.map_legend_markerscale, ) fig.savefig( args.map_plot_file, format=args.map_plot_format, bbox_inches="tight", pad_inches=0.08, ) plt.close(fig) print(f"# wrote {args.map_plot_file}") print(f"# map bounds = {min_lon:g}/{max_lon:g}/{min_lat:g}/{max_lat:g}") print(f"# map overlays: seismicity={n_seis} repeating={n_rep} slab_samples={n_slab}") if args.map_plot_slip and colorbar_extend == "max": print(f"# map colorbar saturated above {args.slip_cmax:g} cm; plotted maximum = {plotted_slip_max:g} cm") def legend_symbol_scale(label, args): """Return the visual-area multiplier for a semantic legend category.""" text = "" if label is None else str(label) target_cross = str(args.target_event_cross_label) target_map = str(args.map_target_event_label) if text == "HINV": return args.legend_scale_hinv if text == "DDRT-3.0": return args.legend_scale_ddrt30 if text == "repeating HINV": return args.legend_scale_repeating_hinv if text == "repeating DDRT-3.0": return args.legend_scale_repeating_ddrt30 if ( text in {target_cross, "target hypocenter", target_map} or (target_map and text.startswith(target_map + " (")) ): return args.legend_scale_target return 1.0 def legend_scatter_proxy(handle, raw_label, base_size, args): """Copy a scatter artist and assign a category-specific legend area.""" from copy import copy from matplotlib.collections import PathCollection if not isinstance(handle, PathCollection): return handle proxy = copy(handle) proxy.set_sizes([ float(base_size) * float(legend_symbol_scale(raw_label, args)) ]) return proxy def scatter_legend_handler_map(marker_size): """Use one absolute size for every PathCollection legend symbol.""" from matplotlib.collections import PathCollection from matplotlib.legend_handler import HandlerPathCollection return { PathCollection: HandlerPathCollection( numpoints=1, sizes=[float(marker_size)], ) } def cross_section_legend_label_and_priority(label, args): """Return configurable display text and the requested semantic order.""" text = "" if label is None else str(label) target_internal = str(args.target_event_cross_label) slab_internal = str(args.slab_label) if text == "HINV": return args.legend_label_hinv, 10 if text == "DDRT-3.0": return args.legend_label_ddrt30, 20 if text == "repeating HINV": return args.legend_label_repeating_hinv, 30 if text == "repeating DDRT-3.0": return args.legend_label_repeating_ddrt30, 40 if text in {target_internal, "target hypocenter"}: return args.legend_label_target, 50 # Slab labels may include an ID suffix, for example "Slab2 interface (cas)". if text in {slab_internal, "Slab2 interface"}: return args.legend_label_slab, 60 if slab_internal and text.startswith(slab_internal + " ("): return args.legend_label_slab + text[len(slab_internal):], 60 if text.startswith("Slab2 interface ("): return args.legend_label_slab + text[len("Slab2 interface"):], 60 # Preserve any future or user-added legend item after the known categories. return text, 100 def cross_section_display_label(label, args): """Backward-compatible display-label helper.""" return cross_section_legend_label_and_priority(label, args)[0] def _ordered_unique_cross_section_entries(handles, labels, args): """Rewrite, sort, and de-duplicate one set of legend entries.""" entries = [] seen = set() for appearance_index, (handle, label) in enumerate(zip(handles, labels)): display, priority = cross_section_legend_label_and_priority(label, args) if not display or display == "_nolegend_" or display in seen: continue seen.add(display) entries.append((priority, appearance_index, handle, display)) entries.sort(key=lambda item: (item[0], item[1])) return [item[2] for item in entries], [item[3] for item in entries] def update_cross_section_legends(axes, args): """Rewrite, order, scale, and optionally combine cross-section legends.""" axes_list = list(np.asarray(axes, dtype=object).ravel()) axes_list = [axis for axis in axes_list if axis is not None] if not axes_list: return def ordered_entries(handles, labels): entries = [] seen = set() for appearance_index, (handle, raw_label) in enumerate(zip(handles, labels)): display, priority = cross_section_legend_label_and_priority( raw_label, args ) if not display or display == "_nolegend_" or display in seen: continue seen.add(display) proxy = legend_scatter_proxy( handle, raw_label, args.seis_legend_marker_size, args, ) entries.append((priority, appearance_index, proxy, display)) entries.sort(key=lambda item: (item[0], item[1])) return [item[2] for item in entries], [item[3] for item in entries] if args.cross_section_legend_mode == "per_segment": for axis in axes_list: handles, labels = axis.get_legend_handles_labels() old_legend = axis.get_legend() if old_legend is not None: old_legend.remove() ordered_handles, ordered_labels = ordered_entries(handles, labels) if ordered_handles: axis.legend( ordered_handles, ordered_labels, loc=args.seis_legend_loc, fontsize=args.seis_legend_fontsize, markerscale=args.seis_legend_markerscale, ) print("# cross-section legend: mode=per_segment; ordered=True") return all_entries = [] seen = set() appearance_index = 0 for axis in axes_list: handles, labels = axis.get_legend_handles_labels() old_legend = axis.get_legend() if old_legend is not None: old_legend.remove() for handle, raw_label in zip(handles, labels): display, priority = cross_section_legend_label_and_priority( raw_label, args ) if not display or display == "_nolegend_" or display in seen: appearance_index += 1 continue seen.add(display) proxy = legend_scatter_proxy( handle, raw_label, args.seis_legend_marker_size, args, ) all_entries.append( (priority, appearance_index, proxy, display) ) appearance_index += 1 all_entries.sort(key=lambda item: (item[0], item[1])) combined_handles = [item[2] for item in all_entries] combined_labels = [item[3] for item in all_entries] if combined_handles: axes_list[0].legend( combined_handles, combined_labels, loc=args.seis_legend_loc, fontsize=args.seis_legend_fontsize, markerscale=args.seis_legend_markerscale, ) print( f"# cross-section legend: mode=segment1_combined " f"ordered=True entries={len(combined_labels)} labels={combined_labels}" ) def overlay_slip_contours(ax, seg, args): """Draw configurable slip-magnitude contours on one fault panel.""" if not args.plot_slip_contours: return [] x_index = pd.to_numeric(seg["x_index"], errors="coerce") z_index = pd.to_numeric(seg["z_index"], errors="coerce") slip = pd.to_numeric(seg["slip_mag"], errors="coerce") valid = x_index.notna() & z_index.notna() & slip.notna() if not valid.any(): return [] work = seg.loc[valid, ["x_index", "z_index", "x_km", "z_km", "slip_mag"]].copy() work["x_index"] = pd.to_numeric(work["x_index"], errors="raise").astype(int) work["z_index"] = pd.to_numeric(work["z_index"], errors="raise").astype(int) x_values = np.sort(work["x_km"].unique()) z_values = np.sort(work["z_km"].unique()) grid = np.full((len(z_values), len(x_values)), np.nan, dtype=float) x_lookup = {float(value): i for i, value in enumerate(x_values)} z_lookup = {float(value): i for i, value in enumerate(z_values)} for row in work.itertuples(index=False): grid[z_lookup[float(row.z_km)], x_lookup[float(row.x_km)]] = float(row.slip_mag) finite = grid[np.isfinite(grid)] if finite.size == 0: return [] data_min, data_max = float(finite.min()), float(finite.max()) interval = float(args.slip_contour_interval_cm) first_level = math.ceil(max(data_min, 0.0) / interval) * interval if first_level <= 0.0: first_level = interval levels = np.arange(first_level, data_max + 0.5 * interval, interval) levels = levels[(levels > data_min) & (levels < data_max)] if levels.size == 0: return [] contour = ax.contour( x_values, z_values, np.ma.masked_invalid(grid), levels=levels, colors=args.slip_contour_color, linewidths=args.slip_contour_linewidth, linestyles=args.slip_contour_linestyle, alpha=args.slip_contour_alpha, zorder=3, ) if args.slip_contour_labels: ax.clabel( contour, contour.levels, inline=True, inline_spacing=2, fmt=lambda value: f"{value:g}", fontsize=args.slip_contour_label_size, colors=args.slip_contour_color, ) levels_list = [float(value) for value in contour.levels] print( f"# plotted slip contours on segment {int(seg['segment_id'].iloc[0])}: " + ", ".join(f"{value:g}" for value in levels_list) + " cm" ) return levels_list def overlay_frw_patch_boundaries_on_axis(ax, seg, args): """Draw selected FRW patch cell boundaries on an along-strike/down-dip panel.""" summaries = getattr(args, "frw_patch_summaries", None) or [] if not summaries: return 0 seg_id = int(seg["segment_id"].iloc[0]) segment_patches = [p for p in summaries if int(p["segment_id"]) == seg_id] if not segment_patches: return 0 if "z_origin_depth_km" in seg.columns: z_origin = float(seg["z_origin_depth_km"].iloc[0]) else: # Combined rake outputs do not retain z_origin_depth_km. Recover it from # z_km = z_origin + z_index * dz using the segment grid spacing. dz_segment = float(seg["z_grid_space_km"].iloc[0]) z_origin_values = ( pd.to_numeric(seg["z_km"], errors="coerce") - pd.to_numeric(seg["z_index"], errors="coerce") * dz_segment ) finite_origin = z_origin_values[np.isfinite(z_origin_values)] if finite_origin.empty: fail(f"cannot recover z origin for FRW patch overlay on segment {seg_id}") z_origin = float(finite_origin.median()) count = 0 for patch in segment_patches: dx = float(patch["dx_km"]) dz = float(patch["dz_km"]) # The FRW mask is selected at cell centers. Plot half-cell-expanded edges so # the dashed rectangle encloses exactly the retained subfault cells. x0 = float(patch["realized_xmin_km"]) - 0.5 * dx x1 = float(patch["realized_xmax_km"]) + 0.5 * dx z0 = z_origin + float(patch["realized_zmin_km"]) - 0.5 * dz z1 = z_origin + float(patch["realized_zmax_km"]) + 0.5 * dz # Keep the dashed patch boundary and P-number annotation on the panel, # but intentionally exclude patch boundaries from the legend. label = "_nolegend_" ax.plot( [x0, x1, x1, x0, x0], [z0, z0, z1, z1, z0], color="0.35", linestyle="--", linewidth=1.25, zorder=6, label=label, ) ax.text( x0 + 0.25 * dx, z0 + 0.75 * dz, f"P{int(patch['patch_id'])}", color="0.25", fontsize=8, ha="left", va="top", zorder=7, bbox={"facecolor": "white", "edgecolor": "none", "alpha": 0.7, "pad": 1.0}, ) print( f"# plotted FRW patch {int(patch['patch_id'])} boundary on segment {seg_id}: " f"x_edges={x0:.3f}..{x1:.3f} km " f"z_edges_local={float(patch['realized_zmin_km']) - 0.5 * dz:.3f}.." f"{float(patch['realized_zmax_km']) + 0.5 * dz:.3f} km" ) count += 1 return count def plot_two_segment_slip_dist(outputs, args, metrics=None, seismicity=None, repeating=None, slab_profile=None, target_projection=None): combined = outputs["combined_slip_dist_latlon"] segment_ids = sorted(combined["segment_id"].dropna().astype(int).unique()) if len(segment_ids) == 0: fail("no segment_id values found for plotting") width_ratios = [] for seg_id in segment_ids: seg = combined[combined["segment_id"].astype(int) == seg_id] xspan, zspan, _, _ = data_span_km(seg) width_ratios.append(xspan / zspan) panel_height = 5.6 if args.panel_width_scale <= 0.0: fail("--panel_width_scale must be positive") fig_width_auto = panel_height * sum(width_ratios) + 0.35 fig_width = max(6.0, fig_width_auto * args.panel_width_scale) fig, axes = plt.subplots( 1, len(segment_ids), figsize=(fig_width, panel_height), squeeze=False, sharey=True, gridspec_kw={ "wspace": args.panel_wspace, "width_ratios": width_ratios, } ) axes = axes[0] # Equal-aspect axes can leave unused horizontal space inside their subplot # allocations. Anchor plotting boxes toward their shared inner boundary. if args.panels_anchor_inward and len(axes) > 1: for ipanel, ax in enumerate(axes): if ipanel == 0: ax.set_anchor("E") elif ipanel == len(axes) - 1: ax.set_anchor("W") else: ax.set_anchor("C") if metrics is not None: plot_top = 0.76 if args.title else 0.78 else: plot_top = 0.90 fig.subplots_adjust( top=plot_top, wspace=args.panel_wspace, ) norm = Normalize(vmin=args.slip_cmin, vmax=args.slip_cmax) scatter_handle = None for ipanel, (ax, seg_id) in enumerate(zip(axes, segment_ids)): seg = combined[combined["segment_id"].astype(int) == seg_id].copy() scatter_handle = ax.scatter(seg["x_km"], seg["z_km"], c=seg["slip_mag"], cmap="hot_r", norm=norm, marker="s", s=args.marker_size, edgecolors="none", linewidths=0.0, zorder=2) overlay_slip_contours(ax, seg, args) for row in seg.itertuples(index=False): slip = float(row.slip_mag) angle = float(row.rake_angle_deg) if not pd.isna(row.rake_angle_deg) else float("nan") if slip <= 0.0 or math.isnan(angle): continue length = slip * args.vec_scale if length <= 0.0: continue dx = length * math.cos(math.radians(angle)) dz = length * math.sin(math.radians(angle)) shaft_width = max(args.arrow_shaft_width_min, min(length * args.arrow_shaft_width_scale, args.arrow_shaft_width_max)) head_width = max(args.arrow_head_width_min, min(length * args.arrow_head_width_scale, args.arrow_head_width_max)) head_length = max(args.arrow_head_length_min, min(length * args.arrow_head_length_scale, args.arrow_head_length_max)) linewidth = max(args.arrow_linewidth_min, min(length * args.arrow_linewidth_scale, args.arrow_linewidth_max)) head_length = min(head_length, length * args.arrow_head_length_fraction_max) ax.arrow(row.x_km, row.z_km, dx, dz, width=shaft_width, head_width=head_width, head_length=head_length, length_includes_head=True, color="black", linewidth=linewidth, zorder=4) if args.plot_seismicity and seismicity is not None: _, n_seis = overlay_seismicity_on_axis(ax, seg, seismicity, args) if n_seis > 0: print(f"# plotted seismicity on segment {seg_id}: {n_seis} events") if args.plot_repeating and repeating is not None: _, n_rep = overlay_repeating_on_axis(ax, seg, repeating, args) if n_rep > 0: print(f"# plotted repeating earthquakes on segment {seg_id}: {n_rep} events") if args.target_event_cross_section and target_projection is not None: _, n_target = overlay_target_event_on_axis( ax, seg, target_projection, args ) if n_target: print(f"# plotted target hypocenter on segment {seg_id}") if args.plot_slab and slab_profile is not None: _, n_slab = overlay_slab_on_axis(ax, seg, slab_profile, args) if n_slab > 0: print(f"# plotted slab profile on segment {seg_id}: {n_slab} samples") n_frw_patches = overlay_frw_patch_boundaries_on_axis(ax, seg, args) if n_frw_patches: print(f"# plotted {n_frw_patches} FRW patch boundary/boundaries on segment {seg_id}") _, _, xpad, zpad = data_span_km(seg) xmin = float(seg["x_km"].min()) - xpad xmax = float(seg["x_km"].max()) + xpad ax.set_xlim(xmax, xmin) if x_axis_flip_for_segment(seg_id, args) else ax.set_xlim(xmin, xmax) ax.set_ylim(float(seg["z_km"].min()) - zpad, float(seg["z_km"].max()) + zpad) ax.invert_yaxis() ax.set_aspect("equal", adjustable="box") strike_deg = strike_for_segment(seg_id, args) dip_deg = dip_for_segment(seg_id, args) dx_grid = float(seg["x_grid_space_km"].median()) dz_grid = float(seg["z_grid_space_km"].median()) ax.set_title( f"Segment {seg_id} - grid dx/dz = {dx_grid:g}/{dz_grid:g} km\n" f"strike = {strike_deg:g} deg; dip = {dip_deg:g} deg", fontsize=9, pad=5, ) ax.set_xlabel("Along-strike distance (km)") if ipanel == 0: ax.set_ylabel("Along-dip distance (km)") else: ax.set_ylabel("") ax.tick_params(labelleft=False) # Major labels remain uncluttered; minor ticks and faint guides show # every local FRW subfault-center spacing for easy patch selection. ax.xaxis.set_minor_locator(MultipleLocator(dx_grid)) ax.yaxis.set_minor_locator(MultipleLocator(dz_grid)) ax.set_axisbelow(True) ax.grid(False, which="major") ax.grid(which="minor", color="0.78", linestyle=":", linewidth=0.45, alpha=0.75) ax.tick_params(axis="both", which="major", direction="out", top=True, right=True) ax.tick_params(axis="both", which="minor", direction="out", top=True, right=True, length=2.5) show_target_legend = ( args.target_event_cross_section and target_projection is not None and not target_projection.empty ) if ( ((args.plot_seismicity or args.plot_repeating) and args.seis_legend) or (args.plot_slab and args.slab_legend) or show_target_legend ): ax.legend(loc=args.seis_legend_loc, fontsize=args.seis_legend_fontsize) # Apply legend text overrides after every segment has been plotted. # Default behavior puts all used entries on segment 1 only. update_cross_section_legends(axes, args) for legend_axis in axes: legend = legend_axis.get_legend() if legend is not None: legend.set_zorder(1000) legend.get_frame().set_alpha(1.0) legend.get_frame().set_facecolor("white") if scatter_handle is not None: # Retain slip_cmax but mark saturation with a maximum-side arrow. plotted_slip = pd.to_numeric( combined["slip_mag"], errors="coerce" ).replace([np.inf, -np.inf], np.nan) plotted_slip_max = ( float(plotted_slip.max()) if plotted_slip.notna().any() else float("nan") ) colorbar_extend = ( "max" if plotted_slip_max > args.slip_cmax else "neither" ) if colorbar_extend == "max": print( f"# colorbar saturated above {args.slip_cmax:g} cm; " f"plotted maximum = {plotted_slip_max:g} cm" ) cbar = fig.colorbar( scatter_handle, ax=list(axes), orientation="horizontal", pad=0.12, fraction=0.05, extend=colorbar_extend, ) cbar.set_label("Slip (cm)") cbar.set_ticks( range(int(args.slip_cmin), int(args.slip_cmax) + 1, 100) ) # Use one unified figure title instead of a separate summary text box. title_lines = [] if args.title: title_lines.append(args.title) if metrics is not None: title_lines.append(summary_text(metrics)) if title_lines: fig.suptitle( "\n".join(title_lines), fontsize=9, y=0.985, va="top", linespacing=1.05, ) fig.savefig(args.plot_file, format=args.plot_format, bbox_inches="tight", pad_inches=0.08) plt.close(fig) print(f"# wrote {args.plot_file}") def plot_mu_vs_depth(args, blocks=None, metrics=None): if args.mu_depth_from_data and blocks: depth_vals = pd.concat([b["depth"] for b in blocks.values()], ignore_index=True) positive_depth = depth_vals[depth_vals > 0.0] if positive_depth.empty: fail("no positive depth values available for mu-depth plot") min_depth = max(float(positive_depth.min()), args.mu_depth_min_km if args.mu_depth_min_km > 0 else 0.001) max_depth = float(positive_depth.max()) else: min_depth = args.mu_depth_min_km if args.mu_depth_min_km > 0 else 0.001 max_depth = args.mu_depth_max_km depth = np.linspace(min_depth, max_depth, args.mu_depth_samples) mu_gpa = mu_pa_from_depth(depth) / 1.0e9 fig, ax = plt.subplots(figsize=(6.2, 7.0)) ax.plot(mu_gpa, depth, linewidth=2) if args.mu_ref_gpa is not None: ax.axvline(args.mu_ref_gpa, color="black", linestyle="--", linewidth=1.4, label=f"{args.mu_ref_gpa:g} GPa reference") ax.legend(loc="lower right") ax.set_ylim(args.mu_depth_axis_min_km, max_depth) ax.invert_yaxis() ax.set_xlabel("mu from depth equation (GPa)") ax.set_ylabel("Depth (km)") ax.set_title(args.mu_plot_title) ax.grid(True, alpha=0.3) if metrics is not None: add_summary_box(fig, metrics) fig.tight_layout() fig.savefig(args.mu_plot_file, format=args.mu_plot_format, bbox_inches="tight") plt.close(fig) out_csv = Path(args.mu_plot_file).with_suffix(".csv") pd.DataFrame({"depth_km": depth, "mu_gpa": mu_gpa}).to_csv(out_csv, index=False) print(f"# wrote {args.mu_plot_file}") print(f"# wrote {out_csv}") def plot_mu_compare(outputs, args, metrics=None): combined = outputs["combined_slip_dist_latlon"] rows = [] for rake in [1, 2]: depth_col = f"mu_gpa_from_depth_rake{rake}" input_slip_mu_col = f"mu_gpa_rake{rake}" if depth_col not in combined.columns or input_slip_mu_col not in combined.columns: continue sub = combined[["segment_id", "depth", depth_col, input_slip_mu_col]].copy() sub = sub.rename(columns={depth_col: "mu_gpa_from_depth", input_slip_mu_col: "mu_gpa_from_input_moment_slip"}) sub["rake"] = rake rows.append(sub) if not rows: fail("no mu comparison columns found") comp = pd.concat(rows, ignore_index=True).replace([np.inf, -np.inf], np.nan) comp = comp.dropna(subset=["depth", "mu_gpa_from_depth", "mu_gpa_from_input_moment_slip"]) if comp.empty: fail("no finite mu comparison values found") comp["mu_gpa_difference_input_minus_depth"] = comp["mu_gpa_from_input_moment_slip"] - comp["mu_gpa_from_depth"] comp["mu_gpa_ratio_input_over_depth"] = comp["mu_gpa_from_input_moment_slip"] / comp["mu_gpa_from_depth"] fig, ax = plt.subplots(figsize=(6.4, 7.0)) depth_curve = comp.groupby("depth", as_index=False)["mu_gpa_from_depth"].mean().sort_values("depth") ax.plot(depth_curve["mu_gpa_from_depth"], depth_curve["depth"], color="black", linewidth=2.0, label="depth equation", zorder=3) for (seg_id, rake), sub in comp.groupby(["segment_id", "rake"]): ax.scatter(sub["mu_gpa_from_input_moment_slip"], sub["depth"], s=args.mu_compare_marker_size, alpha=args.mu_compare_alpha, label=f"seg{int(seg_id)} rake{int(rake)} input moment/slip", zorder=2) if args.mu_ref_gpa is not None: ax.axvline(args.mu_ref_gpa, color="gray", linestyle="--", linewidth=1.2, label=f"{args.mu_ref_gpa:g} GPa reference", zorder=1) ax.set_xlabel("mu (GPa)") ax.set_ylabel("Depth (km)") ax.set_title(args.mu_compare_title) ax.grid(True, alpha=0.3) ax.invert_yaxis() ax.legend(fontsize=8) if metrics is not None: add_summary_box(fig, metrics) fig.tight_layout() fig.savefig(args.mu_compare_file, format=args.mu_compare_format, bbox_inches="tight") plt.close(fig) csv_name = Path(args.mu_compare_file).with_suffix(".csv") comp.to_csv(csv_name, index=False) print(f"# wrote {args.mu_compare_file}") print(f"# wrote {csv_name}") def write_blocks(blocks, output_prefix, write_csv=True, write_txt=True): out_cols = ["lon", "lat", "depth", "slip", "slip_input", "slip_cm_from_moment", "slip_source", "mu_gpa_from_depth", "mu_gpa_from_input_slip", "mu_gpa_from_active_slip", "moment", "segment_id", "rake_id", "local_index", "x_index", "z_index", "x_km", "z_km", "x_grid_space_km", "z_grid_space_km", "subfault_area_km2", "subfault_area_cm2"] for name, block in blocks.items(): out = block[out_cols] if write_csv: csv_name = f"{output_prefix}.block.{name}.csv" out.to_csv(csv_name, index=False) print(f"# wrote {csv_name}") if write_txt: txt_name = f"{output_prefix}.block.{name}.txt" out.to_csv(txt_name, sep=" ", index=False, header=True, float_format="%.6f", na_rep="NaN") print(f"# wrote {txt_name}") def write_outputs(outputs, output_prefix, write_csv=True, write_txt=True): for name, df in outputs.items(): if write_csv: csv_name = f"{output_prefix}.{name}.csv" df.to_csv(csv_name, index=False) print(f"# wrote {csv_name}") if write_txt: txt_name = f"{output_prefix}.{name}.txt" df.to_csv(txt_name, sep=" ", index=False, header=True, float_format="%.6f", na_rep="NaN") print(f"# wrote {txt_name}") def write_summary_metrics(metrics, output_prefix): csv_name = f"{output_prefix}.summary.csv" txt_name = f"{output_prefix}.summary.txt" pd.DataFrame([metrics]).to_csv(csv_name, index=False) with open(txt_name, "w", encoding="utf-8") as fout: for key, value in metrics.items(): fout.write(f"{key} = {value}\n") print(f"# wrote {csv_name}") print(f"# wrote {txt_name}") def print_summary(df, blocks, outputs, args, metrics): print(f"# combined_file = {args.combined_file}") print(f"# rupture_config = {args.rupture_config}") print(f"# total_rows = {len(df)}") print(f"# use_moment_for_slip = {args.use_moment_for_slip}") print(f"# equal_xy_scale = enabled (1 km along strike = 1 km along dip)") print(f"# param_M0 = {args.param_M0:g} # M0(dyne-cm) = moment * param_M0") print(f"# strike1_deg = {args.strike1_deg:g}") print(f"# strike2_deg = {args.strike2_deg:g}") print("# subfault area = (x_length_km/xgrid_num) * (z_length_km/zgrid_num) * (100000^2) cm^2") for name, block in blocks.items(): print(f"# block {name}: dx={block.x_grid_space_km.iloc[0]:.6g} km dz={block.z_grid_space_km.iloc[0]:.6g} km area={block.subfault_area_km2.iloc[0]:.6g} km^2 area_cm2={block.subfault_area_cm2.iloc[0]:.6e}") print("# moment/slip summary") for k in ["total_moment_input_units", "total_M0_dyne_cm", "total_M0_N_m", "Mw", "max_slip_cm", "min_slip_thresh", "min_slip_for_stats_cm", "mean_slip_cm", "median_slip_cm"]: print(f"# {k:28s} = {metrics[k]}") print(f"# n_slip_used/total = {metrics['n_slip_used']}/{metrics['n_slip_total']}") for name, out in outputs.items(): print(f"# output {name}: rows={len(out)} slip_mag_max={out.slip_mag.max():.6f} slip_mag_nonzero={(out.slip_mag > 0).sum()}") def main(): parser = argparse.ArgumentParser(description="Read combined reg_inv lat/lon output, compute slip magnitude/rake angle, and optionally plot.") parser.add_argument("--combined_file", required=True) parser.add_argument("--frw_file", default=None, help="Optional FRW file containing individual moment grids ordered by segment, rake, and time window.") parser.add_argument("--num_time_windows", type=int, default=None, help="Time windows per segment/rake. Defaults to multimeN in --rupture_config, otherwise 1.") parser.add_argument("--frw_rtol", type=float, default=5.0e-4, help="Relative tolerance for cumulative FRW versus combined moment validation.") parser.add_argument("--frw_atol", type=float, default=5.0e-2, help="Absolute tolerance for cumulative FRW versus combined moment validation.") for patch_id in range(1, 4): parser.add_argument( f"--frw_patch{patch_id}", nargs=5, default=None, metavar=("SEGMENT", "XMIN_KM", "XMAX_KM", "ZMIN_KM", "ZMAX_KM"), help=("Retain original FRW values only in this inclusive local x/z rectangle; " "all other cells and segments are written as zero."), ) parser.add_argument( f"--frw_patch{patch_id}_output", default=f"fault_patch{patch_id}.frw", help=f"Output FRW filename for patch {patch_id}." ) parser.add_argument("--frw_patch_summary", default="fault_patch_summary.csv") parser.add_argument( "--frw_patch_beachball_file", default="frw_patch_focal_mechanisms.pdf", help="Separate PDF showing one focal mechanism per configured FRW patch." ) parser.add_argument("--rupture_config", required=True) parser.add_argument("--xgrid_num", required=True, type=int) parser.add_argument("--zgrid_num", required=True, type=int) parser.add_argument("--xgrid_num2", type=int, default=None) parser.add_argument("--zgrid_num2", type=int, default=None) parser.add_argument("--x_length_km", required=True, type=float) parser.add_argument("--z_length_km", required=True, type=float) parser.add_argument("--x_length_km2", type=float, default=None) parser.add_argument("--z_length_km2", type=float, default=None) parser.add_argument("--output_prefix", default="ff_slip_dist_latlon") parser.add_argument("--theta1_deg", type=float, default=-135.0) parser.add_argument("--theta2_deg", type=float, default=-225.0) parser.add_argument("--theta_sign", type=float, default=-1.0) parser.add_argument("--strike1_deg", type=float, default=275.0, help="Strike angle for segment 1, used only for plot title. Default: 275 deg.") parser.add_argument("--strike2_deg", type=float, default=280.0, help="Strike angle for segment 2, used only for plot title. Default: 280 deg.") parser.add_argument("--geometry_tol", type=float, default=1.0e-5) parser.add_argument("--use_moment_for_slip", action="store_true") parser.add_argument("--param_M0", type=float, default=PARAM_M0_DEFAULT, help="Scale factor converting input moment column to dyne-cm: M0 = moment * param_M0. Default: 1e20.") parser.add_argument("--min_slip_thresh", "--min_slip_threh", dest="min_slip_thresh", type=float, default=0.1, help="Fraction of max slip used as minimum slip cutoff for mean/median slip statistics. Default: 0.1.") parser.add_argument("--make_plot", action="store_true") parser.add_argument("--plot_file", default="ff_slip_dist_latlon.pdf") parser.add_argument("--plot_format", default="pdf") parser.add_argument("--title", default="") parser.add_argument("--slip_cmin", type=float, default=0.0) parser.add_argument("--slip_cmax", type=float, default=600.0) add_bool_flag( parser, "plot_slip_contours", default=True, help_text="Draw slip-magnitude contours on fault cross sections." ) parser.add_argument("--slip_contour_interval_cm", type=float, default=100.0) parser.add_argument("--slip_contour_color", default="0.20") parser.add_argument("--slip_contour_linewidth", type=float, default=0.75) parser.add_argument("--slip_contour_linestyle", default="solid") parser.add_argument("--slip_contour_alpha", type=float, default=0.90) add_bool_flag( parser, "slip_contour_labels", default=True, help_text="Label slip contours with values in centimeters." ) parser.add_argument("--slip_contour_label_size", type=float, default=6.5) parser.add_argument("--vec_scale", type=float, default=0.006) parser.add_argument( "--panel_width_scale", type=float, default=1.0, help=( "Scale factor applied to the automatically calculated figure width. " "Values below 1 place equal-aspect panels closer together. " "Try 0.85 to 0.90 for two square segments. Default: 1.0." ), ) parser.add_argument("--marker_size", type=float, default=170.0) parser.add_argument( "--panel_wspace", type=float, default=0.0, help=( "Horizontal spacing between segment subplot allocations, as a fraction " "of average subplot width. Smaller values place panels closer. Default: 0.0." ), ) add_bool_flag( parser, "panels_anchor_inward", default=True, help_text=( "Anchor equal-aspect segment panels toward their shared inner boundary " "to reduce the visible gap." ), ) parser.add_argument("--arrow_shaft_width_scale", type=float, default=0.020) parser.add_argument("--arrow_head_width_scale", type=float, default=0.075) parser.add_argument("--arrow_head_length_scale", type=float, default=0.120) parser.add_argument("--arrow_linewidth_scale", type=float, default=0.080) parser.add_argument("--arrow_shaft_width_min", type=float, default=0.004) parser.add_argument("--arrow_shaft_width_max", type=float, default=0.045) parser.add_argument("--arrow_head_width_min", type=float, default=0.040) parser.add_argument("--arrow_head_width_max", type=float, default=0.250) parser.add_argument("--arrow_head_length_min", type=float, default=0.060) parser.add_argument("--arrow_head_length_max", type=float, default=0.380) parser.add_argument("--arrow_linewidth_min", type=float, default=0.10) parser.add_argument("--arrow_linewidth_max", type=float, default=0.50) parser.add_argument("--arrow_head_length_fraction_max", type=float, default=0.45) add_bool_flag(parser, "x_axis_flip", default=True) add_bool_flag(parser, "x_axis_flip2", default=True) # Latitude/longitude map-view plot (Phase 1) parser.add_argument("--make_map_plot", action="store_true", help="Create an additional longitude/latitude slip map.") parser.add_argument("--map_plot_file", default="ff_slip_dist_latlon_map.pdf") parser.add_argument("--map_plot_format", default="pdf") parser.add_argument("--map_title", default="") add_bool_flag( parser, "map_plot_target_event", default=True, help_text="Plot the target-earthquake epicenter as a star on the map.", ) parser.add_argument("--target_event_lat", type=float, default=40.374, help="Target-earthquake latitude in degrees. Default: 40.374.") parser.add_argument("--target_event_lon", type=float, default=-125.022, help="Target-earthquake longitude in degrees. Default: -125.022.") parser.add_argument( "--target_event_model_depth_km", "--target_event_depth_km", dest="target_event_model_depth_km", type=float, default=9.0, help=( "Target-earthquake model depth in km. " "--target_event_depth_km is retained as a backward-compatible alias. " "Default: 9.0." ), ) parser.add_argument( "--target_event_depth_mode", choices=["model_z", "vertical_depth"], default="model_z", help=( "Interpret --target_event_model_depth_km as fault-plane model z " "or as vertical hypocentral depth. Default: model_z." ), ) parser.add_argument("--map_target_event_marker", default="*") parser.add_argument("--map_target_event_size", type=float, default=180.0) parser.add_argument("--map_target_event_facecolor", default="gold") parser.add_argument("--map_target_event_edgecolor", default="black") parser.add_argument("--map_target_event_linewidth", type=float, default=1.0) parser.add_argument("--map_target_event_alpha", type=float, default=1.0) parser.add_argument("--map_target_event_zorder", type=float, default=25.0) parser.add_argument("--map_target_event_label", default="target earthquake") add_bool_flag( parser, "target_event_cross_section", default=True, help_text="Project and plot the target hypocenter on fault cross-sections.", ) parser.add_argument( "--target_event_segment_mode", choices=["nearest_surface", "nearest_plane", "all", "segment1", "segment2"], default="nearest_surface", help=( "Cross-section target assignment. nearest_surface chooses the " "smallest constant-model-depth geographic residual. " "nearest_plane is retained as an alias. Default: nearest_surface." ), ) parser.add_argument("--target_event_cross_marker", default="*") parser.add_argument("--target_event_cross_size", type=float, default=180.0) parser.add_argument("--target_event_cross_facecolor", default="gold") parser.add_argument("--target_event_cross_edgecolor", default="black") parser.add_argument("--target_event_cross_linewidth", type=float, default=1.0) parser.add_argument("--target_event_cross_alpha", type=float, default=1.0) parser.add_argument("--target_event_cross_zorder", type=float, default=25.0) parser.add_argument("--target_event_cross_label", default="target hypocenter") add_bool_flag( parser, "target_event_cross_clip", default=True, help_text="Clip the projected target to the displayed segment panel bounds.", ) add_bool_flag( parser, "map_target_event_label_depth", default=True, help_text="Include target geoid depth in the map legend label.", ) parser.add_argument("--map_min_lon", type=float, default=None) parser.add_argument("--map_max_lon", type=float, default=None) parser.add_argument("--map_min_lat", type=float, default=None) parser.add_argument("--map_max_lat", type=float, default=None) parser.add_argument("--map_padding_deg", type=float, default=0.05) parser.add_argument("--map_fig_width", type=float, default=8.0) parser.add_argument("--map_fig_height", type=float, default=7.0) parser.add_argument("--map_marker_size", type=float, default=80.0) add_bool_flag( parser, "map_plot_slip", default=True, help_text=( "Plot slip-colored subfault cells and the map slip colorbar. " "Use --no_map_plot_slip for a geometry/EQ-only map." ), ) parser.add_argument( "--map_fault_render", choices=["midpoint_cells", "dip_rectangles", "polygons", "centers", "outline", "outline_centers"], default="midpoint_cells", help="Render physical subfault polygons or square center markers. Default: polygons.", ) parser.add_argument("--map_dip_direction1_deg", type=float, default=None, help="Segment-1 dip-direction azimuth. Default: strike1+90 deg.") parser.add_argument("--map_dip_direction2_deg", type=float, default=None, help="Segment-2 dip-direction azimuth. Default: strike2+90 deg.") parser.add_argument("--map_polygon_alpha", type=float, default=1.0) parser.add_argument("--map_dip_rectangle_edgecolor", default="none") parser.add_argument("--map_dip_rectangle_linewidth", type=float, default=0.0) parser.add_argument("--map_dip_rectangle_alpha", type=float, default=1.0) parser.add_argument("--map_dip_rectangle_zorder", type=float, default=4.0) parser.add_argument("--map_midpoint_cell_edgecolor", default="none") parser.add_argument("--map_midpoint_cell_linewidth", type=float, default=0.0) parser.add_argument("--map_midpoint_cell_alpha", type=float, default=1.0) parser.add_argument("--map_midpoint_cell_zorder", type=float, default=4.0) add_bool_flag( parser, "map_midpoint_draw_outlines", default=True, help_text="Draw segment outlines around midpoint-derived slip cells.", ) parser.add_argument("--map_segment1_outline_color", default="blue") parser.add_argument("--map_segment2_outline_color", default="red") parser.add_argument("--map_segment_outline_color", default="black", help="Fallback outline color for segment IDs above 2.") parser.add_argument("--map_segment_outline_linewidth", type=float, default=2.0) parser.add_argument("--map_segment_outline_linestyle", default="solid") parser.add_argument( "--map_segment1_outline_linewidth", type=float, default=None, help="Segment-1 outline width; default inherits --map_segment_outline_linewidth.", ) parser.add_argument( "--map_segment2_outline_linewidth", type=float, default=None, help="Segment-2 outline width; default inherits --map_segment_outline_linewidth.", ) parser.add_argument( "--map_segment1_outline_linestyle", default=None, help="Segment-1 outline style; default inherits --map_segment_outline_linestyle.", ) parser.add_argument( "--map_segment2_outline_linestyle", default=None, help="Segment-2 outline style; default inherits --map_segment_outline_linestyle.", ) parser.add_argument("--map_segment_outline_zorder", type=float, default=14.0) parser.add_argument("--map_outline_center_marker", default="s") parser.add_argument("--map_outline_center_size", type=float, default=18.0, help="Slip-center marker size for outline_centers mode.") parser.add_argument("--map_outline_center_alpha", type=float, default=1.0) parser.add_argument( "--map_coastline_source", choices=["none", "cartopy", "file"], default="none", help="Coastline source: none, Cartopy Natural Earth, or local line file.", ) parser.add_argument( "--map_coastline_resolution", choices=["110m", "50m", "10m"], default="10m", help="Cartopy Natural Earth coastline resolution. Default: 10m.", ) parser.add_argument("--map_coastline_file", default="", help="Optional GMT multisegment lon/lat line file; '>' separates segments.") parser.add_argument("--map_coastline_color", default="0.25") parser.add_argument("--map_coastline_linewidth", type=float, default=0.8) parser.add_argument("--map_coastline_zorder", type=float, default=1.0) parser.add_argument( "--map_scale_bar_loc", choices=[ "custom", "lower_left", "lower_right", "upper_left", "upper_right", ], default="lower_left", help=( "Named scale-bar position. custom uses --map_scale_bar_x/y " "directly. For named positions, x/y specify the horizontal and " "vertical margin fractions. Default: lower_left." ), ) parser.add_argument("--map_scale_bar_km", type=float, default=20.0, help="Scale-bar length in km; set 0 to disable.") parser.add_argument("--map_scale_bar_x", type=float, default=0.08, help="Scale-bar left position as fraction of map width.") parser.add_argument("--map_scale_bar_y", type=float, default=0.06, help="Scale-bar vertical position as fraction of map height.") parser.add_argument("--map_scale_bar_linewidth", type=float, default=3.0) parser.add_argument("--map_scale_bar_color", default="black") parser.add_argument("--map_scale_bar_fontsize", type=float, default=9.0) parser.add_argument("--map_cmap", default="hot_r") parser.add_argument("--map_fault_edgecolor", default="none") parser.add_argument("--map_fault_linewidth", type=float, default=0.0) add_bool_flag(parser, "map_draw_top_trace", default=True, help_text="Draw the shallowest subfault row for each segment.") parser.add_argument("--map_top_trace_color", default="black") parser.add_argument("--map_top_trace_linewidth", type=float, default=1.2) parser.add_argument("--map_top_trace_linestyle", default="--", help="Segment top-trace line style. Default: --.") parser.add_argument("--map_top_trace_zorder", type=float, default=15.0, help="Top-trace drawing order. Default: 15.") add_bool_flag( parser, "map_plot_selection_polygons", default=False, help_text=( "Draw the exact +/- across-fault event-selection strips used for " "cross-section catalog filtering on the map." ), ) parser.add_argument("--map_selection_segment1_color", default="tab:blue") parser.add_argument("--map_selection_segment2_color", default="tab:orange") parser.add_argument("--map_selection_other_color", default="0.5") parser.add_argument("--map_selection_polygon_alpha", type=float, default=0.12) parser.add_argument("--map_selection_polygon_linewidth", type=float, default=1.2) parser.add_argument("--map_selection_polygon_linestyle", default="dashed") parser.add_argument("--map_selection_polygon_zorder", type=float, default=1.5) add_bool_flag(parser, "map_selection_draw_center", default=True) parser.add_argument("--map_selection_center_marker", default="+") parser.add_argument("--map_selection_center_size", type=float, default=70.0) parser.add_argument( "--map_legend_label_selection_segment1", default="segment 1 selection corridor", ) parser.add_argument( "--map_legend_label_selection_segment2", default="segment 2 selection corridor", ) parser.add_argument("--map_plot_seismicity", action="store_true", help="Plot filtered catalog events at their original lon/lat.") parser.add_argument("--map_plot_repeating", action="store_true", help="Plot repeating earthquakes at their original lon/lat.") parser.add_argument("--map_plot_slab_trace", action="store_true", help="Plot the sampled Slab2 trace at its sampled lon/lat.") parser.add_argument("--map_event_marker_scale", type=float, default=0.65) add_bool_flag(parser, "map_legend", default=True) parser.add_argument("--map_legend_loc", default="best") # Longitude/latitude map legend text. parser.add_argument("--map_legend_label_segment1", default="segment 1 outline") parser.add_argument("--map_legend_label_segment2", default="segment 2 outline") parser.add_argument("--map_legend_label_hinv", default="HINV") parser.add_argument("--map_legend_label_ddrt30", default="DDRT-3.0") parser.add_argument( "--map_legend_label_repeating_hinv", default="repeating HINV" ) parser.add_argument( "--map_legend_label_repeating_ddrt30", default="repeating DDRT-3.0" ) parser.add_argument("--map_legend_label_target", default="target earthquake") parser.add_argument("--map_legend_label_slab", default="Slab2 interface") parser.add_argument("--map_legend_fontsize", type=float, default=8.0) parser.add_argument( "--map_legend_marker_size", type=float, default=48.0, help="Equal scatter-symbol area in the map legend, in pt^2.", ) parser.add_argument( "--map_legend_markerscale", type=float, default=0.60, help=( "Marker scale used only in the longitude/latitude map legend. " "Does not change plotted map symbols. Default: 0.60." ), ) parser.add_argument("--map_lon_tick_interval", type=float, default=None, help="Longitude tick interval in degrees; default is automatic.") parser.add_argument("--map_lat_tick_interval", type=float, default=None, help="Latitude tick interval in degrees; default is automatic.") parser.add_argument("--map_tick_label_size", type=float, default=9.0) add_bool_flag(parser, "map_gridlines", default=True, help_text="Draw light longitude/latitude guide lines.") parser.add_argument("--map_gridline_color", default="0.75") parser.add_argument("--map_gridline_linewidth", type=float, default=0.5) parser.add_argument("--map_gridline_alpha", type=float, default=0.6) parser.add_argument("--map_gridline_linestyle", default=":") parser.add_argument("--plot_mu_depth", action="store_true") parser.add_argument("--mu_plot_file", default="mu_gpa_vs_depth.pdf") parser.add_argument("--mu_plot_format", default="pdf") parser.add_argument("--mu_plot_title", default="mu vs. depth") parser.add_argument("--mu_ref_gpa", type=float, default=30.0) parser.add_argument("--mu_depth_min_km", type=float, default=0.001) parser.add_argument("--mu_depth_axis_min_km", type=float, default=0.0) parser.add_argument("--mu_depth_max_km", type=float, default=40.0) parser.add_argument("--mu_depth_samples", type=int, default=1000) parser.add_argument("--mu_depth_from_data", action="store_true") parser.add_argument("--plot_mu_compare", action="store_true") parser.add_argument("--mu_compare_file", default="mu_input_moment_slip_vs_depth.pdf") parser.add_argument("--mu_compare_format", default="pdf") parser.add_argument("--mu_compare_title", default="mu comparison as a function of depth") parser.add_argument("--mu_compare_marker_size", type=float, default=18.0) parser.add_argument("--mu_compare_alpha", type=float, default=0.65) parser.add_argument("--plot_seismicity", action="store_true", help="Overlay selected earthquakes on the slip plot.") parser.add_argument("--seis_catalog", default="", help="NCAeqDDRT-style earthquake catalog file.") parser.add_argument("--seis_output", default="", help="Optional CSV output of selected seismicity after all filters, including across-distance filter.") parser.add_argument("--seis_candidate_output", default="", help="CSV preserving every seismicity event x segment diagnostic candidate.") parser.add_argument("--seis_min_lat", type=float, default=None) parser.add_argument("--seis_max_lat", type=float, default=None) parser.add_argument("--seis_min_lon", type=float, default=None) parser.add_argument("--seis_max_lon", type=float, default=None) parser.add_argument("--seis_min_depth", type=float, default=None) parser.add_argument("--seis_max_depth", type=float, default=None) parser.add_argument("--seis_min_mag", type=float, default=None) parser.add_argument("--seis_max_mag", type=float, default=None) parser.add_argument("--seis_method", default=None) parser.add_argument("--seis_start_time", default=None) parser.add_argument("--seis_end_time", default=None) parser.add_argument("--seis_polygon", default=None) parser.add_argument("--seis_polygon_radius", type=float, default=0.0) parser.add_argument("--seis_max_across_km", type=float, default=5.0) add_bool_flag( parser, "seis_across_filter", default=True, help_text=( "Apply the +/- --seis_max_across_km corridor before segment " "assignment. With --no_seis_across_filter, every catalog event is " "an assignment candidate for every segment before panel clipping." ), ) parser.add_argument( "--seis_segment_assignment_mode", choices=["nearest", "all_candidates"], default="nearest", help=( "nearest keeps one closest segment per event; all_candidates keeps " "one copy on every segment whose corridor test passes. Default: nearest." ), ) parser.add_argument("--seis_depth_offset_km", type=float, default=0.293) parser.add_argument("--seis_strike1_deg", type=float, default=275.0) parser.add_argument("--seis_strike2_deg", type=float, default=280.0) parser.add_argument("--dip1_deg", "--dip_deg", dest="dip1_deg", type=float, default=84.0) parser.add_argument("--dip2_deg", type=float, default=84.0) add_bool_flag(parser, "ignore_across_strike_distance_depth_projection", default=True, help_text="Collapse EQ/slab depth using depth only. Use --no_ignore_across_strike_distance_depth_projection for full projection.") parser.add_argument("--seis_marker_ddrt30", default="o") parser.add_argument("--seis_marker_hinv", default="D") parser.add_argument("--seis_marker_other", default="x") parser.add_argument("--seis_label", default="earthquakes") parser.add_argument("--seis_label_by_method", action="store_true", default=True) parser.add_argument("--seis_facecolor", default="none") parser.add_argument("--seis_edgecolor", default="cyan") parser.add_argument("--seis_linewidth", type=float, default=0.8) parser.add_argument("--seis_alpha", type=float, default=0.85) parser.add_argument("--seis_zorder", type=float, default=6.0) parser.add_argument("--seis_marker_size", type=float, default=26.0) parser.add_argument("--seis_marker_size_min", type=float, default=8.0) parser.add_argument("--seis_marker_size_max", type=float, default=180.0) parser.add_argument("--seis_scale_by_mag", action="store_true") parser.add_argument("--seis_marker_mag_ref", type=float, default=3.0) parser.add_argument("--seis_marker_mag_scale", type=float, default=1.8) parser.add_argument("--seis_plot_buffer_km", type=float, default=2.0) add_bool_flag(parser, "seis_legend", default=True) parser.add_argument("--seis_legend_loc", default="lower left") parser.add_argument("--seis_legend_fontsize", type=float, default=8.0) # Per-symbol visual compensation in legends. Each value multiplies # the applicable absolute legend marker area. parser.add_argument("--legend_scale_hinv", type=float, default=1.0) parser.add_argument("--legend_scale_ddrt30", type=float, default=1.0) parser.add_argument( "--legend_scale_repeating_hinv", type=float, default=1.0 ) parser.add_argument( "--legend_scale_repeating_ddrt30", type=float, default=1.0 ) parser.add_argument( "--legend_scale_target", type=float, default=2.0, help="Target-star legend area multiplier. Default: 2.0.", ) parser.add_argument( "--seis_legend_marker_size", type=float, default=48.0, help="Equal scatter-symbol area in cross-section legends, in pt^2.", ) parser.add_argument( "--seis_legend_markerscale", type=float, default=0.60, help=( "Marker scale used only in cross-section legends. " "Does not change plotted event or target symbols. Default: 0.60." ), ) parser.add_argument("--plot_repeating", action="store_true") parser.add_argument("--repeating_catalog", default="") parser.add_argument("--repeating_output", default="") parser.add_argument("--repeating_candidate_output", default="", help="CSV preserving every repeating-event x segment diagnostic candidate.") parser.add_argument("--repeating_use_seis_filters", action="store_true", default=True) parser.add_argument("--repeating_depth_offset_km", type=float, default=None) parser.add_argument("--repeating_label_prefix", default="repeating ") parser.add_argument("--repeating_facecolor", default="red") parser.add_argument("--repeating_edgecolor", default="black") parser.add_argument("--repeating_linewidth", type=float, default=0.9) parser.add_argument("--repeating_alpha", type=float, default=0.95) parser.add_argument("--repeating_zorder", type=float, default=8.0) parser.add_argument("--repeating_marker_size", type=float, default=None) parser.add_argument("--repeating_marker_size_min", type=float, default=None) parser.add_argument("--repeating_marker_size_max", type=float, default=None) group_rep_scale = parser.add_mutually_exclusive_group() group_rep_scale.add_argument("--repeating_scale_by_mag", dest="repeating_scale_by_mag", action="store_true") group_rep_scale.add_argument("--no_repeating_scale_by_mag", dest="repeating_scale_by_mag", action="store_false") parser.set_defaults(repeating_scale_by_mag=None) parser.add_argument("--repeating_marker_mag_ref", type=float, default=None) parser.add_argument("--repeating_marker_mag_scale", type=float, default=None) parser.add_argument("--plot_slab", action="store_true") parser.add_argument( "--require_slab_profile", action="store_true", help=( "Treat missing Slab2 geographic coverage as fatal. By default, backend " "execution errors remain fatal, but segment/model combinations with all-NaN " "samples are skipped and plotting continues without unavailable Slab2 curves." ), ) parser.add_argument("--slabdir", default="") parser.add_argument("--slab_id", default="") parser.add_argument("--slab_sample_backend", choices=["auto", "xarray", "gmt", "pygmt"], default="auto") parser.add_argument("--gmt_cmd", default="gmt") parser.add_argument("--slab_sample_dx_km", type=float, default=1.0) parser.add_argument("--slab_depth_offset_km", type=float, default=None) parser.add_argument("--slab_output", default="") parser.add_argument("--slab_color", default="blue") parser.add_argument("--slab_linewidth", type=float, default=2.0) parser.add_argument("--slab_linestyle", default="-") parser.add_argument("--slab_alpha", type=float, default=0.95) parser.add_argument("--slab_zorder", type=float, default=7.0) parser.add_argument("--slab_label", default="Slab2 interface") add_bool_flag(parser, "slab_legend", default=True) parser.add_argument("--no_blocks", action="store_true") parser.add_argument("--no_csv", action="store_true") parser.add_argument("--no_txt", action="store_true") parser.add_argument("--allow_extra_rows", action="store_true") # Cross-section legend text and placement. parser.add_argument( "--cross_section_legend_mode", choices=["segment1_combined", "per_segment"], default="segment1_combined", help=( "segment1_combined collects all used entries from all segment " "panels and places one legend on segment 1; per_segment keeps " "one legend on each panel." ), ) parser.add_argument("--legend_label_hinv", default="HINV") parser.add_argument("--legend_label_ddrt30", default="DDRT-3.0") parser.add_argument( "--legend_label_repeating_hinv", default="repeating HINV" ) parser.add_argument( "--legend_label_repeating_ddrt30", default="repeating DDRT-3.0" ) parser.add_argument("--legend_label_slab", default="Slab2 interface") parser.add_argument("--legend_label_target", default="target hypocenter") args = parser.parse_args() if (args.plot_slab or (args.make_map_plot and args.map_plot_slab_trace)): if args.slab_sample_backend == "auto": fail( "Slab plotting requires an explicit --slab_sample_backend " "(xarray, gmt, or pygmt); auto fallback is disabled." ) if args.slab_depth_offset_km is None: fail( "Slab plotting requires explicit --slab_depth_offset_km; " "implicit inheritance from --seis_depth_offset_km is disabled." ) for option_name, option_value in [ ("--map_segment1_outline_linewidth", args.map_segment1_outline_linewidth), ("--map_segment2_outline_linewidth", args.map_segment2_outline_linewidth), ]: if option_value is not None and option_value <= 0.0: fail(f"{option_name} must be positive") legend_scale_values = { "--legend_scale_hinv": args.legend_scale_hinv, "--legend_scale_ddrt30": args.legend_scale_ddrt30, "--legend_scale_repeating_hinv": args.legend_scale_repeating_hinv, "--legend_scale_repeating_ddrt30": args.legend_scale_repeating_ddrt30, "--legend_scale_target": args.legend_scale_target, } for option_name, option_value in legend_scale_values.items(): if option_value <= 0.0: fail(f"{option_name} must be positive") if args.slip_contour_interval_cm <= 0.0: fail("--slip_contour_interval_cm must be positive") if args.slip_contour_linewidth <= 0.0: fail("--slip_contour_linewidth must be positive") if not 0.0 <= args.slip_contour_alpha <= 1.0: fail("--slip_contour_alpha must be between 0 and 1") if args.slip_contour_label_size <= 0.0: fail("--slip_contour_label_size must be positive") if args.seis_legend_marker_size <= 0.0: fail("--seis_legend_marker_size must be positive") if args.map_legend_marker_size <= 0.0: fail("--map_legend_marker_size must be positive") if args.seis_legend_markerscale <= 0.0: fail("--seis_legend_markerscale must be positive") if args.map_legend_markerscale <= 0.0: fail("--map_legend_markerscale must be positive") df = read_combined_file(args.combined_file) blocks = split_blocks(df, args) num_time_windows = resolve_num_time_windows(args) args.num_time_windows_resolved = num_time_windows args.frw_validation = None if args.frw_file: layout = build_frw_layout(args, num_time_windows) global frw_windows_global frw_windows_global = read_frw_file(args.frw_file, layout) frw_totals = sum_frw_time_windows(frw_windows_global) args.frw_validation = validate_frw_against_combined(frw_totals, blocks, args) print(f"# FRW validation passed: {args.frw_file}") for block_name, summary in args.frw_validation.items(): print( f"# FRW {block_name}: windows={summary['window_count']} " f"max_abs_diff={summary['max_abs_diff']:.6g} " f"max_rel_diff={summary['max_rel_diff']:.6g} " f"zero_mask_mismatches={summary['zero_mask_mismatches']}" ) args.frw_patch_summaries = [] else: args.frw_patch_summaries = [] if any([args.frw_patch1, args.frw_patch2, args.frw_patch3]): fail("FRW patch generation requires --frw_file") outputs = compute_outputs(blocks, args) if args.frw_file: args.frw_patch_summaries = generate_numbered_frw_patches( args, layout, frw_windows_global, outputs ) metrics = compute_summary_metrics(outputs, args) print_summary(df, blocks, outputs, args, metrics) seismicity_map = None seismicity = None need_seismicity = args.plot_seismicity or (args.make_map_plot and args.map_plot_seismicity) if need_seismicity: if not args.seis_catalog: fail("seismicity plotting requires --seis_catalog") seismicity_map = read_and_filter_seismicity(args) write_event_segment_diagnostics( seismicity_map, outputs, args, args.seis_candidate_output, "seismicity" ) if args.plot_seismicity: seismicity = filter_seismicity_by_across_distance(seismicity_map, outputs, args) if args.seis_output: seismicity.to_csv(args.seis_output, index=False) print(f"# wrote {args.seis_output}") repeating_map = None repeating = None need_repeating = args.plot_repeating or (args.make_map_plot and args.map_plot_repeating) if need_repeating: if not args.repeating_catalog: fail("repeating-earthquake plotting requires --repeating_catalog") if not args.seis_catalog: fail("repeating-earthquake plotting requires --seis_catalog because locations are taken from the earthquake catalog") repeating_map = read_and_filter_repeating(args) write_event_segment_diagnostics( repeating_map, outputs, args, args.repeating_candidate_output, "repeating earthquakes", ) if args.plot_repeating: repeating = filter_repeating_by_across_distance(repeating_map, outputs, args) if args.repeating_output: repeating.to_csv(args.repeating_output, index=False) print(f"# wrote {args.repeating_output}") slab_profile = None need_slab = args.plot_slab or (args.make_map_plot and args.map_plot_slab_trace) if need_slab: if not args.slabdir: fail("Slab plotting requires --slabdir") slab_profile = build_slab_profiles(outputs, args) if args.slab_output and slab_profile is not None and not slab_profile.empty: slab_profile.to_csv(args.slab_output, index=False) print(f"# wrote {args.slab_output}") elif args.slab_output and slab_profile is not None and slab_profile.empty: print(f"# WARNING: Slab2 output not written because no valid profile rows exist: {args.slab_output}", file=sys.stderr) # Always calculate and report target-to-plane distance when enabled. # Plot creation is not required for these projection diagnostics. target_projection = None if args.target_event_cross_section: target_projection = target_event_projection_by_segment(outputs, args) if not args.no_blocks: write_blocks(blocks, args.output_prefix, write_csv=not args.no_csv, write_txt=not args.no_txt) write_outputs(outputs, args.output_prefix, write_csv=not args.no_csv, write_txt=not args.no_txt) write_summary_metrics(metrics, args.output_prefix) if args.make_plot: plot_two_segment_slip_dist( outputs, args, metrics=metrics, seismicity=seismicity, repeating=repeating, slab_profile=slab_profile, target_projection=target_projection, ) if args.make_map_plot: plot_latlon_slip_map( outputs, args, metrics=metrics, seismicity_map=seismicity_map, repeating_map=repeating_map, slab_profile=slab_profile, ) if args.plot_mu_depth: plot_mu_vs_depth(args, blocks=blocks, metrics=metrics) if args.plot_mu_compare: plot_mu_compare(outputs, args, metrics=metrics) if __name__ == "__main__": try: main() except Exception as exc: print("", file=sys.stderr) print("============================================================", file=sys.stderr) print("FF SLIP DIST LATLON ERROR", file=sys.stderr) print("============================================================", file=sys.stderr) print(f"Exception type: {type(exc).__name__}", file=sys.stderr) print(f"Message: {exc}", file=sys.stderr) print("-------------------- traceback -----------------------------", file=sys.stderr) traceback.print_exc(file=sys.stderr) print("============================================================", file=sys.stderr) sys.exit(1)