У меня есть сегмент линии (большая часть круга) на земле. Сегмент линии определяется координатами ее концов. Очевидно, что две точки определяют два отрезка линии, поэтому предположим, что меня интересует более короткая.
Мне дан третий пункт, и я ищу (кратчайшее) расстояние между линией и точкой.
Все координаты заданы по долготе\широте (WGS 84).
Как рассчитать расстояние?
Решение на любом разумном языке программирования.
Здесь мое собственное решение, основанное на идее спросить доктора Мат.. Я был бы рад видеть ваши отзывы.
Отказ от ответственности в первую очередь. Это решение верно для сфер. Земля не является сферой, а система координат (WGS 84) не считает ее сферой. Так что это просто приближение, и я не могу оценить ошибку. Кроме того, на очень малых расстояниях, вероятно, также можно получить хорошее приближение, считая, что все является просто копланарным. Опять же, я не знаю, как "маленькие" расстояния должны быть.
Теперь к делу. Я буду называть концы линий A, B и третьей точки C. В основном алгоритм:
Вычислить T, точку на прямой AB, ближайшую к C, используя следующие 3 векторных произведения:
G = A x B
F = C x G
T = G x F
Нормализовать Т и умножить на радиус земли.
Эти шаги достаточно, если вы ищете расстояние между C и большим кругом, определенным A и B. Если вы меня интересуете расстояние между C и более коротким сегментом линии, вам нужно сделать дополнительный шаг проверки того, что 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};
Ближайший node: 50.17493121381319,19.05846668493702
Но у меня проблема с этими данными:
double [] a = {52.00118, 17.53933};
double [] b = {52.00278, 17.54008};
//point
double [] c = {52.008308, 17.542927};
Ближайшим node является: 52.00834987257176,17.542691313436357, что неверно.
Я думаю, что строка, заданная двумя точками, не является замкнутым сегментом.
Если кому-то это нужно, это loleksy ответ портирован на С#
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
Наконец, я только что увидел, что следующее связанное в разделе "Родственные":
Расстояние от функции отличного круга от точки до линии не работает правильно.