prolog = """ **PROGRAM** eq_station_score.py **PURPOSE** | Pure S-score scoring engine for GNSS station vs. earthquake influence-zone | testing. This module implements Part 3 of the project outline: | | "Uses the s-score method to list what stations are within the | coseismic and post-seismic radii of an earthquake (or list of | earthquakes)." | | It contains ONLY the scoring math (Level-1 isotropic + Level-2 Okada), | factored out of eq_station_flag_L1.py and eq_station_flag_L2.py so that | it can be called directly by the Part 4 overhead/orchestrator program | without dragging along CLI parsing, station-file I/O, or USGS fetching. | | This module deliberately does NOT contain: | - argparse / CLI wrapper code (stays in the flag_L*.py scripts | or moves to the Part 4 driver) | - load_stations() file I/O (Part 1/4 concern) | - USGS querying (usgs_earthquake_query.py) | - mask/KMZ generation (Part 2 concern, eq_mask_l1 / | future eq_mask_l2) | | --- Level 1 (isotropic) --- | sc = a * Mw - log10(dist_km) + b (> 0 -> coseismic) | sp = a * Mw - log10(dist_km) + b + log10(k) (> 0 -> postseismic) | where a = 0.5261, b = -1.1478, k = POST_SEISMIC_SCALE_FACTOR = 1.5 | | --- Level 2 (Okada, anisotropic) --- | Builds a boolean displacement mask (|U| >= 1 mm) on a 2D grid centred | on the epicenter using the Okada (1985) elastic half-space model, once | per nodal plane, at both surface and hypocentral depth. The mask is | rescaled so its maximum flagged radius equals dmax = 10^(a*Mw + b) km | (the L1 radius), then queried via KDTree for each station's | azimuthal-equidistant projected position. | | If no nodal-plane data are available for an event, score_stations() | automatically falls back to the L1 isotropic circle, mirroring | geode/pyOkada.py Score.__init__ behaviour. **USAGE** # As a library — typical call from the Part 4 overhead/orchestrator script from eq_station_score import score_stations result = score_stations(event, stations, method='L2') # result['coseismic_stations'], result['postseismic_stations'] # Force L1 isotropic only (fast, no moment tensor needed) result = score_stations(event, stations, method='L1') # Auto: use L2 if nodal planes are present on the event, else L1 result = score_stations(event, stations, method='auto') # Score a single station/event pair directly (no station-list looping) sc, sp = l1_score(magnitude=6.4, dist_km=42.3) eq_score = Score(event_lat=-31.2, event_lon=-69.0, depth_km=10.0, magnitude=6.4, strikes=[12.0], dips=[55.0], rakes=[90.0]) sc, sp = eq_score.score(station_lat=-31.5, station_lon=-69.3) """ epilog = """ **OUTPUTS** score_stations() returns a dict: score_method : 'L1_isotropic' | 'L2_okada' | 'L1_isotropic_fallback' nodal_planes_used : bool strikes, dips, rakes : nodal plane values used (empty list if L1/fallback) n_coseismic : int n_postseismic : int coseismic_stations : list[dict] name, lat, lon, country, city, distance_km, azimuth_deg, sc, sp postseismic_stations : list[dict] same structure Station entries are sorted by distance_km, closest first. No event-identity metadata (event_id, place, time_utc, ...) is attached here — that wrapping is the responsibility of the calling/overhead program (Part 4), which already has the EarthquakeEvent in hand. **METHODOLOGY** S-score formula, coefficients (a, b), and POST_SEISMIC_SCALE_FACTOR from: 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 Score class logic, Okada implementation, fault-dimension scaling, and KDTree-based point-in-mask query adapted from: geode/pyOkada.py by Demian Gomez https://github.com/demiangomez/geode Okada (1985) elastic half-space model: 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 & Coppersmith (1994) empirical scaling relations: 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 **COPYRIGHT** | Copyright 2025, by the California Institute of Technology | United States Government Sponsorship acknowledged | All rights reserved Dependencies: numpy scipy (scipy.spatial.KDTree — only required for method='L2'/'auto') Keywords: gnss, gipsy, earthquake, coseismic, postseismic, s-score, l1, l2, okada, station scoring """ import numpy as np from scipy.spatial import KDTree # --------------------------------------------------------------------------- # Constants — source: Gomez et al. (2024), GPS Solutions 28:214, Eq. (1)-(3) # Coefficients fitted to 809 global seismic events (see paper Section 2). # Identical to the values in eq_station_flag_L1.py and eq_station_flag_L2.py. # --------------------------------------------------------------------------- A = 0.5261 # magnitude scaling coefficient B = -1.1478 # distance offset POST_SEISMIC_SCALE_FACTOR = 1.5 # postseismic zone = coseismic zone x 1.5 # Displacement threshold below which a point is considered unaffected (metres) # Source: geode/pyOkada.py, Score.compute_disp_field(), line 399 DISPLACEMENT_THRESHOLD_M = 1e-3 # 1 mm # Grid density (points per axis of the square mesh centred on the epicenter). # geode default is 250. Higher values improve mask resolution but cost more # compute time (cost scales as density^2 per nodal plane per depth). GRID_DENSITY = 250 # Far-field scale factor for the grid extent. # Source: geode/pyOkada.py Score.__init__, line 330 FAR_FIELD_SCALE = 25 # --------------------------------------------------------------------------- # Degree-trig helpers # Source: geode/pyOkada.py, lines 74-78 # --------------------------------------------------------------------------- cosd = lambda x: np.cos(np.deg2rad(x)) sind = lambda x: np.sin(np.deg2rad(x)) acosd = lambda x: np.rad2deg(np.arccos(x)) asind = lambda x: np.rad2deg(np.arcsin(x)) atand = lambda x: np.rad2deg(np.arctan(x)) # --------------------------------------------------------------------------- # Geometry helpers # Adapted verbatim from geode/pyOkada.py by Demian Gomez # (https://github.com/demiangomez/geode), BSD 3-Clause License # --------------------------------------------------------------------------- def distance_km(lon1: float, lat1: float, lon2: float, lat2: float) -> float: """ Haversine great-circle distance between two points on Earth. Returns distance in kilometres. Source: geode/pyOkada.py, lines 110-127. """ lon1 = lon1 * np.pi / 180 lat1 = lat1 * np.pi / 180 lon2 = lon2 * np.pi / 180 lat2 = lat2 * np.pi / 180 dlon = lon2 - lon1 dlat = lat2 - lat1 d = np.sin(dlat / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2 c = 2 * np.arcsin(np.sqrt(d)) return float((6371.0 * c).item()) def azimuth_deg(lon1: float, lat1: float, lon2: float, lat2: float) -> float: """ Initial bearing (azimuth) in degrees [0, 360) from point 1 to point 2. Source: geode/pyOkada.py, lines 87-108. """ lon1 = lon1 * np.pi / 180 lat1 = lat1 * np.pi / 180 lon2 = lon2 * np.pi / 180 lat2 = lat2 * np.pi / 180 delta_lon = lon2 - lon1 y = np.sin(delta_lon) * np.cos(lat2) x = (np.cos(lat1) * np.sin(lat2) - np.sin(lat1) * np.cos(lat2) * np.cos(delta_lon)) az = np.rad2deg(np.arctan2(y, x)) return float(((az + 360.0) % 360.0).item()) # --------------------------------------------------------------------------- # Level 1 — isotropic S-score # --------------------------------------------------------------------------- def l1_score(magnitude: float, dist_km: float) -> tuple[float, float]: """ Compute the Level-1 (isotropic) coseismic and postseismic S-scores. Returns (sc, sp) where: sc > 0 -> station is within the coseismic influence zone sp > 0 -> station is within the postseismic influence zone Formula source: Gomez et al. (2024), Eq. (1)-(3) and surrounding text. Also visible in geode/pyOkada.py, EarthquakeTable.__init__, lines 170-171. Parameters ---------- magnitude : float Moment magnitude (Mw) of the earthquake. dist_km : float Great-circle distance from station to epicenter, in km. Notes ----- dist_km must be > 0. If the station is at the epicenter (dist_km == 0), log10 is undefined; we return +inf to ensure it is always flagged. """ if dist_km <= 0.0: return (float('inf'), float('inf')) sc = A * magnitude - np.log10(dist_km) + B sp = A * magnitude - np.log10(dist_km) + B + np.log10(POST_SEISMIC_SCALE_FACTOR) return float(np.float64(sc).item()), float(np.float64(sp).item()) # --------------------------------------------------------------------------- # Level 2 — Okada elastic half-space model # Translated from MATLAB by Demian Gomez, based on Okada (1985) BSSA. # Source: geode/pyOkada.py, lines 641-766. # --------------------------------------------------------------------------- def okadakernel(alp, xi, et, q, sd, cd, disl1, disl2, disl3): """ Indefinite integral kernel for the Okada (1985) surface displacement model. Source: geode/pyOkada.py, lines 696-766. """ pi2 = 2 * np.pi f2 = 2 * np.ones_like(xi) xi2 = xi ** 2 et2 = et ** 2 q2 = q ** 2 r2 = xi2 + et2 + q2 r = np.sqrt(r2) d = et * sd - q * cd y = et * cd + q * sd ret = r + et ret[ret < 0] = 0 rd = r + d tt = np.arctan(xi * et / (q * r)) re = np.zeros_like(ret) re[ret != 0] = 1 / ret[ret != 0] dle = np.zeros_like(ret) dle[ret == 0] = -np.log(r[ret == 0] - et[ret == 0]) dle[ret != 0] = np.log(ret[ret != 0]) rrx = 1 / (r * (r + xi)) rre = re / r u1 = np.zeros_like(xi) u2 = np.zeros_like(xi) u3 = np.zeros_like(xi) if abs(cd) > 1e-6: # non-vertical fault td = sd / cd x = np.sqrt(xi2 + q2) a5 = np.zeros_like(xi) a5[xi != 0] = ( alp * f2[xi != 0] / cd * np.arctan( (et[xi != 0] * (x[xi != 0] + q[xi != 0] * cd) + x[xi != 0] * (r[xi != 0] + x[xi != 0]) * sd) / (xi[xi != 0] * (r[xi != 0] + x[xi != 0]) * cd) ) ) a4 = alp / cd * (np.log(rd) - sd * dle) a3 = alp * (y / rd / cd - dle) + td * a4 a1 = -alp / cd * xi / rd - td * a5 else: # vertical fault a1 = -alp / f2 * xi * q / rd ** 2 a3 = alp / f2 * (et / rd + y * q / rd ** 2 - dle) if disl1 != 0: a4 = -alp * q / rd a5 = -alp * xi * sd / rd if disl1 != 0: a2 = -alp * dle - a3 if disl1 != 0: un = disl1 / pi2 req = rre * q u1 -= un * (req * xi + tt + a1 * sd) u2 -= un * (req * y + q * cd * re + a2 * sd) u3 -= un * (req * d + q * sd * re + a4 * sd) if disl2 != 0: un = disl2 / pi2 sdcd = sd * cd u1 -= un * (q / r - a3 * sdcd) u2 -= un * (y * q * rrx + cd * tt - a1 * sdcd) u3 -= un * (d * q * rrx + sd * tt - a5 * sdcd) if disl3 != 0: un = disl3 / pi2 sdsd = sd * sd u1 += un * (q2 * rre - a3 * sdsd) u2 += un * (-d * q * rrx - sd * (xi * q * rre - tt) - a1 * sdsd) u3 += un * (y * q * rrx + cd * (xi * q * rre - tt) - a5 * sdsd) return u1, u2, u3 def okada(alpha, x, y, d, L1, L2, W1, W2, snd, csd, B1, B2, B3): """ Surface displacements (ux, uy, uz) due to a rectangular dislocation in an elastic half-space (Okada 1985). Parameters ---------- alpha : float Elastic constant: mu/(lambda+mu). 0.5 for Poisson solid. x, y : ndarray Station coordinates in Okada's local system (metres). d : float Depth of the dislocation plane intersection with Z-axis (m). L1,L2 : float Along-strike extent of dislocation (m). W1,W2 : float Down-dip extent of dislocation (m). snd : float sin(dip angle) csd : float cos(dip angle) B1 : float Strike-slip Burgers vector component (m). B2 : float Dip-slip Burgers vector component (m). B3 : float Tensile Burgers vector component (m). Returns ------- ux, uy, uz : ndarray Surface displacement components (m). Source: geode/pyOkada.py, lines 641-693. """ if not isinstance(x, np.ndarray): x = np.array([x]) y = np.array([y]) p = y * csd + d * snd q = y * snd - d * csd et = p - W1; xi = x - L1 ux, uy, uz = okadakernel(alpha, xi, et, q, snd, csd, B1, B2, B3) xi = x - L2 u1, u2, u3 = okadakernel(alpha, xi, et, q, snd, csd, B1, B2, B3) ux -= u1; uy -= u2; uz -= u3 et = p - W2; xi = x - L1 u1, u2, u3 = okadakernel(alpha, xi, et, q, snd, csd, B1, B2, B3) ux -= u1; uy -= u2; uz -= u3 xi = x - L2 u1, u2, u3 = okadakernel(alpha, xi, et, q, snd, csd, B1, B2, B3) ux += u1; uy += u2; uz += u3 return ux, uy, uz # --------------------------------------------------------------------------- # Fault dimension scaling # Source: Wells & Coppersmith (1994); geode/pyOkada.py Score.fault_dims(), lines 351-358 # --------------------------------------------------------------------------- def fault_dims(magnitude: float) -> dict: """ Empirical fault dimension scaling from Wells & Coppersmith (1994), Appendix A, Tables 2A and 2B (all-slip-type regression). Returns dimensions in metres and m^2, matching geode/pyOkada.py units. """ along_strike_m = 10.0 ** (-3.22 + 0.69 * magnitude) * 1000.0 downdip_m = 10.0 ** (-1.01 + 0.32 * magnitude) * 1000.0 avg_disp_m = 10.0 ** (-4.80 + 0.69 * magnitude) rupture_area = 10.0 ** (-3.49 + 0.91 * magnitude) # km^2 maximum_disp_m = 10.0 ** (-5.46 + 0.82 * magnitude) return { 'along_strike_m': along_strike_m, 'downdip_m' : downdip_m, 'avg_disp_m' : avg_disp_m, 'rupture_area' : rupture_area, 'maximum_disp_m': maximum_disp_m, } # --------------------------------------------------------------------------- # L2 displacement mask computation # Logic adapted from geode/pyOkada.py Score.__init__ and Score.compute_disp_field() # --------------------------------------------------------------------------- def compute_l2_mask( magnitude: float, depth_km: float, strikes: list[float], dips: list[float], rakes: list[float], scale_factor: float = 1.0, density: int = GRID_DENSITY, threshold_m: float = DISPLACEMENT_THRESHOLD_M, ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """ Compute the boolean L2 displacement mask on a square grid centred at the epicenter in azimuthal-equidistant coordinates (km). The mask is True at grid points where predicted surface displacement magnitude |U| >= threshold_m for ANY nodal plane at ANY of the two reference depths (surface = 0 and hypocentral depth). The raw Okada grid is then scaled so that its maximum non-zero radius equals dmax * scale_factor, where dmax = 10^(A*Mw + B) km. Parameters ---------- magnitude : Mw depth_km : hypocentral depth in km strikes : list of nodal plane strike angles (degrees) dips : list of nodal plane dip angles (degrees) rakes : list of nodal plane rake angles (degrees) scale_factor: 1.0 for coseismic, POST_SEISMIC_SCALE_FACTOR for postseismic density : grid resolution (points per axis) threshold_m: displacement threshold in metres Returns ------- mx, my : ndarray (density x density) Scaled grid coordinates in km (azimuthal equidistant, origin = epicenter). mask : ndarray[bool] (density x density) True where predicted displacement >= threshold. Source logic: geode/pyOkada.py Score.__init__ lines 330-346, Score.compute_disp_field() lines 360-412. """ dims = fault_dims(magnitude) L1_val = -dims['along_strike_m'] / 2.0 L2_val = dims['along_strike_m'] / 2.0 W1_val = -dims['downdip_m'] / 2.0 W2_val = dims['downdip_m'] / 2.0 ad = dims['avg_disp_m'] # L1 dmax in km (used for grid extent and final scaling) dmax = 10.0 ** (A * magnitude + B) # km # Build grid in metres (centred on epicenter) xmax_m = np.ceil(dims['along_strike_m']) * FAR_FIELD_SCALE gx, gy = np.meshgrid( np.linspace(-xmax_m, xmax_m, density), np.linspace(-xmax_m, xmax_m, density) ) # Two reference depths: surface (0) and hypocentral depth (metres) depths_m = [0.0, depth_km * 1000.0] ref_scale = [] for depth_m in depths_m: U = np.zeros_like(gx, dtype=bool) for strike, dip, rake in zip(strikes, dips, rakes): # Ensure fault top doesn't breach the surface # Source: geode/pyOkada.py compute_disp_field() lines 377-381 W2_adj = W2_val depth_adj = depth_m d2 = depth_m - (W2_adj * sind(dip) + 500.0) if d2 < 0: depth_adj = depth_m - d2 # Rotate station coordinates into Okada's strike-aligned frame # Clockwise rotation by (90 - strike) degrees # Source: geode/pyOkada.py compute_disp_field() lines 383-391 angle = 90.0 - strike R = np.array([ [ cosd(angle), sind(angle)], [-sind(angle), cosd(angle)], ]) flat = R @ np.array([gx.flatten(), gy.flatten()]) tx, ty = flat[0, :], flat[1, :] n, e, u = okada( 0.5, # Poisson solid (lambda == mu) tx, ty, depth_adj, L1_val, L2_val, W1_val, W2_adj, sind(dip), cosd(dip), ad * cosd(rake), # B1: strike-slip component ad * sind(rake), # B2: dip-slip component 0.0, # B3: tensile component (zero for shear fault) ) n = np.reshape(n, gx.shape) e = np.reshape(e, gx.shape) u = np.reshape(u, gx.shape) U = np.logical_or( np.sqrt(np.square(n) + np.square(e) + np.square(u)) >= threshold_m, U ) # Record maximum distance of flagged points from epicenter try: ref_scale.append( np.max(np.sqrt(np.square(gx[U]) + np.square(gy[U]))) ) except ValueError: ref_scale.append(0.0) # no points above threshold # Scale the grid so that the maximum flagged radius = dmax * scale_factor # Convert gx/gy from metres to km for the output coordinate system max_ref = np.max(ref_scale) if np.max(ref_scale) > 0 else 1.0 mx = gx / max_ref * dmax * scale_factor # km my = gy / max_ref * dmax * scale_factor # km # Recompute mask at the final (surface) depth for the output U_final = np.zeros_like(gx, dtype=bool) depth_m = depth_km * 1000.0 for strike, dip, rake in zip(strikes, dips, rakes): d2 = depth_m - (W2_val * sind(dip) + 500.0) depth_adj = depth_m - d2 if d2 < 0 else depth_m angle = 90.0 - strike R = np.array([ [ cosd(angle), sind(angle)], [-sind(angle), cosd(angle)], ]) flat = R @ np.array([gx.flatten(), gy.flatten()]) tx, ty = flat[0, :], flat[1, :] n, e, u = okada( 0.5, tx, ty, depth_adj, L1_val, L2_val, W1_val, W2_val, sind(dip), cosd(dip), ad * cosd(rake), ad * sind(rake), 0.0 ) n = np.reshape(n, gx.shape) e = np.reshape(e, gx.shape) u = np.reshape(u, gx.shape) U_final = np.logical_or( np.sqrt(np.square(n) + np.square(e) + np.square(u)) >= threshold_m, U_final ) return mx, my, U_final # --------------------------------------------------------------------------- # Score class — standalone port of geode/pyOkada.py Score # All database and KML dependencies removed. # --------------------------------------------------------------------------- class Score: """ Computes the L2 (or L1 fallback) seismic score for a single earthquake. After construction, call score(lat, lon) to test whether a station falls within the coseismic or postseismic mask. If no nodal-plane data are provided (empty strikes/dips/rakes), the class falls back to the L1 isotropic circular mask, matching the behaviour of geode/pyOkada.py Score.__init__ lines 334-345. """ def __init__( self, event_lat : float, event_lon : float, depth_km : float, magnitude : float, strikes : list = (), dips : list = (), rakes : list = (), density : int = GRID_DENSITY, ): self.lat = float(event_lat) self.lon = float(event_lon) self.mag = float(magnitude) self.depth = [0.0, float(depth_km) * 1000.0] # [surface_m, hypocentral_m] self.strikes = list(strikes) self.dips = list(dips) self.rakes = list(rakes) # L1 dmax (km) — maximum influence radius self.dmax = 10.0 ** (A * self.mag + B) # Fault dimensions dims = fault_dims(magnitude) self.along_strike_m = dims['along_strike_m'] self.downdip_m = dims['downdip_m'] self.avg_disp_m = dims['avg_disp_m'] # Grid extent xmax_m = np.ceil(self.along_strike_m) * FAR_FIELD_SCALE gx_raw, gy_raw = np.meshgrid( np.linspace(-xmax_m, xmax_m, density), np.linspace(-xmax_m, xmax_m, density) ) self.has_nodal_planes = len(self.strikes) > 0 if self.has_nodal_planes: # L2 path: run Okada, build scaled mask self.c_mx, self.c_my, self.c_mask = compute_l2_mask( magnitude, depth_km, self.strikes, self.dips, self.rakes, scale_factor=1.0, density=density, ) self.p_mx, self.p_my, self.p_mask = compute_l2_mask( magnitude, depth_km, self.strikes, self.dips, self.rakes, scale_factor=POST_SEISMIC_SCALE_FACTOR, density=density, ) else: # L1 fallback path: circular mask in km # Grid is in km for direct radius comparison with dmax gx_km = gx_raw / 1000.0 gy_km = gy_raw / 1000.0 self.c_mask = np.sqrt(np.square(gx_km) + np.square(gy_km)) < self.dmax self.p_mask = (np.sqrt(np.square(gx_km) + np.square(gy_km)) < POST_SEISMIC_SCALE_FACTOR * self.dmax) self.c_mx = gx_km self.c_my = gy_km self.p_mx = gx_km self.p_my = gy_km # KDTree for fast point-in-mask lookup # Source: geode/pyOkada.py Score.__init__ lines 348-349 self.kd_c = KDTree(np.column_stack((self.c_mx.flatten(), self.c_my.flatten()))) self.kd_p = KDTree(np.column_stack((self.p_mx.flatten(), self.p_my.flatten()))) def score(self, station_lat: float, station_lon: float) -> tuple[float, float]: """ Return (sc, sp) for a station: sc > 0 -> within coseismic mask sp > 0 -> within postseismic mask Converts station lat/lon to the azimuthal equidistant coordinate system centred on the epicenter, then queries the KDTree. Source: geode/pyOkada.py Score.score(), lines 414-437. """ lat = np.atleast_1d(np.float64(station_lat)) lon = np.atleast_1d(np.float64(station_lon)) c = np.arccos( sind(self.lat) * sind(lat) + cosd(self.lat) * cosd(lat) * cosd(lon - self.lon) ) k = np.zeros_like(c) k[c != 0] = c[c != 0] / np.sin(c[c != 0]) * 6371.0 x = k * cosd(lat) * sind(lon - self.lon) y = k * (cosd(self.lat) * sind(lat) - sind(self.lat) * cosd(lat) * cosd(lon - self.lon)) qp = np.column_stack((x, y)) _, i = self.kd_c.query(qp) sc = float(self.c_mask.flatten()[i.item()]) + 0.0 _, i = self.kd_p.query(qp) sp = float(self.p_mask.flatten()[i.item()]) + 0.0 return sc, sp # --------------------------------------------------------------------------- # Unified scoring entry point # # This is the function the Part 4 overhead/orchestrator script should call. # It replaces flag_stations_l1() / flag_stations_l2() from the old combined # scripts, but returns ONLY the score results — no event-identity metadata, # since the caller already holds the EarthquakeEvent and can attach that # itself (or just pass the event through alongside this function's output). # --------------------------------------------------------------------------- def score_stations( event, stations: list[dict], method: str = 'auto', include_postseismic: bool = True, density: int = GRID_DENSITY, ) -> dict: """ Score every station in `stations` against a single earthquake `event` using the L1 isotropic S-score, the L2 Okada S-score, or 'auto' (use L2 if nodal planes are available on the event, else fall back to L1). Parameters ---------- event : usgs_earthquake_query.EarthquakeEvent Must have .latitude, .longitude, .magnitude set. For method='L2' or 'auto' with nodal planes, .nodal_plane_1 / .nodal_plane_2 and .depth_km are also used. stations : list[dict] Station dicts as produced by load_stations() in eq_station_flag_L1.py / eq_station_flag_L2.py (or any dict with at minimum 'name', 'lat', 'lon' keys — 'country'/'city' default to '' if absent). method : {'auto', 'L1', 'L2'} 'L1' — isotropic epicentral-distance score only. No moment tensor required. Fast. 'L2' — Okada anisotropic score. Requires nodal_plane_1/2 on the event; if absent, raises ValueError (use 'auto' to fall back silently instead). 'auto' — use L2 if the event has at least one usable nodal plane, otherwise fall back to L1 isotropic. This matches the behaviour already implemented in eq_station_flag_L2.py's Score class (geode/pyOkada.py Score.__init__ lines 334-345). include_postseismic : bool If True, stations with sp > 0 (but sc <= 0) are included in postseismic_stations. density : int Okada grid resolution per axis (only relevant for L2/auto-with-NPs). Returns ------- dict with keys: score_method, nodal_planes_used, strikes, dips, rakes, n_coseismic, n_postseismic, coseismic_stations, postseismic_stations Each station entry contains: name, lat, lon, country, city, distance_km, azimuth_deg, sc, sp Notes ----- No event-identity fields (event_id, place, time_utc, magnitude, ...) are included in the return value — attach those in the calling/overhead program, which already holds the EarthquakeEvent object. This keeps the scoring engine decoupled from any particular event-metadata schema. """ if method not in ('auto', 'L1', 'L2'): raise ValueError(f"method must be 'auto', 'L1', or 'L2' — got {method!r}") if event.latitude is None or event.longitude is None or event.magnitude is None: raise ValueError( f"Event {getattr(event, 'event_id', '?')} is missing latitude, " "longitude, or magnitude — cannot score." ) eq_lat = event.latitude eq_lon = event.longitude mag = event.magnitude # --- Extract nodal planes (needed for L2 / auto) --- strikes, dips, rakes = [], [], [] np1 = getattr(event, 'nodal_plane_1', None) np2 = getattr(event, 'nodal_plane_2', None) if np1 is not None and 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 is not None and 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 if method == 'L2' and not nodal_planes_used: raise ValueError( f"method='L2' requested but event {getattr(event, 'event_id', '?')} " "has no usable nodal-plane data. Use method='auto' to fall back to " "L1 isotropic automatically, or method='L1' to score isotropically " "on purpose." ) use_l2 = nodal_planes_used and method in ('L2', 'auto') if use_l2: score_method = 'L2_okada' depth_km = event.depth_km if event.depth_km is not None else 10.0 eq_score = Score( event_lat = eq_lat, event_lon = eq_lon, depth_km = depth_km, magnitude = mag, strikes = strikes, dips = dips, rakes = rakes, density = density, ) else: score_method = 'L1_isotropic' if method == 'L1' else 'L1_isotropic_fallback' strikes, dips, rakes = [], [], [] # not used in L1 path; keep output clean coseismic_stations = [] postseismic_stations = [] for stn in stations: dist = distance_km(eq_lon, eq_lat, stn['lon'], stn['lat']) az = azimuth_deg(eq_lon, eq_lat, stn['lon'], stn['lat']) if use_l2: # Cheap L1 pre-filter to skip the KDTree query for stations that # are clearly out of range. Mirrors geode/pyOkada.py # EarthquakeTable lines 168-171. _, sp_quick = l1_score(mag, max(dist, 1e-3)) if sp_quick <= 0: continue sc, sp = eq_score.score(stn['lat'], stn['lon']) else: sc, sp = l1_score(mag, dist) stn_entry = { 'name' : stn['name'], 'lat' : stn['lat'], 'lon' : stn['lon'], 'country' : stn.get('country', ''), 'city' : stn.get('city', ''), 'distance_km' : round(dist, 3), 'azimuth_deg' : round(az, 2), 'sc' : round(float(sc), 4), 'sp' : round(float(sp), 4), } if sc > 0: coseismic_stations.append(stn_entry) elif sp > 0 and include_postseismic: postseismic_stations.append(stn_entry) coseismic_stations.sort(key=lambda x: x['distance_km']) postseismic_stations.sort(key=lambda x: x['distance_km']) return { 'score_method' : score_method, 'nodal_planes_used' : nodal_planes_used, 'strikes' : strikes, 'dips' : dips, 'rakes' : rakes, 'n_coseismic' : len(coseismic_stations), 'n_postseismic' : len(postseismic_stations), 'coseismic_stations' : coseismic_stations, 'postseismic_stations' : postseismic_stations, }