#!/usr/bin/env python3 """ Create synthetic HypoNet_Hikurangi event files from a 3-D velocity-model grid. Default behaviour ----------------- 3D mode: - Build synthetic hypocenters on a regular Model X / Model Y / Depth grid. - Use the velocity-model table to map Model X/Y nodes to lon/lat. - Apply an edge buffer equal to the requested horizontal/vertical increments. - Reject events deeper than MAX_DEPTH_KM. - Reject events above topography/seafloor using a GMT grid sampled with grdtrack. 2D mode: - Build synthetic hypocenters along PROFILE_START -> PROFILE_END at regular horizontal spacing, and regular depth spacing. - Snap each profile point to the nearest horizontal node in the velocity-model table, so output lon/lat are valid model-node lon/lat values. - Apply the same seafloor/topography and max-depth checks. For each event: - Read station locations from a SIMUL used_stations.txt file. - Select all stations within MAX_OFFSET_KM. - If fewer than STA_MIN stations are within MAX_OFFSET_KM, add nearest extra stations until STA_MIN is satisfied. - Calculate P and S travel times using HypoNet_DeeperPINN. - Add Gaussian random noise to P and S picks. - Write one HypoNet event file per synthetic event. - Write a summary CSV. Run from the HypoNet_Hikurangi/Scripts directory, or adjust paths below. """ from __future__ import annotations import csv import math import os import re import shutil import subprocess import tempfile from dataclasses import dataclass from pathlib import Path from typing import Iterable, Optional import numpy as np import torch import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from HypoNet_Hikurangi.HypoNet_DeeperPINN import HypoNet_DeeperPINN try: from scipy.spatial import cKDTree except Exception: cKDTree = None # ============================================================================= # CONFIGURATION # ============================================================================= # Main mode: "3D" or "2D" MODE = "2D" VELOCITY_MODEL_FILE = Path("../Velocity_model_files/Hikurangi_3D_EQ_model_for_Agata.txt") BATHY_GRID = Path("../Velocity_model_files/Bathy_cut.grd") STATION_FILE = Path( "/scratch/dbassett/Earthquake_relocation/Simul2023/Working_relocation/" "Transect_SIMUL/used_stations.txt" ) OUTPUT_DIR = Path("../Synthetic_Event_Files") SUMMARY_FILE = OUTPUT_DIR / "synthetic_event_summary.csv" # Synthetic hypocenter spacing HORIZONTAL_INCREMENT_KM = 10.0 VERTICAL_INCREMENT_KM = 5.0 # Depth limits. Depth is positive down from sea level. MIN_DEPTH_KM = 0.0 MAX_DEPTH_KM = 70.0 # 2-D profile endpoints are (lat, lon) #PROFILE_START = (-40.2, 174.5) #PROFILE_END = (-41.3, 176.3) PROFILE_START = (-40.2, 174.5) PROFILE_END = (-41.8, 177) # Station selection STA_MIN = 5 MAX_OFFSET_KM = 150.0 # Pick uncertainty written to HypoNet event file, and random noise added to picks P_PICK_ERROR_S = 0.1 S_PICK_ERROR_S = 0.1 GAUSSIAN_NOISE_STD_S = 0.2 RANDOM_SEED = 12345 # HypoNet model input directories INPUT_P = "./input_p/" INPUT_S = "./input_s/" DEVICE = 0 FORCE_CPU_AFTER_MODEL_LOAD = True # If True, remove old event.*.txt files before writing new ones. CLEAN_OUTPUT_DIR_FIRST = True # If True, events where HypoNet returns non-zero errcodes are still written. # Set False if you want to skip phases/stations with non-zero errcodes. KEEP_NONZERO_ERRCODES = True # Progress printing PRINT_EVERY_N_EVENTS = 50 # ============================================================================= # DATA CLASSES # ============================================================================= @dataclass(frozen=True) class HorizontalNode: model_x_km: float model_y_km: float lon: float lat: float @dataclass(frozen=True) class SyntheticEvent: event_index: int model_x_km: float model_y_km: float depth_km: float lon: float lat: float surface_z_m: float source_mode: str profile_distance_km: Optional[float] = None @dataclass(frozen=True) class Station: code: str lat: float lon: float elev_m: float # ============================================================================= # SMALL HELPERS # ============================================================================= def normalize_lon_180(lon: float) -> float: """Return longitude in [-180, 180).""" return ((lon + 180.0) % 360.0) - 180.0 def normalize_lon_360(lon: float) -> float: """Return longitude in [0, 360).""" return lon % 360.0 def normalize_lon_for_hyponet(lon: float) -> float: """ Return longitude in the convention safest for HypoNet_Hikurangi. The Hikurangi grids span the dateline and the bathy grid is 170/185, so points just east of the dateline should be passed as 180-185 rather than -180 to -175. Passing negative longitudes here can make HypoNet's internal geodetic interpolator index outside its EGM array. """ return normalize_lon_360(lon) def angular_lon_delta_deg(lon1: float, lon2: float) -> float: """Shortest signed longitude difference in degrees.""" return ((lon2 - lon1 + 180.0) % 360.0) - 180.0 def haversine_km(lat1: float, lon1: float, lat2: float, lon2: float) -> float: """Great-circle distance in km, dateline-safe.""" r_earth_km = 6371.0 phi1 = math.radians(lat1) phi2 = math.radians(lat2) dphi = math.radians(lat2 - lat1) dlambda = math.radians(angular_lon_delta_deg(lon1, lon2)) a = ( math.sin(dphi / 2.0) ** 2 + math.cos(phi1) * math.cos(phi2) * math.sin(dlambda / 2.0) ** 2 ) return 2.0 * r_earth_km * math.asin(min(1.0, math.sqrt(a))) def interpolate_profile_point( start_lat: float, start_lon: float, end_lat: float, end_lon: float, fraction: float, ) -> tuple[float, float]: """ Approximate point along a short profile. The longitude interpolation is dateline-safe. For your profile near Hikurangi this is effectively identical to a geodesic interpolation at 10 km spacing. """ lat = start_lat + fraction * (end_lat - start_lat) dlon = angular_lon_delta_deg(start_lon, end_lon) lon = normalize_lon_180(start_lon + fraction * dlon) return lat, lon def frange_inclusive(start: float, stop: float, step: float) -> list[float]: vals = [] x = start eps = abs(step) * 1.0e-9 while x <= stop + eps: vals.append(round(x, 8)) x += step return vals def snap_up_to_increment(value: float, inc: float) -> float: return math.ceil((value - 1.0e-9) / inc) * inc def snap_down_to_increment(value: float, inc: float) -> float: return math.floor((value + 1.0e-9) / inc) * inc # ============================================================================= # VELOCITY MODEL READING # ============================================================================= def read_velocity_model_horizontal_nodes(path: Path) -> list[HorizontalNode]: """ Read unique horizontal nodes from the velocity model. The input table repeats lon/lat for every depth at each Model X/Y node. We keep the first occurrence for each X/Y pair. """ if not path.exists(): raise FileNotFoundError(f"Velocity model file not found: {path}") nodes_by_xy: dict[tuple[float, float], HorizontalNode] = {} with path.open("r", encoding="utf-8") as f: for line_number, line in enumerate(f, start=1): line = line.strip() if not line or line.startswith("(") or line.startswith("#"): continue parts = line.split() if len(parts) < 7: continue try: model_x = float(parts[0]) model_y = float(parts[1]) lon = float(parts[5]) lat = float(parts[6]) except ValueError: continue key = (model_x, model_y) if key not in nodes_by_xy: nodes_by_xy[key] = HorizontalNode( model_x_km=model_x, model_y_km=model_y, lon=normalize_lon_for_hyponet(lon), lat=lat, ) nodes = list(nodes_by_xy.values()) if not nodes: raise ValueError(f"No horizontal nodes were read from {path}") print(f"Read unique horizontal velocity-model nodes: {len(nodes)}") return nodes def build_xy_lookup(nodes: list[HorizontalNode]) -> dict[tuple[float, float], HorizontalNode]: return {(n.model_x_km, n.model_y_km): n for n in nodes} def nearest_horizontal_nodes( nodes: list[HorizontalNode], query_lats: np.ndarray, query_lons: np.ndarray, ) -> list[HorizontalNode]: """Find nearest velocity-model horizontal node for each lon/lat query.""" # Equirectangular km coordinates around the mean latitude are fine for nearest # node matching at this scale. mean_lat = float(np.nanmean(query_lats)) coslat = math.cos(math.radians(mean_lat)) node_x = np.array([normalize_lon_360(n.lon) * 111.32 * coslat for n in nodes]) node_y = np.array([n.lat * 111.32 for n in nodes]) query_x = np.array([normalize_lon_360(lon) * 111.32 * coslat for lon in query_lons]) query_y = np.array([lat * 111.32 for lat in query_lats]) node_xy = np.column_stack([node_x, node_y]) query_xy = np.column_stack([query_x, query_y]) if cKDTree is not None: tree = cKDTree(node_xy) _, idx = tree.query(query_xy, k=1) return [nodes[int(i)] for i in idx] # Fallback without scipy. This is slower but okay for modest 2-D profiles. out = [] for qx, qy in query_xy: dist2 = (node_x - qx) ** 2 + (node_y - qy) ** 2 out.append(nodes[int(np.argmin(dist2))]) return out def lonlat_to_local_xy_km(lon: float, lat: float, mean_lat: float) -> tuple[float, float]: """Approximate lon/lat as local Cartesian km, dateline-safe using 0-360 lon.""" x = normalize_lon_360(lon) * 111.32 * math.cos(math.radians(mean_lat)) y = lat * 111.32 return x, y def point_in_polygon(x: float, y: float, polygon_xy: list[tuple[float, float]]) -> bool: """ Ray-casting point-in-polygon test. Points on or very close to the boundary are treated as inside. """ inside = False n = len(polygon_xy) eps = 1.0e-9 for i in range(n): x1, y1 = polygon_xy[i] x2, y2 = polygon_xy[(i + 1) % n] # Boundary check dx = x2 - x1 dy = y2 - y1 cross = (x - x1) * dy - (y - y1) * dx if abs(cross) < eps: dot = (x - x1) * dx + (y - y1) * dy if -eps <= dot <= dx * dx + dy * dy + eps: return True if (y1 > y) != (y2 > y): xinters = (x2 - x1) * (y - y1) / (y2 - y1 + eps) + x1 if x <= xinters: inside = not inside return inside def build_velocity_model_domain_polygon(nodes: list[HorizontalNode]) -> list[tuple[float, float]]: """ Build a polygon around the trained/velocity-model domain. The velocity model is a rotated X/Y grid, so a simple lon/lat min/max box can include stations that are outside the actual trained model. This uses the boundary horizontal nodes: min/max Model X or min/max Model Y. """ xs = [n.model_x_km for n in nodes] ys = [n.model_y_km for n in nodes] x_min, x_max = min(xs), max(xs) y_min, y_max = min(ys), max(ys) boundary_nodes = [ n for n in nodes if ( math.isclose(n.model_x_km, x_min) or math.isclose(n.model_x_km, x_max) or math.isclose(n.model_y_km, y_min) or math.isclose(n.model_y_km, y_max) ) ] if len(boundary_nodes) < 4: raise ValueError("Could not build velocity-model domain polygon from boundary nodes") mean_lat = float(np.mean([n.lat for n in boundary_nodes])) boundary_xy = [lonlat_to_local_xy_km(n.lon, n.lat, mean_lat) for n in boundary_nodes] cx = float(np.mean([p[0] for p in boundary_xy])) cy = float(np.mean([p[1] for p in boundary_xy])) # Sorting all boundary points by angle gives a closed outline of the convex # model footprint, which is appropriate for this rotated rectangular grid. boundary_xy_sorted = sorted(boundary_xy, key=lambda p: math.atan2(p[1] - cy, p[0] - cx)) return boundary_xy_sorted def station_is_inside_velocity_model_domain( station: Station, domain_polygon_xy: list[tuple[float, float]], mean_lat: float, ) -> bool: sx, sy = lonlat_to_local_xy_km(station.lon, station.lat, mean_lat) return point_in_polygon(sx, sy, domain_polygon_xy) def filter_stations_to_velocity_model_domain( stations: list[Station], nodes: list[HorizontalNode], ) -> list[Station]: """Remove stations outside the trained model footprint.""" domain_polygon_xy = build_velocity_model_domain_polygon(nodes) mean_lat = float(np.mean([n.lat for n in nodes])) kept = [] removed = [] for sta in stations: if station_is_inside_velocity_model_domain(sta, domain_polygon_xy, mean_lat): kept.append(sta) else: removed.append(sta) print("Station trained-domain filter:") print(f" Stations inside velocity-model domain: {len(kept)}") print(f" Stations outside velocity-model domain: {len(removed)}") if removed: preview = ", ".join(s.code for s in removed[:20]) suffix = " ..." if len(removed) > 20 else "" print(f" Removed stations: {preview}{suffix}") if not kept: raise ValueError("All stations were excluded by the velocity-model domain filter") return kept # ============================================================================= # BATHYMETRY / TOPOGRAPHY SAMPLING # ============================================================================= def sample_grid_with_gmt_grdtrack( grid_path: Path, lons: Iterable[float], lats: Iterable[float], ) -> np.ndarray: """ Sample a GMT grid using `gmt grdtrack`. Input longitudes are written as 0-360 because Bathy_cut.grd has x_min=170, x_max=185 and the velocity model spans the dateline. """ if not grid_path.exists(): raise FileNotFoundError(f"Bathymetry/topography grid not found: {grid_path}") if shutil.which("gmt") is None: raise RuntimeError( "Could not find GMT executable. Install/load GMT or replace " "sample_grid_with_gmt_grdtrack() with an xarray/netCDF sampler." ) lons_list = list(lons) lats_list = list(lats) if len(lons_list) != len(lats_list): raise ValueError("lons and lats have different lengths") with tempfile.TemporaryDirectory() as tmpdir: in_path = Path(tmpdir) / "points.txt" out_path = Path(tmpdir) / "points_with_z.txt" with in_path.open("w", encoding="utf-8") as f: for lon, lat in zip(lons_list, lats_list): f.write(f"{normalize_lon_360(lon):.8f} {lat:.8f}\n") cmd = ["gmt", "grdtrack", str(in_path), f"-G{grid_path}"] with out_path.open("w", encoding="utf-8") as fout: subprocess.run(cmd, check=True, stdout=fout, text=True) z_vals = [] with out_path.open("r", encoding="utf-8") as f: for line in f: parts = line.split() if len(parts) < 3: z_vals.append(np.nan) continue try: z_vals.append(float(parts[2])) except ValueError: z_vals.append(np.nan) return np.asarray(z_vals, dtype=float) def depth_is_below_surface(depth_km: float, surface_z_m: float) -> bool: """ Return True if hypocenter is below the topography/seafloor. GMT topo/bathy z is metres relative to sea level: z > 0 on land, z < 0 offshore. Depth is km positive down from sea level. The surface depth is therefore -z/1000. Offshore at z=-3000 m -> surface_depth=3 km, so an event must be deeper than 3 km. On land at z=500 m -> surface_depth=-0.5 km, and the MIN_DEPTH_KM setting normally controls the shallowest earthquake depth. """ surface_depth_km = -surface_z_m / 1000.0 min_allowed_depth = max(MIN_DEPTH_KM, surface_depth_km) return depth_km >= min_allowed_depth # ============================================================================= # SYNTHETIC EVENT GENERATION # ============================================================================= def generate_3d_candidate_events(nodes: list[HorizontalNode]) -> list[SyntheticEvent]: xs = sorted({n.model_x_km for n in nodes}) ys = sorted({n.model_y_km for n in nodes}) lookup = build_xy_lookup(nodes) x_min = min(xs) + HORIZONTAL_INCREMENT_KM x_max = max(xs) - HORIZONTAL_INCREMENT_KM y_min = min(ys) + HORIZONTAL_INCREMENT_KM y_max = max(ys) - HORIZONTAL_INCREMENT_KM x_start = snap_up_to_increment(x_min, HORIZONTAL_INCREMENT_KM) x_stop = snap_down_to_increment(x_max, HORIZONTAL_INCREMENT_KM) y_start = snap_up_to_increment(y_min, HORIZONTAL_INCREMENT_KM) y_stop = snap_down_to_increment(y_max, HORIZONTAL_INCREMENT_KM) z_start = max(MIN_DEPTH_KM, VERTICAL_INCREMENT_KM) z_stop = min(MAX_DEPTH_KM, 75.0 - VERTICAL_INCREMENT_KM) x_vals = frange_inclusive(x_start, x_stop, HORIZONTAL_INCREMENT_KM) y_vals = frange_inclusive(y_start, y_stop, HORIZONTAL_INCREMENT_KM) z_vals = frange_inclusive(z_start, z_stop, VERTICAL_INCREMENT_KM) events = [] for x in x_vals: for y in y_vals: node = lookup.get((x, y)) if node is None: # This should not normally happen if increments match model nodes. continue for z in z_vals: events.append( SyntheticEvent( event_index=-1, model_x_km=node.model_x_km, model_y_km=node.model_y_km, depth_km=z, lon=node.lon, lat=node.lat, surface_z_m=np.nan, source_mode="3D", ) ) print( "3D candidates before bathy/topo filter: " f"{len(events)} ({len(x_vals)} X nodes, {len(y_vals)} Y nodes, {len(z_vals)} depths)" ) return events def generate_2d_candidate_events(nodes: list[HorizontalNode]) -> list[SyntheticEvent]: start_lat, start_lon = PROFILE_START end_lat, end_lon = PROFILE_END profile_length_km = haversine_km(start_lat, start_lon, end_lat, end_lon) distances = frange_inclusive(0.0, profile_length_km, HORIZONTAL_INCREMENT_KM) if distances and distances[-1] < profile_length_km: # Do not force the exact endpoint; regular spacing is preferred. pass query_lats = [] query_lons = [] for d in distances: frac = 0.0 if profile_length_km == 0.0 else d / profile_length_km lat, lon = interpolate_profile_point(start_lat, start_lon, end_lat, end_lon, frac) query_lats.append(lat) query_lons.append(lon) snapped_nodes = nearest_horizontal_nodes( nodes, np.asarray(query_lats, dtype=float), np.asarray(query_lons, dtype=float), ) z_start = max(MIN_DEPTH_KM, VERTICAL_INCREMENT_KM) z_stop = min(MAX_DEPTH_KM, 75.0 - VERTICAL_INCREMENT_KM) z_vals = frange_inclusive(z_start, z_stop, VERTICAL_INCREMENT_KM) events = [] seen = set() for d, node in zip(distances, snapped_nodes): # Avoid duplicate profile points snapping to the same model node at very # fine spacing. horiz_key = (node.model_x_km, node.model_y_km, round(d, 6)) if horiz_key in seen: continue seen.add(horiz_key) for z in z_vals: events.append( SyntheticEvent( event_index=-1, model_x_km=node.model_x_km, model_y_km=node.model_y_km, depth_km=z, lon=node.lon, lat=node.lat, surface_z_m=np.nan, source_mode="2D", profile_distance_km=d, ) ) print( "2D candidates before bathy/topo filter: " f"{len(events)} ({len(distances)} profile positions, {len(z_vals)} depths, " f"profile length {profile_length_km:.2f} km)" ) return events def apply_bathy_topo_filter(events: list[SyntheticEvent]) -> list[SyntheticEvent]: if not events: return [] # Only sample each unique horizontal location once. unique_horiz = {} for ev in events: unique_horiz[(ev.lon, ev.lat)] = (ev.lon, ev.lat) unique_lons = [v[0] for v in unique_horiz.values()] unique_lats = [v[1] for v in unique_horiz.values()] unique_z = sample_grid_with_gmt_grdtrack(BATHY_GRID, unique_lons, unique_lats) z_by_lonlat = { key: z for key, z in zip(unique_horiz.keys(), unique_z) } kept = [] skipped_nan_surface = 0 skipped_above_surface = 0 for ev in events: z_m = z_by_lonlat[(ev.lon, ev.lat)] if not np.isfinite(z_m): skipped_nan_surface += 1 continue if not depth_is_below_surface(ev.depth_km, z_m): skipped_above_surface += 1 continue kept.append( SyntheticEvent( event_index=-1, model_x_km=ev.model_x_km, model_y_km=ev.model_y_km, depth_km=ev.depth_km, lon=ev.lon, lat=ev.lat, surface_z_m=float(z_m), source_mode=ev.source_mode, profile_distance_km=ev.profile_distance_km, ) ) print(f"Events after bathy/topo filter: {len(kept)}") print(f" Skipped above topography/seafloor: {skipped_above_surface}") print(f" Skipped with NaN surface samples: {skipped_nan_surface}") return kept # ============================================================================= # STATION FILE PARSING # ============================================================================= def parse_lat_token(token: str) -> Optional[tuple[str, float]]: """ Parse a latitude token. Supported examples: AKTO40S32.39 -> station prefix AKTO, latitude -40.539833 BFZ -> not a latitude token by itself 40S40.78 -> no station prefix, latitude -40.679667 Important: the station-code prefix must be non-greedy. Otherwise a token like AKTO40S32.39 can be misread as prefix=AKTO4, deg=0, or similar. """ m = re.match(r"^([A-Za-z0-9]*?)(\d{2})([NS])(\d+(?:\.\d*)?)$", token) if not m: return None prefix, deg_s, hemi, min_s = m.groups() deg = float(deg_s) minutes = float(min_s) if not (0.0 <= minutes < 60.0): return None val = deg + minutes / 60.0 if hemi.upper() == "S": val *= -1.0 return prefix, val def parse_lon_token(token: str) -> Optional[float]: """Parse a token like 176E27.67, 173W05.00, or decimal lon.""" try: return normalize_lon_for_hyponet(float(token)) except ValueError: pass m = re.match(r"^(\d{1,3})([EW])(\d+(?:\.\d*)?)$", token) if not m: return None deg_s, hemi, min_s = m.groups() val = float(deg_s) + float(min_s) / 60.0 if hemi.upper() == "W": val *= -1.0 return normalize_lon_for_hyponet(val) def read_simul_used_stations(path: Path) -> list[Station]: if not path.exists(): raise FileNotFoundError(f"Station file not found: {path}") stations = [] skipped = 0 with path.open("r", encoding="utf-8") as f: for raw in f: line = raw.strip() if not line: continue parts = line.split() if len(parts) < 3: skipped += 1 continue code = None lat = None lon = None elev_token_index = None # Case 1: station and latitude are combined: # AKTO40S32.39 176E27.67 -432 ... lat_info = parse_lat_token(parts[0]) if lat_info is not None and len(parts) >= 3: prefix, parsed_lat = lat_info parsed_lon = parse_lon_token(parts[1]) if prefix and parsed_lon is not None: code = prefix lat = parsed_lat lon = parsed_lon elev_token_index = 2 # Case 2: station and latitude are separate: # BFZ 40S40.78 176E14.77 -283 ... if code is None and len(parts) >= 4: lat_info = parse_lat_token(parts[1]) parsed_lon = parse_lon_token(parts[2]) if lat_info is not None else None if lat_info is not None and parsed_lon is not None: prefix, parsed_lat = lat_info if prefix == "": code = parts[0] lat = parsed_lat lon = parsed_lon elev_token_index = 3 if code is None or lat is None or lon is None or elev_token_index is None: skipped += 1 continue # Guard against accidental parsing of headers or malformed lines. if not (-50.0 <= lat <= -30.0 and 160.0 <= normalize_lon_360(lon) <= 190.0): print(f"WARNING: skipping suspicious station parse: {line}") print(f" parsed as code={code}, lat={lat}, lon={lon}") skipped += 1 continue try: elev_m = float(parts[elev_token_index]) except ValueError: elev_m = 0.0 stations.append(Station(code=code, lat=lat, lon=normalize_lon_for_hyponet(lon), elev_m=elev_m)) if not stations: raise ValueError(f"No stations parsed from {path}") print(f"Read stations: {len(stations)}") print(f"Skipped non-station/header lines: {skipped}") print( "Station coordinate ranges after parsing:\n" f" lat: {min(s.lat for s in stations):.5f} to {max(s.lat for s in stations):.5f}\n" f" lon: {min(s.lon for s in stations):.5f} to {max(s.lon for s in stations):.5f}" ) return stations # ============================================================================= # HYPONET TRAVEL TIMES # ============================================================================= def tensor(values) -> torch.Tensor: return torch.tensor(values, dtype=torch.float32, device="cpu") def run_inference(model, hypo_lon, hypo_lat, hypo_elev, station_lons, station_lats): # Defensive coordinate check before HypoNet's internal geodetic interpolator. # This gives a useful error instead of an opaque EGM IndexError. for name, arr, lo, hi in [ ("hypo_lon", hypo_lon, 160.0, 190.0), ("station_lons", station_lons, 160.0, 190.0), ("hypo_lat", hypo_lat, -50.0, -30.0), ("station_lats", station_lats, -50.0, -30.0), ]: vals = arr.detach().cpu().numpy() bad = ~np.isfinite(vals) | (vals < lo) | (vals > hi) if np.any(bad): raise ValueError( f"Bad coordinate values passed to HypoNet for {name}: " f"min={np.nanmin(vals):.6f}, max={np.nanmax(vals):.6f}, " f"allowed approx range={lo}..{hi}. First bad values: {vals[bad][:10]}" ) tt, errcodes = model.inference_torch( hypo_lon, hypo_lat, hypo_elev, station_lons, station_lats, ) tt_list = tt.detach().cpu().reshape(-1).tolist() if hasattr(errcodes, "detach"): err_list = errcodes.detach().cpu().reshape(-1).tolist() else: err_list = list(errcodes) return tt_list, err_list def select_stations_for_event(event: SyntheticEvent, stations: list[Station]) -> list[tuple[Station, float]]: station_offsets = [ (sta, haversine_km(event.lat, event.lon, sta.lat, sta.lon)) for sta in stations ] station_offsets.sort(key=lambda x: x[1]) selected = [(sta, dist) for sta, dist in station_offsets if dist <= MAX_OFFSET_KM] if len(selected) < STA_MIN: selected_codes = {sta.code for sta, _ in selected} for sta, dist in station_offsets: if sta.code in selected_codes: continue selected.append((sta, dist)) selected_codes.add(sta.code) if len(selected) >= STA_MIN: break return selected def calculate_event_times( event: SyntheticEvent, selected: list[tuple[Station, float]], model_p, model_s, rng: np.random.Generator, ) -> list[dict]: n = len(selected) if n == 0: return [] stas = [x[0] for x in selected] offsets = [x[1] for x in selected] station_lons = tensor([normalize_lon_for_hyponet(s.lon) for s in stas]) station_lats = tensor([s.lat for s in stas]) hypo_lon = tensor([normalize_lon_for_hyponet(event.lon)] * n) hypo_lat = tensor([event.lat] * n) hypo_elev = tensor([-event.depth_km] * n) # HypoNet uses elevation, km; depth is positive down. p_tt, p_err = run_inference(model_p, hypo_lon, hypo_lat, hypo_elev, station_lons, station_lats) s_tt, s_err = run_inference(model_s, hypo_lon, hypo_lat, hypo_elev, station_lons, station_lats) rows = [] for sta, offset_km, p, pe, s, se in zip(stas, offsets, p_tt, p_err, s_tt, s_err): p_ok = np.isfinite(p) and (KEEP_NONZERO_ERRCODES or int(pe) == 0) s_ok = np.isfinite(s) and (KEEP_NONZERO_ERRCODES or int(se) == 0) p_obs = None s_obs = None if p_ok: p_obs = float(p) + float(rng.normal(0.0, GAUSSIAN_NOISE_STD_S)) if s_ok: s_obs = float(s) + float(rng.normal(0.0, GAUSSIAN_NOISE_STD_S)) rows.append( { "station": sta, "offset_km": offset_km, "p_tt_s": p_obs, "s_tt_s": s_obs, "p_errcode": int(pe), "s_errcode": int(se), } ) return rows # ============================================================================= # OUTPUT WRITING # ============================================================================= def format_pick(value: Optional[float], err_s: float, flag_when_present: str = "0") -> tuple[str, str, str]: if value is None or not np.isfinite(value): return "-", "-", "-" return f"{value:.3f}", f"{err_s:.1f}", flag_when_present def write_event_file(path: Path, phase_rows: list[dict]) -> tuple[int, int]: n_p = 0 n_s = 0 with path.open("w", encoding="utf-8") as f: for row in phase_rows: sta: Station = row["station"] p_val, p_err, p_flag = format_pick(row["p_tt_s"], P_PICK_ERROR_S) s_val, s_err, s_flag = format_pick(row["s_tt_s"], S_PICK_ERROR_S) if p_val != "-": n_p += 1 if s_val != "-": n_s += 1 # Match attached HypoNet event-file format: # station Ptt Perr Stt Serr station_lat station_lon Pflag Sflag f.write( f"{sta.code:<5s} " f"{p_val:>8s} {p_err:>3s} " f"{s_val:>8s} {s_err:>3s} " f"{sta.lat:.4f} {sta.lon:.4f} " f"{p_flag:>1s} {s_flag:>1s}\n" ) return n_p, n_s def prepare_output_dir(path: Path) -> None: path.mkdir(parents=True, exist_ok=True) if CLEAN_OUTPUT_DIR_FIRST: for old in path.glob("event.*.txt"): old.unlink() def write_summary_header(writer: csv.DictWriter) -> None: writer.writeheader() # ============================================================================= # PLOTTING # ============================================================================= def sample_bathy_points(lons: list[float], lats: list[float]) -> np.ndarray: """Convenience wrapper for sampling the bathy/topo grid.""" return sample_grid_with_gmt_grdtrack(BATHY_GRID, lons, lats) def plot_synthetic_events( events: list[SyntheticEvent], stations: list[Station], nodes: list[HorizontalNode], ) -> None: """Create a diagnostic plot for the synthetic event distribution.""" if not events: print("Skipping plot: no synthetic events") return mode = MODE.strip().upper() if mode == "2D": output_png = OUTPUT_DIR / "synthetic_events_2D_profile_and_map.png" plot_2d_synthetic_events(events, stations, output_png) elif mode == "3D": output_png = OUTPUT_DIR / "synthetic_events_3D_slices.png" plot_3d_synthetic_events(events, stations, nodes, output_png) else: print(f"Skipping plot: unknown MODE={MODE}") def plot_2d_synthetic_events( events: list[SyntheticEvent], stations: list[Station], output_png: Path, ) -> None: """ 2-D diagnostic plot: 1. Cross-section with bathymetry/topography and synthetic EQs. 2. Map with profile line, synthetic EQs, and stations included in event files. """ start_lat, start_lon = PROFILE_START end_lat, end_lon = PROFILE_END profile_length_km = haversine_km(start_lat, start_lon, end_lat, end_lon) # Smooth bathy/topo profile. n_profile_samples = max(200, int(profile_length_km / 0.5) + 1) profile_dist = np.linspace(0.0, profile_length_km, n_profile_samples) profile_lats = [] profile_lons = [] for d in profile_dist: frac = 0.0 if profile_length_km == 0.0 else d / profile_length_km lat, lon = interpolate_profile_point(start_lat, start_lon, end_lat, end_lon, frac) profile_lats.append(lat) profile_lons.append(lon) surface_z_m = sample_bathy_points(profile_lons, profile_lats) surface_depth_km = -surface_z_m / 1000.0 ev_profile_dist = np.array([ ev.profile_distance_km if ev.profile_distance_km is not None else np.nan for ev in events ]) ev_depth = np.array([ev.depth_km for ev in events]) ev_lons = np.array([normalize_lon_360(ev.lon) for ev in events]) ev_lats = np.array([ev.lat for ev in events]) station_lons = np.array([normalize_lon_360(s.lon) for s in stations]) station_lats = np.array([s.lat for s in stations]) fig, axes = plt.subplots(2, 1, figsize=(11, 10), constrained_layout=True) ax = axes[0] ax.plot(profile_dist, surface_depth_km, linewidth=1.5, label="Bathymetry/topography") ax.scatter(ev_profile_dist, ev_depth, s=18, edgecolors="k", linewidths=0.3, label="Synthetic EQs") ax.axhline(0.0, linewidth=0.8, linestyle="--") ax.set_xlabel("Distance along profile (km)") ax.set_ylabel("Depth below sea level (km)") ax.set_title("Synthetic events along 2-D profile") max_depth = max(MAX_DEPTH_KM, float(np.nanmax(ev_depth)) if len(ev_depth) else MAX_DEPTH_KM) ax.set_ylim(max_depth + 5.0, min(-2.0, float(np.nanmin(surface_depth_km)) - 1.0)) ax.grid(True, alpha=0.3) ax.legend(loc="best") ax = axes[1] ax.plot([normalize_lon_360(start_lon), normalize_lon_360(end_lon)], [start_lat, end_lat], linewidth=2.0, label="Profile") ax.scatter(station_lons, station_lats, marker="^", s=35, edgecolors="k", linewidths=0.3, label="Stations used") ax.scatter(ev_lons, ev_lats, s=18, edgecolors="k", linewidths=0.3, label="Synthetic EQ map positions") ax.set_xlabel("Longitude (0-360)") ax.set_ylabel("Latitude") ax.set_title("2-D profile, synthetic events, and stations") all_lons = np.concatenate([ev_lons, station_lons, np.array([normalize_lon_360(start_lon), normalize_lon_360(end_lon)])]) all_lats = np.concatenate([ev_lats, station_lats, np.array([start_lat, end_lat])]) ax.set_xlim(float(np.nanmin(all_lons)) - 0.25, float(np.nanmax(all_lons)) + 0.25) ax.set_ylim(float(np.nanmin(all_lats)) - 0.25, float(np.nanmax(all_lats)) + 0.25) ax.grid(True, alpha=0.3) ax.legend(loc="best") output_png.parent.mkdir(parents=True, exist_ok=True) fig.savefig(output_png, dpi=300) plt.close(fig) print(f"Saved 2D synthetic-event plot to: {output_png}") def unique_sorted_values(vals: Iterable[float]) -> list[float]: return sorted({float(v) for v in vals}) def nearest_value(values: list[float], target: float) -> float: return min(values, key=lambda v: abs(v - target)) def plot_3d_synthetic_events( events: list[SyntheticEvent], stations: list[Station], nodes: list[HorizontalNode], output_png: Path, ) -> None: """ 3-D diagnostic plot: 1. Map/depth slice near 10 km showing coastline, stations, EQs. 2. X slice through the nearest middle X node showing topography/bathymetry and EQs. 3. Y slice through the nearest middle Y node showing topography/bathymetry and EQs. """ ev_x = np.array([ev.model_x_km for ev in events]) ev_y = np.array([ev.model_y_km for ev in events]) ev_z = np.array([ev.depth_km for ev in events]) ev_lons = np.array([normalize_lon_360(ev.lon) for ev in events]) ev_lats = np.array([ev.lat for ev in events]) xs = unique_sorted_values(n.model_x_km for n in nodes) ys = unique_sorted_values(n.model_y_km for n in nodes) zs = unique_sorted_values(ev_z) depth_slice_km = nearest_value(zs, 10.0) x_slice_km = nearest_value(xs, 0.5 * (min(xs) + max(xs))) y_slice_km = nearest_value(ys, 0.5 * (min(ys) + max(ys))) depth_tol = max(0.01, VERTICAL_INCREMENT_KM * 0.5) x_tol = max(0.01, HORIZONTAL_INCREMENT_KM * 0.5) y_tol = max(0.01, HORIZONTAL_INCREMENT_KM * 0.5) depth_mask = np.abs(ev_z - depth_slice_km) <= depth_tol x_mask = np.abs(ev_x - x_slice_km) <= x_tol y_mask = np.abs(ev_y - y_slice_km) <= y_tol station_lons = np.array([normalize_lon_360(s.lon) for s in stations]) station_lats = np.array([s.lat for s in stations]) # Build model-domain bathy/topo samples at horizontal nodes. This is used for # coastline on the depth slice and surface lines on X/Y slices. node_lons = [normalize_lon_360(n.lon) for n in nodes] node_lats = [n.lat for n in nodes] node_z_m = sample_bathy_points(node_lons, node_lats) node_x = np.array([n.model_x_km for n in nodes]) node_y = np.array([n.model_y_km for n in nodes]) node_lon = np.array(node_lons) node_lat = np.array(node_lats) node_surface_depth = -node_z_m / 1000.0 fig, axes = plt.subplots(1, 3, figsize=(18, 6), constrained_layout=True) # Depth/map slice with approximate coastline from z=0 contour. ax = axes[0] try: # Contour scattered node bathymetry/topography by interpolating to a regular grid. from scipy.interpolate import griddata grid_lon = np.linspace(float(np.nanmin(node_lon)), float(np.nanmax(node_lon)), 250) grid_lat = np.linspace(float(np.nanmin(node_lat)), float(np.nanmax(node_lat)), 250) glon, glat = np.meshgrid(grid_lon, grid_lat) gz = griddata((node_lon, node_lat), node_z_m, (glon, glat), method="linear") ax.contour(glon, glat, gz, levels=[0.0], linewidths=1.0) except Exception as exc: print(f"WARNING: could not draw coastline contour from bathy grid: {exc}") ax.scatter(station_lons, station_lats, marker="^", s=35, edgecolors="k", linewidths=0.3, label="Stations used") ax.scatter(ev_lons[depth_mask], ev_lats[depth_mask], s=18, edgecolors="k", linewidths=0.3, label=f"EQs at {depth_slice_km:.1f} km") ax.set_xlabel("Longitude (0-360)") ax.set_ylabel("Latitude") ax.set_title(f"Depth slice: {depth_slice_km:.1f} km") ax.set_xlim(float(np.nanmin(node_lon)) - 0.1, float(np.nanmax(node_lon)) + 0.1) ax.set_ylim(float(np.nanmin(node_lat)) - 0.1, float(np.nanmax(node_lat)) + 0.1) ax.grid(True, alpha=0.3) ax.legend(loc="best") # X slice: y versus depth, with surface line along nearest X column. ax = axes[1] x_surface_mask = np.abs(node_x - x_slice_km) <= x_tol order = np.argsort(node_y[x_surface_mask]) if np.any(x_surface_mask): ax.plot(node_y[x_surface_mask][order], node_surface_depth[x_surface_mask][order], linewidth=1.5, label="Bathy/topo") ax.scatter(ev_y[x_mask], ev_z[x_mask], s=18, edgecolors="k", linewidths=0.3, label="Synthetic EQs") ax.axhline(0.0, linewidth=0.8, linestyle="--") ax.set_xlabel("Model Y (km)") ax.set_ylabel("Depth below sea level (km)") ax.set_title(f"X slice: X={x_slice_km:.1f} km") ax.set_ylim(MAX_DEPTH_KM + 5.0, min(-2.0, float(np.nanmin(node_surface_depth[x_surface_mask])) - 1.0 if np.any(x_surface_mask) else -2.0)) ax.grid(True, alpha=0.3) ax.legend(loc="best") # Y slice: x versus depth, with surface line along nearest Y row. ax = axes[2] y_surface_mask = np.abs(node_y - y_slice_km) <= y_tol order = np.argsort(node_x[y_surface_mask]) if np.any(y_surface_mask): ax.plot(node_x[y_surface_mask][order], node_surface_depth[y_surface_mask][order], linewidth=1.5, label="Bathy/topo") ax.scatter(ev_x[y_mask], ev_z[y_mask], s=18, edgecolors="k", linewidths=0.3, label="Synthetic EQs") ax.axhline(0.0, linewidth=0.8, linestyle="--") ax.set_xlabel("Model X (km)") ax.set_ylabel("Depth below sea level (km)") ax.set_title(f"Y slice: Y={y_slice_km:.1f} km") ax.set_ylim(MAX_DEPTH_KM + 5.0, min(-2.0, float(np.nanmin(node_surface_depth[y_surface_mask])) - 1.0 if np.any(y_surface_mask) else -2.0)) ax.grid(True, alpha=0.3) ax.legend(loc="best") output_png.parent.mkdir(parents=True, exist_ok=True) fig.savefig(output_png, dpi=300) plt.close(fig) print(f"Saved 3D synthetic-event slice plot to: {output_png}") # ============================================================================= # MAIN # ============================================================================= def main() -> int: mode = MODE.strip().upper() if mode not in {"2D", "3D"}: raise ValueError("MODE must be '2D' or '3D'") if HORIZONTAL_INCREMENT_KM <= 0.0: raise ValueError("HORIZONTAL_INCREMENT_KM must be > 0") if VERTICAL_INCREMENT_KM <= 0.0: raise ValueError("VERTICAL_INCREMENT_KM must be > 0") if STA_MIN <= 0: raise ValueError("STA_MIN must be > 0") print("Synthetic HypoNet event generation") print(f" MODE: {mode}") print(f" Horizontal increment km: {HORIZONTAL_INCREMENT_KM}") print(f" Vertical increment km: {VERTICAL_INCREMENT_KM}") print(f" Max depth km: {MAX_DEPTH_KM}") print(f" STA_MIN: {STA_MIN}") print(f" MAX_OFFSET_KM: {MAX_OFFSET_KM}") print(f" Gaussian noise std s: {GAUSSIAN_NOISE_STD_S}") nodes = read_velocity_model_horizontal_nodes(VELOCITY_MODEL_FILE) if mode == "3D": candidate_events = generate_3d_candidate_events(nodes) else: candidate_events = generate_2d_candidate_events(nodes) filtered_events = apply_bathy_topo_filter(candidate_events) events = [ SyntheticEvent( event_index=i, model_x_km=ev.model_x_km, model_y_km=ev.model_y_km, depth_km=ev.depth_km, lon=ev.lon, lat=ev.lat, surface_z_m=ev.surface_z_m, source_mode=ev.source_mode, profile_distance_km=ev.profile_distance_km, ) for i, ev in enumerate(filtered_events) ] stations = read_simul_used_stations(STATION_FILE) stations = filter_stations_to_velocity_model_domain(stations, nodes) prepare_output_dir(OUTPUT_DIR) print("Loading HypoNet models...") model_p = HypoNet_DeeperPINN.model(inputdir=INPUT_P, device=DEVICE) model_s = HypoNet_DeeperPINN.model(inputdir=INPUT_S, device=DEVICE) if FORCE_CPU_AFTER_MODEL_LOAD: # Matches the pattern used in Calculate_all_arrival_times.py, because # inference_torch internally calls numpy in some HypoNet versions. model_p.device = torch.device("cpu") model_s.device = torch.device("cpu") rng = np.random.default_rng(RANDOM_SEED) summary_fields = [ "event_index", "event_file", "mode", "model_x_km", "model_y_km", "profile_distance_km", "lon", "lat", "depth_km", "surface_z_m", "surface_depth_km", "stations_selected", "stations_within_max_offset", "nearest_station_offset_km", "farthest_selected_station_offset_km", "n_p_arrivals", "n_s_arrivals", ] total_p = 0 total_s = 0 skipped_no_stations = 0 with SUMMARY_FILE.open("w", newline="", encoding="utf-8") as sf: writer = csv.DictWriter(sf, fieldnames=summary_fields) writer.writeheader() for ev in events: selected = select_stations_for_event(ev, stations) if not selected: skipped_no_stations += 1 continue phase_rows = calculate_event_times(ev, selected, model_p, model_s, rng) event_file = OUTPUT_DIR / f"event.{ev.event_index:08d}.txt" n_p, n_s = write_event_file(event_file, phase_rows) offsets = [d for _, d in selected] n_within = sum(1 for d in offsets if d <= MAX_OFFSET_KM) total_p += n_p total_s += n_s writer.writerow( { "event_index": ev.event_index, "event_file": event_file.name, "mode": ev.source_mode, "model_x_km": f"{ev.model_x_km:.3f}", "model_y_km": f"{ev.model_y_km:.3f}", "profile_distance_km": ( "" if ev.profile_distance_km is None else f"{ev.profile_distance_km:.3f}" ), "lon": f"{normalize_lon_for_hyponet(ev.lon):.6f}", "lat": f"{ev.lat:.6f}", "depth_km": f"{ev.depth_km:.3f}", "surface_z_m": f"{ev.surface_z_m:.3f}", "surface_depth_km": f"{-ev.surface_z_m / 1000.0:.3f}", "stations_selected": len(selected), "stations_within_max_offset": n_within, "nearest_station_offset_km": f"{min(offsets):.3f}", "farthest_selected_station_offset_km": f"{max(offsets):.3f}", "n_p_arrivals": n_p, "n_s_arrivals": n_s, } ) if PRINT_EVERY_N_EVENTS and (ev.event_index + 1) % PRINT_EVERY_N_EVENTS == 0: print(f" Wrote {ev.event_index + 1} / {len(events)} events") print("Done.") print(f" Synthetic events written: {len(events) - skipped_no_stations}") print(f" Skipped with no stations: {skipped_no_stations}") print(f" Total P arrivals: {total_p}") print(f" Total S arrivals: {total_s}") print(f" Event directory: {OUTPUT_DIR}") print(f" Summary file: {SUMMARY_FILE}") plot_synthetic_events(events, stations, nodes) return 0 if __name__ == "__main__": raise SystemExit(main())