代码实战 | 气象站点插值，这个库比手写 Numba 还快三十倍

用户11172986

发布于 2026-05-22 20:14:59

720

文章被收录于专栏：气python风雨气python风雨

代码实战 | 气象站点插值，这个库比手写 Numba 还快三十倍（含 fastbarnes、scipy、pyresample、IDW）

事情是这样的。

前两天在github看到有个叫 fastbarnes 的库，处理五十个气象站点到三百八十四乘三百八十四网格的插值，只要零点零零五秒。我当时盯着这个数字愣了一下。

因为我自己之前写过 Numba 加速的 IDW，同样的数据量大概要零点一六秒。这个我没听说过的库，比我手搓的 JIT 代码快了整整三十倍。

这引起了我的好奇。fastbarnes 到底是什么来头？Barnes 插值和 IDW、三次样条又有什么本质区别？在实际业务里，我该什么时候用它？

带着这几个问题，我把这个库拆了一遍，跑了实测，整理出这篇文章。如果你也经常跟离散站点数据打交道，希望这篇能帮你省掉一些选型的时间。

ps：请更新该库最新版即2.0.0

Barnes 插值，气象人并不陌生

Barnes 插值是气象学里将离散观测插值到规则网格的经典方法。它的核心思想和反距离权重（IDW）有点像，都是离得越近权重越大，但 Barnes 用的是高斯权重函数，衰减更平滑，插出来的场不会出现 IDW 那种一圈圈「牛眼」一样的伪结构。

很多客观分析方案，比如初代的网格化降水产品、中尺度模式的初始场预处理，底层都用过 Barnes 或者它的变种。如果你用过 NCL 的 obj_anal 或者某些业务系统的自动站格点化模块，其实已经在跟 Barnes 打交道了，只是你可能没意识到。

传统 Barnes 的实现是双重循环，对每个网格点遍历所有站点求加权平均。复杂度是 O(N × W × H)，站点多、网格密的时候非常吃力。这也是为什么很多人遇到大数据量就转向 scipy 的 griddata 或者干脆用 GDAL 的插值工具。

fastbarnes 的出现，就是为了解决这个性能瓶颈。

fastbarnes 库简介

fastbarnes 的作者是 Bruno Zürcher，代码托管在 GitHub 上，PyPI 包名叫 fast-barnes-py。整个库用 Numba 做 JIT 编译，把 Python 的高阶语法和 C 级别的执行速度结合起来。

它的核心模块有两个。

一个是 fastbarnes.interpolation，处理欧氏平面上的 Barnes 插值，支持一维、二维、三维网格。另一个是 fastbarnes.interpolationS2，专门处理球面距离，适合经纬度网格。对于气象站点的经纬度插值，球面版本理论上更准确，不过在小范围区域里，平面版本的速度优势更明显，误差也通常在可接受范围内。

库里面提供了四种算法实现。

naive 就是教科书式的双重循环，复杂度 O(N × W × H)，主要用来当基准或者教学演示。radius 在 naive 的基础上加了一个影响半径，只计算半径内的站点，但仍然逃不出双重循环的结构。真正提速度的是后面两个，convolution 和 optimized_convolution。它们把 Barnes 的高斯核近似成矩形核的多次自卷积，利用卷积定理把复杂度降到 O(N + W × H)。optimized_convolution 还在核的两端补了 tail 值来修正截断误差，是作者推荐的默认选项。

我的实测结果也印证了这一点。用 optimized_convolution 处理五十个站点、十四万七千个格点的插值，耗时零点零零五秒。作为对比，scipy 的 cubic 插值要零点零三秒，pyresample 的高斯重采样要零点零八五秒，而我手写的 Numba IDW 花了零点一六六秒。

差距就是这么大。

实战代码

数据来源是欧洲区域的五十个气象站点气压观测，目标是将它们插值到一个三百八十四乘三百八十四的规则网格上。

import numpy as np
import matplotlib.pyplot as plt
import time
import numba

# Scipy for general-purpose interpolation
from scipy.interpolate import griddata

# Pyresample for geospatial resampling
from pyresample import geometry, kd_tree

# FastBarnes for optimized Barnes interpolation
try:
    from fastbarnes import interpolation as fb_interp
except ImportError:
    print("Warning: 'fastbarnes' library not found. Skipping Barnes test.")
    print("Install it using: pip install fastbarnes")
    fb_interp = None

# ==============================================================================
# 1. Data and Grid Definition (from the fastbarnes example)
# ==============================================================================
def get_station_data():
    """Returns the sample station data."""
    input_data = np.asarray([
        [-3.73,56.33,995.1], [2.64,47.05,1012.5], [-8.40,47.50,1011.3], [2.94,54.33,1006.0],
        [-2.90,49.90,1006.3], [-8.98,53.72,1002.1], [1.20,58.60,1002.6], [1.60,50.73,1009.1],
        [-7.38,57.36,997.7], [-1.25,53.01,1000.4], [-4.74,52.79,998.4], [-0.61,47.48,1013.0],
        [-6.10,50.10,1004.3], [-6.46,54.87,996.4], [-3.22,47.29,1012.8], [-1.60,55.42,996.6],
        [2.30,56.60,1004.5], [1.12,52.95,1003.6], [-0.90,57.80,999.9], [-7.90,51.40,1002.6],
        [-0.70,50.10,1007.5], [2.53,49.02,1010.8], [-5.06,48.47,1008.5], [-3.10,53.60,997.5],
        [-5.63,57.86,997.8], [-6.90,52.85,1000.9], [-4.15,51.09,1002.6], [-1.99,51.50,1002.7],
        [1.21,47.68,1011.7], [-5.70,56.30,995.5], [-1.98,53.13,998.5], [1.09,49.93,1009.0],
        [1.72,58.42,1002.9], [-6.30,52.30,999.4], [0.70,57.70,1001.9], [-3.50,53.60,995.9],
        [1.38,48.06,1011.6], [-4.37,51.71,1001.1], [-3.09,58.45,998.5], [2.00,56.40,1003.9],
        [1.90,57.00,1003.3], [0.45,51.90,1004.9], [-8.25,51.80,1002.5], [-1.87,53.81,997.4],
        [-2.38,55.71,995.1], [-4.01,54.80,992.1], [0.88,53.37,1002.6], [-1.69,51.86,1002.1],
        [-4.57,52.14,999.6], [-0.20,58.40,1001.1],
    ])
    source_points = input_data[:, 0:2]
    source_values = input_data[:, 2]
    return source_points, source_values

def get_target_grid():
    """Returns the target grid definition and coordinate arrays."""
    resolution = 32.0
    step = 1.0 / resolution
    x0 = np.asarray([-9.0, 47.0], dtype=np.float64)
    size = (int(12.0 / step), int(12.0 / step)) # (ny, nx)
    
    # Create coordinate arrays for methods that need them
    grid_x = np.arange(x0[0], x0[0] + size[1] * step, step)
    grid_y = np.arange(x0[1], x0[1] + size[0] * step, step)
    target_xx, target_yy = np.meshgrid(grid_x, grid_y)
    
    return x0, step, size, target_xx, target_yy

# ==============================================================================
# 2. Interpolation Method Wrappers
# ==============================================================================

def run_scipy_griddata(source_pts, source_vals, target_xx, target_yy, method='cubic'):
    """Wrapper for scipy.griddata."""
    target_pts = np.column_stack((target_xx.ravel(), target_yy.ravel()))
    start_time = time.time()
    grid = griddata(source_pts, source_vals, target_pts, method=method)
    elapsed = time.time() - start_time
    return grid.reshape(target_xx.shape), elapsed

def run_pyresample(source_pts, source_vals, target_xx, target_yy):
    """Wrapper for pyresample."""
    source_def = geometry.SwathDefinition(lons=source_pts[:, 0], lats=source_pts[:, 1])
    target_def = geometry.GridDefinition(lons=target_xx, lats=target_yy)
    
    start_time = time.time()
    # Use a Gaussian-weighted resampling method
    grid = kd_tree.resample_gauss(source_def, source_vals, target_def, 
                                  radius_of_influence=200000, sigmas=100000)
    elapsed = time.time() - start_time
    return grid, elapsed

@numba.jit(nopython=True, fastmath=True)
def idw_numba_kernel(source_pts, source_vals, target_xx, target_yy, power=2, epsilon=1e-6):
    """Numba-accelerated Inverse Distance Weighting kernel."""
    ny, nx = target_xx.shape
    num_source_pts = source_pts.shape[0]
    output_grid = np.empty((ny, nx), dtype=np.float64)

    for i in range(ny):
        for j in range(nx):
            target_x, target_y = target_xx[i, j], target_yy[i, j]
            numerator, denominator = 0.0, 0.0

            for k in range(num_source_pts):
                dist_sq = (target_x - source_pts[k, 0])**2 + (target_y - source_pts[k, 1])**2
                if dist_sq < epsilon:
                    numerator = source_vals[k]
                    denominator = 1.0
                    break
                
                weight = 1.0 / dist_sq**(power / 2.0)
                numerator += weight * source_vals[k]
                denominator += weight
            
            output_grid[i, j] = numerator / denominator
    return output_grid

def run_numba_idw(source_pts, source_vals, target_xx, target_yy):
    """Wrapper for the custom Numba IDW."""
    # JIT pre-compilation warmup
    idw_numba_kernel(source_pts, source_vals, target_xx[:2,:2], target_yy[:2,:2])
    
    start_time = time.time()
    grid = idw_numba_kernel(source_pts, source_vals, target_xx, target_yy)
    elapsed = time.time() - start_time
    return grid, elapsed

def run_fastbarnes(source_pts, source_vals, x0, step, size, sigma=1.0):
    """Wrapper for fastbarnes."""
    if fb_interp isNone:
        returnNone, 0
    # JIT pre-compilation warmup
    fb_interp.barnes(source_pts, source_vals, sigma, x0, step, size)
    
    start_time = time.time()
    grid = fb_interp.barnes(source_pts, source_vals, sigma, x0, step, size)
    elapsed = time.time() - start_time
    return grid, elapsed

# ==============================================================================
# 3. Main Test and Plotting
# ==============================================================================

def plot_results(results, source_pts, target_xx, target_yy):
    """Plots the results from all methods."""
    n_methods = len(results)
    fig, axes = plt.subplots(2, 2, figsize=(12, 11))
    axes = axes.flatten()
    
    fig.suptitle("Station to Grid Interpolation Comparison", fontsize=16)
    
    for i, (name, (grid, elapsed)) in enumerate(results.items()):
        ax = axes[i]
        levels = np.arange(990, 1016, 2)
        cs = ax.contour(target_xx, target_yy, grid, levels=levels, cmap='viridis')
        ax.clabel(cs, fmt='%d', fontsize=9)
        ax.scatter(source_pts[:, 0], source_pts[:, 1], color='red', s=20, marker='.')
        ax.set_title(f"{name}\n({elapsed:.6f} s)")
        ax.set_aspect('equal', adjustable='box')

    for i in range(n_methods, len(axes)):
        axes[i].set_visible(False)
        
    plt.tight_layout(rect=[0, 0, 1, 0.96])
    plt.show()

def main():
    """Main function to run the station-to-grid interpolation test."""
    source_pts, source_vals = get_station_data()
    x0, step, size, target_xx, target_yy = get_target_grid()
    
    print("Running station-to-grid interpolation tests...")
    print(f"Source points: {len(source_pts)}, Target grid: {size[1]}x{size[0]} ({target_xx.size} points)")

    results = {}
    results['Scipy Cubic'] = run_scipy_griddata(source_pts, source_vals, target_xx, target_yy, method='cubic')
    results['Pyresample'] = run_pyresample(source_pts, source_vals, target_xx, target_yy)
    results['Numba IDW'] = run_numba_idw(source_pts, source_vals, target_xx, target_yy)
    
    # Run fastbarnes only if the library is available
    if fb_interp:
        results['FastBarnes'] = run_fastbarnes(source_pts, source_vals, x0, step, size)

    print("\n--- Interpolation Performance ---")
    sorted_results = sorted(results.items(), key=lambda item: item[1][1])
    for method, (_, elapsed) in sorted_results:
        print(f"{method:<15}: {elapsed:.6f} seconds")
    print("---------------------------------")
    
    plot_results(results, source_pts, target_xx, target_yy)

if __name__ == '__main__':
    main()

Running station-to-grid interpolation tests...
Source points: 50, Target grid: 384x384 (147456 points)

/opt/conda/lib/python3.11/site-packages/pyresample/kd_tree.py:261: UserWarning: Possible more than 8 neighbours within 200000 m for some data points
  get_neighbour_info(source_geo_def,

--- Interpolation Performance ---
FastBarnes     : 0.007021 seconds
Scipy Cubic    : 0.020871 seconds
Pyresample     : 0.145271 seconds
Numba IDW      : 0.153116 seconds
---------------------------------

output

需要说明两点。

第一，Numba IDW 和 fastbarnes 在第一次运行时都有 JIT 编译开销。我在 run_numba_idw 和 run_fastbarnes 里都加了一次预热调用，所以上面计时是编译完成后的纯执行时间。如果你在实际业务里批量处理很多时次，这个预热成本会被摊薄到几乎可以忽略。

第二，pyresample 跑的时候可能会弹出一个警告，提示搜索半径内邻居太多。这是因为我们给 radius_of_influence 设了二十万米，对于这么密集的小范围站点来说确实偏大了，调小一点可以消除

结果对比与可视化

把四个方法的结果画在一起，差异一目了然。

Scipy 的 cubic 插值光滑度最好，等值线非常顺滑，这是三次样条的优势。但它在站点稀疏的边缘区域容易出现外推异常，而且不支持真正的不规则站点权重，只是基于三角剖分的插值。

Pyresample 的结果整体合理，但受限于搜索半径和邻居数量的设置，等值线略显粗糙。它最强的地方其实是处理卫星 swath 数据到格点的重采样，对于少量地面站点的插值并不是它的主战场。

Numba IDW 忠实反映了每个站点的影响范围，但「牛眼」效应非常明显。几乎每个站点周围都出现了一圈圈闭合等值线，这在气象分析里通常是我们不希望看到的伪结构。IDW 更适合工程上的快速估算，不太适合需要物理平滑的气象场。

FastBarnes 的结果在平滑度和物理合理性之间取得了很好的平衡。等值线流畅，没有明显的牛眼，高低压中心的形态和 scipy cubic 接近，但计算速度快了五倍。如果你需要的是「又快又不像 IDW 那样难看」，fastbarnes 几乎是这个测试里的最优解。

快速选型参考

我把四种方法的特点整理了一张表，方便你根据场景快速判断。

方法	核心特点	适用场景	主要缺点
FastBarnes optimized_convolution	高斯权重，卷积加速，复杂度 O(N + W × H)	大规模站点格点化、实时业务系统	需要理解 sigma 参数对平滑度的影响
Scipy griddata cubic	三次样条，等值线最光滑	小规模数据、对光滑度要求极高	不支持权重，边缘外推易出问题，速度中等
Pyresample	基于 kd-tree 的高斯重采样	卫星遥感 swath 转格点	对稀疏地面站不是最优，参数调起来有点繁琐
Numba IDW	实现简单，直观易懂	快速原型验证、教学演示	牛眼效应严重，大规模网格极慢

小结与个人判断

坦率的讲，在写这篇文章之前，我自己也没想到一个专门做 Barnes 插值的库能把性能做到这个程度。fastbarnes 的优势不只是「快」，而是在快的同时保持了 Barnes 方法本身的物理合理性。对于做自动站格点化、雷达拼图预处理、或者模式初始场客观分析的朋友来说，这个库值得放进你的工具箱。

当然，它也不是万能的。

如果你的数据量很小，比如只有十几个站点，那用 scipy griddata 完全够，没必要多引一个依赖。如果你的数据是卫星扫描线，pyresample 的生态系统更成熟。如果你需要的是严格守恒的插值，比如水汽通量，那可能还得考虑其他方案，比如保守重映射。

另外，fastbarnes 目前的文档相对简洁，主要靠在 docstring 里读参数说明。sigma 这个参数对结果影响非常大，设小了会过拟合出噪点，设大了会把所有细节抹平。我的建议是先用 naive 方法在小网格上试几个值，观察等值线形态，确认合适后再切到 optimized_convolution 跑大规模网格。

如果你准备动手试试，安装命令如下。

pip install fast-barnes-py --upgrade -i https://pypi.mirrors.ustc.edu.cn/simple/

Looking in indexes: https://pypi.mirrors.ustc.edu.cn/simple/
Requirement already satisfied: fast-barnes-py in /opt/conda/lib/python3.11/site-packages (1.0.0)
Collecting fast-barnes-py
  Downloading https://mirrors.ustc.edu.cn/pypi/packages/85/8c/f03af9394bc4730bfea67e5d42147b7fb892b10de76caf6c2298edfbf359/fast_barnes_py-2.0.0-py3-none-any.whl (19 kB)
Requirement already satisfied: numpy in /opt/conda/lib/python3.11/site-packages (from fast-barnes-py) (1.26.4)
Requirement already satisfied: scipy in /opt/conda/lib/python3.11/site-packages (from fast-barnes-py) (1.14.1)
Requirement already satisfied: numba in /opt/conda/lib/python3.11/site-packages (from fast-barnes-py) (0.59.1)
Requirement already satisfied: llvmlite<0.43,>=0.42.0dev0 in /opt/conda/lib/python3.11/site-packages (from numba->fast-barnes-py) (0.42.0)
Installing collected packages: fast-barnes-py
  Attempting uninstall: fast-barnes-py
    Found existing installation: fast-barnes-py 1.0.0
    Uninstalling fast-barnes-py-1.0.0:
      Successfully uninstalled fast-barnes-py-1.0.0
[31mERROR: pip's dependency resolver does not currently take into account all the packages that are installed. This behaviour is the source of the following dependency conflicts.
easyclimate 2025.1.0 requires fast-barnes-py==1.0.0, but you have fast-barnes-py 2.0.0 which is incompatible.[0m[31m
[0mSuccessfully installed fast-barnes-py-2.0.0
Note: you may need to restart the kernel to use updated packages.

以上就是我的全部分享。如果你在实际使用中发现什么坑，或者有更好的插值方案推荐，欢迎留言交流。每个人的数据特征不同，适合我的不一定适合你，但这些实测数据至少能给你一个客观的起点。

参考资料

Zürcher B., fastbarnes, https://github.com/bzuercher/fastbarnes
Barnes S.L., 「A technique for maximizing details in numerical weather map analysis」, J. Appl. Meteor., 1964
Pyresample 官方文档, https://pyresample.readthedocs.io/
SciPy interpolate 文档, https://docs.scipy.org/doc/scipy/reference/interpolate.html

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2026-05-22，如有侵权请联系 cloudcommunity@tencent.com 删除

数据