大经纬度网格与小经纬度网格匹配对齐的提速方法问询

现有一组表示地形经纬度的大型等间距数组：

• 经度数组lons：元素以0.005度为间隔分布，示例片段：

lons[0:10] = [-130.0, -129.995, -129.99, -129.985, -129.98, -129.975, -129.97, -129.965, -129.96, -129.955]

• 纬度数组lats：元素同样以0.005度为间隔分布，示例片段：

lats[0:10] = [55.0, 54.995, 54.99, 54.985, 54.98, 54.975, 54.97, 54.965, 54.96, 54.955]

另有一个[m,n]维度的不规则经纬度网格数据集（实际间距约25米），其范围完全处于上述大网格内，需要将小网格的每个点匹配到大网格中最近邻的地形值。目前采用双重循环遍历并通过argmin查找的方法实现，但速度极慢，寻求提速的工具或方法。

方法1：利用等间距网格特性直接计算索引 #

因为大网格是固定间隔0.005度的规则网格，无需遍历或搜索，直接通过坐标差计算对应索引即可，这是效率最高的方案：

import numpy as np

# 假设小网格经度数组为small_lons，纬度为small_lats（形状[m,n]），地形数组为topo（与大网格维度匹配）
lon_step = 0.005
lat_step = 0.005

# 计算经度索引：四舍五入到最近的整数索引
lon_indices = np.round((small_lons - lons[0]) / lon_step).astype(int)
# 计算纬度索引
lat_indices = np.round((small_lats - lats[0]) / lat_step).astype(int)

# 直接通过索引提取匹配的地形值
matched_topo = topo[lat_indices, lon_indices]

该方法时间复杂度为O(m*n)，完全规避Python循环带来的开销。

方法2：使用NumPy向量化替代循环 #

若需兼容非规则间隔场景，可利用NumPy的广播机制替代Python双重循环，依托底层C优化提升速度：

import numpy as np

# 将大网格经纬度转为二维网格
lon_grid, lat_grid = np.meshgrid(lons, lats)
# 计算小网格每个点到大网格所有点的距离平方（省略开根号，不影响argmin结果）
dist_sq = (small_lons[..., np.newaxis, np.newaxis] - lon_grid)**2 + (small_lats[..., np.newaxis, np.newaxis] - lat_grid)**2

# 找到最近邻索引并提取地形值
min_indices = np.argmin(dist_sq, axis=(-2, -1))
matched_topo = topo.flat[min_indices].reshape(small_lons.shape)

方法3：使用SciPy的KDTree进行快速最近邻搜索 #

针对超大规模网格场景，可通过KDTree构建索引实现批量快速查询：

import numpy as np
from scipy.spatial import KDTree

# 构建大网格坐标点集
lon_grid, lat_grid = np.meshgrid(lons, lats)
grid_points = np.column_stack((lon_grid.flatten(), lat_grid.flatten()))
kdtree = KDTree(grid_points)

# 将小网格坐标转为二维点集并查询最近邻
small_points = np.column_stack((small_lons.flatten(), small_lats.flatten()))
_, indices = kdtree.query(small_points, k=1)

# 提取地形值并恢复原形状
matched_topo = topo.flat[indices].reshape(small_lons.shape)

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