如何使用MATLAB找到离给定坐标最近的点？

Question

我需要使用Matlab解决一个最小化问题，想知道哪种方法最简单。我考虑的所有解决方案都需要大量编程工作。

假设我有一个经度/纬度坐标点（A,B），我需要在一个经度/纬度坐标地图中搜索距离这个点最近的点。

特别是，纬度和经度数组是两个2030x1354元素（1公里距离）的矩阵，而要做的就是找到这些矩阵中最小化到坐标（A,B）的距离的唯一索引，即找到最接近给定坐标（A,B）的值。

非常感谢您的任何帮助。

谢谢！

- cardogar

---

你有两个由一维值组成的二维矩阵，对吧？ - macduff

只是确认一下...你是否有第三个矩阵，它在每个纬度/经度处都给出一个值，并且你正在尝试查找该矩阵中任意纬度/经度对应的正确值？如果是这样，那么你需要使用interp2或griddata，具体取决于你的纬度/经度矩阵的特定情况。 - Peter

是的，Macduff，我有带有1D值的2D矩阵。 - cardogar

不，Peter，我没有任何带有额外值的第三个矩阵。但还是谢谢！ - cardogar

Answer 1

设Lat和Long分别表示纬度和经度矩阵，则有：

dist2=sum(bsxfun(@minus, cat(3,A,B), cat(3,Lat,Long)).^2,3);
[I,J]=find(dist2==min(dist2(:)));

I 和 J 包含与最近点对应的A和B中的索引。请注意，如果有多个答案，I和J将不是标量值，而是向量。

- Mohsen Nosratinia

---

谢谢您的回答！我需要补充之前的信息，即所需的（A，B）坐标将位于整个地球上，包括极点。因为它们对应于卫星轨道的跟踪点。根据Rody Oldenhuis的回答，这种方法对此目的将不起作用。 - cardogar

我已经测试过了，它运行得非常好。我将检查杆的情况，但到目前为止，您提出的方法正好满足我的需求。谢谢。 - cardogar

Answer 2

这总是一个有趣的问题 :)

首先，Mohsen Nosratinia的回答可以，只要您不需要知道实际距离，可以绝对确定您永远不会靠近极地地区，也永远不会靠近±180°子午线。对于给定的纬度，-180°和+180°经度实际上是同一点，因此仅查看角度差异是不足够的。这在极地地区将是更大的问题，因为那里的大经度差异对实际距离的影响将更小。

球面坐标在导航、制图等方面非常有用和实用。然而，在表面距离等空间计算方面，球面坐标实际上非常麻烦。虽然可能使用角度直接进行此类计算，但我个人认为这不是很实用：您通常需要具有球面三角学的强背景，并且需要具有相当的经验才能了解其许多缺陷——通常有不稳定性或“特殊点”需要解决（例如极点），需要考虑四象限歧义性，因为您引入了三角函数等。

例如，如果您将纬度和经度转换为3D笛卡尔X、Y、Z坐标，则您的距离问题就非常简单，然后通过简单的公式找到距离：

距离（a，b）= R·arccos(a/|a|·b/|b|) 其中a和b是球体上的两个笛卡尔向量。注意| a |= | b |= R，其中R = 6371是地球的半径。

在MATLAB代码中：

% Some example coordinates (degrees are assumed)
lon = 360*rand(2030, 1354);
lat = 180*rand(2030, 1354) - 90;

% Your point of interest
P = [4, 54];

% Radius of Earth
RE = 6371;

% Convert the array of lat/lon coordinates to Cartesian vectors
% NOTE: sph2cart expects radians
% NOTE: use radius 1, so we don't have to normalize the vectors
[X,Y,Z] = sph2cart( lon*pi/180,  lat*pi/180, 1);

% Same for your point of interest    
[xP,yP,zP] = sph2cart(P(1)*pi/180, P(2)*pi/180, 1);

% The minimum distance, and the linear index where that distance was found
% NOTE: force the dot product into the interval [-1 +1]. This prevents 
% slight overshoots due to numerical artifacts
dotProd = xP*X(:) + yP*Y(:) + zP*Z(:);
[minDist, index] = min( RE*acos( min(max(-1,dotProd),1) ) );

% Convert that linear index to 2D subscripts
[ii,jj] = ind2sub(size(lon), index)

如果您坚持跳过转换为笛卡尔坐标并直接使用纬度/经度，您将不得不使用Haversine公式，例如在此网站中所述的方式，这也是映射工具箱中distance()使用的方法。

现在，所有这些都适用于整个地球，前提是您认为平滑的球形地球足够准确。如果您想包括地球的扁率或一些更高阶的形状模型（或者天哪，包括地形的距离），则需要进行更复杂的操作。但我不认为这是您的目标：）

PS-如果您把我写的一切都写出来，我不会感到惊讶，您可能会重新发现Haversine公式。我只是更喜欢能够仅从第一原理计算沿着球体的简单距离，而不是从您很久以前植入头脑中的某个黑匣子公式计算。

- Rody Oldenhuis

---

尝试使用这两种方法，你几乎总是会得到相同的距离。但并不总是如此，例如在我的特定问题中，Mohsen的方法结果更好。尝试多次迭代这两种方法，你会发现有时Mohsen的模型提供了更好的结果。为什么呢？我也不知道。 - cardogar

如果你尝试几次，你就会明白我的意思。for i=1:15 lon=360rand(2030,1354);lat=180rand(2030,1354)-90;P=[4,54];[X,Y,Z]=sph2cart(latpi/180,lonpi/180,1);[xP,yP,zP]=sph2cart(P(1)pi/180,P(2)pi/180,1);dotProd=xPX(:)+yPY(:)+zPZ(:);[~,index]=min(6371acos(min(max(-1,dotProd),1)));[I,J]=ind2sub(size(lon),index);distdim(distance(P(1),P(2),lat(I,J),lon(I,J)),'deg','kilometers')
dist2=sum(bsxfun(@minus,cat(3,P(1),P(2)),cat(3,lat,lon)).^2,3);[I,J]=find(dist2==min(dist2(:)));distdim(distance(P(1),P(2),lat(I,J),lon(I,J)),'deg','kilometers') end - cardogar

@cardogar：那是因为你在distance()和Mohsen的方法中都把P的纬度/经度顺序搞错了 :) 当我纠正后，发现我的方法产生的距离与distance()相同，而Mohsen的方法产生的距离通常与distance()找到的距离相差甚远，此外，他的最小距离通常比我找到的还要大。 - Rody Oldenhuis

我完全没有意识到这些点在一个球体上。 - Mohsen Nosratinia

解决了。我为找错而疯狂，不明白为什么你没有看到同样的问题。我使用相同的值、相同的公式……但结果不同！问题与变量的精度有关。我的P值被读取为单精度数字，如果你使用[xP,yP,zP]=sph2cart(single(P(2)*pi/180),single(P(1)*pi/180),1);，那么你可能会看到我的结果。但好吧，最终只需使用[xP,yP,zP]=sph2cart(double(P(2)*pi/180),double(P(1)*pi/180),1)就解决了。非常感谢你的所有帮助和耐心。我真的很感激你在这个问题上花费的时间。 - cardogar