预计算三角剖分后LinearNDInterpolator首次调用异常缓慢问题

问题详情 #

我使用Scipy将大型不规则网格（17000×14000，含大量NaN）上的多个场插值到规则网格，已预计算Delaunay三角剖分并通过pickle保存。发现：

• 构造LinearNDInterpolator速度很快，但首次调用插值的耗时是后续调用的10倍以上
• 该耗时不受插值场和目标网格大小影响（调试用目标网格仅17×14，首次调用耗时约10分钟）
• 仅第一次调用regrid函数中的第一个插值器时速度极慢，后续所有插值调用（包括其他regrid执行）都恢复正常

运行环境：14核SGE集群

代码实现 #

import xarray as xr
import os
import matplotlib.pyplot as plt
import gsw
from multiprocessing.pool import ThreadPool 
import dask

import numpy as np
from dask.diagnostics import ProgressBar
from scipy.interpolate import griddata
from scipy.spatial import Delaunay
from scipy.interpolate import LinearNDInterpolator

import matplotlib.tri as mtri
import time
import pickle


def regrid(tri,fnames):
        ds1 = xr.open_dataset(fnames[0],chunks={"xc":10,"yc":10})
        ds2 = xr.open_dataset(fnames[1],chunks={"xc":100})
        ds3 = xr.open_dataset(fnames[2],chunks={"xc":100})
        ds4 = xr.open_dataset(fnames[3],chunks={"xc":100})
        fullds = xr.combine_nested([ds1,ds2,ds3,ds4], concat_dim=["x"])
        im = fullds.sel(zc=52,drop=True)
        maskflat = np.logical_and(im.yc.values>-70,-im.depth.values<ds1.rf.values[52]).flatten()
        del(fullds)
        del(ds1)
        del(ds2)
        del(ds3)
        del(ds4)
        xmin,xmax = np.nanmin(im.xc),np.nanmax(im.xc)
        ymin,ymax = np.nanmin(im.yc),np.nanmax(im.yc)
        xs=np.linspace(xmin,xmax,17)
        ys=np.linspace(-70,-60,14)
        newxs , newys = np.meshgrid(xs,ys)
        uflat = im.u.values.flatten()
        vflat = im.v.values.flatten()
        xy = (im.xc.values.flatten()[maskflat],im.yc.values.flatten()[maskflat])

        print("starting v")
        vflat = vflat[~np.isnan(vflat)]
        start = time.time()
        interpolator = LinearNDInterpolator(tri,vflat)
        newv = interpolator(newxs,newys)
        end = time.time()
        print("newv done: ", end-start)

        print("starting u")
        start = time.time()
        interpolator = LinearNDInterpolator(tri, uflat[~np.isnan(vflat)])
        print("interpolator constructed")
        newu = interpolator(newxs,newys)
        end = time.time()
        print("newu done: ", end - start)
        print("starting mask")
        start = time.time()
        interpolator = LinearNDInterpolator(tri, maskflat[~np.isnan(vflat)])
        maskflat = interpolator(newxs,newys)
        end = time.time()
        print("mask done: ", end-start)
        print("mask")

        ds = xr.Dataset(\
            data_vars=dict(\
                u=("y", "x", newu),\
                v=("y", "x", newv),\
                mask=("y", "x", maskflat),\
            ),\
            coords=dict(\
                x=xs,\
                y=ys,\
            ),\
            attrs=dict(description=" by the slice"),\
        )
        outname = os.path.basename(fnames[0])
        ds.to_netcdf(outname.replace(".nc",'-slice.nc'))

with open("data/tri.pickle","rb") as f:
    tri = pickle.load(f)

#File names redacted
fnames = ["1","2","3","4"]
regrid(tri,fnames)
#File names redacted
fnames = ["a","b","c","d"]
regrid(tri,fnames)

执行输出 #

starting v
newv done:  570.6315619945526
starting u
interpolator constructed
newu done:  0.6614155769348145
starting mask
mask done:  0.4872896671295166
mask
starting v
newv done:  0.1620798110961914
starting u
interpolator constructed
newu done:  0.5560829639434814
starting mask
mask done:  0.47842955589294434
mask

原因分析 #

核心问题在于Scipy的Delaunay三角剖分对象在pickle序列化时，不会保存内部的搜索加速结构（如邻域搜索树）。当你第一次用该tri对象构建LinearNDInterpolator并执行插值时，Scipy会延迟初始化这个加速结构——对于大型三角剖分，这个初始化过程需要遍历大量三角面，因此耗时极长。

一旦加速结构被初始化，它会被缓存到传入的tri对象中，后续所有基于该tri的插值调用都会直接复用这个结构，所以速度会立刻恢复正常。

解决方案 #

在正式执行插值任务前，提前触发一次“预热”插值调用，让Scipy完成加速结构的初始化，后续所有插值操作都会受益。

修改后的关键代码 #

with open("data/tri.pickle","rb") as f:
    tri = pickle.load(f)

# 预热：用一个极小的测试点触发加速结构初始化
test_point = np.array([[tri.points[0, 0], tri.points[0, 1]]])
dummy_interp = LinearNDInterpolator(tri, np.zeros(len(tri.points)))
dummy_interp(test_point)
del dummy_interp  # 清理临时对象

# 后续正常执行regrid调用
fnames = ["1","2","3","4"]
regrid(tri,fnames)
fnames = ["a","b","c","d"]
regrid(tri,fnames)

额外优化建议 #

1. 复用插值器逻辑：如果多个场使用同一三角剖分，可考虑批量构建插值器，避免重复初始化（注意LinearNDInterpolator不支持动态修改values，需重新实例化，但加速结构已缓存）。
2. 主动内存回收：在regrid函数的del操作后，可添加import gc; gc.collect()，避免内存碎片影响集群环境下的性能。
3. 验证pickle内容：如果需要长期复用tri对象，可考虑在保存pickle前先完成预热，再序列化tri对象——这样下次加载时加速结构已存在，但需注意Scipy版本兼容性（不同版本的内部结构可能变化）。

内容的提问来源于stack exchange，提问作者Garrett Dreyfus