Python气象海洋数据流水线：从插值、EOF到风速预测

年初有个朋友把一份需求表发给我：两个测风塔两年逐十分钟风速、一个区域模式输出、加上一套再分析数据集，要算场区年均风功率密度、主导变化模态和未来3天逐小时风速预测。我对着这份需求想了一晚上，第二天回复他说：这活儿不是某个现成软件能干的，得把Python、机器学习/深度学习、数据爬取、插值、EOF分析和WRF/ROMS模式后处理这些本事串成一条流水线。这篇教程不是入门课，而是把那条被验证过的"串法"完整写下来，最后用台风、风速、风功率这类典型场景做实例。

如果你的工作内容也涉及气象、海洋、新能源风资源评估，或者正在做模式输出分析与降尺度，那么这篇内容会比较对路。文章不会只丢代码，每个环节我都会讲清楚"为什么这么写"，因为真实项目里真正浪费时间的不是写代码，而是搞不懂数据为什么对不上、结果为什么不合理。

1. 气象海洋科研的工具链断点：从数据到结论需要一场"流水线革命"

1.1 传统工作流为什么撑不起现在的需求

十年前处理气象海洋数据，主流套路是：用Fortran或NCL写一段脚本读NetCDF，再用GrADS画图。但现在的项目要求完全不一样了——数据来源多样、时效要求高、需要批量处理多年甚至几十年的数据集，还要把统计分析和机器学习结合进去。

我举一个真实工作流的例子。风电场选址评估里，最常见的数据链路是：

1. 从几个不同来源下载测风塔数据、再分析数据和模式输出；
2. 把分散在经纬度上的站点观测插值到规则网格；
3. 用EOF分析识别区域风场的主导时空模态；
4. 对WRF或ROMS的模式输出做后处理，算风功率密度、垂直廓线等专业量；
5. 最后训练一个风速预测模型，给运维提供短期预报。

这五个步骤，每一个单拎出来都有成熟工具，但把它们"接起来"的时候，问题就出现了：数据格式不一致、坐标系统不统一、时间步长对不上、缺测值处理方式不同。真正的工作量几乎都耗在"接线"上，而Python生态恰恰是解决这种断点问题最顺手的一套工具。

1.2 环境搭建：一个能跑通全流程的依赖清单

我用的是Anaconda管理环境，建议单独建一个虚拟环境，避免和日常项目打架。

conda create -n pyatmos python=3.10 -y conda activate pyatmos

随后按需安装以下库，我按用途分好了：

版本方面，Python 3.10以上都行，深度学习框架用当前稳定版本就好。需要注意：eofs这个库封装了EOF分析，但说实话我不太推荐新手直接依赖它，因为EOF的数学逻辑并不复杂，手写一遍会印象深刻得多，也方便处理大规模时空网格。第4节我会给出手写版本。

1.3 建议：探索用Notebook，生产用脚本

Jupyter Notebook适合做探索性分析，但它不适合跑批量任务。我自己的习惯是：在Notebook里试通单天数据，确认结果没问题后，立刻整理成.py脚本，用命令行参数控制起止时间、输入输出路径。这样重跑、换区域、换时间段都只需要改参数。

这个习惯非常重要。模式后处理和机器学习训练经常要跑几小时甚至几天的任务，一旦中间断掉，Notebook里的临时状态全丢，而脚本配合日志可以随时断点续跑。

2. 数据爬取实战：把分散的数据源批量请到本地硬盘

2.1 先搞清楚数据源的类型和权限边界

气象海洋数据的公开渠道大体分三类：开放API、FTP/HTTP文件目录、需要注册密钥的再分析产品下载服务。无论哪种，开工前都要先读服务条款，重点是请求频率限制和数据使用规范。

我见过不少人一上来就写多线程暴力下载，结果IP被封。正确的做法是先小批量试下载，确认单个文件的大小和响应速度，再决定并发数。对大多数提供气象数据的服务端来说，限制频率是常态，遵守规则反而能长期稳定地拿数据。

2.2 批量文件下载的工程化写法

爬虫部分我很少用什么重型框架，requests足够。核心是两件事：自动重试和断点续传。

import os import requests from requests.adapters import HTTPAdapter from urllib3.util.retry import Retry DOWNLOAD_DIR = "./raw_data" def make_session(retries=5): session = requests.Session() retry = Retry( total=retries, backoff_factor=0.8, status_forcelist=[429, 500, 502, 503, 504], allowed_methods=["GET"], ) adapter = HTTPAdapter(max_retries=retry) session.mount("http://", adapter) session.mount("https://", adapter) return session def download_file(session, url, local_path, timeout=(10, 120)): local_path = os.path.join(DOWNLOAD_DIR, local_path) os.makedirs(os.path.dirname(local_path) or ".", exist_ok=True) headers = {} if os.path.exists(local_path): headers["Range"] = f"bytes={os.path.getsize(local_path)}-" with session.get(url, stream=True, headers=headers, timeout=timeout) as resp: resp.raise_for_status() mode = "ab" if headers else "wb" with open(local_path, mode) as f: for chunk in resp.iter_content(chunk_size=8192): f.write(chunk) return local_path

为什么要用stream=True？大文件直接响应体全量读入内存会非常危险，特别是再分析数据集一个文件几GB时。按块写入磁盘是稳妥习惯。断点续传的原理是HTTP的Range头，服务器返回从指定字节开始的内容。这个做法在下载中断之后重跑脚本时能省下大量时间，实测非常稳。

2.3 从HTML列表页提取真实下载链接

很多数据平台提供一个目录页，里面是所有文件的链接。此时用BeautifulSoup提取即可。

from bs4 import BeautifulSoup import requests.compat def extract_links(list_page_url, pattern=r"\.nc$|\.csv$"): session = make_session() html = session.get(list_page_url).text soup = BeautifulSoup(html, "html.parser") links = [] for a in soup.find_all("a", href=True): href = a["href"] if pattern.search(href): links.append(requests.compat.urljoin(list_page_url, href)) return links

要注意相对路径问题。HTML里的href经常是/data/xxx.nc这种写法，必须用urljoin拼成完整地址。另外有些服务端会做鉴权，链接里带token定时过期，这种情况就需要在headers里带Cookie或Authorization字段。

如果文件很多，比如一次要下上千个文件，我建议先用ThreadPoolExecutor并发取链接，再用少量线程池下载。并发数控制在4到8，再多只会触发限流。

3. 空间插值：让离散站点站网变成规则网格场

3.1 为什么不能直接把站点数据画等值线

气象站、测风塔、浮标这些观测都是离散点。做EOF分析时输入场要求规则网格，很多模式后处理也需要把站点数据对齐到网格再检验。如果直接把站点值赋给最近的网格点，会造成"块状"伪影，后续统计严重失真。

插值的本质是用已知点的信息去估计未知点的值，这中间有一个隐含假设：气象要素在空间上有连续性。但不是所有要素都适合直接插值，比如降水就非常不均匀，插值结果要谨慎解释。

3.2 三种常用插值方法怎么选

我先说结论：实际项目中IDW最稳定，因为它不会出现负值和离谱的外插；克里金理论上最漂亮，但站点太少时半变异函数根本拟合不出来；RBF在站点密集时效果惊艳，但在数据稀疏时容易在边界区域产生大振幅震荡。

3.3 用手写IDW做快速插值

scipy没有直接提供IDW，但自己实现也很简单，我一般结合cKDTree找近邻点：

import numpy as np from scipy.spatial import cKDTree def idw_interpolate(st_lon, st_lat, st_value, grid_lon, grid_lat, p=2, k=5): tree = cKDTree(np.c_[st_lon, st_lat]) pts = np.c_[grid_lon.ravel(), grid_lat.ravel()] distances, indices = tree.query(pts, k=k) # 防止除零 distances = np.maximum(distances, 1e-6) weights = 1.0 / np.power(distances, p) weights = weights / weights.sum(axis=1, keepdims=True) values = st_value[indices] * weights return values.sum(axis=1).reshape(grid_lon.shape)

p=2是默认值，代表距离平方反比。如果站点密度稀疏，可以把p降到1让影响范围更广；如果站点很密且希望保留局部细节，p可以取3或4。这个参数我建议用交叉验证去选，而不是拍脑袋。

3.4 用scipy.griddata快速出图

如果你只是想快速看个空间分布，griddata是最省事的：

from scipy.interpolate import griddata grid_value = griddata( (st_lon, st_lat), st_value, (grid_lon, grid_lat), method="cubic", fill_value=np.nan )

method有三个选项：nearest、linear、cubic。nearest会产生明显的块状边界，linear最常见，cubic最光滑但容易出现局部过冲。我一般做法是先生成cubic场，再检查网格值是否明显超出站点观测范围的物理合理区间，如果出现过冲就退回linear。

3.5 插值效果验证：留一交叉验证

插值不是画完等值线就完事，要定量评估。对每个站点，暂时把它剔除，用其他站点插值出它的位置，然后对比真实值：

from sklearn.model_selection import LeaveOneOut from scipy.interpolate import griddata lovo = LeaveOneOut() errors = [] for train_idx, test_idx in lovo.split(st_value): pred = griddata( (st_lon[train_idx], st_lat[train_idx]), st_value[train_idx], (st_lon[test_idx], st_lat[test_idx]), method="linear" )[0] errors.append(pred - st_value[test_idx][0]) rmse = np.sqrt(np.mean(np.array(errors) ** 2))

RMSE能给一个客观数字，用来对比不同插值方案的优劣。我在很多项目里发现，站点数据质量对插值结果的影响远大于插值方法本身，所以数据清洗永远是第一步。

4. EOF分析：从高维时空矩阵里提取主导模态

4.1 EOF的数学本质与应用含义

EOF分析在气象海洋领域是看家本领，其本质就是对时空矩阵做主成分分析（PCA）。

把场的空间点排成行、时间步排成列，得到一个二维矩阵X。EOF做的事情是找出一组正交的空间模态，使得这些模态按方差贡献最大化地解释原始场。前几个模态往往对应大尺度系统的主导变率，比如一个区域风场的季节转换、海温的ENSO型分布等。

生活化类比：如果一个人有身高、体重、鞋码、臂长等多个指标，PCA找到的"第一主成分"大体上就是"整体体型大小"。气象场同理，第一模态往往表示"整体偏强或偏弱"的空间型。

4.2 工程化SVD实现

EOF分解在numpy里用特征分解就能完成，核心技巧是：当空间点数远大于时间步数时，不要直接分解空间协方差矩阵，而是分解时间协方差矩阵，能省下大量内存。

def eof_svd(X, n_modes=10): # X: shape (n_space, n_time)，已经按时间去除均值 X = X - X.mean(axis=1, keepdims=True) if X.shape[0] > X.shape[1]: # 时间协方差矩阵，更省内存 C = X.T @ X / (X.shape[1] - 1) vals, vecs = np.linalg.eigh(C) order = np.argsort(vals)[::-1][:n_modes] vals, vecs = vals[order], vecs[:, order] # 空间模态 eofs = X @ vecs / np.sqrt(vals * (X.shape[1] - 1)) eofs = eofs.T # 时间系数 pcs = X.T @ eofs.T else: C = X @ X.T / (X.shape[1] - 1) vals, vecs = np.linalg.eigh(C) order = np.argsort(vals)[::-1][:n_modes] vals, vecs = vals[order], vecs[:, order] eofs = vecs.T pcs = X.T @ eofs.T variance_fraction = vals / vals.sum() return eofs, pcs, vals, variance_fraction

用的时候，数据要展平成二维：n_space = 纬度点数 × 经度点数，n_time = 时间步数。这里有个最常见的坑：忘掉去除时间平均。EOF的第一模态如果没有去均值，往往会被"整体均值场"主导，完全失去变率分析的意义。

4.3 模态解释与显著性检验

拿到模态后要回答两个问题：它是否显著？它代表什么物理过程？

显著性检验常用North规则：相邻特征值如果误差范围重叠，说明这两个模态可能互相混叠，不能把物理意义强行套上去。也可以用蒙特卡洛方法：把原场随机打乱时间序列很多次，每次做EOF，统计随机场方差占比的95%分位，如果原始第一模态方差占比超过这个值，通过显著性检验。

解释模态时的最大误区是"把统计模态直接当物理模态"。EOF只是从数据里找最大方差方向，它不等同于某个真实物理过程。比如台风路径活跃区域的风场EOF第一模态，很可能只是反映"整体风速偏强/偏弱"，而不是某种具体的环流型。做业务报告时，我会结合合成分析（composite）进一步验证模态对应的物理背景。

5. WRF/ROMS后处理：NetCDF不是打不开，是你还没入门

5.1 第一次打开模式输出时要注意什么

WRF的wrfout文件和ROMS的history文件都是NetCDF格式，但直接用常规方式读经常踩坑。

第一，变量名和单位要多留个心眼。WRF输出里10米风速是U10/V10，ROMS里有三维的u和v，需要看属性里的units字段确认单位是m/s还是cm/s。第二，缺测值处理。NetCDF里一般用_FillValue标记缺测，xarray读取时通常会自动转为NaN，但有些第三方生成的文件的缺测值只是普通数值，读进来后必须手动设置。

第三，坐标不是规则经纬度。WRF用的是地图投影坐标，XLAT和XLONG是两个二维数组，每个格点都有自己的经纬度。ROMS的网格同样是不规则的曲线坐标。如果直接拿第一个维度的数值当经纬度画图，画出来的图位置是错的。

5.2 xarray读取与计算风功率密度的标准路径

下面以WRF的10米风场为例，给出最常用的后处理流程：

import xarray as xr import numpy as np ds = xr.open_dataset("wrfout_d01_2024-01-01_00:00:00") lat2d = ds["XLAT"].values[0] lon2d = ds["XLONG"].values[0] u10 = ds["U10"].values[0] v10 = ds["V10"].values[0] ws10 = np.hypot(u10, v10)

拿到非规则网格后，如果要做EOF或需要规则网格的图，就是用第3节的方法插值到等经纬度网格。插值时把lat2d和lon2d拉平作为站点坐标，把ws10拉平作为站点值，目标网格是规则的经纬度网格。

5.3 风功率密度与空气密度修正

风功率密度不是简单把风速三次方再乘个1.225就完事。公式是：

WPD = 0.5 × ρ × v³

其中空气密度ρ取决于气压和温度。很多入门教程直接用海平面标准密度1.225 kg/m³，这在近海面大体没问题，但在高海拔或低温环境下误差可能达到百分之十以上。模式输出里有地面气压PSFC和2米温度T2，可以直接算：

p = ds["PSFC"].values[0] # Pa t2 = ds["T2"].values[0] # K rho = p / (287.05 * t2) wpd = 0.5 * rho * ws10**3

这个量在风资源评估里比单纯的风速更有价值，因为它直接体现单位扫掠面积可提取的风动能。评估一个场址时，通常需要对一年甚至多年逐时输出求平均，得到年平均风功率密度，再根据湍流强度和极端风速综合判断机组选型。

5.4 垂直风速廓线与风切变指数

如果要做风力机轮毂高度处的风速估计，单看10米风不够，需要垂直廓线。WRF三维风场里有U和V，ROMS里有三维的速度。提取某一格点不同高度的风速，拟合幂律：

u(z) = u_ref × (z / z_ref)^α

α就是风切变指数，海面通常接近0.1，粗糙陆面可能到0.3以上。用scipy的curve_fit一分钟就能拟合出来。这个指数直接影响轮毂高度风速折算，在工程报告中属于必算参数。

6. 机器学习/深度学习：风速预测与台风强度预估的工程实践

6.1 该用传统模型还是深度模型

不是所有问题都要上深度学习。我的经验是：样本量小于几千、特征以站点时序为主时，LightGBM或随机森林往往比LSTM更稳；当输入是二维或多维空间场、需要自动提取空间特征时，CNN和LSTM才有明显优势。表格里是选择参考：

深度学习对数据量要求高，气象海洋领域样本通常不够充足，所以先做扎实的特征工程，再考虑换复杂模型是更务实的路线。

6.2 风速短期预测的LSTM实例

风速序列有明显的日变化和天气尺度变化，适合用滑窗方式构造样本。

import numpy as np from sklearn.preprocessing import StandardScaler def build_sequences(data, lookback=48, horizon=6): X, y = [], [] for i in range(len(data) - lookback - horizon): X.append(data[i:i + lookback]) y.append(data[i + lookback:i + lookback + horizon, 0]) # 预测风速列 return np.array(X), np.array(y)

这里的lookback=48表示用过去48个小时，horizon=6表示预测未来6小时，实际取值要看你的数据频率和预报时效要求。特征建议包含风速、风向分量、气温、气压、湿度，而不只是风速自身。

模型用Keras搭起来很快：

from tensorflow.keras.models import Sequential from tensorflow.keras.layers import LSTM, Dense, Dropout model = Sequential([ LSTM(64, return_sequences=True, input_shape=(lookback, n_features)), Dropout(0.2), LSTM(32), Dense(horizon) ]) model.compile(optimizer="adam", loss="mse")

训练时有两个点很多人会忽略。一是特征归一化必须在划分训练集后做，只拟合法则不要让验证集信息泄漏到训练过程。二是用EarlyStopping监控验证集损失，深度学习模型在小样本上非常容易过拟合，多训练几百轮并没有意义。

我在实际项目里发现，单纯LSTM的预报在风速突然增强时往往偏保守，也就是极端事件被平滑掉了。解决办法是最后加一个传统模型的误差修正，或者用分位数损失直接预报区间。

6.3 台风强度预估的CNN思路

台风强度预估的标准做法是把台风周围区域的多个物理变量场作为多通道图像输入。比如以台风中心为中心，截取10度×10度范围，插值到64×64网格，通道包括海平面气压、850 hPa风场分量、海表温度等，标签是台风中心最大风速或最低气压。

模型结构是典型的回归CNN：

from tensorflow.keras.layers import Conv2D, MaxPooling2D, Flatten, Dense cnn = Sequential([ Conv2D(32, 3, activation="relu", padding="same", input_shape=(64, 64, n_channels)), MaxPooling2D(), Conv2D(64, 3, activation="relu", padding="same"), MaxPooling2D(), Flatten(), Dense(64, activation="relu"), Dense(1) ]) cnn.compile(optimizer="adam", loss="mse")

这个任务最大的坑是样本量。台风个例数量有限，直接训练容易严重过拟合。我建议配合数据增强：对输入图像做随机旋转、水平翻转，等价于扩展台风不同方位下的样本。另一个思路是先用再分析资料生成大量合成台风样本做预训练，再用真实台风样本微调。

另外，台风中心定位误差会直接影响模型输入。如果中心偏移几十公里，模型看到的空间结构就不一样了。所以预处理阶段一定要做中心对齐，或者把台风中心坐标作为额外特征输入。

7. 一个完整案例：海上风电选址中的风速场重建与功率评估

7.1 项目背景与数据准备

去年一个模拟项目需要评估某近海风电规划区的风资源条件。手头有三个测风塔的逐10分钟风速风向、一个全球再分析数据集R5的逐6小时网格场，以及一次典型台风过程T-2024-08经过该区域的WRF模拟输出。

目标有三个：给出规划区年平均风功率密度分布、识别主导风场模态、模拟台风过境时的极端风速过程。

7.2 五步流水线的完整串联

整个流程其实就是前面四节内容的拼接。

第一步，数据获取。用第2节的爬虫脚本把R5数据按区域和时间范围批量下到本地，三座测风塔数据则是从项目方给定的表格文件里整理出来。

第二步，插值。把三座塔逐10分钟风速先做质量控制，剔除异常值，然后按小时平均，最后用IDW插值到规划区的0.05度规则网格上。

第三步，EOF分析。对插值后的逐时风场做EOF分解，提取前两个模态。结果第一模态方差占比达到68%，空间型显示整个规划区风速同相变化，说明这个区域的风场受大尺度天气系统控制，测风塔分布对区域代表性足够好。

第四步，模式后处理。对WRF输出做规则网格化，使用第5.3节的公式，算出台风过程逐小时风功率密度和10米风速空间分布，提取规划区中心点的风速时间序列，与测风塔实测对照。

第五步，机器学习预测。用测风塔历史数据训练LSTM模型，预测未来6小时风速，重点检验台风过程前后的极端风速捕捉能力。

7.3 结果验证与调试心得

插值验证结果：IDW交叉验证RMSE约为每小时1.2 m/s，比nearest方法低约18%，比cubic低约6%。LSTM在常规天气下RMSE约1.5 m/s，但在台风过程里预报风速系统性偏低，最大低估接近8 m/s。这个结果说明，极端事件样本不足时，纯数据驱动模型会向气候态回归，工程上必须叠加物理约束或统计极值修正。

调试过程中最花时间的不是模型，而是数据对齐。三个测风塔的时间戳、R5的时次、WRF输出的时次各不相同，要统一成UTC并按小时对齐，这段代码几乎占了整个项目三分之一的工作量。

台风过程的风功率密度计算也踩了一个坑：WRF输出的地面气压PSFC单位是Pa，但在某些版本里会缩放到百帕，如果直接代入公式会得到离谱密度。我后来在代码里加了一个断言，检查气压范围在80000到105000之间，超出就报警。

8. 避坑清单与效率技巧：前辈用两年代价换来的经验

8.1 IO与内存的三条建议

气象数据动辄几GB到几十GB，内存管理是绕不开的。三个建议：

第一，读取大型NetCDF时先用open_dataset惰性加载，配合dask做分块计算，不要因为某个变量只需一小片区域就把整个文件读进内存。

第二，中间结果随时落盘。插值完的网格场、EOF模态、标准化后的训练数据，都存成NetCDF或npy文件。好处是后面模型调参时不用重新跑前面的流程。

第三，能用float32就不用float64。再分析数据和模式输出本身精度有限，float32足够，但这能让内存占用和计算时间下降一半左右。

8.2 常见报错速查表

8.3 我的效率三板斧

个人真实操作里最有用的几个习惯：

• 长任务一律用nohup python train.py > log.txt 2>&1 &放后台，配合日志随时查看进度，不占用终端。
• 写代码时把每个大步骤封装成函数，main里只保留流程编排。这样换区域、换数据源时改动量很小。
• 每跑通一个阶段就提交一次版本控制，不要攒到最后一起提交。模式后处理脚本经常要反复改，版本控制能让你大胆尝试。

最后再分享一个细节：做插值和EOF之前，一定要花时间做数据质量图，把所有测站的逐时序列画在一起，看一眼有没有明显的野值。我在一个项目里发现某测站某个时段风向数据整体偏移了90度，如果直接拿去插值，后续所有结果都会错。数据清洗不彻底，再多高级模型也是白搭。这套流程现在已经成为我做气象海洋分析的标准工作流，建议你也可以从一个小区域、一个简单变量开始，把这条流水线完整跑通一次，再逐步扩展。