简介面向地震数据分析初学者、地学相关课程作业及对USGS开放数据接口感兴趣的开发者该资源演示了从USGS地震数据平台获取JSON并转换为CSV的完整链路。核心思路是在终端用wget拉取USGS API返回的原始JSON数据再通过Jupyter Notebook完成字段筛选、精简与CSV导出适合希望快速上手地震数据接口调用和格式转换的用户。压缩包共4个文件涵盖ipynb、json、csv、md四种类型包括可直接运行的Notebook示例、原始地震JSON数据、精简后的CSV结果和配套说明文档整体仅1.44MB轻便易用。说明文档还提供了macOS下安装wget、调用查询接口时设置时间范围、最小震级等参数的方法并给出用JSONView检查返回字段的技巧能帮助读者复现从数据获取、清洗到落盘的完整流程也为后续可视化或统计分析打好基础。目前已有650人学习下载。1. 用 USGS earthquake 接口搭实时地震数据管道先跑通再谈分析接过不少地震数据需求之后我发现多数人卡在第一步USGS 官网那个实时地震列表看起来很丰富却不知道它背后是一条公开的 FDSN 查询接口——你给它起止时间和最小震级它就返回成批的全球 earthquake 事件数据而且是 GeoJSON能直接用。要做的事很明确把这条接口做成能每天依赖的数据管道参数怎么传、GeoJSON 怎么清洗成表、成图后又怎么识别假异常。适合要做地震可视化或常态化震情纪录的开发者也适合需要全球地震目录做统计、又不想反复手工导出 CSV 的研究人员。前半程看参数后半程看数据纪律踩坑的部分都写在靠后的章节里。2. 用 USGS earthquake API 的最小请求一条 URL 就是一个数据通道2.1 端点与返回格式为什么我选 GeoJSON 而不是 CSVUSGS earthquake 的 FDSN 事件查询端点是https://earthquake.usgs.gov/fdsnws/event/1/query你不需要爬任何动态页面也不需要登录态直接 HTTP GET 就能拿数据。format参数决定返回格式geojson、csv、xml、kml 都支持。我平时默认用 geojson原因有三个。第一properties 里带着url、detail、magType、status这些元信息CSV 会丢掉相当一部分第二几何数据一次给全geometry.coordinates就是[经度, 纬度, 深度]不用像 CSV 那样自己拼字段第三后续不管进 PostGIS 还是转成 shapefileGeoJSON 都是通用格式。如果只是临时看一眼数据CSV 完全够用但要做数据管道建议从一开始就把 GeoJSON 作为主格式。endpoint 本身不用换变的只是format参数这个切换成本几乎为零。2.2 最小可运行请求requests 传参而不是拼 URL我见过不少同事手工把 URL 拼出来参数一多就开始漏转义时间里的冒号和日期里的 T 尤其容易出问题。用 requests 的params参数可以让框架帮你做 URL 编码这也是一条更不容易翻车的习惯。import requests BASE_URL https://earthquake.usgs.gov/fdsnws/event/1/query def query_usgs(params: dict) - dict: resp requests.get(BASE_URL, paramsparams, timeout30) resp.raise_for_status() # 4xx/5xx 直接抛异常避免拿到空数据还继续跑 return resp.json() params { format: geojson, starttime: 2025-01-01T00:00:00Z, endtime: 2025-01-02T00:00:00Z, minmagnitude: 4.5, } data query_usgs(params) print(data[metadata][count]) # 先看事件总数再决定要不要继续处理requests.get的params参数会自动处理 URL 编码时间字符串里的冒号、日期里的T都不用手动操心。raise_for_status()这行容易被忽略但它很关键如果服务端返回 400 或 500直接抛异常能让你马上发现参数问题而不是拿着一个空字典排查半天。打印metadata.count是排查问题的第一步——先看数量对不对再谈解析和可视化。注意starttime和endtime都要求是 UTC 时间字符串以Z结尾。如果写成2025-01-01T00:00:00不带ZUSGS 后端按 UTC 解析和东八区本地时间会整整差 8 小时时间窗整体错位这是新手最容易踩的第一个坑。2.3 常用参数组速查把时间窗、震级、地理范围写进一次查询USGS earthquake 查询接口的参数不算多但组合起来能覆盖绝大多数场景。下面是我常用的参数表参数作用我常用的取值format返回格式geojson管道用csv临时看starttime/endtime事件发生时间窗口UTC 的 ISO8601带Zminmagnitude/maxmagnitude按震级过滤4.5 或 5.0 起步避免数据太杂eventtype事件类型earthquake过滤塌陷、爆破等非天然事件status复核状态reviewed只取人工复核过的事件latitude/longitude/maxradiuskm圆形地理范围关心某个点周边时用比如maxradiuskm2000minlatitude/maxlatitude/minlongitude/maxlongitude矩形地理范围做大区域切片时用注意跨 180 度经线时要拆开orderby排序time-asc取时间顺序magnitude-desc取最大震级优先limit/offset分页单次 10005000不要拉到上限updatedafter按数据更新时间过滤增量抓取的核心参数第 6 章展开几个常见的组合用法。抓“全球 4.5 级以上、最近一天”的事件starttime和endtime用滚动时间生成minmagnitude4.5就行。抓“某个点周边 100 公里”用latitude、longitude、maxradiuskm三元组省得自己算经纬度边界。抓一个大区域时我习惯用矩形边界但要注意东经为正、西经为负跨 180 度经线的区域需要拆成两个矩形分别请求否则会漏掉太平洋中间那条线附近的事件。相比之下圆形半径虽然没有矩形灵活但逻辑简单跨线问题也少容错率更高。limit的默认值是 200单次请求能设置的上限虽然很高但我不建议真拉到上限。一次拿几万条会让服务端扫描时间变长响应容易超时分页多跑几次更稳妥。综合下来我一般控制在 10005000 条不够就按时间窗切片把大窗口拆成多个小窗口分别请求。这个策略在后面避坑章节里还会再讲。3. 把 GeoJSON 拍平成 DataFrameUSGS earthquake 数据清洗三步3.1 第一次看返回结构先打印 keys 再动手拿到数据后第一件事不是跑模型而是看清结构。GeoJSON 的外层是一个FeatureCollection最顶上有四个 keytype、metadata、features、bbox。metadata里是请求时间和返回数量features才是真正的事件列表。print(data[type]) # FeatureCollection print(data[metadata][count]) # 事件总数 print(data[features][0][properties].keys()) # 第一个事件有哪些属性 print(data[features][0][geometry][coordinates]) # [经度, 纬度, 深度]properties里常见字段有mag、place、time、updated、url、detail、status、magType、tsunami等。geometry.coordinates是三个元素的列表分别是经度、纬度、深度。值得注意的一点是properties里的字段并不是每个事件都完整部分小震级事件可能缺mag人工复核前的事件status是automatic。所以后面解析时要用p.get(字段名)而不是p[字段名]否则一个缺失字段就能让整个管道崩掉。3.2 features 解析函数一行一个地震事件把features拍平成 DataFrame 是这个环节的核心。我一般写一个通用函数把嵌套字典压成扁平的表格同时顺手把时间戳转成可读的时间格式。import pandas as pd def geojson_to_frame(data: dict) - pd.DataFrame: rows [] for f in data[features]: p f[properties] lon, lat, depth f[geometry][coordinates] rows.append({ id: p.get(id), time: p.get(time), # epoch 毫秒 updated: p.get(updated), # 数据修订时间 mag: p.get(mag), mag_type: p.get(magType), place: p.get(place), depth_km: depth, lon: lon, lat: lat, status: p.get(status), url: p.get(url), }) df pd.DataFrame(rows) df[time] pd.to_datetime(df[time], unitms, utcTrue) df[updated] pd.to_datetime(df[updated], unitms, utcTrue) return df.sort_values(time, ascendingFalse).reset_index(dropTrue) df geojson_to_frame(data) print(df.head()) # 最快验证管道的方式这段代码的逻辑是遍历features把每个事件的属性和坐标拆成一行geometry.coordinates按顺序解包成经度、纬度、深度三个变量p.get()保证某个字段缺失时返回None而不是抛KeyError。时间字段从 epoch 毫秒转成datetime类型时unitms必须和utcTrue一起用前者说明输入精度后者避免 pandas 自作主张转成本地时区。最后按时间倒序排方便直接看最新事件。空响应的情况也要考虑如果metadata.count为 0rows是空列表pd.DataFrame(rows)建出来的表没有列后续取df[mag]会直接报错。所以调用这个函数之前先判断一下data[metadata][count]是否为 0是就直接返回空表别硬往下走。3.3 时间、深度、震级的三个数据纪律清洗 USGS earthquake 数据时有三条纪律我踩过不止一次。时间字段必须当 UTC 处理。time列已经是 UTC 的datetime类型生成日报告时不要用df[time].replace(hour...)这种原地替换应该用df[time].dt.tz_convert(Asia/Shanghai)转成目标时区。replace改的是墙钟时间而不是时区转换完的数值是错的。深度可能是负值。USGS 的depth_km默认向下为正但人工事件比如矿震、水库触发的地震偶尔会出现负数画图时会在海平面以上冒出几个点。处理方式很简单df[depth_km] df[depth_km].clip(lower0)负值归零或者单独打个标记不要直接删行。震级可能缺失或为 0。部分事件没有mag字段也有一部分小事件的震级是 0.0 或 None。统计频次、计算 b 值之前必须过滤掉这些记录否则画出来的震级分布图左端会莫名其妙地堆起一座山峰。比较全局的做法是请求阶段就加eventtypeearthquakestatusreviewed清洗阶段再补一道df df.dropna(subset[mag])两道防线更稳妥。4. 全球震中图与 b 值曲线USGS earthquake 数据能直接回答的两个问题4.1 用 plotly 画震中分布大小、深浅一次表达拿到干净的 DataFrame 之后最直观的输出就是全球震中分布图。我推荐 plotly 的scatter_geo它不需要本地底图文件一行代码就能出交互图鼠标悬停能看到事件详情。相比 matplotlib 配 Basemap 的老路少了很多环境依赖的麻烦。import plotly.express as px df_plot df.dropna(subset[mag, lat, lon]) fig px.scatter_geo( df_plot, lonlon, latlat, sizemag, colormag, hover_nameplace, hover_data{time: True, depth_km: True, mag: True}, projectionnatural earth, titleGlobal earthquakes, ) fig.update_traces(markerdict(sizemin1, sizeref2)) fig.show()sizemag让点的大小随震级变化colormag用连续色标表达震级强弱两者叠加后一眼就能看出哪里有大震。sizemin1保证小震级事件在地图上也能看到sizeref2控制点的大小比例数值越大点越小可以根据屏幕分辨率微调。hover_data里放time和depth_km鼠标悬停时不用再查表。如果部署环境没有外网或者不允许加载 plotly 的底图资源退回方案是 matplotlib 直接画经纬度散点图配合经纬度网格线也能看出空间分布只是少了交互能力和地理边界。遇到过内网环境只有 matplotlib 的读者用散点图加南海诸岛范围标注做一个简化版本完全够用。4.2 计算 Gutenberg-Richter b 值小样本也能估算地震学里有一个经典的经验关系Gutenberg-Richter 定律即震级大于等于 M 的事件数量服从log10(N) a - bM。b 值描述大小地震的比例关系通常在 0.81.2 之间。用 USGS earthquake 拉到的数据可以直接估算。import numpy as np import matplotlib.pyplot as plt mags df_plot[mag].dropna().values mc 4.5 # 最小完备震级要根据 minmagnitude 调 def estimate_b(mags, mc): mags mags[mags mc] # 只保留完备震级以上的事件 if len(mags) 20: return np.nan # 样本太少估出来没有统计意义 b np.log10(np.e) / (mags.mean() - mc) return b b_value estimate_b(mags, mc) print(fb {b_value:.2f})这里用的是极大似然估计公式是b log10(e) / (mean(M) - Mc)其中Mc是最小完备震级。为何不用线性回归拟合因为小震级端往往不完备线性拟合会把台网漏检的那段当成真实分布导致 b 值虚高。极大似然只对M Mc的样本求均值对截断更稳健。len(mags) 20的判断很重要少于 20 个样本的估计值可信度很低不如直接返回nan。画频次图的时候把累计频次画出来再看拟合效果M_range np.arange(mc, mags.max() 0.1, 0.1) N_values [(mags m).sum() for m in M_range] plt.plot(M_range, np.log10(N_values), o-) plt.xlabel(magnitude) plt.ylabel(log10(N)) plt.title(Gutenberg-Richter curve) plt.show()频次曲线左端如果明显向下弯说明这个震级段有事件没被检测到右端如果掉得太快说明样本太少。这两个特征都在提醒你先检查数据再谈物理解释。4.3 结果异常时先查数据再查算法跑出 b 值明显偏离 1 的时候我现在的第一反应不是怀疑公式而是回查数据。三个高频数据问题会直接污染结果小震级缺失导致 b 值虚高台网分布不均匀导致某一区域事件明显稀疏事件修订后同一事件在目录里出现两次计数虚增。处理的办法也简单。画图前按id去重df.drop_duplicates(subsetid)能解决修订导致的重复算 b 值前把minmagnitude往上抬一档比如从 4.5 抬到 5.0可以减少台网检测能力不足带来的偏差。空间分布如果明显呈团状不代表地震真的只发生在那些点更多时候只是台网覆盖不均匀。可视化可以做但解读的时候要把这个前提挂在嘴边。5. USGS earthquake 请求的避坑清单五个现场级故障5.1 请求返回空数据count 为 0现象代码看不出问题参数也传了但metadata.count是 0features是空列表。查下来最常见的原因是时间参数出了问题。有人把endtime写成本地时间2025-01-02 00:00:00中间用空格而不是TUSGS 解析失败也有人把starttime和endtime写反了或者结束时间是未来的某个时间点。解决统一用 Python 生成时间参数不要手工敲。datetime.now(timezone.utc)拿到的是 UTC 当前时间strftime(%Y-%m-%dT%H:%M:%S)拼上Z后缀格式和时区都错不了。如果非要在浏览器里先验证参数也先确认 URL 里时间后面带Z再复制到代码里。5.2 一拉长历史时间窗就超时或 504现象想抓 1990 年以来的全球 4.5 级以上事件一次性把参数写成一个 35 年的时间窗limit拉满结果请求直接 504 网关超时或者卡到 timeout 被requests主动掐断。原因USGS 后端要在海量事件里扫描、排序、组装 GeoJSON时间跨度越长响应越慢。把limit拉满是一个误区服务端不是不能处理而是单次请求的耗时不可控。解决按时间窗切片一年一个请求或者半年一个请求。单次limit控制在 10005000 条请求之间time.sleep(0.5)以上给服务端留出余量。失败后重试 3 次每次等待时间递增。把下面这个循环模板用起来for year in range(1990, 2025): params[starttime] f{year}-01-01T00:00:00Z params[endtime] f{year}-12-31T23:59:59Z try: df geojson_to_frame(query_usgs(params)) time.sleep(0.8) except requests.HTTPError: time.sleep(5) # 等 5 秒后重试超时通常一两次就缓过来5.3 地图上出现大量震级为 0 或 None 的点现象画出来的地图边缘区域堆着一批小点hover 一看震级是 0.0 或者显示为 None密集程度明显不符合常理。原因接口默认返回所有事件类型包括爆破、塌陷、矿震等非天然事件这些事件不一定有震级。另外未复核的自动事件里也有一部分mag字段为空。直接把这些数据画图自然会出现大量异常点。解决请求阶段加两个参数eventtypeearthquake过滤非天然事件statusreviewed只保留人工复核过的记录。清洗阶段再加一道dropna(subset[mag])和df[mag] 0的过滤。两层防线下来图面会干净很多。5.4 orderbymagnitude 并没有拿到最大地震现象想取“当前时间段最大震级事件”写了orderbymagnitude结果返回列表第一条是 M2.8而不是 M6.5。原因USGS 的orderby参数比较反直觉orderbymagnitude实际是升序震级从小到大排。需要显式指定方向才能得到预期结果。解决取最大震级优先写orderbymagnitude-desc如果要按时间从新到旧排写orderbytime-desc。这个细节我建议直接写死在代码注释里不然过两个月回来看自己的脚本还会栽一次。5.5 增量脚本每天少数据现象每天跑一次增量脚本抓取前一天的新事件某天突然发现遗漏了一个 M5.5 的事件日志里也没有报错。原因USGS 会对已发布的事件做修订震级、位置、深度都可能被更新。如果只按事件发生时间窗口抓取修订发生在窗口之外的事件就不会被拉回来增量自然就漏了。这不是请求失败而是查询条件本身有盲区。解决增量脚本不要只按starttime过滤要结合updatedafter参数把“这个时间之后有更新”的事件也一并拉回合并时按id去重用updated字段判断新旧。完整实现放在最后一章。6. 让抓取脚本跑成习惯增量更新与每日校验增量更新是数据管道从“能跑”变“能依赖”的关键一步。做法很简单把updatedafter设为上次运行时间starttime往前回溯 7 天这样既覆盖了窗口内新发生的事件也覆盖了窗口外被修订过的事件。from datetime import datetime, timezone, timedelta from pathlib import Path def daily_incremental(since: datetime, save_path: str) - pd.DataFrame: params { format: geojson, starttime: (since - timedelta(days7)).strftime(%Y-%m-%dT%H:%M:%S), endtime: datetime.now(timezone.utc).strftime(%Y-%m-%dT%H:%M:%S), updatedafter: since.strftime(%Y-%m-%dT%H:%M:%S), eventtype: earthquake, status: reviewed, orderby: time-asc, limit: 2000, } df_new geojson_to_frame(query_usgs(params)) if Path(save_path).exists(): df_old pd.read_csv(save_path, parse_dates[time, updated]) df_merged pd.concat([df_old, df_new], ignore_indexTrue) df_merged df_merged.sort_values(updated).drop_duplicates(subsetid, keeplast) else: df_merged df_new df_merged.to_csv(save_path, indexFalse) return df_mergedstarttime往前回溯 7 天是为了兜住“时间窗之外、更新落在窗内”的修订事件代价是多拉一部分重复数据靠drop_duplicates清洗掉。合并前先按updated排序再按id去重并保留最后一条这样旧版本会被新版本覆盖。keeplast这个参数必须配合排序使用否则去重后保留哪一行是不确定的。每日校验我习惯用一个很轻的方式从合并后的数据里取出昨天震级最大的一条打印place、mag、time三个字段扫一眼是否合理。管道如果出了问题最先暴露的往往是这三个字段对不上。这个习惯帮我提前发现过两次数据源字段变更和一次本地区域性明显偏差省下的排查时间远超写这个脚本的投入。希望这个思路对你的定时任务也有用希望帮到你。本文还有配套的精品资源点击获取