#!/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 sys import subprocess 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.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: return np.full(len(lons), np.nan, dtype=float) vals = [] for line in result.stdout.splitlines(): try: vals.append(float(line.split()[0])) except Exception: vals.append(np.nan) out = np.full(len(lons), np.nan, dtype=float) out[:min(len(vals), len(out))] = vals[:min(len(vals), len(out))] return out 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): errors = [] for backend in ([args.slab_sample_backend] if args.slab_sample_backend != "auto" else ["xarray", "gmt", "pygmt"]): 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: errors.append(f"{backend}: {exc}") if args.slab_sample_backend == backend: raise fail("could not sample slab grid with requested backend(s): " + "; ".join(errors)) 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 = args.seis_depth_offset_km if args.slab_depth_offset_km is None else 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: print("# SKIP: no depth grid file found") continue 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: print("# WARNING: all depth samples are NaN/invalid for this slab ID") 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: print("# WARNING: no valid slab profile samples found") 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 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: xnum = args.xgrid_num2 if args.xgrid_num2 is not None else args.xgrid_num znum = args.zgrid_num2 if args.zgrid_num2 is not None else args.zgrid_num xlen = args.x_length_km2 if args.x_length_km2 is not None else args.x_length_km zlen = args.z_length_km2 if args.z_length_km2 is not None else args.z_length_km 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 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) 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)}") if not tmp.empty: candidates.append(tmp) if not candidates: print(f"# WARNING: no {label} events within +/- {args.seis_max_across_km:g} km across any segment") return eq_df.iloc[0:0].copy() near = pd.concat(candidates, ignore_index=True) if "id" in near.columns and near["id"].notna().any(): key_cols = ["id"] elif "repeating_evid" in near.columns: key_cols = ["csid", "repeating_evid"] else: key_cols = ["lat", "lon", "dep", "mag"] near = near.sort_values("seis_abs_across_km").drop_duplicates(subset=key_cols, keep="first").reset_index(drop=True) print(f"# {label} across-fault filter") print(f"# max_across_km = {args.seis_max_across_km:g}") print(f"# after across filter = {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 lon/lat at an already defined fault-plane model depth.""" target_z = float(target["dep_model"].iloc[0]) zmin = float(seg["z_km"].min()) zmax = float(seg["z_km"].max()) z_tolerance = max(float(seg["z_grid_space_km"].median()) * 1.0e-6, 1.0e-8) inside_z = (target_z >= zmin - z_tolerance) and (target_z <= zmax + z_tolerance) line_rows = [] for x_index, column in seg.groupby("x_index", sort=True): column = column.sort_values("z_km") z = column["z_km"].to_numpy(float) if target_z < z.min() - z_tolerance or target_z > z.max() + z_tolerance: continue line_rows.append({ "x_index": int(x_index), "x_km": float(np.interp(target_z, z, column["x_km"].to_numpy(float))), "lon": float(np.interp(target_z, z, column["lon"].to_numpy(float))), "lat": float(np.interp(target_z, z, column["lat"].to_numpy(float))), }) if not line_rows: return { "x_km": float("nan"), "z_km": target_z, "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 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" ) print("# model-depth interpretation= fault-plane z; no dip conversion") 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 == "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 == "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_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", ) 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_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 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) 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") _, _, 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) ax.set_title( f"Segment {seg_id}\nstrike = {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) ax.grid(False) ax.tick_params(axis="both", which="both", direction="out", top=True, right=True) 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) 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("--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) 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("--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.") 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_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) 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_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("--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() 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.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) outputs = compute_outputs(blocks, args) 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) 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) 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: slab_profile.to_csv(args.slab_output, index=False) print(f"# wrote {args.slab_output}") # 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("") print("============================================================") print("FF SLIP DIST LATLON ERROR") print("============================================================") print(exc) print("============================================================") sys.exit(1)