Как рассчитать расстояние от точки до отрезка прямой, на сфере?
У меня есть отрезок линии (часть большого круга) на земле. Отрезок линии определяется координатами ее концов. Очевидно, что две точки определяют два отрезка прямой, поэтому предположим, что меня интересует более короткий отрезок.
Мне дана третья точка, и я ищу (кратчайшее) расстояние между линией и точкой.
Все координаты указаны в долготешироте (WGS 84).
Как рассчитать расстояние?
Решение в любом разумном язык программирования подойдет.
7 ответов:
Вот мое собственное решение, основанное на идее в спросите доктора математики. Я был бы рад увидеть ваши отзывы.
Сначала оговорка. Это решение верно для сфер. Земля не является сферой, и система координат (WGS 84) не предполагает, что это сфера. Так что это всего лишь аппроксимация, и я не могу реально оценить ошибку. Кроме того, для очень малых расстояний, вероятно, также можно получить хорошее приближение, предположив, что все является просто копланаром. И снова я не знаю как. "маленькими" должны быть расстояния.
Теперь перейдем к делу. Я буду называть концы линий A, B и третью точку C. В основном алгоритм заключается в следующем:Этих шагов достаточно, если вы ищете расстояние между C и большим кругом, определенным A и B. Если вас, как и меня, интересует расстояние между C и более коротким отрезком прямой, вам нужно сделать дополнительный шаг, чтобы убедиться, что T действительно находится на этом отрезке. Если это не так, то обязательно ближайшей точкой является один из концов A или B - самый простой способ проверить, какой из них. В общих чертах идея, лежащая в основе трех векторных произведений, заключается в следующем. Первая из них (G) дает нам плоскость большого круга A и B (то есть плоскость,содержащую A, B и начало координат). Вторая (F) дает нам большую окружность, проходящую через C и перпендикулярную G. тогда T-пересечение больших окружностей, определенных F и G, приведенное к правильная длина путем нормализации и умножения на R.
- преобразуем координаты сначала в декартовы координаты (с началом координат в центре Земли) - например, здесь.
Вычислите T, точку на прямой AB, ближайшую к C, используя следующие 3 векторных произведения:
G = A x B
F = C x G
T = G x F
Нормализуем T и умножаем на радиус Земли.
- преобразования t вернуться к долготы\широты.
- вычислите расстояние между T и C - , например, здесь.
Вот некоторый частичный 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); }Нахождение ближайшей точки На отрезке:
Это простой метод проверки, находится ли точка T, которая, как мы знаем, находится на той же большой окружности, что и A и B, на более коротком отрезке этой большой окружности. Однако там есть более эффективные методы для этого: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; }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; }
Попробуйте расстояние от точки до Большого круга, от Спросите доктора математики. Вам все еще нужно преобразовать долготу/широту в сферические координаты и масштаб для радиуса Земли, но это кажется хорошим направлением.
Это полный код для принятого ответа как 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, что это неправильно.
Я думаю, что линия, указанная двумя точками, не является замкнутым отрезком.
Если кому-то это нужно, это ответ лолексы, портированный на 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-плоскости (с помощью функции line_locate_point PostGIS), а затем (2) вычисление расстояния между этой и третьей точками на сфероиде. Я понятия не имею, является ли математически Шаг (1) правильным, но я был бы удивлен.
Http://postgis.refractions.net/pipermail/postgis-users/2009-July/023903.html
Наконец, я только что увидел, что следующее связано под "Related":
Расстояние от точки до Функция линии большого круга не работает правильно.