首先，你的全局欧氏距离计算思路是对的，但面对1200×1500这样的大规模数组，确实可以通过空间索引或分块策略进一步提速。下面我会先针对你提到的四等分缩小范围思路给出实现指导，再推荐一种更成熟的高效方案。

一、分块缩小范围的实现思路 #

你提到的分块思路核心是逐步缩小搜索范围，但难点在于判断目标点所在的非矩形区域。其实我们可以换个角度：不用严格判断点是否在多边形子区域内，而是先基于坐标的数值范围做初步分块过滤——因为你的坐标数组看起来是规则网格排列的（纬度从54.496到54.531递增，经度从21.77到21.70递减），本身具备可利用的网格结构！

具体步骤： #

1. 拆解坐标矩阵：
   先把coords数组拆分为纬度矩阵和经度矩阵，方便快速获取任意子块的坐标区间：

   lat = coords[:, :, 0]
   lon = coords[:, :, 1]
   

2. 迭代分块筛选：
   以整个网格为初始范围，每次将当前行、列范围各拆分为两半，得到4个子块。计算每个子块的坐标包围盒（最小/最大纬度、经度），判断目标点pt是否落在该包围盒内，保留所有包含目标点的子块（若两个相邻子块都包含，合并范围避免漏点）。
3. 精确计算最近点：
   当子块的行/列范围缩小到16以内时，取出该子块的所有坐标，计算与pt的欧氏距离，找到最近点后映射回原数组的全局索引。

代码示例（简化版）： #

import numpy as np

def find_nearest_block(coords, pt, target_size=16):
    rows, cols = coords.shape[0], coords.shape[1]
    lat, lon = coords[:, :, 0], coords[:, :, 1]
    
    # 初始搜索范围
    row_range = [0, rows-1]
    col_range = [0, cols-1]
    
    while (row_range[1] - row_range[0] + 1) > target_size or (col_range[1] - col_range[0] + 1) > target_size:
        # 拆分行范围为两半
        mid_row = (row_range[0] + row_range[1]) // 2
        row_blocks = [[row_range[0], mid_row], [mid_row+1, row_range[1]]]
        # 拆分列范围为两半
        mid_col = (col_range[0] + col_range[1]) // 2
        col_blocks = [[col_range[0], mid_col], [mid_col+1, col_range[1]]]
        
        # 筛选包含目标点纬度范围的行块
        valid_row = []
        for r_start, r_end in row_blocks:
            block_lat_min = lat[r_start:r_end+1, :].min()
            block_lat_max = lat[r_start:r_end+1, :].max()
            if block_lat_min <= pt[0] <= block_lat_max:
                valid_row.append([r_start, r_end])
        # 若两个行块都包含，保留完整范围避免漏点
        row_range = valid_row[0] if len(valid_row)==1 else [row_blocks[0][0], row_blocks[1][1]]
        
        # 同理筛选列块
        valid_col = []
        for c_start, c_end in col_blocks:
            block_lon_min = lon[:, c_start:c_end+1].min()
            block_lon_max = lon[:, c_start:c_end+1].max()
            if block_lon_min <= pt[1] <= block_lon_max:
                valid_col.append([c_start, c_end])
        col_range = valid_col[0] if len(valid_col)==1 else [col_blocks[0][0], col_blocks[1][1]]
    
    # 在缩小后的范围内计算最近点
    sub_coords = coords[row_range[0]:row_range[1]+1, col_range[0]:col_range[1]+1]
    dist_sq = ((sub_coords - pt)**2).sum(axis=2)
    sub_r, sub_c = np.unravel_index(dist_sq.argmin(), dist_sq.shape)
    # 映射回原数组索引
    return row_range[0]+sub_r, col_range[0]+sub_c

二、更成熟的高效方案：KD-Tree最近邻搜索 #

如果你的坐标数组存在局部不规则情况，或者想省去手动分块的麻烦，推荐使用scipy.spatial.KDTree——它专门为高维空间的最近邻搜索做了优化，构建树的时间复杂度为O(n log n)，单次查询仅需O(log n)，处理1.8e6个点的效率极高。

代码示例： #

from scipy.spatial import KDTree
import numpy as np

# 将二维网格展平为N×2的数组
coords_flat = coords.reshape(-1, 2)
# 构建KD-Tree
tree = KDTree(coords_flat)
# 查询最近点的距离和索引
dist, idx = tree.query(pt, k=1)
# 映射回原网格的行列索引
r1 = idx // coords.shape[1]
c1 = idx % coords.shape[1]

优势说明： #

• 无需手动实现分块逻辑，底层算法已做极致优化
• 即使坐标不是严格规则网格，也能高效工作
• 支持批量查询多个目标点，扩展性更强

三、原始方法的小优化 #

如果不想引入scipy依赖，也可以对原始numpy计算做一点简化——利用numpy的广播机制和np.unravel_index直接获取行列索引，比手动计算索引更简洁：

dist_sq = np.sum((coords - pt)**2, axis=2)
r1, c1 = np.unravel_index(dist_sq.argmin(), dist_sq.shape)

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