Как рассчитать (наименьшее) расстояние между набором линий и точкой? - PullRequest
1 голос
/ 01 июня 2019

Так что в основном я пытаюсь пометить, если некоторые точки GPS находятся очень близко к дороге.

У меня есть векторы широковых координат, которые при соединении с последующей / предшествующей строкой в ​​векторе образуют линию, представляющую дорогу (изображенную в виде графиков matplotlib в цветах) У меня также есть простые географические точки (лат-лон) (изображены как диаграмма рассеяния matplotlib в черном) two muppets

Я хотел бы отметить, находится ли точка рядом с дорогой (например, в пределах расстояния 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]

Как можно решить эту проблему? Я не могу придумать, как решить эту проблему, но, глядя на сюжет, это кажется простым ... Есть ли библиотека, которая могла бы помочь?

1 Ответ

2 голосов
/ 01 июня 2019

Я предполагаю, что вы работаете со сферической моделью Земли.Тогда я полагаю, что то, что вы называете «линиями», на самом деле являются дуговыми сегментами больших кругов (прямые сферической геометрии).Другими словами, у вас есть одна точка 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 она ближе всего к одному из концов интервала дуги.

Редактировать. Теперь, если кто-то хочет использовать евклидово (плоское) приближение, игнорируя кривизну Земли, может потребоваться преобразовать широтные координаты в евклидовы координаты плоского приближения.Чтобы избежать каких-либо специфических для карты координат, можно попытаться построить график широты относительно координат долготы.Это может быть хорошо для экватора, но чем ближе к полюсам, тем более неточными эти координаты становятся при представлении данных о расстоянии.Это связано с тем, что ближе к полюсам расстояние вдоль фиксированной широты намного меньше расстояния вдоль фиксированной долготы.Вот почему мы должны исправить это несоответствие.Это делается с помощью римановой метрики на сфере в координатах широты или долготы, или просто путем просмотра трехмерной геометрии кругов широты и долготы вблизи заданной точки на сфере.

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 должно работать точно так же.

Вы выбираете, какую функцию хотите, сферическую или приблизительно евклидову.Последнее, вероятно, быстрее.Вы можете удалить квадратный корень и рассчитать квадрат расстояния, чтобы сделать вещи еще быстрее.

Я написал это на скорую руку, поэтому могут быть некоторые неточности.

...