如何计算一组线和一个点之间的（最小）距离？

Posted 2023-03-24

【中文标题】如何计算一组线和一个点之间的（最小）距离？

【英文标题】：How to calculate the (smallest) distance between a set of lines and a point?

【发布时间】：2019-10-17 13:16:33

【问题描述】：

所以基本上我会尝试标记一些 GPS 点是否非常靠近道路。

我有 lat-lon 坐标的向量，当连接到向量中的后续/先前行时，形成一条代表道路的线（描绘为 matplotlib 彩色图） 我也有简单的地理点（经纬度）（用黑色表示为 matplotlib 散点图）

我想标记一个点是否在道路附近（例如，在 0.001 弧度距离内） - 为此，我猜需要计算一个点到这组向量的最近距离。

#example vector 1
[[-84.3146272, 33.7741084], [-84.3145183, 33.7741196]]
#example vector 2
[[-84.4043106, 33.7700542], [-84.4045421, 33.770055]]

#example point to predict wether it will be near one of these two lines
[-84.31106, 33.73887]

如何解决这个问题？我想不出解决这个问题的方法，但是看情节时似乎很简单……有没有可以提供帮助的库？

【问题讨论】：

检查 sicpy cdist 或 nnumpy.linalg.norm

你有十进制度坐标。 0.001 弧度不能转换为平面距离，除非纬度差异很小，因为经度间距会随着向极地移动而减小。如果您的北纬和南纬距离很近，您会变得很接近

我更新了我的帖子，并在算法中添加了欧几里得近似值，这可能会运行得更快并提供可比较的结果。选择您喜欢的方法。我还没有测试过第二个。

【参考方案1】：

我假设您正在使用地球的球形模型。那么我想你所说的“线”实际上是大圆的弧段（球面几何的直线）。换句话说，您在球体表面上有一个点p1，经纬度坐标为[lat1, lon1]，球体上有另一个点p2，经纬度坐标为[lat2, lon2]。通过p1 和p2 的球面上被认为是“直线”的是球面与通过球心的平面相交得到的唯一圆（称为大圆）。两点p1 和p2。然后，你所说的“线”是这个大圆的两条圆弧中较小的一条，以p1 和p2 两点为界。

您要计算的距离（以弧度为单位）可以是从第三个点p 与经纬坐标[lat, lon] 到由两个点p1 和p2 确定的弧段的距离.所述距离应为从大圆经过p并垂直于p1和p2的大圆的弧长。这个垂直的大圆是由球面与垂直于p1和p2的大圆平面并通过点p和球心的平面的交点确定的。如果垂直大圆与大圆弧p1 p2的交点是圆弧段p1 p2内的一点h，那么大圆弧p h的长度就是所寻求的距离。但是，如果h 在弧线p1 p2 之外，则寻求的距离是p p1 或p p2，以较小者为准。

这是一些计算点和弧间隔之间最短距离的 Matlab 代码：

lat_lon = [lat, lon];
lat_lon1 = [lat1, lon1];
lat_lon2 = [lat2, lon2];

function dist = dist_point_2_road(lat_lon, lat_lon1, lat_lon2)

   % you may have to convert lat-long angles from degrees to radians 

   % First, convert from lat-long coordinates to 3D coordinates of points on the unit    
   % sphere. Since Earth's radius cancels out in our computations, we simply assume it 
   % is R = 1 
   lat = lat_lon(1);
   lon = lat_lon(2);

   lat1 = lat_lon1(1);
   lon1 = lat_lon1(2);

   lat2 = lat_lon2(1);
   lon2 = lat_lon2(2);

   p1 = [ cosd(lat1)*cosd(lon1),  cosd(lat1)*sind(lon1),  sind(lat1) ]; %cosd = cos(degrees)
   p2 = [ cosd(lat2)*cosd(lon2),  cosd(lat2)*sind(lon2),  sind(lat2) ]; %sind = sin(degrees)
   p = [ cosd(lat)*cosd(lon),  cosd(lat)*sind(lon),  sind(lat) ];

   % n12 is the unit vector perpendicular to the plane of the great circle 
   % determined by the points p1 and p2  
   n12 = cross(p1, p2);
   n12 = n12 / sqrt(dot(n12, n12));
   sin_of_dist = dot(p, n12); % sine of the angle that equals arc-distance 
                              % from point p to the great arc p1 p2  

   dist = pi/2 - acos(abs(sin_of_dist)); % acos = arccos, abs() = absolute value
          % dist is the shortest distance in radians from p to the 
          % great circle determined by the points p1 and p2

   n1 = cross(p1, p); 
   n1 = n1 / sqrt(dot(n1, n1));
   % unit normal vector perpendicular to the great-arc determined by p and p1

   n2 = cross(p, p2);
   n2 = n1 / sqrt(dot(n2, n2));
   % unit normal vector perpendicular to the great-arc determined by p and p2

   if dot(n12, n1) < 0 % if the angle of spherical triangle p p1 p2 at vertex p1 is obtuse 
      dist = acos(dot(p, p1)); % the shortest distance is p p1 
   elseif dot(n12, n2) < 0 % if the angle of spherical triangle p p1 p2 at vertex p2 is obtuse 
      dist = acos(dot(p, p2)); % the shortest distance is p p2 
   end

   % the function returns the appropriate dist as output 

end

您可以对形成道路的弧间隔序列进行迭代，并选择到弧间隔的最小距离。

根据这个计算，点到第一个“向量1”的距离为 0.0000615970599633145 弧度，到第二个“向量 2”的距离是 0.00162840939265068 弧度。该点最接近向量 1 内的点，但对于向量 2，它最接近弧间隔的端点之一。

编辑。 现在，如果要使用欧几里得（平面）近似，忽略地球曲率，可能需要将经纬坐标转换为欧几里得平面近似坐标。为避免出现任何地图特定坐标，可能会尝试将纬度与经度坐标进行对比。在赤道附近可能没问题，但是越靠近两极，这些坐标在表示距离数据时就越不准确。这是因为越靠近两极，沿固定纬度的距离远小于沿固定经度的距离。这就是为什么我们需要纠正这种差异。这是通过在经纬度坐标中使用球体上的黎曼度量来完成的，或者只是通过查看球体上给定点附近的纬度和经度圆的 3D 几何图形来完成。

lat_lon = [lat, lon];
lat_lon1 = [lat1, lon1];
lat_lon2 = [lat2, lon2];

%center of approximate Euclidean coordinate system is point p 
% with lat_long coordinates and the scaling coefficient of longitude, 
% which equalizes longitude and latitude distance at point p, is

a = cosd(lat_long(1));  

function  x = convert_2_Eucl(lat_long1, lat_long, a)
   x = [lat_long1(1) - lat_long(1),  a*(lat_long1(2) - lat_long(2))];
end

% convert from lat-long to approximate Euclidean coordinates
x1 = convert_2_Eucl(lat_long1, lat_long, a);
x2 = convert_2_Eucl(lat_long2, lat_long, a);

function dist = dist_point_2_road(x1, x2)

   dist = dot(x1, x1) * dot(x2 - x1, x2 - x1) - dot(x1, x2 - x1)^2 ;
   dist = sqrt( dist / ( dot(x2 - x1, x2 - x1)^2) );
   % dist is the distance from the point p, which has Eucl coordinates [0,0] 
   % to the straight Euclidean interval x1 x2 representing the interval p1 p2

   if dot(x1, x2 - x1) > 0
      dist = sqrt( dot(x1, x1) );
   elseif dot(x2, x1 - x2) > 0
      dist = sqrt( dot(x2, x2) );
   end

end

备注：后一个函数计算距离，但只计算dist^2可能同样方便，避免计算平方根sqrt以加快性能。关于 dist^2 的测量也应该同样有效。

您可以选择所需的函数，球面函数或近似欧几里得函数。后者可能更快。您可以选择删除平方根并计算距离平方，以使事情变得更快。

我写的很匆忙，所以可能有一些不准确的地方。

【讨论】：

hmm.. 虽然有些道路可能一直通向罗马，但在这种情况下，它们足够短且局部性足以忽略地球的曲率，以确定哪条道路最接近.

@user2297550 是的，我知道，但您的数据是球坐标。要使用欧几里得平面近似，您需要对其进行转换。最重要的是，地图上的直线取决于地图，如果地图发生变化，曾经的直线不再是直线。我写的是地图独立。你有球坐标到欧几里得的转换器吗？您不能只根据纬度绘制经度，因为在北欧国家和阿拉斯加，经度的变化比纬度要快得多。

感谢您的澄清。我不是提问者；只是一个路人，但你的澄清会对每个人都有帮助。

@2297550 抱歉，现在我注意到了不同之处。我实际上添加了一个欧几里得近似值，它可能会更快。感谢您的评论。

你真好。继续添加很棒的答案:-)