python - 如何计算一组线和一个点之间的(最小)距离？

Question

所以基本上我试图标记一些 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]

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

Answer 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是一个点h弧段内p1 p2 , 那么大圆弧的长度 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，它最接近弧间隔的端点之一。

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

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

写的比较匆忙，可能有一些不准确的地方。