这次我们来看一个很有意思的模型。NASA 最新发布的 3D 模型，不是某个 AI 绘画工具，也不是本地视频生成模型，而是一个把地球真实形状“变形放大”之后的三维可视化：地球表面并不是完美球体，而是一个凹凸不平、看起来像“土豆”的大块头。因为这个特征，外媒几乎都用同一个词来形容它——Lumpy Mess。
这个模型的本质是“大地水准面”可视化，英文叫 Geoid。它是从 GRACE、GOCE 等卫星重力观测数据反演出来的一张全球等势面，代表地球重力场在“平均海面”意义下的真实起伏。NASA 把它做成 3D 模型后，地球看起来不再是光滑的椭球，而是一个有隆起、有凹陷、低处坑洼、高处鼓包的不规则表面。
这篇文章会围绕这个模型做一次技术拆解：数据从哪来、用什么格式、怎么在本地读取、如何在 Python 里生成 3D 可视化、能做哪些验证和批处理、以及实际工作时最容易踩的坑。如果你平时处理遥感数据、做地理可视化，或者想研究卫星重力数据，这篇可以直接收藏。
1. 核心能力速览
先把这个模型的关键信息整理成一张表。需要说明的是，NASA 发布页面的具体版本、分辨率会随时间更新，下面的描述以常见公开版本为基础，具体参数请以数据源页面为准。
| 能力项 | 说明 |
|---|---|
| 模型类型 | 地球重力场 / 大地水准面三维模型 |
| 核心数据来源 | GRACE、GRACE-FO、GOCE 等卫星重力观测，结合地面和海洋测量数据 |
| 主要功能 | 展示大地水准面起伏、全球重力异常、地球三维不规则形状 |
| 数据格式 | netCDF、GeoTIFF、txt / csv xyz 网格，部分版本提供 OBJ/GLTF 网格 |
| 常见分辨率 | 1 度、0.5 度、0.25 度、0.1 度，视发布版本而定 |
| 推荐处理工具 | Python、GMT、QGIS、ParaView、PyVista |
| 三维可视化方式 | 将高度或重力异常值做垂直夸张后，映射到椭球面上 |
| 是否支持批量任务 | 支持，可批量处理多时相重力场数据并合成动画 |
| 是否支持接口调用 | 可自行封装 REST API，官方主要以数据下载为主 |
| 适合场景 | 地球科学展示、测绘基准、海洋面研究、卫星轨道计算、地理教学 |
从材料看，这个项目不是“双击启动”的本地软件，而是一整套公开数据和处理流程。它的价值在于，你不需要自己发射卫星，也能用公开重力场模型在本地重现 NASA 的可视化效果。
2. 适用场景与使用边界
这类模型适合谁？
第一类是做地球物理和大地测量的人。重力场模型是确定全球高程基准、研究地壳结构、反演地下密度分布的基础数据。第二类是搞遥感影像或 GIS 开发的人。Geoid 数据可以用来做高程转换，比如 GPS 测到的是椭球高，真实海拔需要减去大地水准面差距。第三类是科普和技术展示。NASA 这种“土豆地球”可视化适合做互动界面、教学课件或数据大屏。
但它也有非常明确的使用边界。最重要的边界是：这个凹凸不平的表面不是真实地形。
真实地球表面有山脉、海沟、大陆，那是地形学概念。而 Geoid 是一个重力等势面，可以理解为“全球静止海平面穿越大陆延伸后形成的形状”。它和参考椭球之间的差异通常在正负一百米左右，肉眼根本看不出来。NASA 的 3D 模型为了让观众“看见”这种差异，会把高度放大几百到几千倍。所以你看到的是夸张后的重力场示意，不是地球真实比例。
另一个边界是数据合规。NASA 公开的重力场模型大多是开放共享的，但它依赖的局部精细数据、商业数据可能受许可约束。使用时要注意数据源页面给出的引用规范和授权范围。如果涉及本国高精度局部重力数据，更要遵守数据管理政策，不能随意传播或使用。
3. 关键技术原理：卫星重力数据是怎么变成“土豆地球”的
先把原理讲清楚，后面的代码才有意义。
地球的重力场可以用一组球谐系数来描述，类似傅里叶展开。球谐展开的阶数越高，能表达的空间细节越丰富。低阶系数描述地球的整体扁率和大尺度质量分布，高阶系数刻画山脉、海沟、密度异常引起的重力变化。
GRACE 和 GRACE-FO 卫星主要是通过测量两颗卫星之间的距离变化，反演地球重力场的长期变化，特别适合监测地表水、冰盖、地下水等质量迁移。GOCE 卫星则通过测量重力梯度，得到更高分辨率的静态重力场。把这些数据融合，再加上地面重力测量、海洋测高数据，就能构建出全球大地水准面模型。
模型里面真正展示的表面，通常是一个规则经纬度网格。每个网格点上有三个关键数值：
- 经度 longitude
- 纬度 latitude
- 大地水准面高度 geoid height，单位一般是米
Geoid height 表示相对于参考椭球面的起伏。正值表示该处重力等势面高于参考椭球，负值表示低于参考椭球。在 3D 可视化时，把每个经纬度点放到对应椭球位置，再沿径向向外或向内偏移一个高度值，同时乘以一个夸张系数，就得到了我们看到的高低不平地球。
NASA 的新模型之所以叫“Lumpy Mess”，就是因为它把这种起伏显示得非常明显。从视觉上看，赤道附近会鼓出来，南亚、东南亚有较明显的隆起，印度洋部分区域有凹陷，整体看起来确实不太“圆”。
4. 环境准备与数据获取
在本地复现这套可视化，不需要很夸张的硬件。CPU 就能跑，内存和显卡取决于数据分辨率和渲染方式。如果只处理 1 度网格，4GB 内存足够；如果处理 0.1 度全球网格，内存需求会明显上升。
建议环境如下：
- 操作系统：Windows 10/11、Ubuntu 20.04 以上、macOS 均可
- Python 版本：3.9 以上
- 核心库：numpy、matplotlib、xarray、netCDF4、pyvista、cartopy
- 可选工具：GMT、QGIS、ParaView
- 磁盘空间：至少预留 10GB，用于存放模型文件和处理结果
数据获取推荐几个渠道：
- NASA Earthdata：需要注册账号，下载 GRACE/GRACE-FO 数据
- ICGEM：International Centre for Global Earth Models，可直接在线计算 Geoid 格网
- GFZ、CSR、JPL：三家机构发布不同版本的 GRACE 时变重力场
如果你只想快速验证，可以找一个公开的 Geoid 模型文件，比如 EGM2008、GOCO 系列、XGM2019 的网格数据。文件扩展名常见为.nc、.grd、.txt、.xyz。
安装 Python 依赖：
pip install numpy matplotlib xarray netCDF4 pyvista cartopy如果读取 netCDF 时遇到时间相关依赖缺失，可以补装：
pip install cftime5. 本地数据处理与三维可视化
下面给出一套通用的处理和可视化流程。实际数据字段名、单位、坐标变量可能因文件而异，代码里我会标注需要替换的位置。
5.1 读取 Geoid 格网数据
假设你拿到的是 netCDF 格式，用 xarray 打开：
import xarray as xr import numpy as np # 打开 Geoid 数据，实际路径需要替换 ds = xr.open_dataset("geoid_model.nc") print(ds) # 取出经纬度和高度数组 # 常见变量名是 "geoid"、"h"、"z"，建议先打印变量列表确认 lat = np.array(ds["lat"]) lon = np.array(ds["lon"]) geoid = np.array(ds["geoid"])如果你的数据是文本网格，也可以用 pandas 读取：
import pandas as pd df = pd.read_csv("geoid_xyz.txt", delim_whitespace=True, header=None, names=["lon", "lat", "geoid"]) lon = df["lon"].to_numpy() lat = df["lat"].to_numpy() geoid = df["geoid"].to_numpy()这里最容易出错的是经纬度顺序。很多地理数据文件遵循“经度、纬度、值”，但也有的文件是“纬度、经度、值”。处理前务必确认。
5.2 生成二维平面图
先画一张二维伪彩图，确认数据范围是否正常：
import matplotlib.pyplot as plt import cartopy.crs as ccrs fig = plt.figure(figsize=(12, 6)) ax = plt.axes(projection=ccrs.PlateCarree()) pcm = ax.pcolormesh(lon, lat, geoid, cmap="RdBu_r", shading="auto") ax.coastlines(linewidth=0.5) ax.set_title("Geoid Height (m)") plt.colorbar(pcm, ax=ax, orientation="horizontal", pad=0.05) plt.savefig("geoid_2d.png", dpi=150)这一步如果正常，你会在海面区域看到明显的重力异常带，比如印度洋低值区、西太平洋高值区。
5.3 生成“土豆地球”三维模型
三维可视化的核心思路：先把经纬度网格映射到三维球面坐标，然后在径向叠加 Geoid 高度值。为了看到形变，夸张系数建议从 1000 倍开始。
import numpy as np import pyvista as pv # 将经纬度转成弧度 lon_r = np.deg2rad(lon) lat_r = np.deg2rad(lat) # 使用 WGS84 椭球参数，这里简化为平均半径 R = 6371000.0 # 网格化经纬度 lon_grid, lat_grid = np.meshgrid(lon_r, lat_r) geoid_grid = geoid # 确保 geoid 维度与 lon_grid 一致 # 夸张系数，可视情况调整 exaggeration = 1000.0 # 径向距离 = 地球半径 + 夸张后的 geoid 高度 r = R + geoid_grid * exaggeration # 转笛卡尔坐标 x = r * np.cos(lat_grid) * np.cos(lon_grid) y = r * np.cos(lat_grid) * np.sin(lon_grid) z = r * np.sin(lat_grid) # 构造 PyVista 网格 points = np.column_stack([x.ravel(), y.ravel(), z.ravel()]) cloud = pv.PolyData(points) # 通过 Delaunay 三维采样构建表面，或用结构化网格更高效 grid = pv.StructuredGrid() grid.points = points grid.dimensions = lon_grid.shape # 添加高度值作为标量，用于着色 grid.point_data["geoid"] = geoid_grid.ravel() # 绘制 plotter = pv.Plotter() plotter.add_mesh(grid, scalars="geoid", cmap="RdBu_r", show_edges=False) plotter.add_axes() plotter.show()这个代码是通用模板。实际使用中，结构化网格要求 lon 和 lat 是严格递增的一维数组，并且维度顺序与 reshape 一致。如果遇到索引错位，建议先用二维伪彩图比对数据方向，再做三维映射。
5.4 判断可视化是否成功
判断标准有三点：
- 三维模型整体保持球形，但表面有明显起伏。
- 颜色分布与二维伪彩图一致，高值区在图上呈现红色隆起，低值区呈蓝色凹陷。
- 旋转模型时，能看到赤道附近略微膨胀，两极方向偏平，且局部区域有明显的鼓包。
如果模型出现撕裂、飞点，先检查经纬度网格是否单调，再检查 geoid 值里是否有 NaN 或异常大值。
6. 功能测试与效果验证
只做出一张图还不够，建议按下面几个维度验证数据和处理流程。
6.1 与参考椭球对比
生成一组基准球面，半径等于地球平均半径，颜色设成灰色，再叠加 Geoid 模型。这样可以直观看到形变大小：
sphere = pv.Sphere(radius=R, theta_resolution=200, phi_resolution=200) sphere.point_data["geoid_ref"] = np.zeros(sphere.n_points) plotter = pv.Plotter() plotter.add_mesh(sphere, color="lightgray", opacity=0.3) plotter.add_mesh(grid, scalars="geoid", cmap="RdBu_r") plotter.show()如果 Geoid 模型确实是夸张后的，你会看到红色部分向外鼓出，蓝色部分向里塌陷，灰色球体夹在中间。
6.2 沿特定纬度做剖面
为了验证数据不是“乱画出来的”，可以沿某一条经度或纬度画剖面曲线。例如取北纬 30 度上的 Geoid 高度，横轴是经度，纵轴是高度：
import matplotlib.pyplot as plt # 假设 lat 和 geoid 都是二维数组 # 找到最接近 30°N 的纬度索引 target_lat = 30.0 lat_index = np.argmin(np.abs(lat - target_lat)) plt.figure(figsize=(10, 4)) plt.plot(lon, geoid[lat_index, :]) plt.xlabel("Longitude") plt.ylabel("Geoid Height (m)") plt.title(f"Geoid Profile at Latitude {target_lat}N") plt.grid(True) plt.savefig("geoid_profile.png", dpi=150)剖面图应该呈现连续、有物理意义的起伏，而不是随机噪声。
6.3 批量处理多时相数据
GRACE 数据是按月发布的，如果你拿到多个时期的 Geoid 或重力异常网格，可以批量渲染成动画，观察地表质量迁移。批量处理脚本的关键是循环遍历文件，并固定投影参数、色标范围，否则动画会出现颜色跳动。
for f in geoid_*.nc; do python render_frame.py --input "$f" --output "frame_${f%.nc}.png" done生成 PNG 帧后，用 ffmpeg 合成视频：
ffmpeg -framerate 10 -i "frame_%04d.png" -c:v libx264 -pix_fmt yuv420p geoid_animation.mp4批量任务最容易出现的问题是部分数据缺失、NaN 值、坐标范围不一致。建议在脚本里加一个“数据完整性检查”步骤，先统计每个文件的 min、max、mean，再做渲染。
7. 接口 API 与任务服务化
如果你的目标不是自己做一张图，而是把 Geoid 数据接入业务系统，可以考虑把数据处理模块封装成 REST API。
下面是一个 FastAPI 示例骨架，提供“查询指定窗口内 Geoid 统计值”的接口。实际部署时需要替换数据加载逻辑。
from fastapi import FastAPI, Query import xarray as xr import numpy as np app = FastAPI() # 全局加载一次，避免每次请求都读磁盘 ds = xr.open_dataset("geoid_model.nc") @app.get("/geoid/stats") def geoid_stats( lon_min: float = Query(...), lon_max: float = Query(...), lat_min: float = Query(...), lat_max: float = Query(...) ): # 用 xarray 切片 region = ds.sel( lon=slice(lon_min, lon_max), lat=slice(lat_min, lat_max) ) geoid = region["geoid"].values return { "min": float(np.nanmin(geoid)), "max": float(np.nanmax(geoid)), "mean": float(np.nanmean(geoid)), "std": float(np.nanstd(geoid)), "count": int(np.count_nonzero(~np.isnan(geoid))) } if __name__ == "__main__": import uvicorn uvicorn.run(app, host="127.0.0.1", port=8000)启动接口服务：
uvicorn main:app --host 127.0.0.1 --port 8000调用接口：
curl "http://127.0.0.1:8000/geoid/stats?lon_min=80&lon_max=120&lat_min=10&lat_max=50"返回结果示例：
{ "min": -32.5, "max": 48.2, "mean": 12.6, "std": 18.4, "count": 120000 }需要注意，上面的返回数字只是演示，实际数值取决于你加载的模型版本和区域。接口服务化以后，可以继续扩展切片下载、PNG 渲染、批量任务队列等功能。批量任务建议加任务 ID、日志文件和失败重试机制，不要在主进程里跑长时间渲染。
8. 资源占用与性能观察
这个模型和常见的 AI 推理不同，它主要吃内存、磁盘和一定的 GPU 渲染资源。
数据占用的估算公式很简单：全球网格点数等于纬度点数乘以经度点数。比如 1 度网格是 180×360=64800 个点；0.5 度是 360×720=259200 个点；0.25 度是 720×1440=1036800 个点。每个点如果是 float64，占 8 字节，100 万个点大概是 8MB。乍一看不大，但如果有多个时间步、多个变量，加上经纬度坐标和属性，整体数据量会成倍增加。
三维渲染时，PyVista 会把网格点构造成曲面。顶点数到几十万时，一般显卡都能流畅旋转；如果到了几百万点，建议先做下采样，否则旋转会卡顿，显存占用也会明显上升。降低资源占用的常用方法：
- 下载或加载数据前先裁剪目标区域
- 把 float64 转成 float32
- 三维渲染时降低球体分辨率
- 剖面图用二维画，不需要开三维窗口
- 批量渲染时关闭显卡抗锯齿，降低截图分辨率
观察资源占用，Windows 可以用任务管理器，Linux 可以用htop和nvidia-smi。如果显存不足，优先减网格点数，而不是换显卡。
9. 常见问题与排查方法
下面把最容易遇到的问题列成排查表。
| 问题现象 | 可能原因 | 排查方式 | 解决方案 |
|---|---|---|---|
| netCDF 文件打开失败 | 缺少 netCDF4 或 cftime | 打印报错信息 | 安装 netCDF4、cftime |
| 经纬度维度顺序不对 | 数据文件是 lat/lon，代码按 lon/lat 读取 | 查看变量 shape | 转置数组或调整索引 |
| 三维模型出现飞点 | Geoid 值含 NaN 或极大异常值 | 检查 min、max、NaN 数量 | 清洗数据，插值或剔除异常值 |
| 球体表面出现撕裂 | 经纬度网格不是严格递增 | 打印 lat、lon 数组 | 按升序重排网格 |
| 颜色标尺不明显 | Geoid 范围太小，夸张系数不足 | 查看数据 min/max | 增大夸张系数或调整配色范围 |
| 显存不足 | 三维网格点太多 | 查看显卡占用 | 下采样、裁剪区域、减少球体精度 |
| 批量任务卡住 | 单个文件读取太慢或内存溢出 | 打印日志定位到具体文件 | 增加超时控制、任务日志、跳过坏文件 |
| 剖面曲线不连续 | 网格中间有缺测值 | 检查 NaN 分布 | 插值或使用掩膜 |
最值得关注的是“坐标方向错误”。地理数据里lat有时从南到北，有时从北到南。如果生成的三维模型高低分布和官方图相反，大概率是纬度方向反了。先画二维图对比，再进三维，能省很多时间。
10. 最佳实践与使用建议
做完这套可视化流程，后面再碰到类似的重力场、地形、海洋测高数据，都可以沿用同一套思路。
第一，数据版本要记录清楚。模型文件名称、来源机构、发布时间、网格分辨率、单位、夸张系数，最好写进一个 config 文件或 README。别人拿到你的代码和图片时，能知道你到底画的是什么。
model: name: "XGM2019-like geoid" source: "ICGEM / NASA Earthdata" resolution_deg: 0.5 vertical_exaggeration: 1000 unit: "meter" output: format: "png" dpi: 150 colorbar: "RdBu_r"第二，三维可视化和平面的“平面图”要配套使用。平面图适合看全球重力异常分布，三维模型适合看整体形变感受。发布或汇报时，两者结合更有说服力。
第三，不要轻视合规问题。全球公开重力模型作为开放数据可以用于研究和教学，但商用前要确认许可条款。如果你的工作涉及高精度局部重力数据，务必遵守主管部门的使用规定，不要随意扩散。
第四，涉及可视化成果时，要在图注里明确“高度经过放大展示”。否则读者可能会误以为地球真实形状就是那么凹凸不平，这在科普上是误导。
第五，工程化批量任务时，一定要给每个任务加日志。GRACE 数据一个月一个文件，批量处理几十个月的数据很容易中途失败。记录每个文件的处理状态、输出路径、异常原因，是稳定跑完批量任务的关键。
11. 总结与下一步
这个 NASA 新 3D 模型最值得尝试的地方，并不是它本身有多炫酷，而是它把抽象的重力场变成了肉眼可见的形状。只要学会读取 Geoid 网格、映射到三维球面、加一个垂直夸张系数，你完全可以复现类似的“土豆地球”。
建议先从 1 度或 0.5 度公开网格开始，用 matplotlib 画二维图，再用 PyVista 生成三维模型。先验证数据方向和数据范围，再逐步加入夸张系数、剖面曲线、批量动画和 API 服务。
最容易踩的坑有两个：一是把 Geoid 当成真实地形，二是在三维映射时忽略经纬度顺序。这两个问题想清楚，整套流程基本不会卡住。
后续可以继续向两个方向扩展：一个是时变方向，用 GRACE/GRACE-FO 月重力场做时间序列动画，观察地下水、极地冰盖变化；另一个是拟合方向，把 Geoid 模型和全球地形、卫星测高、GNSS 高程数据做联合分析，解决高精度高程转换问题。
