☰
实时地震数据管道实战:USGS接口、GeoJSON清洗与增量更新
2026/9/29 17:04:23 网站建设 项目流程

简介:面向地震数据分析初学者、地学相关课程作业及对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 而不是 CSV

USGS 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 还是转成 shapefile,GeoJSON 都是通用格式。如果只是临时看一眼数据,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, params=params, timeout=30) 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不带Z,USGS 后端按 UTC 解析,和东八区本地时间会整整差 8 小时,时间窗整体错位,这是新手最容易踩的第一个坑。

2.3 常用参数组速查:把时间窗、震级、地理范围写进一次查询

USGS earthquake 查询接口的参数不算多,但组合起来能覆盖绝大多数场景。下面是我常用的参数表:

参数作用我常用的取值
format返回格式geojson(管道用),csv(临时看)
starttime/endtime事件发生时间窗口UTC 的 ISO8601,带Z
minmagnitude/maxmagnitude按震级过滤4.5 或 5.0 起步,避免数据太杂
eventtype事件类型earthquake,过滤塌陷、爆破等非天然事件
status复核状态reviewed,只取人工复核过的事件
latitude/longitude/maxradiuskm圆形地理范围关心某个点周边时用,比如maxradiuskm=2000
minlatitude/maxlatitude/minlongitude/maxlongitude矩形地理范围做大区域切片时用,注意跨 180 度经线时要拆开
orderby排序time-asc取时间顺序,magnitude-desc取最大震级优先
limit/offset分页单次 1000~5000,不要拉到上限
updatedafter按数据更新时间过滤增量抓取的核心参数,第 6 章展开

几个常见的组合用法。抓“全球 4.5 级以上、最近一天”的事件,starttime和endtime用滚动时间生成,minmagnitude=4.5就行。抓“某个点周边 100 公里”,用latitude、longitude、maxradiuskm三元组,省得自己算经纬度边界。抓一个大区域时我习惯用矩形边界,但要注意东经为正、西经为负,跨 180 度经线的区域需要拆成两个矩形分别请求,否则会漏掉太平洋中间那条线附近的事件。相比之下圆形半径虽然没有矩形灵活,但逻辑简单,跨线问题也少,容错率更高。

limit的默认值是 200,单次请求能设置的上限虽然很高,但我不建议真拉到上限。一次拿几万条会让服务端扫描时间变长,响应容易超时,分页多跑几次更稳妥。综合下来我一般控制在 1000~5000 条,不够就按时间窗切片,把大窗口拆成多个小窗口分别请求。这个策略在后面避坑章节里还会再讲。

3. 把 GeoJSON 拍平成 DataFrame:USGS earthquake 数据清洗三步

3.1 第一次看返回结构:先打印 keys 再动手

拿到数据后第一件事不是跑模型,而是看清结构。GeoJSON 的外层是一个FeatureCollection,最顶上有四个 key:type、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"], unit="ms", utc=True) df["updated"] = pd.to_datetime(df["updated"], unit="ms", utc=True) return df.sort_values("time", ascending=False).reset_index(drop=True) df = geojson_to_frame(data) print(df.head()) # 最快验证管道的方式

这段代码的逻辑是:遍历features,把每个事件的属性和坐标拆成一行;geometry.coordinates按顺序解包成经度、纬度、深度三个变量;p.get()保证某个字段缺失时返回None而不是抛KeyError。时间字段从 epoch 毫秒转成datetime类型时,unit="ms"必须和utc=True一起用,前者说明输入精度,后者避免 pandas 自作主张转成本地时区。最后按时间倒序排,方便直接看最新事件。

空响应的情况也要考虑:如果metadata.count为 0,rows是空列表,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(lower=0),负值归零,或者单独打个标记,不要直接删行。

震级可能缺失或为 0。部分事件没有mag字段,也有一部分小事件的震级是 0.0 或 None。统计频次、计算 b 值之前必须过滤掉这些记录,否则画出来的震级分布图左端会莫名其妙地堆起一座山峰。比较全局的做法是请求阶段就加eventtype=earthquake&status=reviewed,清洗阶段再补一道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, lon="lon", lat="lat", size="mag", color="mag", hover_name="place", hover_data={"time": True, "depth_km": True, "mag": True}, projection="natural earth", title="Global earthquakes", ) fig.update_traces(marker=dict(sizemin=1, sizeref=2)) fig.show()

size="mag"让点的大小随震级变化,color="mag"用连续色标表达震级强弱,两者叠加后一眼就能看出哪里有大震。sizemin=1保证小震级事件在地图上也能看到,sizeref=2控制点的大小比例,数值越大点越小,可以根据屏幕分辨率微调。hover_data里放time和depth_km,鼠标悬停时不用再查表。

如果部署环境没有外网或者不允许加载 plotly 的底图资源,退回方案是 matplotlib 直接画经纬度散点图,配合经纬度网格线,也能看出空间分布,只是少了交互能力和地理边界。遇到过内网环境只有 matplotlib 的读者,用散点图加南海诸岛范围标注做一个简化版本,完全够用。

4.2 计算 Gutenberg-Richter b 值:小样本也能估算

地震学里有一个经典的经验关系:Gutenberg-Richter 定律,即震级大于等于 M 的事件数量服从log10(N) = a - bM。b 值描述大小地震的比例关系,通常在 0.8~1.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(f"b = {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(subset="id")能解决修订导致的重复;算 b 值前把minmagnitude往上抬一档,比如从 4.5 抬到 5.0,可以减少台网检测能力不足带来的偏差。空间分布如果明显呈团状,不代表地震真的只发生在那些点,更多时候只是台网覆盖不均匀。可视化可以做,但解读的时候要把这个前提挂在嘴边。

5. USGS earthquake 请求的避坑清单:五个现场级故障

5.1 请求返回空数据,count 为 0

现象:代码看不出问题,参数也传了,但metadata.count是 0,features是空列表。查下来最常见的原因是时间参数出了问题。有人把endtime写成本地时间2025-01-02 00:00:00,中间用空格而不是T,USGS 解析失败;也有人把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控制在 1000~5000 条,请求之间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字段为空。直接把这些数据画图,自然会出现大量异常点。

解决:请求阶段加两个参数,eventtype=earthquake过滤非天然事件,status=reviewed只保留人工复核过的记录。清洗阶段再加一道dropna(subset=["mag"])和df["mag"] > 0的过滤。两层防线下来,图面会干净很多。

5.4 orderby=magnitude 并没有拿到最大地震

现象:想取“当前时间段最大震级事件”,写了orderby=magnitude,结果返回列表第一条是 M2.8,而不是 M6.5。

原因:USGS 的orderby参数比较反直觉,orderby=magnitude实际是升序,震级从小到大排。需要显式指定方向,才能得到预期结果。

解决:取最大震级优先,写orderby=magnitude-desc;如果要按时间从新到旧排,写orderby=time-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(days=7)).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_index=True) df_merged = df_merged.sort_values("updated").drop_duplicates(subset="id", keep="last") else: df_merged = df_new df_merged.to_csv(save_path, index=False) return df_merged

starttime往前回溯 7 天是为了兜住“时间窗之外、更新落在窗内”的修订事件,代价是多拉一部分重复数据,靠drop_duplicates清洗掉。合并前先按updated排序,再按id去重并保留最后一条,这样旧版本会被新版本覆盖。keep="last"这个参数必须配合排序使用,否则去重后保留哪一行是不确定的。

每日校验我习惯用一个很轻的方式:从合并后的数据里取出昨天震级最大的一条,打印place、mag、time三个字段,扫一眼是否合理。管道如果出了问题,最先暴露的往往是这三个字段对不上。这个习惯帮我提前发现过两次数据源字段变更和一次本地区域性明显偏差,省下的排查时间远超写这个脚本的投入。希望这个思路对你的定时任务也有用,希望帮到你。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询