prolog = """
**PROGRAM**
    gps_day_eq_pipeline.py

**PURPOSE**
    | Orchestrator script that combines earthquake querying, station scoring,
    | and mask generation into a single end-to-end pipeline that, for one GPS day:
    |
    |   a. Finds every earthquake globally above a specified magnitude
    |      (default M4.0) that occurred within the GPS day's 30-hour window.
    |      A "GPS day" runs from 21:00 UTC on the day before the target date
    |      through 03:00 UTC on the day after -- i.e. 30 hours centered on
    |      noon UTC of the target date. This definition is verbatim from
    |      Bertiger et al. (2020), Section 6.4.2 (see METHODOLOGY below).
    |   b. Builds the L1/L2 coseismic + postseismic displacement mask for
    |      each event using the Score class from eq_station_score.py.
    |   c. Scores every station in the supplied station list against every
    |      event using score_stations() from eq_station_score.py (Part 3).
    |   d. Writes, per event: a KMZ file with the coseismic/postseismic mask
    |      contours, and a JSON file with the mask's mathematical parameters
    |      (dmax, nodal planes, fault dimensions, grid density, score
    |      results) -- in a stable per-event JSON schema, so downstream
    |      tooling does not need to change. These
    |      per-event files are kept in their own events/ subfolder (named
    |      gps_day_<date>[_<add_filename>]_M<mag>[_postseismic_included]_events/),
    |      separate from the day-level summary files (see OUTPUT LAYOUT below).
    |      This prevents mixing dozens of per-event KMZ/JSON pairs with the
    |      summary files in the main output directory.
    |   e. Pools every station found to be coseismic OR postseismic across
    |      ALL events that day into a single deduplicated 2xN [lat, lon]
    |      array, optionally saving it to disk.
    |
    | This script does not reimplement scoring math; it imports score_stations(),
    | l1_score(), and the Score class from eq_station_score.py, and
    | query_earthquakes()/EarthquakeEvent from usgs_earthquake_query.py.
    | The additional logic implemented in this module includes:
    |   - GPS-day window construction (gps_day_window())
    |   - KMZ/contour export (export_event_mask_kmz()), based on
    |     geode/pyOkada.py, as eq_station_score.py omits KMZ generation.
    |   - The day-level loop, output writing, and station-pooling logic.

**USAGE**
    # Default: today's GPS day, M4+, gcore stations, current directory
    python gps_day_eq_pipeline.py --date 2024-01-17

    # Custom magnitude floor and station catalog (IGS JSON or NGL .nglx format)
    python gps_day_eq_pipeline.py --date 2024-01-17 --min-mag 4.5 \
        --stations_igs_format gdet_1.nglx

    # Use a custom GipsyX station database to generate the catalog
    python gps_day_eq_pipeline.py --date 2024-01-17 \
        --station_list_db stations.txt

    # Control scoring method, displacement threshold, and output location
    python gps_day_eq_pipeline.py --date 2024-01-17 --method L2 \
        --threshold 0.002 --output-dir ./eq_runs

    # Disable specific outputs (KMZ masks or the pooled station array file)
    python gps_day_eq_pipeline.py --date 2024-01-17 --no-kmz --no-array-file

    # Fast pass: skip per-event USGS detail fetches and suppress progress printing
    python gps_day_eq_pipeline.py --date 2024-01-17 --no-details --quiet

    # Adjust Okada grid density for higher resolution L2 masks
    python gps_day_eq_pipeline.py --date 2024-01-17 --density 200

    # As a library
    from gps_day_eq_pipeline import run_gps_day_pipeline
    result = run_gps_day_pipeline("2024-01-17", min_magnitude=4.0,
                                    stations_igs_path="igs_stations.json",
                                    output_dir=".")
"""

epilog = """
**OUTPUT LAYOUT**
    --output-dir is the directory where all output is placed (default: current
    directory). To avoid collisions between runs with different magnitude
    thresholds or settings, files are prefixed with
    gps_day_<date>[_<add_filename>]_M<min_magnitude>[_postseismic_included]:

        <output_dir>/gps_day_<date>[_<add_filename>]_M<min_magnitude>[_postseismic_included]_summary.json
        <output_dir>/gps_day_<date>[_<add_filename>]_M<min_magnitude>[_postseismic_included]_station_array.txt
        <output_dir>/gps_day_<date>[_<add_filename>]_M<min_magnitude>[_postseismic_included]_stations.json
        <output_dir>/gps_day_<date>[_<add_filename>]_M<min_magnitude>[_postseismic_included]_events/
            <event_id>_mask.kmz            <- per event: coseismic
                                                (blue) + postseismic
                                                (orange) contour
                                                polygons, epicenter
                                                marker.
            <event_id>_stations.json       <- per event: event identity +
                                                score_method +
                                                coseismic_stations +
                                                postseismic_stations.

    Rationale: a global M4+ GPS day typically has 50-70+ events, so putting
    per-event files into a specific events/ subfolder prevents the main
    output directory from becoming cluttered, while the prefix ensures that
    different runs for the same day (e.g., M4.0 vs M5.0) do not overwrite
    each other.

    The run_dir and events_dir keys in run_gps_day_pipeline()'s return dict
    give the actual paths used.

**METHODOLOGY**
    GPS-day window (30 hours, centered on noon UTC, 6-hour total overlap with
    adjacent days -- 3 hours on each end):

    S-score formula (L1 isotropic) and Okada-mask rescaling (L2 anisotropic),
    coefficients a=0.5261, b=-1.1478, POST_SEISMIC_SCALE_FACTOR=1.5:
      Gomez, D.D., Caccamise II, D.J. & Bevis, M. (2024). An automated
      approach to identify event-related deformation in continuous GNSS
      networks: an examination of GipsyX coseismic and postseismic
      displacement parameters. GPS Solutions 28:214.
      https://doi.org/10.1007/s10291-024-01758-9

    Okada (1985) elastic half-space model, Score class internals, and the
    contour-extraction / KMZ-export approach (find_contours threshold,
    shapely outer/hole polygon separation, simplekml polygon styling):
      geode/pyOkada.py by Demian Gomez, https://github.com/demiangomez/geode
      (Score.save_masks(), lines 439-581). Ported here rather than imported
      because geode_pyOkada.py has hard relative imports on JPL-internal
      modules (pyDate, pyETM, dbConnection) that are not available outside
      the GeoDE package, and because eq_station_score.py intentionally
      excluded KMZ generation from its scope (see that file's prolog).

    Okada, Y. (1985). Surface deformation due to shear and tensile faults in
      a half-space. Bulletin of the Seismological Society of America,
      75(4), 1135-1154.

    Wells, D.L. & Coppersmith, K.J. (1994). New empirical relationships among
      magnitude, rupture length, rupture width, rupture area, and surface
      displacement. Bulletin of the Seismological Society of America,
      84(4), 974-1002.

**AUTHORS**
    | Hasani Spann
    | Jet Propulsion Laboratory, California Institute of Technology
    | Pasadena, CA, USA
    | Assistance from ChatHPC

**COPYRIGHT**
    | Copyright 2025, by the California Institute of Technology
    | United States Government Sponsorship acknowledged
    | All rights reserved

Dependencies:
  numpy
  scipy                     (KDTree, via eq_station_score.py)
  requests                  (via usgs_earthquake_query.py)
  simplekml, shapely, scikit-image   (KMZ export only -- pipeline still runs
                                       without them, with KMZ output skipped
                                       and a warning printed)
  eq_station_score.py       (must be in the same directory or on PYTHONPATH)
  usgs_earthquake_query.py  (must be in the same directory or on PYTHONPATH)

Keywords: gnss, gipsy, gipsyx, earthquake, coseismic, postseismic, s-score,
          gps day, okada, kmz, orchestrator, overhead, pipeline
"""

import argparse
import json
import os
import shutil
import sys
from datetime import datetime, timedelta, timezone
from typing import Optional

import numpy as np
from gcore.StationDataBase import StationDataBase
from gcore.dateUtils import cal2sec

try:
    import simplekml
    from shapely.geometry import Polygon
    from skimage.measure import find_contours
    _KMZ_AVAILABLE = True
except ImportError:
    _KMZ_AVAILABLE = False

from eq_mask_L2 import compute_l2_mask

# ---------------------------------------------------------------------------
# Import Part 1 (USGS query) and Part 3 (S-score engine)
# ---------------------------------------------------------------------------
try:
    from usgs_earthquake_query import (
        query_earthquakes,
        events_to_dict,
        EarthquakeEvent,
    )
except ImportError as exc:
    sys.exit(
        "ERROR: could not import usgs_earthquake_query.py. It must be in "
        f"the same directory or on PYTHONPATH. Original error: {exc}\n"
        "NOTE: usgs_earthquake_query.py currently imports geopandas at "
        "module level even though it is unused -- if geopandas is not "
        "installed, this import will fail. See this file's epilog, "
        "'KNOWN ISSUE IN A DEPENDENCY'."
    )

try:
    from eq_station_score import (
        score_stations,
        Score,
        l1_score,
        compute_l1_radius,
        distance_km,
        load_stations,
        inv_azimuthal,
        A as L1_A,
        B as L1_B,
        POST_SEISMIC_SCALE_FACTOR,
        GRID_DENSITY,
    )
except ImportError as exc:
    sys.exit(
        "ERROR: could not import eq_station_score.py. It must be in the "
        f"same directory or on PYTHONPATH. Original error: {exc}"
    )

# ---------------------------------------------------------------------------
# Station database exporter for pipeline input
#
# Reads geodetic coordinates from a GipsyX StationDataBase and exports them
# to a JSON file that follows the IGS schema expected by load_stations().
# This ensures a consistent station catalog for a specific epoch is placed
# directly in the run folder.
# ---------------------------------------------------------------------------

def export_stadb_to_pipeline_run(
    date,
    min_magnitude,
    output_dir,
    epoch_date,
    stations=None,
    db_path='$GOA_VAR/sta_info/sta_db_qlflinn',
    staList=None,
    add_filename=None
):
    """
    Retrieve station locations from a StationDataBase and write them to a
    JSON file compatible with the pipeline's load_stations() function.
    
    The output is placed directly in the output_dir using a prefix
    (gps_day_<date>[_<add_filename>]_M<mag>[_postseismic_included]) to
    maintain synchronization with other pipeline outputs.

    Parameters
    ----------
    date : str
        Target calendar date, "YYYY-MM-DD".
    min_magnitude : float
        Minimum event magnitude (used to construct the output filename prefix).
    output_dir : str
        Parent output directory.
    epoch_date : str
        The epoch for the station coordinates, "YYYY-MM-DD".
    stations : list of str, optional
        List of station names to export. If None, all stations in the 
        database are retrieved.
    db_path : str
        Path to the station database file. Defaults to the qlflinn 
        database via the $GOA_VAR environment variable.

    Returns
    -------
    output_file : str
        The full path to the generated JSON stations file.
    """

    # 1. Resolve file path per pipeline convention
    mag_tag = f"{min_magnitude:g}"
    filename_tag = f"_{add_filename}" if add_filename else ""
    output_file = os.path.join(output_dir, f"gps_day_{date}{filename_tag}_M{mag_tag}_stations.json")
    os.makedirs(output_dir, exist_ok=True)
    
    # 2. Resolve database path and convert date string to y2k seconds
    resolved_db_path = os.path.expandvars(db_path)
    y2k_epoch = cal2sec(epoch_date)

    # Parse staList file if provided: extract first 4 chars of each line as the station ID
    parsed_sta_list = None
    if staList:
        resolved_staList_path = os.path.expandvars(staList)
        try:
            with open(resolved_staList_path, 'r') as f:
                parsed_sta_list = [line[:4] for line in f if line.strip()]
        except Exception as e:
            print(f"[StadbExport] WARNING: Could not parse staList file {resolved_staList_path}: {e}")

    db = StationDataBase(dataBase=resolved_db_path, t=y2k_epoch, staList=parsed_sta_list)
    
    # 3. Default to all stations if no specific list is provided
    if stations is None:
        stations = db.getValidStations(y2k_epoch)
    
    # 4. Retrieve geodetic coordinates (Lat, Lon, Height)
    lats, lons, heights = db.dumpLatLonHeights(stations)
    
    # 5. Format as IGS JSON schema: {"name": "STN", "llh": [lat, lon, height]}
    pipeline_stations = []
    for stn, lat, lon, h in zip(stations, lats, lons, heights):
        pipeline_stations.append({
            "name": stn,
            "llh": [lat, lon, h]
        })
    
    # 6. Write to output directory
    with open(output_file, 'w') as f:
        json.dump(pipeline_stations, f, indent=2)
    
    print(f"[StadbExport] Exported {len(pipeline_stations)} stations to {output_file}")
    return output_file

def load_stations(path: str) -> list:
    """
    Load a station list from either a JSON file or a plain-text,
    whitespace-delimited file. Format is detected by file extension:

    - .json          ->  IGS JSON schema (array of objects with 'llh' field)
    - ANYTHING ELSE  ->  plain-text, whitespace-delimited: <id> <lat> <lon>

    Plain-text rules:
      - One station per line; columns separated by any run of whitespace.
      - First 3 tokens are used (id, lat, lon in decimal degrees).
      - Blank lines and lines beginning with '#' are skipped.
      - Malformed lines (fewer than 3 tokens or non-float coordinates) are 
        skipped with a warning.
      - Height, country, and city default to 0.0 or empty strings.

    Returns a list of dicts with keys: name, lat, lon, height, country, city.
    """
    ext = os.path.splitext(path)[1].lower()

    if ext == ".json":
        with open(path, "r") as f:
            raw = json.load(f)
        stations = []
        for stn in raw:
            llh = stn.get("llh")
            if not llh or len(llh) < 2 or llh[0] is None or llh[1] is None:
                continue
            stations.append({
                "name":    stn.get("name", ""),
                "lat":     float(llh[0]),
                "lon":     float(llh[1]),
                "height":  float(llh[2]) if len(llh) > 2 and llh[2] is not None else 0.0,
                "country": stn.get("country", ""),
                "city":    stn.get("city", ""),
            })
        return stations

    else:
        # Any non-.json extension (.txt, .nglx, .dat, none, ...) -> plain text
        stations = []
        with open(path, "r") as f:
            for lineno, line in enumerate(f, 1):
                line = line.strip()
                if not line or line.startswith("#"):
                    continue
                parts = line.split()
                if len(parts) < 3:
                    continue
                try:
                    stations.append({
                        "name":    parts[0],
                        "lat":     float(parts[1]),
                        "lon":     float(parts[2]),
                        "height":  0.0,
                        "country": "",
                        "city":    "",
                    })
                except ValueError:
                    print(f"[Stations] WARNING: skipping malformed line "
                          f"{lineno}: {line!r}", flush=True)
        return stations



# -------------------------------------------------------------------------
# GPS-day window construction
# Source: Bertiger et al. (2020), Section 6.4.2, p. 482 (see epilog).
# ---------------------------------------------------------------------------

def gps_day_window(date_str: str) -> tuple[str, str]:
    """
    Compute the 30-hour GPS-day window for a given calendar date, per
    Bertiger et al. (2020): 30 hours of data centered on noon UTC of the
    target date, giving a 6-hour total overlap with adjacent days (3 hours
    bleeding into the prior day, 3 hours bleeding into the next).

    Parameters
    ----------
    date_str : str
        Target date, "YYYY-MM-DD".

    Returns
    -------
    (start_iso, end_iso) : tuple[str, str]
        ISO-8601 UTC start/end timestamps suitable for
        usgs_earthquake_query.query_earthquakes(start_time=..., end_time=...).

    Example
    -------
    >>> gps_day_window("2024-01-17")
    ('2024-01-16T21:00:00', '2024-01-18T03:00:00')
    """
    target_date = datetime.strptime(date_str.strip(), "%Y-%m-%d").replace(
        tzinfo=timezone.utc
    )
    noon = target_date + timedelta(hours=12)
    start = noon - timedelta(hours=15)
    end   = noon + timedelta(hours=15)
    return start.strftime("%Y-%m-%dT%H:%M:%S"), end.strftime("%Y-%m-%dT%H:%M:%S")


# ---------------------------------------------------------------------------
# Mask parameter export + KMZ contour export
#
# This function adds the KMZ/contour export that eq_station_score.py omits.
# It leverages the boolean mask arrays from the Score class for L2 events,
# or constructs simple circular masks for L1 events.
# ---------------------------------------------------------------------------

def export_event_mask_kmz(
    event: "EarthquakeEvent",
    eq_score: Optional[Score],
    kmz_path: str,
    include_postseismic: bool = True,
    threshold_m: float = 1e-3,
) -> bool:
    """
    Write a KMZ file with the coseismic (blue) and postseismic (orange)
    displacement-mask contour polygons for one event, plus an epicenter
    marker pin.

    Parameters
    ----------
    event : EarthquakeEvent
        Used for identity metadata in the KML description balloon (id,
        time, magnitude, depth) and the epicenter marker location.
    eq_score : Score or None
        A Score object already constructed for this event (see
        score_stations() / build_event_score() below). If None (e.g. the
        event had no usable nodal planes and an L1-only circular mask is
        being exported instead), an L1 circular mask is built on the fly.
    kmz_path : str
        Output .kmz file path.
    include_postseismic : bool
        Whether to also contour and export the postseismic (1.5x) mask.

    Returns
    -------
    bool
        True if the KMZ was written, False if KMZ libraries are unavailable
        or mask grids could not be constructed.

    Source: geode/pyOkada.py Score.save_masks(), lines 439-581 (Demian Gomez,
    https://github.com/demiangomez/geode). The beachball focal-mechanism
    image overlay from that original implementation (which requires obspy)
    is intentionally omitted here to keep this script's dependency list
    light; everything else -- contour extraction, outer/hole polygon
    separation, coseismic/postseismic color coding, epicenter marker -- is
    preserved.
    """
    if not _KMZ_AVAILABLE:
        print(
            "[KMZ] WARNING: simplekml / shapely / scikit-image not "
            "installed -- skipping KMZ export. Install with: "
            "pip install simplekml shapely scikit-image",
            flush=True,
        )
        return False

    # --- Resolve the mask grids to contour ---
    if eq_score is not None:
        c_mx, c_my, c_mask = eq_score.c_mx, eq_score.c_my, eq_score.c_mask
        p_mx, p_my, p_mask = eq_score.p_mx, eq_score.p_my, eq_score.p_mask
        eq_lat, eq_lon = eq_score.lat, eq_score.lon
    else:
        # L1-only fallback: build a simple circular mask grid in km,
        # mirroring eq_station_score.py's Score.__init__ L1 fallback path.
        dmax = compute_l1_radius(event.magnitude, threshold_m=threshold_m)
        xmax = dmax * 2.0
        gx, gy = np.meshgrid(
            np.linspace(-xmax, xmax, GRID_DENSITY),
            np.linspace(-xmax, xmax, GRID_DENSITY),
        )
        c_mask = np.sqrt(gx ** 2 + gy ** 2) < dmax
        p_mask = np.sqrt(gx ** 2 + gy ** 2) < POST_SEISMIC_SCALE_FACTOR * dmax
        c_mx, c_my, p_mx, p_my = gx, gy, gx, gy
        eq_lat, eq_lon = event.latitude, event.longitude

    # --- inv_azimuthal: Cartesian (km) -> lon/lat. Source: geode/pyOkada.py
    #     lines 130-138. No external dependency beyond numpy. ---
    def _inv_azimuthal(x, y, lon0, lat0):
        cosd = lambda v: np.cos(np.deg2rad(v))
        sind = lambda v: np.sin(np.deg2rad(v))
        asind = lambda v: np.rad2deg(np.arcsin(v))
        atand = lambda v: np.rad2deg(np.arctan(v))

        r = np.sqrt(np.square(x) + np.square(y)).flatten()
        c = r / 6371.0
        r_safe = np.where(r == 0, 1e-12, r)  # guard against div-by-zero at origin

        i_lat = asind(
            np.cos(c) * sind(lat0)
            + y.flatten() * np.sin(c) * cosd(lat0) / r_safe
        )
        i_lon = lon0 + atand(
            (x.flatten() * np.sin(c))
            / (r_safe * cosd(lat0) * np.cos(c) - y.flatten() * sind(lat0) * np.sin(c))
        )
        return i_lon, i_lat

    c_lon, c_lat = _inv_azimuthal(c_mx, c_my, eq_lon, eq_lat)
    c_lon = np.reshape(c_lon, c_mask.shape)
    c_lat = np.reshape(c_lat, c_mask.shape)

    if include_postseismic:
        p_lon, p_lat = _inv_azimuthal(p_mx, p_my, eq_lon, eq_lat)
        p_lon = np.reshape(p_lon, p_mask.shape)
        p_lat = np.reshape(p_lat, p_mask.shape)

    # --- Contour extraction. Threshold 1e-16 on a boolean mask isolates the
    #     True/False boundary. Source: geode/pyOkada.py lines 458-460. ---
    coseismic_contours = find_contours(c_mask.astype(float), 0.5)
    postseismic_contours = (
        find_contours(p_mask.astype(float), 0.5) if include_postseismic else []
    )

    kml = simplekml.Kml()

    epicenter = kml.newpoint(name=event.place or event.event_id,
                              coords=[(eq_lon, eq_lat)])
    epicenter.style.iconstyle.icon.href = (
        "https://maps.google.com/mapfiles/kml/shapes/star.png"
    )
    epicenter.style.iconstyle.scale = 1.5
    epicenter.style.iconstyle.color = simplekml.Color.yellow
    epicenter.style.labelstyle.scale = 0

    depth_str = f"{event.depth_km:.0f} km" if event.depth_km is not None else "N/A"
    time_str = (
        event.time_utc.strftime("%Y-%m-%d %H:%M:%S") if event.time_utc else "N/A"
    )
    epicenter.description = (
        '<table style="border-collapse: collapse;">'
        f'<tr><td style="font-weight:bold;">ID:</td><td>{event.event_id}</td></tr>'
        f'<tr><td style="font-weight:bold;">Date:</td><td>{time_str}</td></tr>'
        f'<tr><td style="font-weight:bold;">Mag:</td><td>{event.magnitude:.1f}</td></tr>'
        f'<tr><td style="font-weight:bold;">Depth:</td><td>{depth_str}</td></tr>'
        '<tr><td colspan="2">'
        f'<a href="https://earthquake.usgs.gov/earthquakes/eventpage/{event.event_id}" '
        'target="_blank">USGS event page</a></td></tr>'
        "</table>"
    )

    def _process_contours(contours, lon_grid, lat_grid, color, name):
        """Source: geode/pyOkada.py process_contours(), lines 519-564."""
        polygons = []
        for contour in contours:
            row_idx = contour[:, 0].astype(int)
            col_idx = contour[:, 1].astype(int)
            lon = lon_grid[row_idx, col_idx]
            lat = lat_grid[row_idx, col_idx]
            poly = Polygon(np.column_stack((lon, lat)))
            if poly.is_valid and poly.area > 0.0025:
                # Ignore polygons smaller than ~10 km^2 (expressed in deg^2),
                # same threshold as geode/pyOkada.py.
                polygons.append(poly)

        outer_polys = []
        for poly in polygons:
            is_hole = any(outer.contains(poly) for outer in outer_polys)
            if not is_hole:
                outer_polys.append(poly)

        for outer in outer_polys:
            kml_poly = kml.newpolygon(
                name=name, outerboundaryis=list(outer.exterior.coords)
            )
            kml_poly.style.linestyle.color = color
            kml_poly.style.linestyle.width = 3
            kml_poly.style.polystyle.color = simplekml.Color.changealphaint(
                0, simplekml.Color.white
            )

    _process_contours(
        coseismic_contours, c_lon, c_lat, simplekml.Color.blue, "Coseismic mask"
    )
    if include_postseismic:
        _process_contours(
            postseismic_contours, p_lon, p_lat, simplekml.Color.orange,
            "Postseismic mask",
        )

    kml.savekmz(kmz_path)
    return True


def build_event_score(
    event: "EarthquakeEvent", method: str, density: int, threshold_m: float = 1e-3
) -> tuple:
    """
    Construct a Score object for one event (reusing eq_station_score.py's
    nodal-plane extraction logic so this stays in sync with score_stations()
    rather than duplicating that logic a third time).

    Returns
    -------
    (eq_score, score_method, nodal_planes_used) : tuple[Score or None, str, bool]
        eq_score is None for the pure-L1 case (no Score object needed; the
        caller falls back to the on-the-fly circular mask inside
        export_event_mask_kmz()).
    """
    strikes, dips, rakes = [], [], []
    np1, np2 = event.nodal_plane_1, event.nodal_plane_2

    if np1.strike is not None and np1.dip is not None and np1.rake is not None:
        strikes.append(float(np1.strike))
        dips.append(float(np1.dip))
        rakes.append(float(np1.rake))
    if np2.strike is not None and np2.dip is not None and np2.rake is not None:
        strikes.append(float(np2.strike))
        dips.append(float(np2.dip))
        rakes.append(float(np2.rake))

    nodal_planes_used = len(strikes) > 0
    use_l2 = nodal_planes_used and method in ("L2", "auto")

    if method == "L2" and not nodal_planes_used:
        raise ValueError(
            f"method='L2' requested but event {event.event_id} has no "
            "usable nodal-plane data."
        )

    if not use_l2:
        score_method = "L1_isotropic" if method == "L1" else "L1_isotropic_fallback"
        return None, score_method, nodal_planes_used

    depth_km = event.depth_km if event.depth_km is not None else 10.0
    eq_score = Score(
        event_lat=event.latitude,
        event_lon=event.longitude,
        depth_km=depth_km,
        magnitude=event.magnitude,
        strikes=strikes,
        dips=dips,
        rakes=rakes,
        density=density,
        threshold_m=threshold_m,
    )
    return eq_score, "L2_okada", nodal_planes_used


def build_mask_parameters(
    event: "EarthquakeEvent", eq_score: Optional[Score], score_method: str, threshold_m: float = 1e-3
) -> dict:
    """
    Collect the mathematical mask parameters for one event into a plain
    dict, for the "outputs the mathematical parameters for it" requirement
    of outline Part 2.
    """
    # Base L1 parameters
    dmax = 10.0 ** (L1_A * event.magnitude + L1_B)
    params = {
        "score_method": score_method,
        "post_seismic_scale_factor": POST_SEISMIC_SCALE_FACTOR,
    }

    if score_method == "L2_okada":
        l2_meta = compute_l2_mask(event, threshold_m=threshold_m)
        if 'error' not in l2_meta:
            params.update({
                "coseismic_min_radius_km": l2_meta.get("coseismic_min_radius_km"),
                "coseismic_max_radius_km": l2_meta.get("coseismic_max_radius_km"),
                "postseismic_min_radius_km": l2_meta.get("postseismic_min_radius_km"),
                "postseismic_max_radius_km": l2_meta.get("postseismic_max_radius_km"),
                "threshold_m": l2_meta.get("threshold_m"),
            })
        if eq_score is not None:
            params.update({
                "strikes_deg": eq_score.strikes,
                "dips_deg": eq_score.dips,
                "rakes_deg": eq_score.rakes,
                "along_strike_fault_length_m": round(eq_score.along_strike_m, 1),
                "downdip_fault_width_m": round(eq_score.downdip_m, 1),
                "average_slip_m": round(eq_score.avg_disp_m, 4),
                "grid_density": int(eq_score.c_mx.shape[0]),
            })
    else:
        # L1 parameters
        params.update({
            "l1_dmax_km": round(dmax, 3),
            "l1_dmax_postseismic_km": round(dmax * POST_SEISMIC_SCALE_FACTOR, 3),
            "l1_coefficients": {"a": L1_A, "b": L1_B},
        })

    return params


# ---------------------------------------------------------------------------
# Main day-level orchestration (Part 4 entry point)
# ---------------------------------------------------------------------------

def run_gps_day_pipeline(
    date: str,
    min_magnitude: Optional[float] = 4.0,
    stations_igs_path: Optional[str] = None,
    station_db: Optional[str] = None,
    station_list_db: Optional[str] = None,
    method: Optional[str] = "auto",
    output_dir: Optional[str] = ".",
    density: Optional[int] = GRID_DENSITY,
    threshold_m: Optional[float] = 1e-3,
    write_kmz: Optional[bool] = True,
    write_array_file: Optional[bool] = True,
    fetch_details: Optional[bool] = True,
    verbose: Optional[bool] = True,
    clean: Optional[bool] = False,
    include_postseismic: Optional[bool] = False,
    add_filename: Optional[str] = None,
) -> dict:
    """
    Run the full end-to-end pipeline for one GPS day.

    Parameters
    ----------
    date : str
        Target calendar date, "YYYY-MM-DD". The actual USGS query window is
        the 30-hour GPS day centered on noon UTC of this date (see
        gps_day_window()).
    min_magnitude : float
        Minimum event magnitude (outline Part 4a specifies M4.0 as the
        global default).
    stations_igs_path : str
        Path to a station catalog. If a .json file is provided, it is read
        as an IGS JSON schema; any other extension (.txt, .nglx, .dat, etc.)
        is read as plain-text (<id> <lat> <lon>). Auto-detected via load_stations().
    station_db : str
        Path to the GipsyX StationDataBase file (the 'dataBase' input for gcore.StationDataBase).
        Defaults to the qlflinn database via the $GOA_VAR environment variable.
    station_list_db : str
        Path to a text file containing a list of station names (one per line).
        This is used as the 'staList' input for gcore.StationDataBase to filter
        the database when generating the IGS catalog.
    method : {'auto', 'L1', 'L2'}
        Scoring method passed through to score_stations() and used to
        decide whether build_event_score() builds an Okada Score object.
    output_dir : str
        The directory where all output is placed. All files are prefixed with
        gps_day_<date>_M<min_magnitude> (see OUTPUT LAYOUT in this module's
        epilog for the full tree). Created if it does not exist. Defaults to
        the current directory, per outline Part 4d.
    density : int
        Okada grid density per axis (L2/auto with nodal planes only).
    write_kmz : bool
        If False, skip KMZ generation entirely (outline Part 4d marks KMZ
        output as optional -- "(or KMZ file)" / "(possibly saving...)").
    write_array_file : bool
        If False, skip writing the pooled 2xN station array to disk (it is
        always returned in the result dict regardless).
    fetch_details : bool
        Passed through to usgs_earthquake_query.query_earthquakes(). When
        True (default), every event triggers a SECOND HTTP request to fetch
        its moment-tensor / nodal-plane data -- this is what makes L2/auto
        scoring possible, but on a global M4+ day (~50-70 events over a
        30-hour GPS day) it means ~50-70 sequential detail requests, each
        with up to a REQUEST_TIMEOUT-second budget. If the pipeline appears
        to hang with no new output for a long time, this is almost always
        why. Set False for a much faster pass that scores every event with
        the L1 isotropic method only (no nodal planes are fetched, so
        build_event_score() always falls back to L1_isotropic_fallback
        regardless of --method).
    verbose : bool
        Print per-event progress.

    Returns
    -------
    dict with keys:
        date, gps_day_start_utc, gps_day_end_utc, min_magnitude,
        n_events, events (list of per-event result dicts, each containing
            event identity fields + score_method + mask_parameters +
            coseismic_stations + postseismic_stations),
        station_array (ndarray, shape (N, 2), columns = [lat, lon] for every
            UNIQUE station found coseismic OR postseismic in ANY event that
            day), station_array_names (list[str], parallel to station_array
            rows), run_dir (str, the actual folder everything was written
            into -- see OUTPUT LAYOUT above), events_dir (str, the
            per-event KMZ/JSON subfolder inside run_dir).
    """
    # --- Define output paths using a consistent prefix to avoid collisions ---
    # Prefix encodes the date and magnitude floor to keep runs separate.
    mag_tag = f"{min_magnitude:g}"
    filename_tag = f"_{add_filename}" if add_filename else ""
    prefix = f"gps_day_{date}{filename_tag}_M{mag_tag}"
    if include_postseismic:
        prefix += "_postseismic_included"
    events_dir = os.path.join(output_dir, f"{prefix}_events")
    os.makedirs(events_dir, exist_ok=True)

    start_iso, end_iso = gps_day_window(date)
    if verbose:
        print(f"[GPS day] {date}  ->  window {start_iso}  to  {end_iso}  UTC",
              flush=True)
        print(f"[GPS day] Output directory: {os.path.abspath(output_dir)}", flush=True)

    # --- Part 4a: global earthquake query over the GPS-day window ---
    if verbose and fetch_details:
        print("[GPS day] fetch_details=True -- this fetches a SECOND request "
              "per event for moment-tensor data. A global M4+ day typically "
              "has 50-70+ events over the 30-hour window, so this step can "
              "take a while with one 'Parsing event i/N' line printed per "
              "event as it goes. Pass --no-details for a much faster "
              "L1-only pass if you don't need Okada/L2 scoring.", flush=True)
    events = query_earthquakes(
        start_time=start_iso,
        end_time=end_iso,
        min_magnitude=min_magnitude,
        fetch_details=fetch_details,
    )
    if verbose:
        print(f"[Events] {len(events)} earthquake(s) >= M{min_magnitude} found.",
              flush=True)

    # --- Load station catalog once (Part 1/4 file I/O) ---
    if stations_igs_path is None:
        # Dynamically generate the station file from the database for this specific run
        # We use the pipeline's current 'date' as the epoch for the station coordinates
        stations_igs_path = export_stadb_to_pipeline_run(
            date=date,
            min_magnitude=min_magnitude,
            output_dir=output_dir,
            epoch_date=date,
            db_path=station_db if station_db else '$GOA_VAR/sta_info/sta_db_qlflinn',
            staList=station_list_db,
            add_filename=add_filename
        )


    if not os.path.isfile(stations_igs_path):
        sys.exit(f"ERROR: station file not found: {stations_igs_path}")
    stations = load_stations(stations_igs_path)

    if verbose:
        print(f"[Stations] {len(stations)} station(s) loaded from {stations_igs_path}.",
              flush=True)

    day_results = []
    pooled = {}   # name -> (lat, lon, max_disp), deduplicated across all events

    for i, ev in enumerate(events, 1):
        if ev.latitude is None or ev.longitude is None or ev.magnitude is None:
            print(f"  [SKIP] Event {i}/{len(events)} {ev.event_id} is missing "
                  "location/magnitude -- cannot score.", flush=True)
            continue

        if verbose:
            place_str = str(ev.place) if ev.place else "Unknown"
            try:
                print(f"  [{i}/{len(events)}] {ev.event_id}  M{ev.magnitude}  "
                      f"{place_str}", flush=True)
            except UnicodeEncodeError:
                # Fallback for terminals that cannot handle certain Unicode characters
                encoding = sys.stdout.encoding or 'utf-8'
                sanitized_place = place_str.encode(encoding, errors='replace').decode(encoding)
                print(f"  [{i}/{len(events)}] {ev.event_id}  M{ev.magnitude}  "
                      f"{sanitized_place}", flush=True)

        # --- Part 4b: build the mask (Score object) for this event ---
        try:
            eq_score, score_method, nodal_planes_used = build_event_score(
                ev, method=method, density=density, threshold_m=threshold_m
            )
        except ValueError as exc:
            print(f"    [SKIP] {exc}", flush=True)
            continue

        mask_params = build_mask_parameters(ev, eq_score, score_method, threshold_m=threshold_m)

        # --- Part 4c: score every station against this event ---
        score_result = score_stations(
            ev, stations, method=method, include_postseismic=include_postseismic,
            density=density, threshold_m=threshold_m,
        )

        # --- Part 4d: per-event KMZ + JSON output, into events/ ---
        event_dict = {
            "event_id": ev.event_id,
            "place": ev.place,
            "time_utc": ev.time_utc.isoformat() if ev.time_utc else None,
            "latitude": ev.latitude,
            "longitude": ev.longitude,
            "depth_km": ev.depth_km,
            "magnitude": ev.magnitude,
            "magnitude_type": ev.magnitude_type,
            "mask_parameters": mask_params,
            **score_result,
        }

        json_path = os.path.join(events_dir, f"{ev.event_id}_stations.json")
        with open(json_path, "w") as f:
            json.dump(event_dict, f, indent=2)

        kmz_written = False
        if write_kmz:
            kmz_path = os.path.join(events_dir, f"{ev.event_id}_mask.kmz")
            kmz_written = export_event_mask_kmz(
                ev, eq_score, kmz_path, include_postseismic=include_postseismic, threshold_m=threshold_m
            )

        if verbose:
            print(f"    -> {score_result['n_coseismic']} coseismic, "
                  f"{score_result['n_postseismic']} postseismic "
                  f"[{score_result['score_method']}]"
                  f"{'  (KMZ written)' if kmz_written else ''}", flush=True)

        day_results.append(event_dict)

        # --- Part 4e: pool coseismic + postseismic stations across events ---
        for stn in score_result["coseismic_stations"] + score_result["postseismic_stations"]:
            name = stn["name"]
            disp = stn.get("displacement_m", 0.0)
            if name not in pooled or disp > pooled[name][2]:
                pooled[name] = (stn["lat"], stn["lon"], disp)

    # --- Build the 2xN [lat, lon] array (outline Part 4e) ---
    if pooled:
        names = sorted(list(pooled.keys()))
        station_array = np.array([pooled[n] for n in names], dtype=float)
    else:
        names = []
        station_array = np.empty((0, 3), dtype=float)

    if verbose:
        print(f"[Pooled] {len(names)} unique station(s) flagged across "
              f"{len(day_results)} event(s) on {date}.", flush=True)

    # --- Write day-level summary + pooled array directly in output_dir ---
    # (summary-level outputs only -- per-event files live in events_dir)
    summary_path = os.path.join(output_dir, f"{prefix}_summary.json")
    with open(summary_path, "w") as f:
        json.dump({
            "date": date,
            "gps_day_start_utc": start_iso,
            "gps_day_end_utc": end_iso,
            "min_magnitude": min_magnitude,
            "n_events": len(day_results),
            "events": day_results,
        }, f, indent=2)

    if write_array_file:
        array_path = os.path.join(output_dir, f"{prefix}_station_array.txt")
        with open(array_path, "w") as f:
            f.write("# StationName Lat Lon Displacement_m\n")
            for name in names:
                lat, lon, disp = pooled[name]
                f.write(f"{name.upper()} {lat:.6f} {lon:.6f} {disp:.6f}\n")

    if clean:
        if verbose:
            print(f"[Cleanup] Cleaning output directory {output_dir}, preserving station array file...", flush=True)
        
        array_filename = f"{prefix}_station_array.txt"
        for item in os.listdir(output_dir):
            if item.startswith(prefix) and item != array_filename:
                item_path = os.path.join(output_dir, item)
                try:
                    if os.path.isfile(item_path) or os.path.islink(item_path):
                        os.remove(item_path)
                    elif os.path.isdir(item_path):
                        shutil.rmtree(item_path)
                except Exception as e:
                    if verbose:
                        print(f"    [Cleanup] WARNING: Could not delete {item_path}: {e}", flush=True)

    return {
        "date": date,
        "gps_day_start_utc": start_iso,
        "gps_day_end_utc": end_iso,
        "min_magnitude": min_magnitude,
        "n_events": len(day_results),
        "events": day_results,
        "station_array": station_array,
        "station_array_names": names,
        "run_dir": output_dir,
        "events_dir": events_dir,
        "include_postseismic": include_postseismic,
    }


# ---------------------------------------------------------------------------
# CLI
# ---------------------------------------------------------------------------

def build_parser() -> argparse.ArgumentParser:
    p = argparse.ArgumentParser(
        prog="gps_day_eq_pipeline.py",
        description=prolog,
        epilog=epilog,
        formatter_class=argparse.RawDescriptionHelpFormatter,
    )
    p.add_argument(
        "--date", required=True, metavar="YYYY-MM-DD",
        help="Target calendar date. The actual query window is a 30-hour 'GPS day' "
              "centered on noon UTC of this date (extending 15 hours before and after).",
    )
    p.add_argument(
        "--stations_igs_format", default=None,
        help="Station catalog path (IGS JSON or plain-text). Will be auto-generated "
              "from the database if not provided. Overwrites --station_list_db "
              "and --station_db if provided. (default: None).",
    )
    p.add_argument(
        "--station_db", default=None,
        help="Path to the GipsyX StationDataBase file (the 'dataBase' input for gcore.StationDataBase). "
             "Defaults to the qlflinn database via the $GOA_VAR environment variable. (default: None).",
    )
    p.add_argument(
        "--station_list_db", default=None,
        help="Station list path (StationDataBase friendly formatted): Path to a list of station names "
             "(the 'staList' input for gcore.StationDataBase). Will be overwritten if --stations_igs_format is provided."
             "(default: None).",
    )
    p.add_argument(
        "--min-mag", type=float, default=4.0,
        help="Minimum earthquake magnitude (default: 4.0, per outline Part 4a).",
    )
    p.add_argument(
        "--threshold", type=float, default=1e-3,
        help="Displacement threshold for L2 scoring in metres (default: 0.001). Warning: any earthquake using L1 method will ignore this parameter and behave as if it were 0.001, even on the (auto) L2 fallback path.",
    )
    p.add_argument(
        "--method", choices=("auto", "L1", "L2"), default="auto",
        help="Scoring method (default: auto). Note that when --method=L2 is requested, any event without usable nodal-plane data will be skipped with a warning. Use --method=auto to score those events with L1_isotropic_fallback instead." ,
    )
    p.add_argument(
        "--output-dir", default=".",
        help="PARENT directory for this run's output folder (default: "
              "current directory). The actual files land in "
              "<output-dir>/gps_day_<date>[_<add_filename>]_M<min-mag>[_postseismic_included]/, "
              "with per-event KMZ/JSON in an events/ subfolder of that.",
    )
    p.add_argument(
        "--density", type=int, default=GRID_DENSITY,
        help=f"Okada grid density per axis (default: {GRID_DENSITY}).",
    )
    p.add_argument(
        "--no-kmz", action="store_true",
        help="Skip KMZ mask export (JSON output is always written).",
    )
    p.add_argument(
        "--no-array-file", action="store_true",
        help="Skip writing the pooled station array to disk (still "
             "returned when called as a library).",
    )
    p.add_argument(
        "--no-details", action="store_true",
        help="Skip the per-event USGS moment-tensor detail fetch. A global "
             "M4+ GPS day typically has 50-70+ events, and each one normally "
             "triggers a SECOND HTTP request to get nodal-plane data -- this "
             "is almost always why the pipeline appears to hang with no new "
             "output. Use this flag for a much faster run; every event will "
             "be scored L1_isotropic_fallback only (no L2/Okada), regardless "
             "of --method.",
    )
    p.add_argument(
        "--quiet", action="store_true",
        help="Suppress per-event progress printing.",
    )
    p.add_argument(
        "--clean", action="store_true",
        help="Delete all generated files in the run directory at the end, except for the station array file.",
    )
    p.add_argument(
        "--include-postseismic", action="store_true",
        help="Include postseismic masks and stations in the output.",
    )
    p.add_argument(
        "--add-filename", default=None,
        help="Optional string to add to the output filenames (e.g. --add-filename 1cm).",
    )
    return p


def main():
    args = build_parser().parse_args()
    result = run_gps_day_pipeline(
        date=args.date,
        min_magnitude=args.min_mag,
        stations_igs_path=args.stations_igs_format,
        station_db=args.station_db,
        station_list_db=args.station_list_db,
        method=args.method,
        output_dir=args.output_dir,
        density=args.density,
        threshold_m=args.threshold,
        write_kmz=not args.no_kmz,
        write_array_file=not args.no_array_file,
        fetch_details=not args.no_details,
        verbose=not args.quiet,
        clean=args.clean,
        include_postseismic=args.include_postseismic,
        add_filename=args.add_filename,
    )
    print("\n" + "=" * 60)
    post_str = " (including postseismic)" if result['include_postseismic'] else " (coseismic only)"
    print(f"GPS day {result['date']}: {result['n_events']} event(s), "
          f"{len(result['station_array_names'])} unique flagged station(s){post_str}.")
    print(f"Output directory: {os.path.abspath(result['run_dir'])}")
    print(f"Events folder:    {os.path.abspath(result['events_dir'])}")
    print("=" * 60)


if __name__ == "__main__":
    main()