Как рассчитать расстояние от точки до сегмента линии, на сфере?
https://stackoverflow.com/questions/1299567
-
18-09-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianВопрос
У меня есть линейный сегмент (отличная часть круга) на Земле. Линейный сегмент определяется координатами его концов. Очевидно, что две точки определяют два сегмента линии, поэтому предположим, что я заинтересован в более короткой.
Мне дали третью точку, и я ищу (самое короткое) расстояние между линией и точкой.
Все координаты приведены в долготе Широта (WGS 84).
Как рассчитать расстояние?
Решение на любом разумном языке программирования подойдет.
Решение
Вот мое собственное решение, основанное на идее в Спросите доктора Мат. Анкет Я был бы рад увидеть ваши отзывы.
Отказ от ответственности первым. Это решение верно для сфер. Земля не сфера, и система координат (WGS 84) не предполагает, что это сфера. Так что это просто приближение, и я не могу по -настоящему оценить ошибку. Кроме того, на очень маленьких расстояниях, вероятно, также можно получить хорошее приближение, предполагая, что все является просто копланаром. Опять же, я не знаю, какими «маленькими» должны быть расстояния.
Теперь к бизнесу. Я назову концы строк A, B и третьей точкой C. В основном алгоритм должен:
- Сначала преобразовать координаты в картезианские координаты (с началом происхождения в центре Земли) - Например.
Рассчитайте T, точка на линии AB, которая близка к C, используя следующие 3 векторные продукты:
G = a x b
F = C X G
T = g x f
Нормализуйте t и умножьте на радиус Земли.
- Преобразовать T обратно в долготу широта.
- Рассчитайте расстояние между T и C - Например.
Этих шагов достаточно, если вы ищете расстояние между C и великим кругом, определяемым A и B. Если, как и я T действительно в этом сегменте. Если это не так, то обязательно ближайшая точка - один из целей A или B - самый простой способ - проверить, какой из них.
В общих чертах, идея трех векторных продуктов является следующей. Первый (G) дает нам плоскость великого круга A и B (так плоскость, содержащая A, B и начало). Второй (f) дает нам великий круг, который проходит через C и перпендикулярно G. Тогда t - это пересечение великих кругов, определяемых F и G, приведенным в правильную длину путем нормализации и умножения на R.
Вот какой -то частичный код Java для этого.
Найти ближайшую точку на великом круге. Входы и вывод-массивы длины-2. Промежуточные массивы имеют длину 3.
double[] nearestPointGreatCircle(double[] a, double[] b, double c[])
{
double[] a_ = toCartsian(a);
double[] b_ = toCartsian(b);
double[] c_ = toCartsian(c);
double[] G = vectorProduct(a_, b_);
double[] F = vectorProduct(c_, G);
double[] t = vectorProduct(G, F);
normalize(t);
multiplyByScalar(t, R_EARTH);
return fromCartsian(t);
}
Поиск ближайшей точки на сегменте:
double[] nearestPointSegment (double[] a, double[] b, double[] c)
{
double[] t= nearestPointGreatCircle(a,b,c);
if (onSegment(a,b,t))
return t;
return (distance(a,c) < distance(b,c)) ? a : c;
}
Это простой метод тестирования, если точка T, которая, как мы знаем, находится на том же великом круге, что и A и B, находится на более коротком сегменте этого великого круга. Однако есть более эффективные методы для этого:
boolean onSegment (double[] a, double[] b, double[] t)
{
// should be return distance(a,t)+distance(b,t)==distance(a,b),
// but due to rounding errors, we use:
return Math.abs(distance(a,b)-distance(a,t)-distance(b,t)) < PRECISION;
}
Другие советы
Пытаться Расстояние от точки до большого круга, от Ask Dr. Math. Вам все еще нужно преобразовать долготу/широту в сферические координаты и масштабировать радиус Земли, но это кажется хорошим направлением.
Это полный код для принятого ответа как IdeOne Fiddle (найденный здесь):
import java.util.*;
import java.lang.*;
import java.io.*;
/* Name of the class has to be "Main" only if the class is public. */
class Ideone
{
private static final double _eQuatorialEarthRadius = 6378.1370D;
private static final double _d2r = (Math.PI / 180D);
private static double PRECISION = 0.1;
// Haversine Algorithm
// source: http://stackoverflow.com/questions/365826/calculate-distance-between-2-gps-coordinates
private static double HaversineInM(double lat1, double long1, double lat2, double long2) {
return (1000D * HaversineInKM(lat1, long1, lat2, long2));
}
private static double HaversineInKM(double lat1, double long1, double lat2, double long2) {
double dlong = (long2 - long1) * _d2r;
double dlat = (lat2 - lat1) * _d2r;
double a = Math.pow(Math.sin(dlat / 2D), 2D) + Math.cos(lat1 * _d2r) * Math.cos(lat2 * _d2r)
* Math.pow(Math.sin(dlong / 2D), 2D);
double c = 2D * Math.atan2(Math.sqrt(a), Math.sqrt(1D - a));
double d = _eQuatorialEarthRadius * c;
return d;
}
// Distance between a point and a line
public static void pointLineDistanceTest() {
//line
//double [] a = {50.174315,19.054743};
//double [] b = {50.176019,19.065042};
double [] a = {52.00118, 17.53933};
double [] b = {52.00278, 17.54008};
//point
//double [] c = {50.184373,19.054657};
double [] c = {52.008308, 17.542927};
double[] nearestNode = nearestPointGreatCircle(a, b, c);
System.out.println("nearest node: " + Double.toString(nearestNode[0]) + "," + Double.toString(nearestNode[1]));
double result = HaversineInM(c[0], c[1], nearestNode[0], nearestNode[1]);
System.out.println("result: " + Double.toString(result));
}
// source: http://stackoverflow.com/questions/1299567/how-to-calculate-distance-from-a-point-to-a-line-segment-on-a-sphere
private static double[] nearestPointGreatCircle(double[] a, double[] b, double c[])
{
double[] a_ = toCartsian(a);
double[] b_ = toCartsian(b);
double[] c_ = toCartsian(c);
double[] G = vectorProduct(a_, b_);
double[] F = vectorProduct(c_, G);
double[] t = vectorProduct(G, F);
return fromCartsian(multiplyByScalar(normalize(t), _eQuatorialEarthRadius));
}
@SuppressWarnings("unused")
private static double[] nearestPointSegment (double[] a, double[] b, double[] c)
{
double[] t= nearestPointGreatCircle(a,b,c);
if (onSegment(a,b,t))
return t;
return (HaversineInKM(a[0], a[1], c[0], c[1]) < HaversineInKM(b[0], b[1], c[0], c[1])) ? a : b;
}
private static boolean onSegment (double[] a, double[] b, double[] t)
{
// should be return distance(a,t)+distance(b,t)==distance(a,b),
// but due to rounding errors, we use:
return Math.abs(HaversineInKM(a[0], a[1], b[0], b[1])-HaversineInKM(a[0], a[1], t[0], t[1])-HaversineInKM(b[0], b[1], t[0], t[1])) < PRECISION;
}
// source: http://stackoverflow.com/questions/1185408/converting-from-longitude-latitude-to-cartesian-coordinates
private static double[] toCartsian(double[] coord) {
double[] result = new double[3];
result[0] = _eQuatorialEarthRadius * Math.cos(Math.toRadians(coord[0])) * Math.cos(Math.toRadians(coord[1]));
result[1] = _eQuatorialEarthRadius * Math.cos(Math.toRadians(coord[0])) * Math.sin(Math.toRadians(coord[1]));
result[2] = _eQuatorialEarthRadius * Math.sin(Math.toRadians(coord[0]));
return result;
}
private static double[] fromCartsian(double[] coord){
double[] result = new double[2];
result[0] = Math.toDegrees(Math.asin(coord[2] / _eQuatorialEarthRadius));
result[1] = Math.toDegrees(Math.atan2(coord[1], coord[0]));
return result;
}
// Basic functions
private static double[] vectorProduct (double[] a, double[] b){
double[] result = new double[3];
result[0] = a[1] * b[2] - a[2] * b[1];
result[1] = a[2] * b[0] - a[0] * b[2];
result[2] = a[0] * b[1] - a[1] * b[0];
return result;
}
private static double[] normalize(double[] t) {
double length = Math.sqrt((t[0] * t[0]) + (t[1] * t[1]) + (t[2] * t[2]));
double[] result = new double[3];
result[0] = t[0]/length;
result[1] = t[1]/length;
result[2] = t[2]/length;
return result;
}
private static double[] multiplyByScalar(double[] normalize, double k) {
double[] result = new double[3];
result[0] = normalize[0]*k;
result[1] = normalize[1]*k;
result[2] = normalize[2]*k;
return result;
}
public static void main(String []args){
System.out.println("Hello World");
Ideone.pointLineDistanceTest();
}
}
Он отлично работает для комментариев:
//line
double [] a = {50.174315,19.054743};
double [] b = {50.176019,19.065042};
//point
double [] c = {50.184373,19.054657};
Ближайший узел: 50.17493121381319,19.05846668493702
Но у меня есть проблемы с этими данными:
double [] a = {52.00118, 17.53933};
double [] b = {52.00278, 17.54008};
//point
double [] c = {52.008308, 17.542927};
Ближайший узел: 52.00834987257176,17.542691313436357, что неправильно.
Я думаю, что линия, указанная двумя точками, не является закрытым сегментом.
Если это кому -то нужно, это Loleksy отвечает, портированный на C#
private static double _eQuatorialEarthRadius = 6378.1370D;
private static double _d2r = (Math.PI / 180D);
private static double PRECISION = 0.1;
// Haversine Algorithm
// source: http://stackoverflow.com/questions/365826/calculate-distance-between-2-gps-coordinates
private static double HaversineInM(double lat1, double long1, double lat2, double long2) {
return (1000D * HaversineInKM(lat1, long1, lat2, long2));
}
private static double HaversineInKM(double lat1, double long1, double lat2, double long2) {
double dlong = (long2 - long1) * _d2r;
double dlat = (lat2 - lat1) * _d2r;
double a = Math.Pow(Math.Sin(dlat / 2D), 2D) + Math.Cos(lat1 * _d2r) * Math.Cos(lat2 * _d2r)
* Math.Pow(Math.Sin(dlong / 2D), 2D);
double c = 2D * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1D - a));
double d = _eQuatorialEarthRadius * c;
return d;
}
// Distance between a point and a line
static double pointLineDistanceGEO(double[] a, double[] b, double[] c)
{
double[] nearestNode = nearestPointGreatCircle(a, b, c);
double result = HaversineInKM(c[0], c[1], nearestNode[0], nearestNode[1]);
return result;
}
// source: http://stackoverflow.com/questions/1299567/how-to-calculate-distance-from-a-point-to-a-line-segment-on-a-sphere
private static double[] nearestPointGreatCircle(double[] a, double[] b, double [] c)
{
double[] a_ = toCartsian(a);
double[] b_ = toCartsian(b);
double[] c_ = toCartsian(c);
double[] G = vectorProduct(a_, b_);
double[] F = vectorProduct(c_, G);
double[] t = vectorProduct(G, F);
return fromCartsian(multiplyByScalar(normalize(t), _eQuatorialEarthRadius));
}
private static double[] nearestPointSegment (double[] a, double[] b, double[] c)
{
double[] t= nearestPointGreatCircle(a,b,c);
if (onSegment(a,b,t))
return t;
return (HaversineInKM(a[0], a[1], c[0], c[1]) < HaversineInKM(b[0], b[1], c[0], c[1])) ? a : b;
}
private static bool onSegment (double[] a, double[] b, double[] t)
{
// should be return distance(a,t)+distance(b,t)==distance(a,b),
// but due to rounding errors, we use:
return Math.Abs(HaversineInKM(a[0], a[1], b[0], b[1])-HaversineInKM(a[0], a[1], t[0], t[1])-HaversineInKM(b[0], b[1], t[0], t[1])) < PRECISION;
}
// source: http://stackoverflow.com/questions/1185408/converting-from-longitude-latitude-to-cartesian-coordinates
private static double[] toCartsian(double[] coord) {
double[] result = new double[3];
result[0] = _eQuatorialEarthRadius * Math.Cos(deg2rad(coord[0])) * Math.Cos(deg2rad(coord[1]));
result[1] = _eQuatorialEarthRadius * Math.Cos(deg2rad(coord[0])) * Math.Sin(deg2rad(coord[1]));
result[2] = _eQuatorialEarthRadius * Math.Sin(deg2rad(coord[0]));
return result;
}
private static double[] fromCartsian(double[] coord){
double[] result = new double[2];
result[0] = rad2deg(Math.Asin(coord[2] / _eQuatorialEarthRadius));
result[1] = rad2deg(Math.Atan2(coord[1], coord[0]));
return result;
}
// Basic functions
private static double[] vectorProduct (double[] a, double[] b){
double[] result = new double[3];
result[0] = a[1] * b[2] - a[2] * b[1];
result[1] = a[2] * b[0] - a[0] * b[2];
result[2] = a[0] * b[1] - a[1] * b[0];
return result;
}
private static double[] normalize(double[] t) {
double length = Math.Sqrt((t[0] * t[0]) + (t[1] * t[1]) + (t[2] * t[2]));
double[] result = new double[3];
result[0] = t[0]/length;
result[1] = t[1]/length;
result[2] = t[2]/length;
return result;
}
private static double[] multiplyByScalar(double[] normalize, double k) {
double[] result = new double[3];
result[0] = normalize[0]*k;
result[1] = normalize[1]*k;
result[2] = normalize[2]*k;
return result;
}
Для расстояния до нескольких тысяч метров я бы упростил проблему с сферы до самолета. Затем проблема довольно просто, так как можно использовать простой расчет треугольника:
У нас есть точки A и B, и мы ищем расстояние x до линии AB. Затем:
Location a;
Location b;
Location x;
double ax = a.distanceTo(x);
double alfa = (Math.abs(a.bearingTo(b) - a.bearingTo(x))) / 180
* Math.PI;
double distance = Math.sin(alfa) * ax;
Самое короткое расстояние между двумя точками на сфере - это меньшая сторона большого круга, проходящего через две точки. Я уверен, что вы уже знаете это. Здесь есть аналогичный вопрос http://www.physicsforums.com/archive/index.php/t-178252.html Это может помочь вам математически моделировать его.
Я не уверен, насколько вероятно, что вы получите кодированный пример этого, если честно.
Я в основном ищу то же самое прямо сейчас, за исключением того, что я строго говорю, не заботится о том, чтобы иметь сегмент большого круга, а просто хочу расстояния до любой точки на полном круге.
Две ссылки, которые я сейчас расследую:
Эта страница Упоминает «перекрестное расстояние», которое в основном кажется тем, что вы ищете.
Кроме того, в следующем потоке в списке рассылки PostGIS попытка, похоже, (1) определит ближайшую точку на великом круге с той же формулой, используемой для линейного дистанции на 2D-плоскости (с PostGIS 'LINE_LOCATE_POINT), а затем, а затем, а затем, а затем, а затем, а затем, а затем, а затем и затем (2) Расчет расстояния между этим и третьей точкой на сфероиде. Я понятия не имею, является ли математически шаг (1) правильным, но я был бы удивлен.
http://postgis.refractions.net/pipermail/postgis-users/2009-july/023903.html
Наконец, я только что увидел, что следующее связано с «Связанным»:
Расстояние от точки до линии отличная функция круга не работает правильно.
