fuzzyworld

shadow-zone / files / scripts / fetch-earthquake-data.mjs

fetch-earthquake-data.mjs

11.3 KB · history · blame

1#!/usr/bin/env node
2/**
3 * Fetch real seismic phase arrival data from public sources and write formatted
4 * JSON into src/data/ for use by the validation module.
5 *
6 * Data source: USGS ComCat (QuakeML product)
7 * 1. GET /fdsnws/event/1/query?format=geojson&eventid=<id>
8 * → find products['phase-data'][0].contents['quakeml.xml'].url
9 * 2. Parse QuakeML: match <pick> publicIDs with <arrival> pickIDs to get
10 * station, epicentral distance (degrees), and observed travel time.
11 *
12 * Usage:
13 * node scripts/fetch-earthquake-data.mjs
14 * node scripts/fetch-earthquake-data.mjs --event tohoku-2011
15 * node scripts/fetch-earthquake-data.mjs --max-stations 300
16 *
17 * Writes: src/data/{id}.json
18 */
19
20import { writeFileSync, mkdirSync } from 'fs';
21import { dirname, join } from 'path';
22import { fileURLToPath } from 'url';
23
24const __dirname = dirname(fileURLToPath(import.meta.url));
25const OUT_DIR = join(__dirname, '..', 'src', 'data');
26
27mkdirSync(OUT_DIR, { recursive: true });
28
29// ── CLI args ──────────────────────────────────────────────────────────────────
30const args = process.argv.slice(2);
31const filterEvent = argVal(args, '--event');
32const maxStations = parseInt(argVal(args, '--max-stations') ?? '500', 10);
33
34function argVal(arr, flag) {
35 const i = arr.indexOf(flag);
36 return i !== -1 ? arr[i + 1] : null;
37}
38
39// ── Event registry ────────────────────────────────────────────────────────────
40const EVENTS = [
41 {
42 id: 'tohoku-2011',
43 name: '2011 Tōhoku, Japan (Mw 9.1)',
44 originTime: '2011-03-11T05:46:24Z',
45 latitude: 38.297,
46 longitude: 142.373,
47 depthKm: 29,
48 magnitudeMw: 9.1,
49 usgsId: 'official20110311054624120_30',
50 note: 'Largest recorded earthquake in Japan. Shallow source, excellent global coverage. Clear P shadow zone near 105°–140°.',
51 },
52 {
53 id: 'chile-2010',
54 name: '2010 Maule, Chile (Mw 8.8)',
55 originTime: '2010-02-27T06:34:11Z',
56 latitude: -35.909,
57 longitude: -72.733,
58 depthKm: 22,
59 magnitudeMw: 8.8,
60 usgsId: 'official20100227063411530_30',
61 note: 'Megathrust earthquake, shallow source. Complementary geometry to Tōhoku — stations in shadow zone are different parts of the global network.',
62 },
63 {
64 id: 'okhotsk-2013',
65 name: '2013 Sea of Okhotsk (Mw 8.3, depth 608 km)',
66 originTime: '2013-05-24T05:44:49Z',
67 latitude: 54.874,
68 longitude: 153.281,
69 depthKm: 608,
70 magnitudeMw: 8.3,
71 usgsId: 'usb000h4jh',
72 note: 'Very deep focus earthquake. Rays launch from 608 km depth, altering shadow zone geometry. Good for demonstrating that depth matters in the model.',
73 },
74];
75
76// ── Main ──────────────────────────────────────────────────────────────────────
77const toFetch = filterEvent ? EVENTS.filter(e => e.id === filterEvent) : EVENTS;
78
79if (toFetch.length === 0) {
80 console.error(`Unknown event id: ${filterEvent}`);
81 console.error('Available:', EVENTS.map(e => e.id).join(', '));
82 process.exit(1);
83}
84
85for (const event of toFetch) {
86 console.log(`\n── ${event.id} ──────────────────────────────────`);
87 try {
88 const arrivals = await fetchArrivals(event, maxStations);
89 const dataset = { ...event, arrivals };
90 const outPath = join(OUT_DIR, `${event.id}.json`);
91 writeFileSync(outPath, JSON.stringify(dataset, null, 2));
92 console.log(`✅ Wrote ${arrivals.length} arrivals → ${outPath}`);
93 if (arrivals.length === 0) {
94 console.warn(' ⚠️ No arrivals found — check the usgsId or run with --debug');
95 }
96 } catch (err) {
97 console.error(`❌ Failed: ${err.message}`);
98 if (process.env.DEBUG) console.error(err.stack);
99 }
100}
101
102// ── USGS ComCat + QuakeML parser ──────────────────────────────────────────────
103
104async function fetchArrivals(event, maxStations) {
105 // Step 1: get the QuakeML URL from the event detail
106 const detailUrl =
107 `https://earthquake.usgs.gov/fdsnws/event/1/query?format=geojson&eventid=${event.usgsId}`;
108 console.log(` Fetching event detail from USGS…`);
109 const detail = await fetchJSON(detailUrl);
110
111 // Search all phase-data product variants for a quakeml.xml content.
112 // 'trump-phase-data' is sometimes present but with empty contents, so we
113 // scan all candidates and take the first that has a real quakeml.xml.
114 const allProducts = [
115 ...(detail?.properties?.products?.['phase-data'] ?? []),
116 ...(detail?.properties?.products?.['trump-phase-data'] ?? []),
117 ];
118
119 if (allProducts.length === 0) {
120 throw new Error('No phase-data product found in ComCat response');
121 }
122
123 let quakemlUrl = null;
124 for (const product of allProducts) {
125 const contents = product?.contents;
126 if (!contents || typeof contents !== 'object' || Array.isArray(contents)) continue;
127 const qml = Object.entries(contents).find(([k]) => k.endsWith('quakeml.xml'));
128 if (qml) {
129 quakemlUrl = qml[1].url;
130 break;
131 }
132 }
133 if (!quakemlUrl) {
134 throw new Error('quakeml.xml not found in any phase-data product');
135 }
136
137 console.log(` Downloading QuakeML (${quakemlUrl.split('/').slice(-3).join('/')})…`);
138 const xml = await fetchText(quakemlUrl);
139 console.log(` Parsing (${(xml.length / 1024).toFixed(0)} KB)…`);
140
141 return parseQuakeML(xml, event.originTime, maxStations);
142}
143
144/**
145 * Parse USGS QuakeML to extract P/S travel times.
146 *
147 * QuakeML structure:
148 * <pick publicID=".../<sta>/<phase>">
149 * <time><value>ISO8601</value></time>
150 * <waveformID networkCode="IU" stationCode="ANMO" />
151 * <phaseHint>P</phaseHint>
152 * </pick>
153 *
154 * <arrival publicID="...">
155 * <pickID>.../<sta>/<phase></pickID> ← matches pick publicID
156 * <phase>P</phase>
157 * <distance>28.14</distance> ← degrees
158 * </arrival>
159 *
160 * Travel time = pick.time − origin.time.
161 */
162function parseQuakeML(xml, fallbackOriginTime, maxStations) {
163 // Extract origin time (prefer the first reviewed origin)
164 const originTimeStr = extractOriginTime(xml) ?? fallbackOriginTime;
165 const originMs = new Date(originTimeStr).getTime();
166 if (isNaN(originMs)) throw new Error(`Bad origin time: ${originTimeStr}`);
167 console.log(` Origin time: ${originTimeStr}`);
168
169 // Build pick map: publicID → { time, network, station, phaseHint }
170 const picks = new Map();
171 for (const block of iterXmlElements(xml, 'pick')) {
172 const pid = extractAttr(block, 'publicID');
173 // <time><value>ISO8601</value></time> — nested, so use a two-level match
174 const timeStr = extractNestedTag(block, 'time', 'value');
175 const network = extractAttr(block, 'networkCode') ?? '';
176 const station = extractAttr(block, 'stationCode') ?? '';
177 const phase = extractTag(block, 'phaseHint') ?? '';
178
179 if (!pid || !timeStr) continue;
180 const tMs = new Date(timeStr).getTime();
181 if (isNaN(tMs)) continue;
182
183 picks.set(pid, { tMs, network, station, phase });
184 }
185 console.log(` Parsed ${picks.size} picks`);
186
187 // Process arrivals — each arrival references a pick and adds distance
188 const stationMap = new Map(); // key: net.sta → record
189
190 for (const block of iterXmlElements(xml, 'arrival')) {
191 const pickId = extractTag(block, 'pickID');
192 const phase = extractTag(block, 'phase') ?? '';
193 const distStr = extractTag(block, 'distance');
194
195 if (!pickId || !distStr) continue;
196
197 const distDeg = parseFloat(distStr);
198 if (isNaN(distDeg) || distDeg <= 0 || distDeg > 180) continue;
199
200 const pick = picks.get(pickId);
201 if (!pick) continue;
202
203 // Accept direct P/S phases only (not PP, PKP, PcP, etc.)
204 const isP = /^P[ng]?$/.test(phase);
205 const isS = /^S[ng]?$/.test(phase);
206 if (!isP && !isS) continue;
207
208 const travelSec = (pick.tMs - originMs) / 1000;
209 if (travelSec <= 0 || travelSec > 4000) continue;
210
211 const key = `${pick.network}.${pick.station}`;
212 if (!stationMap.has(key)) {
213 stationMap.set(key, {
214 station: pick.station,
215 network: pick.network,
216 distanceDeg: parseFloat(distDeg.toFixed(2)),
217 pTravelTimeSec: null,
218 sTravelTimeSec: null,
219 });
220 }
221
222 const rec = stationMap.get(key);
223 if (isP && (rec.pTravelTimeSec === null || travelSec < rec.pTravelTimeSec)) {
224 rec.pTravelTimeSec = Math.round(travelSec * 10) / 10;
225 }
226 if (isS && (rec.sTravelTimeSec === null || travelSec < rec.sTravelTimeSec)) {
227 rec.sTravelTimeSec = Math.round(travelSec * 10) / 10;
228 }
229 }
230
231 const result = [...stationMap.values()]
232 .filter(r => r.pTravelTimeSec !== null || r.sTravelTimeSec !== null)
233 .sort((a, b) => a.distanceDeg - b.distanceDeg)
234 .slice(0, maxStations);
235
236 const pCount = result.filter(r => r.pTravelTimeSec !== null).length;
237 const sCount = result.filter(r => r.sTravelTimeSec !== null).length;
238 console.log(` P: ${pCount}, S: ${sCount}, total stations: ${result.length}`);
239 return result;
240}
241
242// ── Minimal XML helpers (no library dependency) ───────────────────────────────
243
244/** Iterate over all <tagName ...>...</tagName> blocks in XML text. */
245function* iterXmlElements(xml, tagName) {
246 const open = new RegExp(`<${tagName}[\\s>]`, 'g');
247 const close = `</${tagName}>`;
248 let match;
249 while ((match = open.exec(xml)) !== null) {
250 // Rewind 1 to include the '<'
251 const start = match.index;
252 const end = xml.indexOf(close, start);
253 if (end === -1) continue;
254 yield xml.slice(start, end + close.length);
255 }
256}
257
258/**
259 * Extract text of the first <outer><inner>TEXT</inner></outer> nested pair.
260 * Handles the common QuakeML pattern where values are doubly-wrapped.
261 */
262function extractNestedTag(xml, outer, inner) {
263 const re = new RegExp(`<${outer}[^>]*>\\s*<${inner}[^>]*>([^<]+)</${inner}>`);
264 const m = xml.match(re);
265 return m ? m[1].trim() : null;
266}
267
268/** Extract text of the first <tag>TEXT</tag> (no nested tags inside). */
269function extractTag(xml, tag) {
270 const m = xml.match(new RegExp(`<${tag}[^>]*>([^<]+)</${tag}>`));
271 return m ? m[1].trim() : null;
272}
273
274/** Extract a named attribute value from anywhere in the string. */
275function extractAttr(xml, attr) {
276 const m = xml.match(new RegExp(`\\b${attr}="([^"]*)"`));
277 return m ? m[1] : null;
278}
279
280/** Extract the preferred origin time from the first <origin> block. */
281function extractOriginTime(xml) {
282 for (const block of iterXmlElements(xml, 'origin')) {
283 const t = extractNestedTag(block, 'time', 'value');
284 if (t) return t;
285 }
286 return null;
287}
288
289// ── HTTP helpers ──────────────────────────────────────────────────────────────
290
291async function fetchJSON(url) {
292 const res = await fetch(url, { headers: { 'User-Agent': 'shadow-zone-edu/1.0' } });
293 if (!res.ok) throw new Error(`HTTP ${res.status} for ${url}`);
294 return res.json();
295}
296
297async function fetchText(url) {
298 const res = await fetch(url, { headers: { 'User-Agent': 'shadow-zone-edu/1.0' } });
299 if (!res.ok) throw new Error(`HTTP ${res.status} for ${url}`);
300 return res.text();
301}