使用MATLAB查找给定点的最近唯一经纬度对问题排查

我有两组大型纬度（Lat）和经度（Lon）向量，想要找到与给定点[lat_deg, lon_deg]距离最短的唯一一对经纬度。当前使用的代码如下：

P = ([lat_deg, lon_deg]);
PQ = [Lat, Lon];
[k,dist] = dsearchn(P,PQ);

运行后得到所有点的距离，且向量k全为1。需要指导该函数是否适用，若适用如何修正；若不适用，应使用什么函数。

示例向量：

Lat                      Lon
39.2591200000000    -85.9394200000000
39.2591300000000    -85.9392000000000
39.2590800000000    -85.9406300000000
39.2593500000000    -85.9406200000000
39.1949800000000    -85.9633400000000
39.1954200000000    -85.9633500000000
39.1954200000000    -85.9633500000000
39.1963300000000    -85.9633600000000
39.1957400000000    -85.9678800000000
39.1959300000000    -85.9682400000000

P=39.2005981000000  -85.9045842000000

1. 修正dsearchn的参数顺序 #

dsearchn的正确调用格式是dsearchn(X,Y)，其中：

• X是已知点集（每行对应一个经纬度点，即你的PQ = [Lat, Lon]）
• Y是待查询的目标点（即你的P）

你当前参数顺序写反，导致函数逻辑错误，才会出现k全为1的结果。修正后的代码：

P = [lat_deg, lon_deg];
PQ = [Lat, Lon];
[k, dist] = dsearchn(PQ, P);

运行后：

• k返回距离目标点最近的点在PQ中的索引
• dist返回该点与目标点的欧氏距离

2. 经纬度距离的准确计算 #

dsearchn默认用欧氏距离，小范围经纬度数据误差不大，但跨越大范围时建议用球面距离计算，两种方案可选：

方案一：转笛卡尔坐标后用dsearchn #

先将经纬度转换为三维笛卡尔坐标（基于地球半径），再用dsearchn计算：

% 转弧度
lat_rad = deg2rad(Lat);
lon_rad = deg2rad(Lon);
p_lat_rad = deg2rad(lat_deg);
p_lon_rad = deg2rad(lon_deg);

% 地球半径（单位：km）
R = 6371;
% 转换为笛卡尔坐标
X = R .* cos(lat_rad) .* cos(lon_rad);
Y = R .* cos(lat_rad) .* sin(lon_rad);
Z = R .* sin(lat_rad);
PQ_cart = [X, Y, Z];

p_X = R * cos(p_lat_rad) * cos(p_lon_rad);
p_Y = R * cos(p_lat_rad) * sin(p_lon_rad);
p_Z = R * sin(p_lat_rad);
P_cart = [p_X, p_Y, p_Z];

% 查找最近点
[k, dist] = dsearchn(PQ_cart, P_cart);

方案二：用Haversine公式直接计算球面距离 #

手动计算每个点与目标点的球面距离，再取最小值：

% 转弧度
lat_rad = deg2rad(Lat);
lon_rad = deg2rad(Lon);
p_lat_rad = deg2rad(lat_deg);
p_lon_rad = deg2rad(lon_deg);

% Haversine公式计算球面距离
d_lat = lat_rad - p_lat_rad;
d_lon = lon_rad - p_lon_rad;
a = sin(d_lat/2).^2 + cos(p_lat_rad) .* cos(lat_rad) .* sin(d_lon/2).^2;
c = 2 * atan2(sqrt(a), sqrt(1-a));
R = 6371; % 地球半径（km）
distances = R * c;

% 找最小距离的索引
[min_dist, k] = min(distances);
% 获取最近点的经纬度
closest_lat = Lat(k);
closest_lon = Lon(k);

3. 处理重复点确保唯一性 #

如果你的经纬度向量中有重复点（比如示例中第6、7行），可以先对PQ去重，保证得到唯一的最近点：

% 去重，保留唯一的经纬度对
[unique_PQ, ~, idx] = unique(PQ, 'rows');
% 计算去重后点的最近点索引
[k_unique, dist] = dsearchn(unique_PQ, P);
% 若需要映射回原数组的索引
original_k = find(idx == k_unique);

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