prolog = """ **PROGRAM** eq_mask_L2.py **PURPOSE** | For one or more earthquake events, computes the Level-2 (Okada-based, | anisotropic) coseismic and postseismic influence-zone masks. | | The L2 method: | 1. Uses fault dimensions (length, width, slip) retrieved from USGS. | 2. Runs the Okada (1985) elastic half-space model over a 2D grid | of points centred on the epicenter, once per nodal plane, for both | surface and hypocentral depths. | 3. The boolean mask is formed where predicted displacement |U| >= threshold. | 4. The coseismic mask is based on the physical displacement threshold; | the postseismic mask is a geometrically scaled version of the | coseismic mask (using POST_SEISMIC_SCALE_FACTOR). | | For each event this program: | 1. Computes the L2 coseismic and postseismic masks | 2. Generates each mask as a list of geographic (lon, lat) vertices | 3. Saves one KMZ file per event (default dir = current directory) | 4. Optionally saves mask parameters to a JSON file (--output) | | If nodal plane data are unavailable, the program falls back to the L1 | isotropic circular mask. | | Primary design: the core functions (compute_l2_mask, save_mask_kmz) are | intended to be imported and called from a parent pipeline script. | | **USAGE — CLI** | python eq_mask_L2.py --event-id us7000abcd | python eq_mask_L2.py --date 2024-01-01 --min-mag 5.0 --threshold 0.001 | python eq_mask_L2.py --eq-json earthquake_events.json | python eq_mask_L2.py --event-id us7000abcd --no-kmz | python eq_mask_L2.py --event-id us7000abcd --no-json | python eq_mask_L2.py --event-id us7000abcd --no-postseismic """ epilog = """ **OUTPUTS** JSON file (--output, default: eq_masks_L2.json) with one entry per earthquake. Each entry contains: event_id, place, time_utc, magnitude, magnitude_type, longitude, latitude, depth_km, threshold_m, mask_type ('L2_okada' or 'L1_isotropic_fallback'), A, B, PSF, coseismic_polygon, postseismic_polygon If --kmz-dir is supplied, one KMZ per event is written there, named: _L2_mask.kmz Each KMZ contains: - A star icon at the epicenter with an HTML description balloon - A blue polygon for the coseismic mask - An orange polygon for the postseismic mask (omitted with --no-postseismic) **METHODOLOGY** 1. L2 Mask Generation: - Input: Fault dimensions (length, width, slip) and nodal planes retrieved from USGS. - Grid: A square mesh is created around the epicenter. - Simulation: Okada (1985) model computes displacements for each nodal plane at both surface and hypocentral depths. - Thresholding: Points where |U| >= threshold (e.g. 1mm) are flagged to create the coseismic mask. - Scaling: The postseismic mask is created by scaling the grid coordinates of the coseismic mask by POST_SEISMIC_SCALE_FACTOR. - Coordinates: Final vertices are mapped to Lon/Lat using an inverse azimuthal projection. 2. Polygon Extraction: -Contours are extracted from the boolean mask using find_contours (skimage). -Points are mapped to Lon/Lat using inverse azimuthal projection. 3. KMZ export: -Style and simplekml usage adapted from geode/pyOkada.py Score.save_masks(). **REFERENCES** - Gomez, D.D. et al. (2024). GPS Solutions 28:214. - Okada, Y. (1985). BSSA 75(4), 1135-1154. - Wells, D.L. & Coppersmith, K.J. (1994). BSSA 84(4), 974-1002. **AUTHORS** | Hasani Spann | Jet Propulsion Laboratory, California Institute of Technology | Pasadena, CA, USA | Assistance from ChatHPC """ import argparse import json import os import sys import numpy as np from datetime import datetime, timezone try: import simplekml from shapely.geometry import Polygon from skimage.measure import find_contours HAS_KMZ_LIBS = True except ImportError: HAS_KMZ_LIBS = False 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 ( inv_azimuthal, distance_km, destination_point, circle_polygon, okadakernel, okada, compute_l2_grid_mask ) # Constants POST_SEISMIC_SCALE_FACTOR = 1.5 # postseismic zone = coseismic zone x 1.5 EARTH_RADIUS_KM = 6371.0 GRID_DENSITY = 250 FAR_FIELD_SCALE = 25 # --- Trig Helpers --- cosd = lambda x: np.cos(np.deg2rad(x)) sind = lambda x: np.sin(np.deg2rad(x)) asind = lambda x: np.rad2deg(np.arcsin(x)) atand = lambda x: np.rad2deg(np.arctan(x)) # --- L2 Mask Core --- def compute_l2_mask_polygon( event: usgs.EarthquakeEvent, threshold_m: float = 1e-3, scale_factor: float = 1.0, density: int = GRID_DENSITY, n_points: int = 360 ) -> list: """ Compute the L2 (Okada) mask for one event and return it as a list of [lon, lat] vertices. Fallback to L1 circular mask if nodal planes are missing. Parameters ---------- event : EarthquakeEvent The earthquake event to process. threshold_m : float Displacement threshold in meters (default: 1mm). scale_factor : float Scale factor for the mask radius/area (default: 1.0 for coseismic). density : int Okada grid density per axis. n_points : int Number of points to use for the fallback L1 circular polygon. Returns ------- list A list of [lon, lat] vertex pairs defining the mask boundary polygon, or an empty list if no mask could be constructed. """ if event.latitude is None or event.longitude is None or event.magnitude is None: return [] # Check for nodal planes strikes, dips, rakes = [], [], [] np1, np2 = event.nodal_plane_1, event.nodal_plane_2 if all(getattr(np1, a, None) is not None for a in ['strike', 'dip', 'rake']): strikes.append(float(np1.strike)); dips.append(float(np1.dip)); rakes.append(float(np1.rake)) if all(getattr(np2, a, None) is not None for a in ['strike', 'dip', 'rake']): strikes.append(float(np2.strike)); dips.append(float(np2.dip)); rakes.append(float(np2.rake)) if not strikes: # L1 Fallback - without A and B, use a default distance based on length if available r_fallback = (event.fault_length_km if event.fault_length_km else 100.0) * scale_factor return circle_polygon(event.longitude, event.latitude, r_fallback, n_points) # Use the consolidated L2 grid mask function from eq_station_score depth_km = (event.depth_km if event.depth_km is not None else 10.0) mx, my, mask, _ = compute_l2_grid_mask( magnitude=event.magnitude, depth_km=depth_km, strikes=strikes, dips=dips, rakes=rakes, scale_factor=scale_factor, density=density, threshold_m=threshold_m ) # Extract contours and convert to lon/lat if not HAS_KMZ_LIBS: return [] contours = find_contours(mask.astype(float), 0.5) polygons = [] for c in contours: row_idx, col_idx = c[:, 0].astype(int), c[:, 1].astype(int) lons, lats = inv_azimuthal(mx[row_idx, col_idx], my[row_idx, col_idx], event.longitude, event.latitude) poly = Polygon(np.column_stack((lons, lats))) if poly.area > 0.0025: polygons.append(poly) if not polygons: return [] # Keep only outer polygons outer_polys = [] for p in polygons: if not any(outer.contains(p) for outer in outer_polys): outer_polys.append(p) main_poly = max(outer_polys, key=lambda p: p.area) return [list(pt) for pt in main_poly.exterior.coords] def compute_l2_mask( event: usgs.EarthquakeEvent, threshold_m: float = 1e-3, n_points: int = 360, density: int = GRID_DENSITY, include_postseismic: bool = False ) -> dict: """ Wrapper to compute L2 mask parameters and polygons for a single earthquake event. Parameters ---------- event : EarthquakeEvent The earthquake event to process. threshold_m : float Displacement threshold in meters (default: 1mm). n_points : int Number of points for fallback circular polygons. density : int Okada grid density per axis. include_postseismic : bool Whether to compute the postseismic mask (scaled by POST_SEISMIC_SCALE_FACTOR). Returns ------- dict A dictionary containing event metadata, mask type, min/max radii, and the vertex lists for the coseismic and postseismic polygons. """ if event.latitude is None or event.longitude is None or event.magnitude is None: return {'event_id': event.event_id, 'error': 'Missing event data.'} c_poly = compute_l2_mask_polygon(event, threshold_m, scale_factor=1.0, n_points=n_points, density=density) p_poly = compute_l2_mask_polygon(event, threshold_m, scale_factor=POST_SEISMIC_SCALE_FACTOR, n_points=n_points, density=density) if include_postseismic else [] def get_radii(poly): if not poly: return None, None dists = [distance_km(event.longitude, event.latitude, pt[0], pt[1]) for pt in poly] return min(dists), max(dists) c_min, c_max = get_radii(c_poly) p_min, p_max = get_radii(p_poly) return { '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, 'threshold_m': threshold_m, 'mask_type': 'L2_okada' if c_poly else 'L1_isotropic_fallback', 'PSF': POST_SEISMIC_SCALE_FACTOR, 'coseismic_min_radius_km': c_min, 'coseismic_max_radius_km': c_max, 'postseismic_min_radius_km': p_min, 'postseismic_max_radius_km': p_max, 'coseismic_polygon': c_poly, 'postseismic_polygon': p_poly, } def save_mask_kmz(mask: dict, kmz_path: str, include_postseismic: bool = True): """ Save a KMZ file containing the epicenter and the coseismic/postseismic displacement mask polygons for one event. Parameters ---------- mask : dict The mask result dictionary containing event data and polygons. kmz_path : str Output path for the .kmz file. include_postseismic : bool Whether to also export the postseismic polygon. """ if not HAS_KMZ_LIBS: print("[KMZ] WARNING: simplekml/shapely/skimage not installed. Skipping."); return if mask.get('error') or not mask.get('coseismic_polygon'): return kml = simplekml.Kml() epicenter = kml.newpoint(name=mask['place'], coords=[(mask['longitude'], mask['latitude'])]) epicenter.style.iconstyle.icon.href = 'https://maps.google.com/mapfiles/kml/shapes/star.png' epicenter.style.iconstyle.color = simplekml.Color.yellow epicenter.style.labelstyle.scale = 0 eid = mask['event_id'] epicenter.description = f"ID: {eid} | Mag: {mask['magnitude']} | Depth: {mask['depth_km']} km | Threshold: {mask['threshold_m']} m" def _add_poly(coords, name, color): if not coords: return poly = kml.newpolygon(name=name, outerboundaryis=[tuple(pt) for pt in coords]) poly.style.linestyle.color = color poly.style.linestyle.width = 3 poly.style.polystyle.color = simplekml.Color.changealphaint(0, simplekml.Color.white) _add_poly(mask['coseismic_polygon'], 'Coseismic mask (L2)', simplekml.Color.blue) if include_postseismic: _add_poly(mask['postseismic_polygon'], 'Postseismic mask (L2)', simplekml.Color.orange) kml.savekmz(kmz_path) print(f"[KMZ] Saved: {kmz_path}", flush=True) def _load_events_from_json(path: str) -> list: """ Load a list of EarthquakeEvent objects from a JSON file. Parameters ---------- path : str Path to the JSON file. Returns ------- list A list of EarthquakeEvent objects. """ with open(path, 'r') as f: records = json.load(f) events = [] for rec in records: t_str = rec.get('time_utc') time_utc = datetime.fromisoformat(t_str).replace(tzinfo=timezone.utc) if t_str else None 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', '') ) ev.nodal_plane_1 = usgs.NodalPlane(rec.get('np1_strike_deg'), rec.get('np1_dip_deg'), rec.get('np1_rake_deg')) ev.nodal_plane_2 = usgs.NodalPlane(rec.get('np2_strike_deg'), rec.get('np2_dip_deg'), rec.get('np2_rake_deg')) events.append(ev) return events def build_parser() -> argparse.ArgumentParser: """Construct the CLI argument parser for eq_mask_L2.py.""" p = argparse.ArgumentParser(prog='eq_mask_L2.py', description=prolog, epilog=epilog, formatter_class=argparse.RawDescriptionHelpFormatter) src = p.add_mutually_exclusive_group(required=True) src.add_argument('--event-id', metavar='USGS_ID') src.add_argument('--date', metavar='YYYY-MM-DD') src.add_argument('--eq-json', metavar='PATH') p.add_argument('--end-date', metavar='YYYY-MM-DD', default=None) p.add_argument('--min-mag', type=float, default=5.0) p.add_argument('--threshold', type=float, default=1e-3, help='Displacement limit in metres (default 1mm)') p.add_argument('--n-points', type=int, default=360) p.add_argument('--include-postseismic', action='store_true') p.add_argument('--output', default='eq_masks_L2.json') p.add_argument('--kmz-dir', default='.') p.add_argument('--no-kmz', action='store_true') p.add_argument('--no-json', action='store_true') p.add_argument('--verbose', '-v', action='store_true') return p def main(): """Main entry point for the eq_mask_L2 CLI.""" parser = build_parser() args = parser.parse_args() include_post = args.include_postseismic if args.event_id: ev = usgs.fetch_event_by_id(args.event_id) events = [ev] if ev else [] elif args.eq_json: events = _load_events_from_json(args.eq_json) else: events = usgs.query_earthquakes(start_time=args.date, end_time=args.end_date, min_magnitude=args.min_mag) if not events: return os.makedirs(args.kmz_dir, exist_ok=True) masks = [] for ev in events: mask = compute_l2_mask(ev, threshold_m=args.threshold, n_points=args.n_points, include_postseismic=include_post) masks.append(mask) if not args.no_kmz and not mask.get('error'): save_mask_kmz(mask, os.path.join(args.kmz_dir, f"{ev.event_id}_L2_mask.kmz"), include_postseismic=include_post) if not args.no_json: with open(args.output, 'w') as f: json.dump(masks, f, indent=2) if __name__ == '__main__': main()