prolog = """ **PROGRAM** eq_mask_l1.py **PURPOSE** | For one or more earthquake events, computes the Level-1 (isotropic, | circular) coseismic and postseismic influence-zone masks using the | S-score boundary formulation from Gómez et al. (2024, GPS Solutions 28:214). | | The L1 mask boundary is the circle at which the coseismic S-score equals | zero. Starting from the L1 S-score formula: | | sc = A * Mw - log10(dist_km) + B | | Setting sc = 0 and solving for dist_km gives the coseismic circle radius: | | r_coseismic = 10^(A * Mw + B) [km] | r_postseismic = 10^(A * Mw + B) * PSF [km] | | where A = 0.5261, B = -1.1478, PSF = 1.5 (POST_SEISMIC_SCALE_FACTOR) | are the empirically fitted coefficients from Gómez et al. (2024). | | For each event this program: | 1. Computes the coseismic and postseismic circle radii | 2. Generates each circle as a list of geographic (lon, lat) vertices | using the spherical destination-point formula | 3. Saves one KMZ file per event (always; default dir = current directory) | 4. Optionally saves mask parameters to a JSON file (--output) | | Primary design: the core functions (compute_l1_mask, save_mask_kmz) are | intended to be imported and called from a parent pipeline script. | The CLI is a convenience wrapper around those same functions. **USAGE — Python module (primary design intent)** import usgs_earthquake_query as usgs import eq_mask_l1 as ml # Compute masks from a list of EarthquakeEvent objects — returns plain dicts events = usgs.query_earthquakes('2024-01-17', min_magnitude=5.0) masks = [ml.compute_l1_mask(ev) for ev in events] # KMZ is always written when calling save_mask_kmz directly import os os.makedirs('./masks', exist_ok=True) for m in masks: ml.save_mask_kmz(m, f'./masks/{m["event_id"]}_L1_mask.kmz') # Optionally convert to JSON yourself import json with open('my_masks.json', 'w') as f: json.dump(masks, f, indent=2) # Or load from a pre-fetched JSON file events = ml._load_events_from_json('earthquake_events.json') masks = [ml.compute_l1_mask(ev) for ev in events] **USAGE — CLI** # KMZ files are written to the current directory by default python eq_mask_l1.py --event-id us7000abcd # Fetch all M>=5 events on a calendar date python eq_mask_l1.py --date 2024-01-01 --min-mag 5.0 # Multi-day window; save KMZ to a specific directory python eq_mask_l1.py --date 2024-01-01 --end-date 2024-01-07 \\ --min-mag 5.0 --kmz-dir ./masks/ # From a pre-fetched earthquake_events.json python eq_mask_l1.py --eq-json earthquake_events.json # Suppress KMZ output (JSON only) python eq_mask_l1.py --event-id us7000abcd --no-kmz # Suppress JSON output (KMZ only) python eq_mask_l1.py --event-id us7000abcd --no-json # Coseismic mask only (no postseismic ring) python eq_mask_l1.py --event-id us7000abcd --no-postseismic """ epilog = """ **OUTPUTS** JSON file (--output, default: eq_masks_L1.json) with one entry per earthquake. Each entry contains: event_id, place, time_utc, magnitude, magnitude_type, longitude, latitude, depth_km — event identity and location mask_type — 'L1_isotropic_circle' A, B, PSF — the empirical coefficients used coseismic_radius_km — r_c = 10^(A*Mw + B) postseismic_radius_km — r_p = r_c * PSF n_polygon_points — number of vertices per circle coseismic_polygon — list of [lon, lat] pairs tracing the coseismic circle postseismic_polygon — same for the postseismic circle (empty if --no-postseismic) If --kmz-dir is supplied, one KMZ per event is written there, named: _L1_mask.kmz Each KMZ contains: - A star icon at the epicenter with an HTML description balloon - A blue polygon for the coseismic circle - An orange polygon for the postseismic circle (omitted with --no-postseismic) **METHODOLOGY** 1. S-score boundary derivation The L1 S-score (Gómez et al. 2024, Eq. 1) is: sc = A * Mw - log10(dist_km) + B Setting sc = 0 and solving algebraically: log10(r_c) = A * Mw + B => r_c = 10^(A * Mw + B) [km] This is the exact formula used in geode/pyOkada.py line 304: self.dmax = 10. ** (a * self.mag + b) The postseismic radius follows the same log-scale shift: r_p = r_c * POST_SEISMIC_SCALE_FACTOR 2. Circle polygon generation For each bearing angle θ ∈ [0°, 360°) the circle vertex is the destination point at distance r from the epicenter, computed using the standard spherical destination-point formula (derived from the spherical law of cosines — see any spherical trigonometry reference): φ2 = arcsin( sin(φ1)*cos(δ) + cos(φ1)*sin(δ)*cos(θ) ) λ2 = λ1 + arctan2( sin(θ)*sin(δ)*cos(φ1), cos(δ) - sin(φ1)*sin(φ2) ) where δ = r / R_earth (angular distance, radians). This formula is exact for a spherical earth (R = 6371 km), which is consistent with the Haversine distance used throughout this codebase. 3. KMZ export Style and simplekml usage adapted from: geode/pyOkada.py Score.save_masks() by Demian Gomez https://github.com/demiangomez/geode (BSD-3-Clause) **REFERENCES** S-score formula, coefficients, and POST_SEISMIC_SCALE_FACTOR: Gómez, D.D. et al. (2024). A method to identify and model co-seismic and post-seismic displacements in a GNSS network. GPS Solutions 28:214. https://doi.org/10.1007/s10291-024-01758-9 S-score implementation cross-referenced from: geode/pyOkada.py by Demian Gomez https://github.com/demiangomez/geode **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 simplekml (optional; required only for KMZ export) Install with: pip install simplekml requests (via usgs_earthquake_query) usgs_earthquake_query.py (must be in the same directory or on PYTHONPATH) Keywords: gnss, gipsy, earthquake, coseismic, postseismic, mask, l1, isotropic, circle, kmz """ import argparse import json import os import sys import numpy as np from datetime import datetime, timezone # --------------------------------------------------------------------------- # Optional KMZ dependency — simplekml # --------------------------------------------------------------------------- try: import simplekml HAS_SIMPLEKML = True except ImportError: HAS_SIMPLEKML = False # --------------------------------------------------------------------------- # Import the USGS query module (must be co-located or on PYTHONPATH) # --------------------------------------------------------------------------- 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." ) from eq_station_score import ( A, B, POST_SEISMIC_SCALE_FACTOR, destination_point, circle_polygon ) # --------------------------------------------------------------------------- # JSON loader # Reads a file previously written by usgs_earthquake_query.events_to_dict() # and reconstructs EarthquakeEvent objects. The flat field names follow # the events_to_dict() schema in usgs_earthquake_query.py lines 834-871. # --------------------------------------------------------------------------- def _load_events_from_json(path: str) -> list: """ Load a list of EarthquakeEvent objects from a JSON file previously written by usgs_earthquake_query.events_to_dict() / json.dump(). The JSON must be a list of flat dicts with keys matching the events_to_dict() output schema (see usgs_earthquake_query.py). Fields not present default to None / empty string. Parameters ---------- path : str Path to the JSON file (e.g. earthquake_events.json). Returns ------- list[usgs.EarthquakeEvent] """ with open(path, 'r') as f: records = json.load(f) events = [] for rec in records: # Parse the ISO-8601 time string back to a timezone-aware datetime time_utc = None t_str = rec.get('time_utc') if t_str: try: time_utc = datetime.fromisoformat(t_str) if time_utc.tzinfo is None: time_utc = time_utc.replace(tzinfo=timezone.utc) except ValueError: pass ev = usgs.EarthquakeEvent( event_id = rec.get('event_id', ''), place = rec.get('place', ''), time_utc = time_utc, longitude = rec.get('longitude'), latitude = rec.get('latitude'), depth_km = rec.get('depth_km'), magnitude = rec.get('magnitude'), magnitude_type = rec.get('magnitude_type', ''), ) # Principal axes ev.p_axis = usgs.PrincipalAxis( azimuth = rec.get('p_axis_azimuth_deg'), plunge = rec.get('p_axis_plunge_deg'), length = rec.get('p_axis_length_Nm'), ) ev.t_axis = usgs.PrincipalAxis( azimuth = rec.get('t_axis_azimuth_deg'), plunge = rec.get('t_axis_plunge_deg'), length = rec.get('t_axis_length_Nm'), ) ev.n_axis = usgs.PrincipalAxis( azimuth = rec.get('n_axis_azimuth_deg'), plunge = rec.get('n_axis_plunge_deg'), length = rec.get('n_axis_length_Nm'), ) # Nodal planes ev.nodal_plane_1 = usgs.NodalPlane( strike = rec.get('np1_strike_deg'), dip = rec.get('np1_dip_deg'), rake = rec.get('np1_rake_deg'), ) ev.nodal_plane_2 = usgs.NodalPlane( strike = rec.get('np2_strike_deg'), dip = rec.get('np2_dip_deg'), rake = rec.get('np2_rake_deg'), ) # Moment tensor provenance ev.moment_tensor_source = rec.get('moment_tensor_source', '') ev.moment_tensor_type = rec.get('moment_tensor_type', '') # Derived fault dimensions ev.fault_length_km = rec.get('fault_length_km') ev.fault_width_km = rec.get('fault_width_km') ev.fault_avg_slip_m = rec.get('fault_avg_slip_m') ev.fault_top_depth_km = rec.get('fault_top_depth_km') events.append(ev) return events # --------------------------------------------------------------------------- # L1 radius formulas # Both are algebraic rearrangements of the L1 S-score formula. # Cross-referenced with geode/pyOkada.py Score.__init__, lines 304 and 340. # --------------------------------------------------------------------------- def l1_coseismic_radius(magnitude: float) -> float: """ Compute the L1 coseismic influence-zone radius in km. Derived by setting sc = 0 in the L1 S-score formula: sc = A * Mw - log10(dist_km) + B (Gómez et al. 2024, Eq. 1) 0 = A * Mw - log10(r_c) + B => r_c = 10^(A * Mw + B) [km] This is 'dmax' in geode/pyOkada.py (line 304): self.dmax = 10. ** (a * self.mag + b) Parameters ---------- magnitude : float Moment magnitude Mw. Returns ------- float Coseismic circle radius in km. """ return 10.0 ** (A * magnitude + B) def l1_postseismic_radius(magnitude: float) -> float: """ Compute the L1 postseismic influence-zone radius in km. Derived by setting sp = 0 in the postseismic S-score formula: sp = A * Mw - log10(dist_km) + B + log10(PSF) => r_p = 10^(A * Mw + B + log10(PSF)) = r_c * PSF Cross-referenced with geode/pyOkada.py Score.__init__, line 340: self.p_mask = ... < (POST_SEISMIC_SCALE_FACTOR * (self.dmax * 1000.)) Parameters ---------- magnitude : float Moment magnitude Mw. Returns ------- float Postseismic circle radius in km. """ return l1_coseismic_radius(magnitude) * POST_SEISMIC_SCALE_FACTOR # --------------------------------------------------------------------------- # Core mask computation # --------------------------------------------------------------------------- def compute_l1_mask( event: usgs.EarthquakeEvent, n_points: int = 360, include_postseismic: bool = True, ) -> dict: """ Compute the L1 (circular) coseismic and postseismic mask for one event. Parameters ---------- event : usgs.EarthquakeEvent Earthquake data object from usgs_earthquake_query. n_points : int Number of vertices (per circle) in the polygon representation. include_postseismic : bool If True, compute and include the postseismic polygon. Returns ------- dict Mask parameters suitable for json.dump(), with keys: event_id, place, time_utc, magnitude, magnitude_type, longitude, latitude, depth_km, mask_type, A, B, PSF, coseismic_radius_km, postseismic_radius_km, n_polygon_points, coseismic_polygon, postseismic_polygon On missing data, returns a minimal dict with an 'error' key instead. """ if event.latitude is None or event.longitude is None or event.magnitude is None: return { 'event_id' : event.event_id, 'error' : 'Missing latitude, longitude, or magnitude — cannot compute mask.', } r_c = l1_coseismic_radius(event.magnitude) r_p = l1_postseismic_radius(event.magnitude) c_polygon = circle_polygon(event.longitude, event.latitude, r_c, n_points) p_polygon = ( circle_polygon(event.longitude, event.latitude, r_p, n_points) if include_postseismic else [] ) return { # --- Event identity --- 'event_id' : event.event_id, 'place' : event.place, 'time_utc' : event.time_utc.isoformat() if event.time_utc else None, 'magnitude' : event.magnitude, 'magnitude_type' : event.magnitude_type, 'longitude' : event.longitude, 'latitude' : event.latitude, 'depth_km' : event.depth_km, # --- Mask parameters --- 'mask_type' : 'L1_isotropic_circle', 'A' : A, 'B' : B, 'PSF' : POST_SEISMIC_SCALE_FACTOR, 'coseismic_radius_km' : round(r_c, 4), 'postseismic_radius_km': round(r_p, 4) if include_postseismic else None, 'n_polygon_points' : n_points, # --- Polygon vertices --- 'coseismic_polygon' : c_polygon, 'postseismic_polygon' : p_polygon, } # --------------------------------------------------------------------------- # KMZ export # --------------------------------------------------------------------------- def save_mask_kmz( mask: dict, kmz_path: str, include_postseismic: bool = True, ) -> None: """ Save one L1 circular mask as a KMZ file viewable in Google Earth. The KMZ contains: - A starred yellow point at the epicenter with an HTML description balloon showing event details and a link to the USGS event page. - A solid blue outline polygon for the coseismic circle. - A solid orange outline polygon for the postseismic circle (omitted if include_postseismic=False). Style and simplekml usage adapted from geode/pyOkada.py Score.save_masks() by Demian Gomez (https://github.com/demiangomez/geode, BSD-3-Clause). Parameters ---------- mask : dict Output dict from compute_l1_mask(). kmz_path : str Destination file path (should end in .kmz). include_postseismic : bool If True and postseismic polygon data is present, include it. """ if not HAS_SIMPLEKML: print( "[KMZ] WARNING: simplekml is not installed; skipping KMZ export.\n" " Install with: pip install simplekml", flush=True, ) return if mask.get('error'): print( f"[KMZ] Skipping {mask.get('event_id', '?')} — no valid mask data.", flush=True, ) return kml = simplekml.Kml() # ------------------------------------------------------------------ # # Epicenter marker # ------------------------------------------------------------------ # place = mask.get('place', mask['event_id']) mag = mask.get('magnitude', 'N/A') depth = mask.get('depth_km', 'N/A') t_utc = mask.get('time_utc', 'N/A') eid = mask.get('event_id', 'N/A') r_c_km = mask.get('coseismic_radius_km', 'N/A') r_p_km = mask.get('postseismic_radius_km', 'N/A') epicenter = kml.newpoint( name = place, coords = [(mask['longitude'], mask['latitude'])], ) 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 epicenter.description = ( "" f"" f"" f"" f"" f"" f"" f"" f"" "
ID:{eid}
Date (UTC):{t_utc}
Magnitude:{mag}
Depth:{depth} km
Method:L1 isotropic circle
Coseismic radius:{r_c_km} km
Postseismic radius:{r_p_km} km
USGS event page" f"
" ) # ------------------------------------------------------------------ # # Helper: add an outline-only polygon ring to the KML # ------------------------------------------------------------------ # def _add_ring(polygon_coords: list, name: str, color) -> None: if not polygon_coords: return # simplekml expects (lon, lat) tuples; our polygon is [[lon, lat], ...] coords = [(pt[0], pt[1]) for pt in polygon_coords] poly = kml.newpolygon(name=name, outerboundaryis=coords) poly.style.linestyle.color = color poly.style.linestyle.width = 3 # Transparent fill — outline only poly.style.polystyle.color = simplekml.Color.changealphaint( 0, simplekml.Color.white ) # ------------------------------------------------------------------ # # Coseismic circle — blue # ------------------------------------------------------------------ # _add_ring( mask.get('coseismic_polygon', []), 'Coseismic mask (L1)', simplekml.Color.blue, ) # ------------------------------------------------------------------ # # Postseismic circle — orange # ------------------------------------------------------------------ # if include_postseismic and mask.get('postseismic_polygon'): _add_ring( mask['postseismic_polygon'], 'Postseismic mask (L1)', simplekml.Color.orange, ) kml.savekmz(kmz_path) print(f"[KMZ] Saved: {kmz_path}", flush=True) # --------------------------------------------------------------------------- # CLI # --------------------------------------------------------------------------- def build_parser() -> argparse.ArgumentParser: p = argparse.ArgumentParser( prog='eq_mask_l1.py', description=prolog, epilog=epilog, formatter_class=argparse.RawDescriptionHelpFormatter, ) # --- Earthquake source (mutually exclusive) --- 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", ' 'treated as 00:00:00 UTC) or a full ISO-8601 timestamp. ' 'Combine with --end-date for multi-day or hour-precise windows.' ) ) src.add_argument( '--eq-json', metavar='PATH', help='Load pre-fetched earthquake_events.json instead of querying USGS.' ) # --- Query filters (only used when --date is active) --- p.add_argument( '--end-date', metavar='YYYY-MM-DD or YYYY-MM-DDTHH:MM:SS', default=None, help=( 'End of the USGS query window. Accepts the same formats as --date. ' 'A bare date is treated as 23:59:59 UTC. ' 'If omitted, defaults to end-of-day on the --date value.' ) ) p.add_argument( '--min-mag', type=float, default=5.0, help='Minimum magnitude filter for USGS queries (default: 5.0).' ) p.add_argument('--min-lat', type=float, default=None, help='Southern boundary of geographic bounding box (decimal degrees).') p.add_argument('--max-lat', type=float, default=None, help='Northern boundary of geographic bounding box (decimal degrees).') p.add_argument('--min-lon', type=float, default=None, help='Western boundary of geographic bounding box (decimal degrees).') p.add_argument('--max-lon', type=float, default=None, help='Eastern boundary of geographic bounding box (decimal degrees).') p.add_argument( '--limit', type=int, default=100, help='Max number of events to fetch from USGS (default: 100, hard cap: 20000).' ) # --- Mask options --- p.add_argument( '--n-points', type=int, default=360, help=( 'Number of vertices in each circle polygon (default: 360, i.e. one ' 'vertex per degree of bearing). Increase for smoother export at the ' 'cost of larger file size.' ) ) p.add_argument( '--no-postseismic', action='store_true', help='Compute and save only the coseismic mask; omit the postseismic circle.' ) # --- Output --- p.add_argument( '--output', metavar='PATH', default='eq_masks_L1.json', help='Output JSON file path (default: eq_masks_L1.json). Ignored if --no-json.' ) p.add_argument( '--kmz-dir', metavar='DIR', default='.', help=( 'Directory in which to save one KMZ file per event ' '(default: current directory "."). ' 'The directory is created if it does not exist. ' 'To suppress KMZ output entirely, use --no-kmz.' ) ) p.add_argument( '--no-kmz', action='store_true', help='Suppress KMZ output entirely. Overrides --kmz-dir.' ) p.add_argument( '--no-json', action='store_true', help='Suppress JSON output. Useful when only KMZ files are needed.' ) p.add_argument( '--verbose', '-v', action='store_true', help='Print a per-event summary table to stdout.' ) return p # --------------------------------------------------------------------------- # Pretty printer # --------------------------------------------------------------------------- def _print_mask_result(mask: dict) -> None: """Print a formatted summary of one compute_l1_mask() result.""" sep = '─' * 60 print(sep) if mask.get('error'): print(f" Event : {mask.get('event_id', '?')} ERROR: {mask['error']}") print(sep) return print(f" Event : {mask['event_id']} M{mask['magnitude']} {mask['place']}") print(f" Time (UTC) : {mask.get('time_utc', 'N/A')}") print(f" Location : {mask['latitude']:.4f}°, {mask['longitude']:.4f}° " f"depth {mask.get('depth_km')} km") print(f" Mask type : {mask['mask_type']}") print(f" Cos. radius : {mask['coseismic_radius_km']:.2f} km") if mask.get('postseismic_radius_km') is not None: print(f" Post. radius : {mask['postseismic_radius_km']:.2f} km") print(f" Polygon pts : {mask['n_polygon_points']} vertices per circle (closed)") print(sep) # --------------------------------------------------------------------------- # Main # --------------------------------------------------------------------------- def main(): parser = build_parser() args = parser.parse_args() include_post = not args.no_postseismic # ------------------------------------------------------------------ # # 1. Load earthquake 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] if ev is not None else [] elif args.eq_json: if not os.path.isfile(args.eq_json): sys.exit(f"ERROR: file not found: {args.eq_json}") print(f"[JSON] Loading events from {args.eq_json} ...", flush=True) events = _load_events_from_json(args.eq_json) print(f"[JSON] {len(events)} event(s) loaded.", flush=True) else: # Date-range query — normalisation handled inside query_earthquakes() print( f"[USGS] Querying {args.date} → " f"{args.end_date or '(end of day)'} " f"(minMag={args.min_mag}) ...", 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, limit = args.limit, ) print(f"[USGS] {len(events)} event(s) returned.", flush=True) if not events: print("[WARNING] No events to process. Exiting.") return [] # ------------------------------------------------------------------ # # 2. Prepare KMZ output directory # ------------------------------------------------------------------ # write_kmz = not args.no_kmz if write_kmz: os.makedirs(args.kmz_dir, exist_ok=True) if not HAS_SIMPLEKML: print( "[WARNING] simplekml is not installed; KMZ files will not be written.\n" " Install with: pip install simplekml\n" " Use --no-kmz to silence this warning.", flush=True, ) write_kmz = False # can't write without the library # ------------------------------------------------------------------ # # 3. Compute masks # ------------------------------------------------------------------ # masks = [] for i, ev in enumerate(events, 1): print( f"[Mask] {i}/{len(events)} {ev.event_id} " f"M{ev.magnitude} {ev.place}", flush=True, ) mask = compute_l1_mask(ev, n_points=args.n_points, include_postseismic=include_post) masks.append(mask) if args.verbose: _print_mask_result(mask) # Write one KMZ per event unless suppressed if write_kmz and not mask.get('error'): kmz_name = f"{ev.event_id}_L1_mask.kmz" kmz_path = os.path.join(args.kmz_dir, kmz_name) save_mask_kmz(mask, kmz_path, include_postseismic=include_post) # ------------------------------------------------------------------ # # 4. Write JSON output (optional) # ------------------------------------------------------------------ # if not args.no_json: with open(args.output, 'w') as f: json.dump(masks, f, indent=2) print(f"\n[Output] {len(masks)} mask(s) written to {args.output}") return masks if __name__ == '__main__': main()