格点数据插值到站点这个话题,气象、海洋、环境领域几乎每周都会碰到:再分析产品是规则的格点场,观测站是零散的点位,两者要放在一起算误差、做检验,首先就得把格点数据插值到站点上。很多新手上来就写双层循环,一边找最近格点一边调用插值函数,数据一大就卡得人想摔键盘,而且最邻近插值和双线性插值到底怎么选,心里也没底。这篇文章我把这两种方法掰开揉碎讲清楚,给出可直接跑通的 Scipy 和 Xarray 两套代码,再把性能优化、边界处理、NaN传播这些实际项目里必然遇到的坑一起说了。适合正在做模式评估、站点检验、再分析资料处理、空气质量预报验证的数据分析人员和研究生,看完基本能直接用到自己的项目里。
1. 站点和格点维度对不上,是一切麻烦的起点
1.1 两种数据形态的本质差异
先理清一个基本事实:格点数据是一张大表。拿一张典型的0.25°再分析温度场来说,纬度维度可能有241个点,经度维度可能有1440个点,数据场本身的shape是(nlat, nlon),它背后隐含的坐标信息完全由lat和lon这两个一维数组决定。比如data[10, 20]这个值,对应的就是lat[10]这条纬线和lon[20]这条经线交叉处的格点。
站点数据则完全是另一回事:它不是矩阵,而是一张点位表。每一行记录一个站点的经纬度、海拔、编号,以及对应的观测值。这些站点是空间上零散的、分布不均匀的,有的地方几十公里一个站,有的地方几百公里没有站。
这两个数据结构之间没有天然的对应关系,你想拿格点值去跟站点观测值做差、算均方根误差、画散点图,就必须先解决一个问题:格点数据如何落到站点位置上去。这就是"格点数据插值到站点"这个需求的来源。
1.2 验证场景里为什么要反着插
很多新手会问:既然格点数据是规则场,站点是散点,那我能不能把站点数据插到格点上,然后直接对比两个格点场?理论上可以,但做验证评估时,强烈建议不要这么干。
原因很简单:观测站点分布极不均匀,把站点插值到格点,等于在站点稀疏区"脑补"了一片观测值,这会人为扭曲观测场原本的真实空间结构。反过来做——把格点插到站点——则不会改变观测数据本身,你只是从模式场里取了站点所在位置的值来做一对一对比,模式场没有被污染,观测值也没有被污染。
所以做检验评估、偏差订正、TS评分这类工作,行业惯例基本一致:格点向站点插值。这篇文章所有代码和讨论也都基于这个方向。
1.3 插值选择会影响你的结论
还有一个容易被低估的点:插值方法不是打印机设置,选错了真的会改变你的业务结论。最邻近插值会把站点值强行取成"距离最近的格点值",当网格分辨率较粗时,这个格点可能离站点还有几十公里,地形差异、海陆差异带来的偏差全都被忽略了。双线性插值会平滑掉一部分极端值,如果你正在做极端降水事件评估,双线性插值后的格点峰值可能被明显拉低,导致你误判模式对极端事件的再现能力。
所以,不要觉得"能出数就是对的"。理解了每种插值在数学上到底做了什么,你才知道它适不适合你的变量和场景。
2. 最邻近和双线性插值的算法本质与直觉类比
2.1 最邻近插值:给每个站点找一个"最顺路"的格点
最邻近插值的逻辑非常简单:对每个站点,找到所有格点中距离它最近的那个,直接把这个格点的值赋给站点。
在规则网格上,这个"找最近格点"的过程不需要遍历所有格点。如果你的网格经度从100°E开始、分辨率0.5°、共60个格点,那么一个116.39°E的站点,落在哪个经度索引附近,本质上就是解一个一元一次方程:
lon_idx = round((116.39 - 100.0) / 0.5)纬度方向同理。算出来的索引对应的格点(lat_idx, lon_idx)就是该站点的最近邻。注意这里有个非常关键的细节:真正的最邻近应该按球面距离来算,但对于常规的区域网格、站点又在格点覆盖范围内的情况下,直接用经纬度网格坐标算"矩形距离"误差很小,工程上都这么用。
打个比方:你叫外卖,平台按配送距离给你匹配"最近的门店"。最邻近插值就是这个逻辑——不计算附近几家店可能有什么差异,直接取最近那家店的出餐结果。
2.2 双线性插值:在一片矩形格子里横竖各插一次
双线性插值比最邻近多走了一步:它不只找最近的一个点,而是先找到环绕站点的那四个格点,然后在经度方向做两次线性插值,再在纬度方向做一次线性插值。
具体拆开看:假设站点坐标为(lat_s, lon_s),它落在经度索引i和i+1、纬度索引j和j+1围成的矩形里。先在经度方向,用lon_i和lon_{i+1}两个位置上的值,按站点经度与两端点的距离比例,插出上下两条边上lon_s处的值;接着在这两个新值之间,再按纬度比例做一次线性插值,最终得到站点处的值。
这个过程的物理解释是:双线性插值假设物理量在网格单元内沿经度方向和纬度方向都近似线性变化。它比最邻近传递了更多空间信息,结果更平滑,更接近"站点周围格点值共同影响站点值"的自然直觉。
类比一下:你站在两块地砖交界附近,想知道脚下这块地面上的温度。最邻近只靠最近的一平方米地砖的温度告诉你,双线性则是把你周围四平方米地砖的温度按距离远近加权平均一下给你一个估计值。
2.3 两种方法各自的"软肋"
最邻近插值的问题在于阶梯效应。它会让插值结果场呈现明显的块状突变,尤其当网格较粗时,一个站点的值完全由某一个格点决定,另一个距离不远的站点可能因为落在了不同格点中心区域而得到差异很大的值,这个跳变没有物理含义。
双线性插值的问题在于平滑过度。它是落在四个格点值之间的凸组合,所以结果永远不会大于四角的最大值、也不会小于四角的最小值——这保证了它不会产生降水负值这类物理上不可能的结果,但同时也意味着它天然抹掉了一部分局地极值。你如果拿双线性插值结果去评估模式对台风中心强度的模拟,那个中心气压可能被平滑得不那么极端了。
正因为各有软肋,实际业务里从来不是"哪个好就永远用哪个",而是根据变量类型和评估目标来选。
3. 动手写代码:Scipy和Xarray两条路线都够用
3.1 准备环境与模拟数据
先准备环境。需要 numpy、pandas、scipy,以及可选但强烈建议安装的 xarray。用 pip 直接装就行:
pip install numpy pandas scipy xarray插值本身只依赖 scipy,xarray 是给那些直接用 NetCDF 文件、懒加载、带坐标标签的工作流准备的。下面这段代码先构造一个模拟的0.5°格点温度场和一组合法站点,后面所有演示都基于这份数据。实际项目中,把模拟部分替换成xarray.open_dataset()或者scipy.io.netcdf_file读取的再分析数据即可。
import numpy as np import pandas as pd # 构造规则格点:北纬20°~40°,东经100°~130°,分辨率0.5° lat = np.arange(20, 40.01, 0.5) lon = np.arange(100, 130.01, 0.5) lon2d, lat2d = np.meshgrid(lon, lat) rng = np.random.default_rng(2024) # 模拟一个带有空间趋势和随机扰动的温度场 temp = ( 20 + 2 * np.sin(np.deg2rad(lat2d - 20)) + 3 * np.cos(np.deg2rad(lon2d - 100)) + rng.normal(0, 0.2, size=lat2d.shape) ) print("网格数据 shape:", temp.shape) # (40, 60),第0维是纬度,第1维是经度 # 模拟站点观测表,每一行是一个站的经纬度 stations = pd.DataFrame({ "name": ["北京", "上海", "广州"], "lat": [39.91, 31.23, 23.13], "lon": [116.39, 121.47, 113.26], }) print(stations)3.2 路线一:Scipy 的 RegularGridInterpolator
Scipy 里有两个常用工具做网格插值:RegularGridInterpolator和griddata。前者要求数据是规则网格、坐标轴是单调数组,性能好、接口也够干净;后者处理的是散乱点,走的是LinearNDInterpolator底层,适合不规则网格数据。这里先说RegularGridInterpolator,因为它最贴近"规则格点插值到站点"这个场景。
from scipy.interpolate import RegularGridInterpolator # 传入参数顺序:points=(lat, lon),对应数据temp的第0维纬度和第1维经度 interp_near = RegularGridInterpolator( (lat, lon), temp, method="nearest", bounds_error=False, # 站点超出网格范围时不报错 fill_value=np.nan # 超范围站点返回NaN ) interp_bilin = RegularGridInterpolator( (lat, lon), temp, method="linear", bounds_error=False, fill_value=np.nan ) # 关键点:查询坐标顺序必须和points顺序一致,也就是(lat, lon),不是常见的(lon, lat) pts = stations[["lat", "lon"]].to_numpy() stations["temp_nearest"] = interp_near(pts) stations["temp_bilinear"] = interp_bilin(pts) print(stations)这个方法里最容易翻车的,就是坐标顺序。很多气象数据习惯把经度放前面、纬度放后面,写(lon, lat),但RegularGridInterpolator的规则是:points里坐标数组的顺序,必须和数据temp的维度顺序完全对应。temp.shape的第0维是lat、第1维是lon,所以points=(lat, lon),查询点也必须是[lat_s, lon_s]。如果你数据维度正好反过来,那points=(lon, lat)也没问题,关键是前后保持一致。
method="nearest"和method="linear"分别对应前文说的两种插值。bounds_error=False配合fill_value=np.nan的意思是:当站点落在网格覆盖范围之外时,不要抛异常,统一给一个 NaN 标记,后续再统一筛查。这个做法在生产环境里很实用,因为你不能保证每一批站点的经纬度都在网格范围内。
3.3 路线二:Xarray 的 interp,一行顶十行
如果你是 NetCDF/NC 数据的老用户,工作流里早就满是xarray.Dataset,那用interp方法是最自然的。它会把坐标标签、维度顺序、单位、缺失值全都帮你处理好,写起来非常简洁。
import xarray as xr # 把模拟数据包装成 DataArray,带上坐标标签 da = xr.DataArray( temp, dims=("lat", "lon"), coords={"lat": lat, "lon": lon}, ) # 单站插值 v_near = da.interp(lat=39.91, lon=116.39, method="nearest").item() v_bilin = da.interp(lat=39.91, lon=116.39, method="linear").item() print("单站最近邻:", v_near, "单站双线性:", v_bilin) # 批量站点插值:把站点经纬度构造成 dims="station" 的 DataArray sta_lat = xr.DataArray(stations["lat"].values, dims="station", name="lat") sta_lon = xr.DataArray(stations["lon"].values, dims="station", name="lon") res = da.interp(lat=sta_lat, lon=sta_lon, method="linear") stations["temp_bilinear_xr"] = res.values print(stations)注意这里sta_lat和sta_lon的长度必须一致,它们分别表示每个站点在lat维和lon维上的插值目标坐标。interp会把两个一维坐标按dims="station"对齐,一次性算出每个站点的插值结果。
xarray 还支持method="cubic",对应三次样条插值,这个后文会提到,但需要谨慎使用,不是所有变量都适合。
3.4 两条路径怎么选
我的建议很直接:数据已经是xr.Dataset或者你习惯用坐标标签管理数据,直接用 xarray 的interp,代码最少、可读性最好。如果你要做的是大批量生产、需要精细控制边界行为、或者希望完全绕开 xarray 的额外依赖,用 Scipy 的RegularGridInterpolator。两条路线的插值数学原理完全一致,结果几乎相同,不存在谁更准的问题。
还有一点值得说:griddata我在这篇文章里刻意不推荐作为首选,因为它走的是 Delaunay 三角剖分道路,对规则网格来说没有利用数据本身的结构优势,速度慢一个数量级,而且会产生边界剖分上的额外麻烦。只有当你面临的是不规则网格(比如站点密不规则的三角网格、观测雷达径向数据)才需要考虑它。
4. 站点数量上来之后:三种性能优化思路
4.1 最忌讳的做法:循环里逐个插值
很多业务数据动辄上万个站点,循环写法是最直观也最容易想到的,但它极慢。原因有两层:第一,Python 层for循环本身就是性能瓶颈;第二,每调用一次插值函数,都要重新做一些参数检查、坐标解析、索引计算,这些固定开销在小批量下可以忽略,但循环一万次就被放大了。
# 这种写法强烈不推荐,数据量一大就肉眼可见地卡 for i in range(len(stations)): v = interp_bilin([(stations["lat"].iloc[i], stations["lon"].iloc[i])])RegularGridInterpolator本身支持传入二维数组作为查询点,一次调用就能批量处理全部站点,完全没必要逐点调用。对于五千、一万个站点,单个批处理调用耗时通常在几十毫秒到几百毫秒之间,比循环快几十倍不止。
4.2 自己写向量化双线性:连函数调用都省了
如果你的站点量级到了十万级,或者插值会被循环调用几千次(比如做蒙特卡洛扰动试验),那连RegularGridInterpolator都嫌重。这时候可以针对规则网格手动实现插值逻辑。规则网格的定位本质上就是searchsorted一趟的事,不需要任何搜索树。
原理很简单:对每个站点,用searchsorted找到它左边那个纬度索引和经度索引,然后取i和i+1两个方向上的四个格点,按距离权重做双线性组合。
def bilinear_interp_on_grid(lat, lon, field, sta_lat, sta_lon): # 定位左侧索引 lat_idx = np.searchsorted(lat, sta_lat, side="right") - 1 lon_idx = np.searchsorted(lon, sta_lon, side="right") - 1 # 边界保护:站点在网格内且不是最右/最上时,保证有右侧/上侧格点可用 lat_idx = np.clip(lat_idx, 0, len(lat) - 2) lon_idx = np.clip(lon_idx, 0, len(lon) - 2) # 计算归一化权重 w_lat = (sta_lat - lat[lat_idx]) / (lat[lat_idx + 1] - lat[lat_idx]) w_lon = (sta_lon - lon[lon_idx]) / (lon[lon_idx + 1] - lon[lon_idx]) # 四个角点的值 f00 = field[lat_idx, lon_idx] f10 = field[lat_idx, lon_idx + 1] f01 = field[lat_idx + 1, lon_idx] f11 = field[lat_idx + 1, lon_idx + 1] return ( f00 * (1 - w_lat) * (1 - w_lon) + f10 * (1 - w_lat) * w_lon + f01 * w_lat * (1 - w_lon) + f11 * w_lat * w_lon )这段代码全部是向量化运算,没有 Python 循环,性能非常好。理解它的关键是:searchsorted给出的是"插入位置",减1之后就指向站点左侧那个索引,右侧索引就是i+1。权重w_lat和w_lon表示站点离左右两侧格点的相对距离,四个格点的贡献权重加起来永远等于1,这正是双线性插值等于凸组合的体现。
4.3 最邻近插值的快速路径
最邻近插值在规则网格上更简单。最直接的做法是用坐标间隔做圆整:
dlat = lat[1] - lat[0] dlon = lon[1] - lon[0] lat_idx = np.clip(np.round((sta_lat - lat[0]) / dlat).astype(int), 0, len(lat) - 1) lon_idx = np.clip(np.round((sta_lon - lon[0]) / dlon).astype(int), 0, len(lon) - 1) values_nearest = temp[lat_idx, lon_idx]注意一个潜在的浮点坑:如果站点经度恰好落在两个格点正中间,round会因为有舍入误差而可能选到错误一侧。更稳妥的做法还是用searchsorted先找到左右索引,再比较站点与左右两侧格点的距离,选近的那一个。代码量也不多,关键是边界安全:
left_lat = np.clip(np.searchsorted(lat, sta_lat, side="left") - 1, 0, len(lat) - 1) right_lat = np.clip(left_lat + 1, 0, len(lat) - 1) near_lat = np.where( np.abs(sta_lat - lat[left_lat]) <= np.abs(sta_lat - lat[right_lat]), left_lat, right_lat )如果网格确实完全规则、站点覆盖范围也都在网格内部,round写法完全够用;只要站点接近边界或者网格坐标存在非均匀间隔,建议换searchsorted版本,鲁棒性不是一个级别。
4.4 不规则网格的兜底方案:cKDTree
上面讲的所有优化都建立在"规则网格"这个前提上。一旦网格不是规则的——比如曲率坐标系、区域加密网格、三角形网格——就不能用searchsorted这套了。这时候做最近邻插值的正确姿势是scipy.spatial.cKDTree:
from scipy.spatial import cKDTree tree = cKDTree(np.column_stack([lon2d.ravel(), lat2d.ravel()])) dist, idx = tree.query(pts, k=1) temp_flat = temp.ravel() values_nearest_unstructured = temp_flat[idx]cKDTree的优势是查询复杂度接近O(logN),网格上万个点也毫无压力。但它默认用欧氏距离,当你处理的区域范围很大、靠近高纬度时,"1°经度的物理长度"和"1°纬度的物理长度"不一致,直接用经纬度坐标算距离会失真。如果要做高纬区域的严格最近邻,建议先把经纬度投影到等距平面上(比如用pyproj转成兰伯特等角投影或极射赤面投影),再建 KDTree。这一点在地面观测站点纬度很高、或区域跨经度很大的场景下尤其要注意。
5. 我踩过的坑:经度范围、NaN传播与边界站点
5.1 经度起点不统一,插值结果全是NaN
这是最坑、也最常见的问题。很多全球模式或再分析产品输出的是0~360°E的经度坐标,比如lon = array([0, 1.5, ..., 358.5]),而你的站点经度表是-180~180°的习惯写法,比如西经120度写成-120.0。两者不统一,插值函数压根找不到-120对应的位置,批量计算后一列 NaN。
解决办法是先把站点经度统一到数据经度的区间。如果数据是0~360,就把站点经度取模:
station_lon_unified = station_lon % 360如果反过来,数据是-180~180、站点经度是240,则:
station_lon_unified = (station_lon + 180) % 360 - 180建议把统一坐标的代码放在插值之前,并且加个打印检查,确认转换前后站点经度没有异常跨越。这个步骤很多人忽略,但它能省下大量后续排查时间。
5.2 纬度从北到南排列,Scipy直接报错
RegularGridInterpolator对坐标数组要求很严格:必须是严格递增的。有些卫星资料、模式输出的纬度维度习惯于从北极向南极排列,也就是lat[0]是90、lat[-1]是-90。你把它直接丢进RegularGridInterpolator,它会抛异常或者给你一个错得离谱的插值结果。
解决办法是检测后用[::-1]翻转纬度轴,同时翻转数据场的第0维:
if lat[1] < lat[0]: lat = lat[::-1] temp = temp[::-1, :]翻转之后,数据和坐标必须同步翻,别只改一个。这种错误不会产生警告,int类型翻转一切正常,只有结果看起来不对劲,所以最好在读取数据后立即打印lat[0], lat[-1]做断言检查。
5.3 NaN传播比你想的更严重
格点数据经常带缺测:海洋上的海表温度、被地形遮蔽的下层大气变量、云遮挡的卫星反演产品,都会产生 NaN。问题是,双线性插值一旦碰到四个角点里混入一个 NaN,由于权重相乘,结果大概率也是 NaN。尤其当网格上有一小片缺口、站点又恰好落在缺口附近时,一个 NaN 能污染一大片站点。
处理策略要看业务目标。如果你只是临时把浓度场插到站点做一张对比散点图,可以直接把 NaN 站点筛掉:
valid = ~np.isnan(stations["temp_bilinear"]) print("有效站点占比:", valid.mean()) stations_valid = stations[valid]但如果你的下游是客观分析、同化系统,缺测站点不能直接丢弃,就得考虑先用有效格点填充 NaN。最简单的做法是先跑一遍最近邻插值,用最近邻结果填充格点场中的 NaN,然后再做双线性插值:
from scipy.interpolate import griddata # 提取有效格点 valid_mask = ~np.isnan(temp) pts_valid = np.column_stack([lon2d[valid_mask].ravel(), lat2d[valid_mask].ravel()]) values_valid = temp[valid_mask].ravel() # 用最近邻方式填充完整网格 lon_full, lat_full = np.meshgrid(lon, lat) temp_filled = griddata( pts_valid, values_valid, (lon_full, lat_full), method="nearest" )注意griddata的method="nearest"会把最近有效值填进缺口,防止 NaN 扩散。这样的填充场再做双线性插值到站点,结果会更稳定。但这个方案只适用于缺口不大、周围有效值密度足够的情况,如果缺测面积太大,填充出来的值已经不具备真实的物理代表性,那时候要回头检查数据源了。
5.4 边界站点:超范围不等于报错
bounds_error=False很贴心,但代价是你可能静默地拿到一堆 NaN。站点如果稍微落在网格边界外,比如网格北界是 40°N,站点是 40.2°N,最邻近插值和双线性插值都会返回 NaN,除非你用fill_value指定其他默认值。很多人在下游算平均值时没检查 NaN,直接得到 0 或者空结果,这种 bug 排查起来非常耗费时间。
我的习惯是每次插值完都做一次完整性审计,至少打印缺失数量和在网格外的站点明细:
nan_mask = stations["temp_bilinear"].isna() if nan_mask.any(): print("以下站点在网格范围内无有效插值结果:") print(stations.loc[nan_mask, ["name", "lat", "lon"]])5.5 跨越180°经线的区域要单独处理
如果研究区域跨越东西经边界(比如包含白令海、南太平洋岛弧),站点经度可能一边是 179.9,另一边是 -179.8,物理距离很近,但数值上相差 359.7°。RegularGridInterpolator不知道地球是圆的,它会试图在 179.9 和 -179.8 之间做线性插值,结果自然是垃圾。
这种场景的处理思路是:要么把数据的经度坐标整体统一成0~360区间,同时把站点经度也取模进去,让原本在 -180 附近的站点变成 180 附近;要么把数据沿经度方向做一次 roll 平移,把 180° 经线挪到数据边界,而不是让插值跨越它。具体选哪种,取决于你的数据范围和站点分布,但核心原则是:不要让插值算法去处理它不理解的不连续边界。
6. 到底怎么选:连续量、分类量和极端值的不同答案
6.1 一张表看清两个选项的适用边界
| 维度 | 最邻近插值 | 双线性插值 |
|---|---|---|
| 数学操作 | 取最近格点原值 | 四角格点按距离加权 |
| 结果平滑性 | 差,呈阶梯状 | 好,平滑过渡 |
| 是否保留格点原值 | 是 | 否,是加权合成新值 |
| 极端值保留 | 好 | 差,容易被平滑 |
| 计算成本 | 极低 | 低 |
| 分类变量适用性 | 完全适用 | 不适用 |
| 连续物理量适用性 | 可用但有系统偏差 | 推荐 |
| 缺测传播 | 受最近格点影响 | 受周围四格点影响,更易扩散 |
这里最值得强的一句话:双线性插值的结果是四个角值的凸组合,所以它不会产生比四角最小值更低、比四角最大值更高的新值。这带来了一个好处——降水、湿度这种非负变量,只要四角都非负,插值结果一定非负,不会出现物理上不可能的负值。但代价是它磨掉了峰值,做极端事件检验时不适合。
6.2 变量类型是最优先的决策依据
碰到土地利用类型、土壤质地分类、天气现象编码、云量类型、冻土状态这类分类变量,别无选择,只能用最邻近插值。做最近邻时把数值编码当成连续量去平均,结果会产生毫无意义的"第3.7类土壤",这种错误常出现在把双线性插值一股脑套用到所有变量的代码里。
连续变量也不是全都适合双线性。气温、气压、位势高度这类空间连续性好、变化平缓的变量,双线性插值表现非常好,几乎是行业默认。降水、对流有效位能、云顶亮温这类局地性极强、空间突变明显的变量,双线性会明显平滑掉强中心,导致模式极值被低估。对这类变量,很多业务团队的方案是"评估平均态用双线性,评估极端值用最邻近",两种结果都保留,分析差异来源。
6.3 高阶插值不是万金油
Scipy 和 xarray 都提供了三次样条插值(method="cubic"),它能给出比双线性更平滑的曲面,但代价是可能在格点之间的区域产生超出原始场值域的过冲,也就是出现不真实的虚假极值。对温度这类平滑变量,过冲幅度小,通常可以接受;对降水、湿度、气溶胶浓度这类非负变量,过冲可能产生负浓度,这在物理上完全不可行。如果你需要更高阶的光滑,请务必在插值后做物理约束检查,至少确认结果最小值不小于0。
还有一种常见的"高阶"做法是:先把格点场投影到更高分辨率的网格,再插到站点。这个中间步骤并不会增加真实信息,反而把插值误差又转了一道手。我的经验是,能一步到位就不要做两级插值。
6.4 插值解决不了代表性问题
最后说一个经常被忽略、但实际工作中影响最大的点:插值只是几何操作,它解决不了站点代表性不足的问题。
山区尤其明显。模式格点代表的是一个格点区域内的平均状态,网格内如果有几百米的高差,格点温度代表的是那个平均高度的温度,而站点可能在山顶也可能在山谷。即使双线性插值在空间位置上精确命中了站点经纬度,它也没有办法把模式平均地形与站点真实海拔之间的差异订正掉。复杂地形区,温度偏差达到几度是常有的事。
所以,你要是看到插值后站点检验偏差系统性偏大,先别急着换插值算法。先检查海拔差、先检查数据经度范围、先检查 NaN 和边界站点,这些常规因素排干净了,再考虑是不是插值方法本身不合适。这件事我在项目里栽过不止一次,写出来希望能帮你少走几小时弯路。
实操层面的最终建议很简单:低纬到中纬、地形平坦区域,连续变量无脑双线性;分类变量、极端事件检验、粗网格数据,用最邻近;两者的结果都应纳入后期的敏感性分析。插值只是管道,不是决策本身,理解了它在哪里会引入误差,你才知道自己评估结论里有多少是信号的贡献、多少是插值的噪声。