prolog = """ **PROGRAM** gipsy_eq_station_flag.py **PURPOSE** | Loads station locations from GipsyX output files, queries USGS for | earthquake events over a specified time window, and flags which stations | fall within the coseismic and/or postseismic influence zone of each event. | | This file is infrastructure-complete. The only section requiring | completion is the GipsyX-specific parser inside parse_gipsy_stations(). | Everything else -- USGS querying, L1/L2 scoring, and output -- is fully | operational and imports directly from the existing flag scripts. | | Two scoring modes are available: | L1 -- isotropic epicentral distance (fast, no moment-tensor needed) | L2 -- Okada-based anisotropic mask (requires nodal planes from USGS) | Falls back to L1 automatically if no moment-tensor is available. **USAGE** # L1 scoring, single event python gipsy_eq_station_flag.py --gipsy-file path/to/gipsy_output \\ --event-id us7000abcd --method L1 # L2 scoring, date range python gipsy_eq_station_flag.py --gipsy-file path/to/gipsy_output \\ --date 2024-01-01 --end-date 2024-01-07 --min-mag 5.5 --method L2 # Pre-fetched earthquake JSON python gipsy_eq_station_flag.py --gipsy-file path/to/gipsy_output \\ --eq-json earthquake_events.json --method L1 # Geographic bounding box python gipsy_eq_station_flag.py --gipsy-file path/to/gipsy_output \\ --date 2024-01-01 --min-mag 5.5 --method L2 \\ --min-lat -40 --max-lat -30 --min-lon -75 --max-lon -65 """ epilog = """ **OUTPUTS** JSON file (--output) with one entry per earthquake, each containing: - event identity fields (id, place, time, magnitude, coordinates) - score_method: 'L1_isotropic', 'L2_okada', or 'L1_isotropic_fallback' - coseismic_stations, postseismic_stations: lists of station dicts with name, lat, lon, distance_km, azimuth_deg, sc, sp **GIPSY PARSING -- ACTION REQUIRED** The function parse_gipsy_stations() contains a clearly marked stub. Fill in the parsing logic for your specific GipsyX output format. The function must return a list of dicts with keys: name, lat, lon, height, country, city lat/lon must be in decimal degrees (North/East); height in metres. An ECEF (XYZ) -> lat/lon/height converter is provided as ecef_to_llh(). **METHODOLOGY** L1 scoring: Gomez et al. (2024), GPS Solutions 28:214. https://doi.org/10.1007/s10291-024-01758-9 L2 scoring: same paper, plus Okada (1985) BSSA 75(4):1135-1154 and Wells & Coppersmith (1994) BSSA 84(4):974-1002. USGS API: https://earthquake.usgs.gov/fdsnws/event/1/ **AUTHORS** | Hasani Spann | Jet Propulsion Laboratory, California Institute of Technology | Pasadena, CA, USA **COPYRIGHT** | Copyright 2025, by the California Institute of Technology | United States Government Sponsorship acknowledged | All rights reserved Dependencies: numpy scipy (L2 only, via eq_station_flag_L2) requests (via usgs_earthquake_query) eq_station_flag_L1.py (must be in the same directory or on PYTHONPATH) eq_station_flag_L2.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, station flagging, okada """ import argparse import json import os import sys from datetime import datetime # --------------------------------------------------------------------------- # Import scoring infrastructure from existing flag scripts # --------------------------------------------------------------------------- # We import the individual functions rather than calling main() so we can # pass GipsyX station dicts directly without touching argparse or disk I/O. try: import usgs_earthquake_query as usgs except ImportError: sys.exit( "ERROR: usgs_earthquake_query.py not found. " "Place it in the same directory as this script or add it to PYTHONPATH." ) # --------------------------------------------------------------------------- # GipsyX station parser # --------------------------------------------------------------------------- def parse_gipsy_stations(gipsy_file_path: str) -> list[dict]: """ Parse a GipsyX output file and return a list of station dicts compatible with the flag_stations_l1() and flag_stations_l2() functions. Each returned dict must contain: name : str — station identifier (e.g. 'JPL1', 'ABMF') lat : float — latitude in decimal degrees North (-90 to 90) lon : float — longitude in decimal degrees East (-180 to 180) height : float — ellipsoidal height in metres country : str — optional, used only in output labels ('' if unknown) city : str — optional, used only in output labels ('' if unknown) Parameters ---------- gipsy_file_path : str Path to the GipsyX file containing station locations. The exact file type depends on your GipsyX workflow — common candidates are: - A .tdp (time-dependent position) file - A .stacov (station covariance) file - A site info / station metadata file - A smoothFinal or tdpFinal output file Returns ------- list[dict] One dict per station with the keys described above. Notes ----- GipsyX internal coordinate convention may use XYZ (ECEF metres) rather than lat/lon. If so, convert using the helper below: lat, lon, height = ecef_to_llh(x, y, z) That helper is provided at the bottom of this file. """ stations = [] # ------------------------------------------------------------------ # >>> ACTION REQUIRED: GipsyX-specific parsing goes here <<< # # The block below is a stub. Replace it with code that: # 1. Opens gipsy_file_path and reads station data # 2. Extracts station name, coordinates (lat/lon/height OR XYZ ECEF) # 3. Converts XYZ to lat/lon/height if needed (use ecef_to_llh below) # 4. Appends each station to the `stations` list in the format above # # Example structure (do NOT leave this raise in production): # ------------------------------------------------------------------ raise NotImplementedError( "parse_gipsy_stations() has not been implemented yet.\n" "Open gipsy_eq_station_flag.py and fill in the parsing logic\n" "inside this function for your specific GipsyX file format.\n" f"File to parse: {gipsy_file_path}" ) # ------------------------------------------------------------------ # When implemented, your code will look something like one of these # patterns depending on your GipsyX file format: # # --- Pattern A: whitespace-delimited columns (e.g. .tdp files) --- # # with open(gipsy_file_path, 'r') as f: # for line in f: # if line.startswith('#') or not line.strip(): # continue # skip comments / blank lines # parts = line.split() # name = parts[?] # column index TBD # x_m = float(parts[?]) # ECEF X in metres # y_m = float(parts[?]) # ECEF Y in metres # z_m = float(parts[?]) # ECEF Z in metres # lat, lon, height = ecef_to_llh(x_m, y_m, z_m) # stations.append({ # 'name': name, 'lat': lat, 'lon': lon, # 'height': height, 'country': '', 'city': '' # }) # # --- Pattern B: JSON or XML metadata file --- # # with open(gipsy_file_path, 'r') as f: # data = json.load(f) # for entry in data['stations']: # stations.append({ # 'name' : entry['id'], # 'lat' : float(entry['latitude']), # 'lon' : float(entry['longitude']), # 'height' : float(entry['height']), # 'country': entry.get('country', ''), # 'city' : entry.get('city', ''), # }) # # --- Pattern C: GipsyX .stacov format (if lat/lon are columns) --- # # Similar to Pattern A but with different column layout. # ------------------------------------------------------------------ return stations # --------------------------------------------------------------------------- # ECEF → lat/lon/height conversion helper # Provided for use inside parse_gipsy_stations() if GipsyX gives XYZ coords. # --------------------------------------------------------------------------- def ecef_to_llh(x_m: float, y_m: float, z_m: float) -> tuple[float, float, float]: """ Convert ECEF (Earth-Centred Earth-Fixed) Cartesian coordinates in metres to geodetic latitude, longitude, and ellipsoidal height. Uses the WGS-84 ellipsoid via Bowring's iterative method, which converges to millimetre accuracy in 2-3 iterations for any point on or near Earth. Parameters ---------- x_m, y_m, z_m : float ECEF coordinates in metres. Returns ------- lat : float — geodetic latitude, decimal degrees North (-90 to 90) lon : float — geodetic longitude, decimal degrees East (-180 to 180) height : float — ellipsoidal height, metres References ---------- Bowring, B.R. (1985). The geodesic line and the normal section. Survey Review 28(218): 233-243. WGS-84 ellipsoid parameters: a=6378137.0 m, f=1/298.257223563 """ import math # WGS-84 ellipsoid constants a = 6378137.0 # semi-major axis (metres) f = 1.0 / 298.257223563 # flattening b = a * (1.0 - f) # semi-minor axis e2 = 1.0 - (b / a) ** 2 # first eccentricity squared ep2 = (a / b) ** 2 - 1.0 # second eccentricity squared lon = math.atan2(y_m, x_m) p = math.sqrt(x_m ** 2 + y_m ** 2) # distance from Z-axis theta = math.atan2(z_m * a, p * b) # parametric latitude (initial) # Bowring's iterative formula lat = math.atan2( z_m + ep2 * b * math.sin(theta) ** 3, p - e2 * a * math.cos(theta) ** 3 ) # One refinement iteration (sufficient for geodetic purposes) for _ in range(3): sin_lat = math.sin(lat) N = a / math.sqrt(1.0 - e2 * sin_lat ** 2) # prime vertical radius lat_new = math.atan2(z_m + e2 * N * sin_lat, p) if abs(lat_new - lat) < 1e-12: break lat = lat_new sin_lat = math.sin(lat) N = a / math.sqrt(1.0 - e2 * sin_lat ** 2) height = p / math.cos(lat) - N if abs(math.cos(lat)) > 1e-10 else abs(z_m) - b return ( math.degrees(lat), math.degrees(lon), height, ) # --------------------------------------------------------------------------- # Earthquake loader (shared logic for all three input modes) # --------------------------------------------------------------------------- def load_events(args) -> list[usgs.EarthquakeEvent]: """ Load earthquake events from whichever source the user specified: --event-id : single USGS event fetch --date : USGS time-window query --eq-json : pre-fetched JSON file Returns a list of EarthquakeEvent objects ready for scoring. The fetch_details flag is set based on --method: L2 needs moment-tensor data; L1 does not. """ need_details = (args.method == 'L2') events = [] if args.event_id: print(f"[USGS] Fetching event {args.event_id} ...", flush=True) ev = usgs.fetch_event_by_id(args.event_id) events = [ev] elif args.date: print( f"[USGS] Querying events {args.date}" + (f" → {args.end_date}" if args.end_date else "") + f" min_mag={args.min_mag} fetch_details={need_details} ...", flush=True ) events = usgs.query_earthquakes( start_time = args.date, end_time = args.end_date, min_magnitude = args.min_mag, min_latitude = args.min_lat, max_latitude = args.max_lat, min_longitude = args.min_lon, max_longitude = args.max_lon, fetch_details = need_details, limit = args.limit, ) print(f"[USGS] {len(events)} event(s) returned.", flush=True) elif args.eq_json: if not os.path.isfile(args.eq_json): sys.exit(f"ERROR: Earthquake JSON not found: {args.eq_json}") print(f"[EQ JSON] Loading from {args.eq_json} ...", flush=True) with open(args.eq_json, 'r') as f: raw_events = json.load(f) for r in raw_events: ev = usgs.EarthquakeEvent( event_id = r.get('event_id', ''), place = r.get('place', ''), longitude = r.get('longitude'), latitude = r.get('latitude'), depth_km = r.get('depth_km'), magnitude = r.get('magnitude'), magnitude_type = r.get('magnitude_type', ''), ) if r.get('time_utc'): try: ev.time_utc = datetime.fromisoformat(r['time_utc']) except ValueError: pass # Reconstruct nodal planes if present (needed for L2) ev.nodal_plane_1 = usgs.NodalPlane( strike = r.get('np1_strike_deg'), dip = r.get('np1_dip_deg'), rake = r.get('np1_rake_deg'), ) ev.nodal_plane_2 = usgs.NodalPlane( strike = r.get('np2_strike_deg'), dip = r.get('np2_dip_deg'), rake = r.get('np2_rake_deg'), ) events.append(ev) print(f"[EQ JSON] {len(events)} event(s) loaded.", flush=True) return events # --------------------------------------------------------------------------- # Main pipeline # --------------------------------------------------------------------------- def main(): parser = build_parser() args = parser.parse_args() # ------------------------------------------------------------------ # 1. Parse GipsyX stations # ------------------------------------------------------------------ if not os.path.isfile(args.gipsy_file): sys.exit(f"ERROR: GipsyX file not found: {args.gipsy_file}") print(f"[GipsyX] Parsing stations from {args.gipsy_file} ...", flush=True) try: stations = parse_gipsy_stations(args.gipsy_file) except NotImplementedError as e: sys.exit(f"\nERROR: {e}") if not stations: sys.exit("ERROR: parse_gipsy_stations() returned no stations. " "Check your GipsyX file path and parsing logic.") print(f"[GipsyX] {len(stations)} station(s) loaded.", flush=True) # Quick sanity check on parsed coordinates for stn in stations[:3]: print(f" Sample: {stn['name']} lat={stn['lat']:.4f} " f"lon={stn['lon']:.4f} h={stn['height']:.1f} m", flush=True) # ------------------------------------------------------------------ # 2. Load / fetch earthquakes # ------------------------------------------------------------------ events = load_events(args) if not events: sys.exit("No earthquake events to process. Exiting.") # ------------------------------------------------------------------ # 3. Score each event against the GipsyX station list # ------------------------------------------------------------------ include_post = not args.no_postseismic results = [] for i, ev in enumerate(events, 1): print( f"[Scoring] {i}/{len(events)} {ev.event_id} " f"M{ev.magnitude} {ev.place} [{args.method}]", flush=True ) if args.method == 'L1': result = flag_stations_l1( ev, stations, include_postseismic=include_post, ) else: # L2 result = flag_stations_l2( ev, stations, include_postseismic=include_post, density=args.density, ) results.append(result) if args.verbose: print(f" → Method : {result.get('score_method', args.method)}") print(f" → Coseismic : {result['n_coseismic']} station(s)") print(f" → Postseismic : {result['n_postseismic']} station(s)") if result['coseismic_stations']: names = [s['name'] for s in result['coseismic_stations'][:10]] print(f" → Coseismic (first 10): {', '.join(names)}") # ------------------------------------------------------------------ # 4. Write output # ------------------------------------------------------------------ with open(args.output, 'w') as f: json.dump(results, f, indent=2) print(f"\n[Output] Written to {args.output} ({len(results)} event(s))") # --------------------------------------------------------------------------- # CLI argument parser # --------------------------------------------------------------------------- def build_parser() -> argparse.ArgumentParser: p = argparse.ArgumentParser( prog='gipsy_eq_station_flag.py', description=prolog, epilog=epilog, formatter_class=argparse.RawDescriptionHelpFormatter, ) # --- GipsyX station file (always required) --- p.add_argument( '--gipsy-file', metavar='PATH', required=True, help='Path to the GipsyX output file containing station locations.' ) # --- Earthquake source (pick one) --- src = p.add_mutually_exclusive_group(required=True) src.add_argument( '--event-id', metavar='USGS_ID', help='Fetch a single event by USGS event ID (e.g. us7000abcd).' ) src.add_argument( '--date', metavar='YYYY-MM-DD or YYYY-MM-DDTHH:MM:SS', help=( 'Start of the USGS query window. Accepts a bare date ("2024-01-01") ' 'or a full ISO-8601 timestamp ("2024-01-01T06:00:00").' ) ) src.add_argument( '--eq-json', metavar='PATH', help='Load pre-fetched earthquake_events.json instead of querying USGS.' ) # --- Time window (used with --date) --- p.add_argument( '--end-date', metavar='YYYY-MM-DD or YYYY-MM-DDTHH:MM:SS', default=None, help=( 'End of the USGS query window. If omitted, defaults to end-of-day ' 'on --date (single day).' ) ) # --- Query filters --- p.add_argument('--min-mag', type=float, default=5.0, help='Minimum magnitude (default: 5.0; use >=5.5 with L2).') p.add_argument('--min-lat', type=float, default=None, help='Bounding box south edge (decimal degrees).') p.add_argument('--max-lat', type=float, default=None, help='Bounding box north edge (decimal degrees).') p.add_argument('--min-lon', type=float, default=None, help='Bounding box west edge (decimal degrees).') p.add_argument('--max-lon', type=float, default=None, help='Bounding box east edge (decimal degrees).') p.add_argument('--limit', type=int, default=100, help='Max events from USGS (default: 100, hard cap: 20000).') # --- Scoring method --- p.add_argument( '--method', choices=['L1', 'L2'], default='L1', help=( 'Scoring method: L1 = isotropic distance-only (fast, no moment-tensor ' 'needed); L2 = Okada-based anisotropic mask (slower, requires nodal ' 'planes, falls back to L1 if unavailable). Default: L1.' ) ) # --- L2-specific --- p.add_argument( '--density', type=int, default=250, help='Okada grid resolution for L2 (default: 250). Higher = more accurate, slower.' ) # --- Output --- p.add_argument( '--output', metavar='PATH', default='gipsy_flagged_stations.json', help='Output JSON file path (default: gipsy_flagged_stations.json).' ) p.add_argument('--no-postseismic', action='store_true', help='Omit postseismic station flagging from output.') p.add_argument('--verbose', '-v', action='store_true', help='Print per-event summaries to stdout.') return p if __name__ == '__main__': main()