import requests from dataclasses import dataclass, field from datetime import datetime, timezone from typing import Optional import json import math # Author: Hasani Spann # Date: June 8, 2025 # Organization: JPL, Caltech prolog=""" **PROGRAM** usgs_earthquake_query.py **PURPOSE** | Queries the USGS FDSNWS Earthquake Catalog API for events over a specified | day-long window and extracts: | - Hypocentral location (longitude, latitude, depth) | - Earthquake magnitude (and magnitude type) | - Principal stress axes from the moment-tensor product: | * P-axis (max compression): azimuth, plunge, length | * T-axis (max extension): azimuth, plunge, length | * N-axis (null/neutral): azimuth, plunge, length | - Nodal plane geometry from the moment-tensor product: | * Nodal plane 1: strike, dip, rake | * Nodal plane 2: strike, dip, rake | - Derived fault dimensions (Wells & Coppersmith 1994 scaling relations): | * Along-strike length (km) | * Down-dip length (km) | * Average slip magnitude (m) | * Shallowest non-outcropping fault depth (km) | | The nodal plane and derived fault dimension fields support the Okada | (1985) elastic half-space model used in Gomez et al. (2024) GPS Solutions | to predict coseismic displacements at GNSS stations. """ epilog=""" **EXAMPLE** EqDataExtractor.py event = fetch_event_by_id("us7000abcd") **COPYRIGHT** | Copyright 2025, by the California Institute of Technology | United States Government Sponsorship acknowledged | All rights reserved **AUTHORS** | Hasani Spann | Jet Propulsion Laboratory | California Institute of Technology | Pasadena, CA, USA Data sources used to build this file: - USGS GeoJSON Summary format: https://earthquake.usgs.gov/earthquakes/feed/v1.0/geojson.php - USGS GeoJSON Detail format: https://earthquake.usgs.gov/earthquakes/feed/v1.0/geojson_detail.php - USGS PDL moment-tensor product property keys: https://ghsc.code-pages.usgs.gov/hazdev/pdl/userguide/products/known-types/moment-tensor.html - USGS FDSNWS API endpoint documentation: https://earthquake.usgs.gov/fdsnws/event/1/ NOTE on GipsyX integration: The EarthquakeEvent dataclass at the bottom is intentionally kept flat so downstream GipsyX helper functions can consume it directly. You will need to verify the exact field names/units required by GipsyX against its own documentation -- this file only handles the USGS side. **CHANGELOG** 2025: Patched two issues found during Part 4 overhead-script integration (gps_day_eq_pipeline.py): 1. Removed unused module-level imports (sys, pandas, geopandas, matplotlib.pyplot, functools.reduce, argparse, io.StringIO, numpy). None of these were referenced anywhere in this file. geopandas in particular is a heavy, sometimes fragile dependency on Windows, and its absence previously blocked importing this module (and anything that imports it) even though it was never used. 2. query_earthquakes() now checks for a None return from fetch_event_list() (which happens on any failed request -- network error, blocked domain, non-2xx HTTP status, timeout) and returns an empty list instead of crashing with "TypeError: 'NoneType' object is not iterable". 3. Added flush=True to the per-event "Parsing event i/N" progress print (plus a few adjacent status prints). A global M4+ query window normally returns 50-70+ events, and with fetch_details=True each one triggers a SECOND sequential HTTP request for moment-tensor data. Without flush=True, progress output could appear to stall on some terminals (observed on Windows PowerShell) even though the script was actively working through the event list -- this looks identical to a genuine hang. The fix makes each line appear the moment that event's processing starts. If the pipeline still feels slow, see gps_day_eq_pipeline.py's --no-details flag, which skips the per-event detail fetch entirely (L1-only scoring) for a much faster pass. 4. Replaced every non-ASCII character (->, --, deg, N-m in place of the original arrow/em-dash/degree/middot glyphs) in both live print statements and the PrincipalAxis/NodalPlane __repr__ methods. Windows consoles frequently default to a legacy codepage (cp1252 or similar) rather than UTF-8 for stdout. Printing a character outside that codepage raises UnicodeEncodeError ("'charmap' codec can't encode character...") at the exact print() call, which can terminate the process with little or no visible traceback in some terminal configurations (observed: VS Code's integrated terminal with PowerShell) -- this presents as the script silently doing nothing and returning to the prompt, which is easily mistaken for "the request hung" or "nothing happened" when the real cause is an encoding crash on the very first non-ASCII print. Verified by reproducing the crash under a forced cp1252 stdout encoder and confirming it no longer occurs after this fix. Dependencies: pip install requests Keywords: usgs, earthquake, seismicity, moment tensor, focal mechanism, p-axis, t-axis, n-axis, fdsnws, geojson """ # --------------------------------------------------------------------------- # Constants # --------------------------------------------------------------------------- USGS_QUERY_URL = "https://earthquake.usgs.gov/fdsnws/event/1/query" # Default request timeout in seconds REQUEST_TIMEOUT = 30 # --------------------------------------------------------------------------- # Data structures # --------------------------------------------------------------------------- @dataclass class PrincipalAxis: """ Represents one principal stress axis (P, T, or N). Azimuth: degrees clockwise from North (0-360) Plunge: degrees below horizontal (0-90) Length: scalar moment contribution in Newton-meters """ azimuth: Optional[float] = None plunge: Optional[float] = None length: Optional[float] = None def is_available(self) -> bool: return all(v is not None for v in (self.azimuth, self.plunge, self.length)) def __repr__(self): if self.is_available(): return (f"PrincipalAxis(az={self.azimuth:.1f}deg, " f"pl={self.plunge:.1f}deg, len={self.length:.3e} N-m)") return "PrincipalAxis(unavailable)" @dataclass class NodalPlane: """ Represents one nodal plane of a focal mechanism solution. Strike: degrees clockwise from North that the fault plane strikes (0-360) Dip: degrees the fault plane dips below horizontal (0-90) Rake: degrees describing the slip direction on the fault plane (-180 to 180) 0 = left-lateral strike-slip 90 = reverse/thrust -90 = normal 180 = right-lateral strike-slip These three parameters are the direct inputs to the Okada (1985) elastic half-space model used to compute surface displacement fields. Property keys sourced from the verified USGS PDL documentation at: https://ghsc.code-pages.usgs.gov/hazdev/pdl/userguide/products/known-types/moment-tensor.html """ strike: Optional[float] = None # degrees, clockwise from North dip: Optional[float] = None # degrees below horizontal rake: Optional[float] = None # degrees, slip direction def is_available(self) -> bool: return all(v is not None for v in (self.strike, self.dip, self.rake)) def __repr__(self): if self.is_available(): return (f"NodalPlane(strike={self.strike:.1f}deg, " f"dip={self.dip:.1f}deg, rake={self.rake:.1f}deg)") return "NodalPlane(unavailable)" @dataclass class EarthquakeEvent: """ Flat container for the earthquake parameters we extract from USGS. Designed to be passed into downstream GipsyX helper functions. Coordinate convention follows USGS ComCat: longitude: degrees East (-180 to 180) latitude: degrees North (-90 to 90) depth: kilometers below surface (positive down) """ # Event identity event_id: str = "" event_url: str = "" place: str = "" time_utc: Optional[datetime] = None # Hypocentral location longitude: Optional[float] = None # degrees East latitude: Optional[float] = None # degrees North depth_km: Optional[float] = None # km, positive downward # Magnitude magnitude: Optional[float] = None magnitude_type: str = "" # Principal axes from moment-tensor product # P-axis = axis of maximum compression p_axis: PrincipalAxis = field(default_factory=PrincipalAxis) # T-axis = axis of maximum extension t_axis: PrincipalAxis = field(default_factory=PrincipalAxis) # N-axis = null (intermediate) axis n_axis: PrincipalAxis = field(default_factory=PrincipalAxis) # Note: The field(default_factory=...) is needed here instead # of a plain default value because p_axis is a mutable object. If # you wrote p_axis: PrincipalAxis = PrincipalAxis() directly, all # instances of EarthquakeEvent would share the same PrincipalAxis object, # which would cause bugs. default_factory tells Python to call # PrincipalAxis() fresh for each new instance. # Moment tensor source metadata (useful for provenance) moment_tensor_source: str = "" # e.g. "us", "gcmt" moment_tensor_type: str = "" # e.g. "Mww", "Mwc", "Mwb" # Nodal planes from moment-tensor product (both planes stored; Okada model # is applied to each and results are combined via logical OR per Gomez 2024) # Property keys: nodal-plane-1-strike, nodal-plane-1-dip, nodal-plane-1-rake # Source: https://ghsc.code-pages.usgs.gov/hazdev/pdl/userguide/products/known-types/moment-tensor.html nodal_plane_1: NodalPlane = field(default_factory=NodalPlane) nodal_plane_2: NodalPlane = field(default_factory=NodalPlane) # ------------------------------------------------------------------------- # Derived fault dimensions -- computed from magnitude via Wells & # Coppersmith (1994) empirical scaling relations (Appendix A, Gomez 2024). # These are inputs to the Okada (1985) elastic half-space model. # Reference: Wells, D.L. & Coppersmith, K.J. (1994). New empirical # relationships among magnitude, rupture length, rupture width, rupture # area, and surface displacement. BSSA, 84(4), 974-1002. # ------------------------------------------------------------------------- fault_length_km: Optional[float] = None # along-strike length (km) fault_width_km: Optional[float] = None # down-dip length (km) fault_avg_slip_m: Optional[float] = None # average slip magnitude (m) fault_top_depth_km: Optional[float] = None # shallowest non-outcropping # depth of fault (km) def has_moment_tensor(self) -> bool: return self.p_axis.is_available() and self.t_axis.is_available() def has_nodal_planes(self) -> bool: return self.nodal_plane_1.is_available() and self.nodal_plane_2.is_available() def has_fault_dimensions(self) -> bool: return all(v is not None for v in ( self.fault_length_km, self.fault_width_km, self.fault_avg_slip_m, self.fault_top_depth_km )) def summary(self) -> str: lines = [ f"Event: {self.event_id} ({self.place})", f"Time UTC: {self.time_utc}", f"Location: lon={self.longitude}, lat={self.latitude}, depth={self.depth_km} km", f"Magnitude: {self.magnitude} {self.magnitude_type}", ] if self.has_moment_tensor(): lines += [ f"P-axis: {self.p_axis}", f"T-axis: {self.t_axis}", f"N-axis: {self.n_axis}", f"MT source: {self.moment_tensor_source} / {self.moment_tensor_type}", ] else: lines.append("P/T axes: not available for this event") if self.has_nodal_planes(): lines += [ f"NP1: {self.nodal_plane_1}", f"NP2: {self.nodal_plane_2}", ] else: lines.append("Nodal planes: not available for this event") if self.has_fault_dimensions(): lines += [ f"Fault length (along-strike): {self.fault_length_km:.2f} km", f"Fault width (down-dip): {self.fault_width_km:.2f} km", f"Avg slip: {self.fault_avg_slip_m:.3f} m", f"Shallowest fault depth: {self.fault_top_depth_km:.2f} km", ] return "\n".join(lines) # --------------------------------------------------------------------------- # API call helpers # --------------------------------------------------------------------------- def _safe_float(value) -> Optional[float]: """Convert a value to float, returning None if it cannot be converted.""" try: return float(value) except (TypeError, ValueError): return None def fetch_event_list( start_time: str, end_time: str, min_magnitude: float = 0.001, min_latitude: Optional[float] = None, max_latitude: Optional[float] = None, min_longitude: Optional[float] = None, max_longitude: Optional[float] = None, limit: int = 100, ) -> list[dict]: """ Step 1 of 2: Query the USGS FDSNWS API for a list of events. Parameters ---------- start_time : str ISO-8601 UTC start time, e.g. "2024-01-17T00:00:00" end_time : str ISO-8601 UTC end time, e.g. "2024-01-18T00:00:00" min_magnitude : float Minimum magnitude threshold (default 0.001) min/max_latitude, min/max_longitude : float, optional Bounding box filter. If omitted, global search is performed. limit : int Maximum number of events to return (USGS hard cap is 20000) Returns ------- list[dict] Raw list of GeoJSON Feature dicts from the USGS summary response. API endpoint documented at: https://earthquake.usgs.gov/fdsnws/event/1/ """ params = { "format": "geojson", "starttime": start_time, "endtime": end_time, "minmagnitude": min_magnitude, "limit": limit, "orderby": "time", # most recent first } # Optional geographic bounding box if min_latitude is not None: params["minlatitude"] = min_latitude if max_latitude is not None: params["maxlatitude"] = max_latitude if min_longitude is not None: params["minlongitude"] = min_longitude if max_longitude is not None: params["maxlongitude"] = max_longitude print(f"[USGS] Querying event list: {start_time} -> {end_time}, " f"minMag={min_magnitude} ...", flush=True) try: response = requests.get(USGS_QUERY_URL, params=params, timeout=REQUEST_TIMEOUT) response.raise_for_status() # raises HTTPError for 4xx/5xx except requests.RequestException as exc: print(f"[ERROR] USGS query failed: {exc}") return None data = response.json() features = data.get("features", []) print(f"[USGS] Found {len(features)} event(s).", flush=True) return features def fetch_event_detail(detail_url: str) -> Optional[dict]: """ Step 2 of 2: Fetch the GeoJSON Detail response for a single event. The detail URL is provided by the summary response in feature["properties"]["detail"]. Returns the parsed JSON dict, or None on failure. Detail format documented at: https://earthquake.usgs.gov/earthquakes/feed/v1.0/geojson_detail.php """ try: response = requests.get(detail_url, timeout=REQUEST_TIMEOUT) response.raise_for_status() return response.json() except requests.RequestException as exc: print(f" [WARNING] Could not fetch detail from {detail_url}: {exc}", flush=True) return None def _extract_axes_from_moment_tensor(mt_product: dict) -> tuple[PrincipalAxis, PrincipalAxis, PrincipalAxis]: """ Parse P, T, and N axes out of a single moment-tensor product dict. The moment-tensor product properties are documented at: https://ghsc.code-pages.usgs.gov/hazdev/pdl/userguide/products/known-types/moment-tensor.html Keys used: p-axis-azimuth, p-axis-plunge, p-axis-length (max compression) t-axis-azimuth, t-axis-plunge, t-axis-length (max extension) n-axis-azimuth, n-axis-plunge, n-axis-length (null axis) All angles in degrees; lengths in Newton-meters. """ props = mt_product.get("properties", {}) p_axis = PrincipalAxis( azimuth = _safe_float(props.get("p-axis-azimuth")), plunge = _safe_float(props.get("p-axis-plunge")), length = _safe_float(props.get("p-axis-length")), ) t_axis = PrincipalAxis( azimuth = _safe_float(props.get("t-axis-azimuth")), plunge = _safe_float(props.get("t-axis-plunge")), length = _safe_float(props.get("t-axis-length")), ) n_axis = PrincipalAxis( azimuth = _safe_float(props.get("n-axis-azimuth")), plunge = _safe_float(props.get("n-axis-plunge")), length = _safe_float(props.get("n-axis-length")), ) return p_axis, t_axis, n_axis def _extract_nodal_planes(mt_product: dict) -> tuple[NodalPlane, NodalPlane]: """ Parse both nodal planes out of a single moment-tensor product dict. Both planes are needed because the Gomez (2024) method applies the Okada model to each plane independently and combines results via logical OR, producing a more conservative displacement mask. Property keys sourced from the verified USGS PDL documentation at: https://ghsc.code-pages.usgs.gov/hazdev/pdl/userguide/products/known-types/moment-tensor.html Keys used: nodal-plane-1-strike, nodal-plane-1-dip, nodal-plane-1-rake nodal-plane-2-strike, nodal-plane-2-dip, nodal-plane-2-rake Note: The older PDL 2.9.1 documentation uses "slip" instead of "rake" for this key name. The current version uses "rake". If values come back as None, inspect the raw API response to confirm the key name in use. """ props = mt_product.get("properties", {}) np1 = NodalPlane( strike = _safe_float(props.get("nodal-plane-1-strike")), dip = _safe_float(props.get("nodal-plane-1-dip")), rake = _safe_float(props.get("nodal-plane-1-rake")), ) np2 = NodalPlane( strike = _safe_float(props.get("nodal-plane-2-strike")), dip = _safe_float(props.get("nodal-plane-2-dip")), rake = _safe_float(props.get("nodal-plane-2-rake")), ) return np1, np2 def _compute_fault_dimensions( magnitude: float, depth_km: float, dip_deg: float, ) -> tuple[Optional[float], Optional[float], Optional[float], Optional[float]]: """ Compute empirical fault dimensions from moment magnitude using the Wells & Coppersmith (1994) scaling relations, as applied in Appendix A and B of Gomez et al. (2024), GPS Solutions 28:214. These are required inputs to the Okada (1985) elastic half-space model for computing surface displacement fields. Scaling relations (Appendix A, Gomez 2024): Along-strike length: log10(L) = -3.22 + 0.69 * Mw Down-dip width: log10(W) = -1.01 + 0.32 * Mw Average slip: log10(s) = -4.80 + 0.69 * Mw Shallowest non-outcropping depth (Appendix B, Gomez 2024): dmin = depth - (W/2 * sin(dip)) If dmin < 0: fault_top_depth = depth - dmin (shift down to avoid outcrop) Else: fault_top_depth = dmin Reference: 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. Parameters ---------- magnitude : float Moment magnitude (Mw). depth_km : float Hypocentral depth in km (positive downward). dip_deg : float Fault dip angle in degrees (from nodal plane 1 or 2). Returns ------- fault_length_km, fault_width_km, fault_avg_slip_m, fault_top_depth_km All float values, or None if inputs are invalid. """ try: mw = float(magnitude) d = float(depth_km) dip = float(dip_deg) except (TypeError, ValueError): return None, None, None, None # Wells & Coppersmith (1994) scaling relations fault_length_km = 10 ** (-3.22 + 0.69 * mw) # along-strike length, km fault_width_km = 10 ** (-1.01 + 0.32 * mw) # down-dip width, km fault_avg_slip_m = 10 ** (-4.80 + 0.69 * mw) # average slip, m # Shallowest non-outcropping depth (Appendix B, Gomez 2024) dip_rad = math.radians(dip) d_min = d - (fault_width_km / 2.0) * math.sin(dip_rad) if d_min < 0: # Fault would outcrop at this depth -- shift it down fault_top_depth_km = d - d_min else: fault_top_depth_km = d_min return fault_length_km, fault_width_km, fault_avg_slip_m, fault_top_depth_km def parse_event(summary_feature: dict, fetch_details: bool = True) -> EarthquakeEvent: """ Convert a raw GeoJSON Feature (from the summary list) into an EarthquakeEvent dataclass, optionally fetching the detail endpoint to populate the moment-tensor / pnt axis fields. Parameters ---------- summary_feature : dict A single Feature dict from the USGS summary GeoJSON response. fetch_details : bool If True (default), follow the detail URL and extract moment-tensor axes. Set to False for a quick pass that only gives location/mag. """ props = summary_feature.get("properties", {}) geometry = summary_feature.get("geometry", {}) coords = geometry.get("coordinates", [None, None, None]) # --- Convert epoch milliseconds to UTC datetime --- time_ms = props.get("time") time_utc = ( datetime.fromtimestamp(time_ms / 1000.0, tz=timezone.utc) if time_ms is not None else None ) event = EarthquakeEvent( event_id = summary_feature.get("id", ""), event_url = props.get("url", ""), place = props.get("place", ""), time_utc = time_utc, longitude = _safe_float(coords[0]), latitude = _safe_float(coords[1]), depth_km = _safe_float(coords[2]), magnitude = _safe_float(props.get("mag")), magnitude_type = props.get("magType", ""), ) if not fetch_details: return event # --- Follow the detail URL to get moment-tensor data --- detail_url = props.get("detail") if not detail_url: return event detail = fetch_event_detail(detail_url) if detail is None: return event # The "products" block is inside detail["properties"]["products"] products = detail.get("properties", {}).get("products", {}) # Prefer "moment-tensor" product; fall back to "focal-mechanism" if absent. # USGS recommends moment-tensor when available (more data, more accurate). # See: https://geohazards.usgs.gov/pipermail/realtime-feed-users/2019-July/000095.html mt_list = products.get("moment-tensor") or products.get("focal-mechanism") if not mt_list: # This event has no focal mechanism data (common for small earthquakes) return event # Use the first (highest-weight / preferred) moment tensor solution mt_product = mt_list[0] mt_props = mt_product.get("properties", {}) p_axis, t_axis, n_axis = _extract_axes_from_moment_tensor(mt_product) np1, np2 = _extract_nodal_planes(mt_product) event.p_axis = p_axis event.t_axis = t_axis event.n_axis = n_axis event.moment_tensor_source = mt_product.get("source", "") event.moment_tensor_type = mt_props.get("derived-magnitude-type", "") event.nodal_plane_1 = np1 event.nodal_plane_2 = np2 # Compute derived fault dimensions using nodal plane 1 dip as primary. # We use NP1 dip preferentially; if unavailable we fall back to NP2. # The calling code (Gomez L2 method) should run Okada for both planes # independently using nodal_plane_1 and nodal_plane_2 directly. dip_for_dims = ( np1.dip if np1.dip is not None else np2.dip if np2.dip is not None else None ) if event.magnitude is not None and event.depth_km is not None and dip_for_dims is not None: (event.fault_length_km, event.fault_width_km, event.fault_avg_slip_m, event.fault_top_depth_km) = _compute_fault_dimensions( event.magnitude, event.depth_km, dip_for_dims ) return event # --------------------------------------------------------------------------- # High-level query function # --------------------------------------------------------------------------- def _normalise_time(time_str: str, end_of_day: bool = False) -> str: """ Accept either a bare date ("YYYY-MM-DD") or a full ISO-8601 timestamp ("YYYY-MM-DDTHH:MM:SS") and always return a full ISO-8601 string. Parameters ---------- time_str : str Date or datetime string. end_of_day : bool When True and a bare date is given, append "T23:59:59" instead of "T00:00:00". Ignored when a full timestamp is supplied. Returns ------- str ISO-8601 datetime string, e.g. "2024-01-17T00:00:00". """ time_str = time_str.strip() if "T" in time_str: # Already a full timestamp -- pass through unchanged return time_str # Bare date -- attach a time component suffix = "T23:59:59" if end_of_day else "T00:00:00" return time_str + suffix def query_earthquakes( start_time: str, end_time: Optional[str] = None, min_magnitude: float = 0.001, min_latitude: Optional[float] = None, max_latitude: Optional[float] = None, min_longitude: Optional[float] = None, max_longitude: Optional[float] = None, fetch_details: bool = True, limit: int = 100, ) -> list[EarthquakeEvent]: """ Main entry point: query USGS for earthquakes over an arbitrary UTC time window and return a list of EarthquakeEvent objects. Parameters ---------- start_time : str Start of the query window. Accepts either: - A bare date: "YYYY-MM-DD" -> treated as 00:00:00 UTC - A full ISO-8601 time: "YYYY-MM-DDTHH:MM:SS" Examples: "2024-01-17" -> 2024-01-17 00:00:00 UTC "2024-01-17T06:00:00" -> 2024-01-17 06:00:00 UTC end_time : str, optional End of the query window. Accepts the same formats as start_time. When a bare date is given it is treated as 23:59:59 UTC on that day. If omitted, defaults to end-of-day on the same date as start_time, preserving the original single-day behaviour. Examples: None -> end of the start_time date (single day) "2024-01-17" -> 2024-01-17 23:59:59 UTC "2024-01-20T12:00:00" -> 2024-01-20 12:00:00 UTC min_magnitude : float Minimum magnitude threshold (default 0.001). Moment-tensor data is rarely available for events below ~M4.5, so pnt axes will mostly be empty for smaller events. min/max_latitude, min/max_longitude : float, optional Optional geographic bounding box. fetch_details : bool Whether to follow detail URLs to fetch moment-tensor / nodal plane data (default True). Set False for location + magnitude only. limit : int Maximum events to return (USGS hard cap is 20000). Returns ------- list[EarthquakeEvent] Examples -------- # Single day (original behaviour preserved) events = query_earthquakes("2024-01-17") # Explicit multi-day range events = query_earthquakes("2024-01-17", "2024-01-20") # Hour-precise window events = query_earthquakes("2024-02-06T00:00:00", "2024-02-06T12:00:00") # Multi-month range with bounding box events = query_earthquakes( start_time = "2023-01-01", end_time = "2023-12-31", min_magnitude = 6.0, min_latitude = 30.0, max_latitude = 45.0, min_longitude = 25.0, max_longitude = 45.0, ) """ # Normalise start time -- always treat bare date as start of day t_start = _normalise_time(start_time, end_of_day=False) # Normalise end time -- default to end of the start date if not given if end_time is None: # Extract just the date portion of the (possibly already full) start t_end = _normalise_time(t_start[:10], end_of_day=True) else: t_end = _normalise_time(end_time, end_of_day=True) print(f"[query_earthquakes] Window: {t_start} -> {t_end}", flush=True) features = fetch_event_list( start_time = t_start, end_time = t_end, min_magnitude = min_magnitude, min_latitude = min_latitude, max_latitude = max_latitude, min_longitude = min_longitude, max_longitude = max_longitude, limit = limit, ) if features is None: # fetch_event_list() returns None (not []) on a failed request -- # network error, blocked domain, non-2xx HTTP status, timeout, etc. # Without this check, the loop below raises # "TypeError: 'NoneType' object is not iterable" instead of failing # gracefully. Returning [] lets callers (e.g. the Part 4 overhead # script) treat "no events" and "request failed" the same way. print("[query_earthquakes] No events returned (request failed) -- " "returning an empty list.") return [] events = [] for i, feature in enumerate(features): event_id = feature.get("id", f"event_{i}") print(f" Parsing event {i+1}/{len(features)}: {event_id}", flush=True) event = parse_event(feature, fetch_details=fetch_details) events.append(event) return events # Query to extract a single event by its USGS ComCat event ID, e.g. "us7000abcd". def fetch_event_by_id(event_id: str) -> Optional[EarthquakeEvent]: """ Fetch a single event directly by its USGS ComCat event ID. This bypasses the summary list and goes straight to the detail endpoint, so moment-tensor data is always included if available. The eventid parameter is documented at: https://earthquake.usgs.gov/fdsnws/event/1/ Parameters ---------- event_id : str USGS ComCat event ID, e.g. "us7000abcd". Visible in the URL of any USGS event page: https://earthquake.usgs.gov/earthquakes/eventpage/us7000abcd Example ------- event = fetch_event_by_id("us7000abcd") """ params = { "format": "geojson", "eventid": event_id, } print(f"[USGS] Fetching event by ID: {event_id}") response = requests.get(USGS_QUERY_URL, params=params, timeout=REQUEST_TIMEOUT) response.raise_for_status() detail = response.json() # The eventid query returns a detail-format response directly, # so we can parse axes from it immediately without a second call. # We construct a minimal summary-style feature to reuse parse_event(), # but pass fetch_details=False since we already have the full response. props = detail.get("properties", {}) geometry = detail.get("geometry", {}) coords = geometry.get("coordinates", [None, None, None]) time_ms = props.get("time") time_utc = ( datetime.fromtimestamp(time_ms / 1000.0, tz=timezone.utc) if time_ms is not None else None ) event = EarthquakeEvent( event_id = detail.get("id", event_id), event_url = props.get("url", ""), place = props.get("place", ""), time_utc = time_utc, longitude = _safe_float(coords[0]), latitude = _safe_float(coords[1]), depth_km = _safe_float(coords[2]), magnitude = _safe_float(props.get("mag")), magnitude_type = props.get("magType", ""), ) # Extract moment tensor from the products block in this same response products = props.get("products", {}) mt_list = products.get("moment-tensor") or products.get("focal-mechanism") if mt_list: mt_product = mt_list[0] mt_props = mt_product.get("properties", {}) p_axis, t_axis, n_axis = _extract_axes_from_moment_tensor(mt_product) np1, np2 = _extract_nodal_planes(mt_product) event.p_axis = p_axis event.t_axis = t_axis event.n_axis = n_axis event.moment_tensor_source = mt_product.get("source", "") event.moment_tensor_type = mt_props.get("derived-magnitude-type", "") event.nodal_plane_1 = np1 event.nodal_plane_2 = np2 dip_for_dims = ( np1.dip if np1.dip is not None else np2.dip if np2.dip is not None else None ) if event.magnitude is not None and event.depth_km is not None and dip_for_dims is not None: (event.fault_length_km, event.fault_width_km, event.fault_avg_slip_m, event.fault_top_depth_km) = _compute_fault_dimensions( event.magnitude, event.depth_km, dip_for_dims ) return event # --------------------------------------------------------------------------- # Optional: export to JSON (useful for GipsyX pipeline hand-off) # --------------------------------------------------------------------------- def events_to_dict(events: list[EarthquakeEvent]) -> list[dict]: """ Serialize a list of EarthquakeEvent objects to a list of plain dicts. Suitable for json.dump() or passing to a GipsyX helper function. """ results = [] for ev in events: results.append({ # Identity "event_id": ev.event_id, "place": ev.place, "time_utc": ev.time_utc.isoformat() if ev.time_utc else None, # Hypocentral location "longitude": ev.longitude, "latitude": ev.latitude, "depth_km": ev.depth_km, # Magnitude "magnitude": ev.magnitude, "magnitude_type": ev.magnitude_type, # Principal axes "p_axis_azimuth_deg": ev.p_axis.azimuth, "p_axis_plunge_deg": ev.p_axis.plunge, "p_axis_length_Nm": ev.p_axis.length, "t_axis_azimuth_deg": ev.t_axis.azimuth, "t_axis_plunge_deg": ev.t_axis.plunge, "t_axis_length_Nm": ev.t_axis.length, "n_axis_azimuth_deg": ev.n_axis.azimuth, "n_axis_plunge_deg": ev.n_axis.plunge, "n_axis_length_Nm": ev.n_axis.length, # Moment tensor provenance "moment_tensor_source": ev.moment_tensor_source, "moment_tensor_type": ev.moment_tensor_type, # Nodal planes (Okada model inputs) "np1_strike_deg": ev.nodal_plane_1.strike, "np1_dip_deg": ev.nodal_plane_1.dip, "np1_rake_deg": ev.nodal_plane_1.rake, "np2_strike_deg": ev.nodal_plane_2.strike, "np2_dip_deg": ev.nodal_plane_2.dip, "np2_rake_deg": ev.nodal_plane_2.rake, # Derived fault dimensions (Wells & Coppersmith 1994) "fault_length_km": ev.fault_length_km, "fault_width_km": ev.fault_width_km, "fault_avg_slip_m": ev.fault_avg_slip_m, "fault_top_depth_km": ev.fault_top_depth_km, }) return results # --------------------------------------------------------------------------- # Example usage / quick test # --------------------------------------------------------------------------- if __name__ == "__main__": # ----------------------------------------------------------------------- # Configuration -- edit these to suit your query # ----------------------------------------------------------------------- # # Time window -- three supported formats: # # Single day (original behaviour): # START_TIME = "2024-01-01" # END_TIME = None # defaults to end of START_TIME date # # Multi-day range: # START_TIME = "2024-01-01" # END_TIME = "2024-01-07" # # Hour-precise window: # START_TIME = "2024-02-06T00:00:00" # END_TIME = "2024-02-06T12:00:00" # START_TIME = "2024-01-01" END_TIME = None # set to a date/datetime to span multiple days MIN_MAGNITUDE = 5.0 # Geographic bounding box (optional -- set to None for global) MIN_LAT, MAX_LAT = None, None MIN_LON, MAX_LON = None, None OUTPUT_JSON = "earthquake_events.json" # set to None to skip file output # ----------------------------------------------------------------------- # Run query # ----------------------------------------------------------------------- events = query_earthquakes( start_time = START_TIME, end_time = END_TIME, min_magnitude = MIN_MAGNITUDE, min_latitude = MIN_LAT, max_latitude = MAX_LAT, min_longitude = MIN_LON, max_longitude = MAX_LON, fetch_details = True, # set False for location/mag only (faster) limit = 50, ) single_event = fetch_event_by_id("us6000m0xj") events.append(single_event) print(f"\nFetched single event by ID: {single_event.event_id} - {single_event.place}") # ----------------------------------------------------------------------- # Print summaries # ----------------------------------------------------------------------- print("\n" + "="*60) print(f"Results: {len(events)} earthquake(s)") print("="*60) for ev in events: print() print(ev.summary()) # ----------------------------------------------------------------------- # Optionally write to JSON for downstream GipsyX functions # ----------------------------------------------------------------------- if OUTPUT_JSON: data = events_to_dict(events) with open(OUTPUT_JSON, "w") as f: json.dump(data, f, indent=2) print(f"\n[Output] Written to {OUTPUT_JSON}")