Если вы измеряете расстояния меньше (возможно) 1 градуса широты / долготы, ищете очень высокую производительность приближение и готовы принять более неточность , чем у Haversine формулу, рассмотрим эти две альтернативы:
(1) "Формула плоской земли с полярными координатами" из Вычислительные расстояния :
a = pi/2 - lat1
b = pi/2 - lat2
c = sqrt( a^2 + b^2 - 2 * a * b * cos(lon2 - lon1) )
d = R * c
(2) Теорема Пифагора с поправкой на широту , как видно в посте Эвана Тодда :
d_ew = (long1 - long0) * cos(average(lat0, lat1))
d_ns = (lat1 - lat0)
d = sqrt(d_ew * d_ew + d_ns * d_ns)
ПРИМЕЧАНИЯ:
По сравнению с постом Эвана я заменил average(lat0, lat1) на lat0 внутри cos( lat0 ).
# 2 неясно, являются ли значения градусами, радианами или километрами; Вам также понадобится код конверсии. Смотрите мой полный код в нижней части этого поста.
# 1 предназначен для эффективной работы даже вблизи полюсов, хотя, если вы измеряете расстояние, конечные точки которого находятся на «противоположных» сторонах полюса (долготы различаются более чем на 90 градусов?), Вместо этого рекомендуется использовать Haversine, даже на небольшие расстояния.
Я не очень тщательно измерил ошибки этих подходов, поэтому вы должны взять репрезентативные баллы для вашего приложения и сравнить результаты с какой-то высококачественной библиотекой, чтобы решить, приемлемы ли эти точности. Для расстояний менее нескольких километров я чувствую, что они находятся в пределах 1% от правильного измерения.
Альтернативный способ получить высокую производительность (если применимо):
Если у вас есть большой набор статических точек, в пределах одного или двух градусов долготы / широты, то вы будете затем рассчитывать расстояния от небольшого числа динамического (перемещение ) точек, рассмотрите возможность преобразования ваших статических точек ОДИН РАЗ в содержащую UTM-зону (или в любую другую локальную декартову систему координат), а затем выполните все свои математические операции в этой декартовой системе координат.
Декартово = плоская земля = теорема Пифагора, поэтому distance = sqrt(dx^2 + dy^2).
Тогда легко получить стоимость точного преобразования немногих подвижных точек в UTM.
CAVEAT для # 1 (полярный): может быть очень неправильным для расстояний менее 0,1 (?) Метра. Даже при использовании математики с двойной точностью следующие координаты, истинное расстояние которых составляет около 0,005 метра, были заданы как «ноль» моей реализацией алгоритма Полярный:
Входы:
lon1Xdeg 16.6564465477996 double
lat1Ydeg 57.7760262271983 double
lon2Xdeg 16.6564466358281 double
lat2Ydeg 57.776026248554 double
Результаты:
Oblate spheroid formula:
0.00575254911118364 double
Haversine:
0.00573422966122257 double
Polar:
0
это было связано с двумя факторами u и v, которые точно взаимно отменяли друг друга:
u 0.632619944868587 double
v -0.632619944868587 double
В другом случае он дал расстояние 0.067129 m, когда сжатый сфероидный ответ был 0.002887 m. Проблема заключалась в том, что cos(lon2 - lon1) было слишком близко к 1, поэтому функция cos вернула точно 1.
Помимо измерения расстояний до метра, максимальные ошибки (по сравнению с сжатой формулой сфероида) я нашел для ограниченных данных о малых расстояниях, которые я до сих пор вводил:
maxHaversineErrorRatio 0.00350976281908381 double
maxPolarErrorRatio 0.0510789996931342 double
где «1» будет означать 100% ошибку в ответе; например когда он возвращал «0», это была ошибка «1» (исключено сверху «maxPolar»). Таким образом, «0,01» будет ошибкой «1 части на 100» или 1%.
Сравнение полярной ошибки с ошибкой Хаверсайна на расстояниях менее 2000 метров , чтобы увидеть, насколько хуже эта более простая формула. Пока что худшее, что я видел, это 51 часть на 1000 для Polar против 4 частей на 1000 для Haversine. Приблизительно на 58 градусах широты.
Теперь реализовано "пифагорейское с настройкой широты".
Он НАМНОГО более последовательный, чем полярный для расстояний <2000 м. <br>
Я изначально думал, что полярные проблемы были только тогда, когда <1 м, <br>
но результат, показанный непосредственно ниже, весьма тревожит.
Когда расстояния приближаются к нулю, пифагорейская / широта приближается к haversine.
Например это измерение ~ 217 метров:
lon1Xdeg 16.6531667510102 double
lat1Ydeg 57.7751705615804 double
lon2Xdeg 16.6564468739869 double
lat2Ydeg 57.7760263007586 double
oblate 217.201200413731
haversine 216.518428601051
polar 226.128616011973
pythag-cos 216.518428631907
havErrRatio 0.00314349925958048
polErrRatio 0.041102054598393
pycErrRatio 0.00314349911751603
Polar имеет гораздо худшую ошибку с этими входами; либо в моем коде есть какая-то ошибка, либо в функции Cos, на которой я работаю, или я должен порекомендовать не использовать Polar, хотя большинство измерений Polar были намного ближе, чем эта.
OTOH, пифагорейский, даже с регулировкой * cos(latitude), имеет ошибку, которая увеличивается быстрее, чем расстояние (отношение max_error / distance увеличивается для больших расстояний), поэтому вам нужно тщательно учитывать максимальное расстояние, которое вы будете измерять, и допустимая ошибка. Кроме того, не рекомендуется сравнивать два почти равных расстояния, используя пифагорейское, чтобы решить, какое из них короче, так как ошибка различна в разных НАПРАВЛЕНИЯХ (доказательства не показаны).
Измерения в наихудшем случае, errorRatio = Abs(error) / distance (Швеция; до 2000 м):
t_maxHaversineErrorRatio 0.00351012021578681 double
t_maxPolarErrorRatio 66.0825360597085 double
t_maxPythagoreanErrorRatio 0.00350976281416454 double
Как упоминалось ранее, экстремальные полярные ошибки относятся к расстояниям менее метра, где он может сообщать о нуле вместо 6 см или сообщать о расстоянии более 0,5 м на расстоянии 1 см (отсюда наихудший случай "66 x", показанный в t_maxPolarErrorRatio), но есть и плохие результаты на больших расстояниях. [Необходимо еще раз проверить с помощью функции косинуса, которая, как известно, очень точна.]
Измерения, сделанные в коде C # в Xamarin. Android работает на Moto E4.
C # код:
// x=longitude, y= latitude. oblate spheroid formula. TODO: From where?
public static double calculateDistanceDD_AED( double lon1Xdeg, double lat1Ydeg, double lon2Xdeg, double lat2Ydeg )
{
double c_dblEarthRadius = 6378.135; // km
double c_dblFlattening = 1.0 / 298.257223563; // WGS84 inverse
// flattening
// Q: Why "-" for longitudes??
double p1x = -degreesToRadians( lon1Xdeg );
double p1y = degreesToRadians( lat1Ydeg );
double p2x = -degreesToRadians( lon2Xdeg );
double p2y = degreesToRadians( lat2Ydeg );
double F = (p1y + p2y) / 2;
double G = (p1y - p2y) / 2;
double L = (p1x - p2x) / 2;
double sing = Math.Sin( G );
double cosl = Math.Cos( L );
double cosf = Math.Cos( F );
double sinl = Math.Sin( L );
double sinf = Math.Sin( F );
double cosg = Math.Cos( G );
double S = sing * sing * cosl * cosl + cosf * cosf * sinl * sinl;
double C = cosg * cosg * cosl * cosl + sinf * sinf * sinl * sinl;
double W = Math.Atan2( Math.Sqrt( S ), Math.Sqrt( C ) );
if (W == 0.0)
return 0.0;
double R = Math.Sqrt( (S * C) ) / W;
double H1 = (3 * R - 1.0) / (2.0 * C);
double H2 = (3 * R + 1.0) / (2.0 * S);
double D = 2 * W * c_dblEarthRadius;
// Apply flattening factor
D = D * (1.0 + c_dblFlattening * H1 * sinf * sinf * cosg * cosg - c_dblFlattening * H2 * cosf * cosf * sing * sing);
// Transform to meters
D = D * 1000.0;
// tmstest
if (true)
{
// Compare Haversine.
double haversine = HaversineApproxDistanceGeo( lon1Xdeg, lat1Ydeg, lon2Xdeg, lat2Ydeg );
double error = haversine - D;
double absError = Math.Abs( error );
double errorRatio = absError / D;
if (errorRatio > t_maxHaversineErrorRatio)
{
if (errorRatio > t_maxHaversineErrorRatio * 1.1)
Helper.test();
t_maxHaversineErrorRatio = errorRatio;
}
// Compare Polar Coordinate Flat Earth.
double polarDistanceGeo = ApproxDistanceGeo_Polar( lon1Xdeg, lat1Ydeg, lon2Xdeg, lat2Ydeg, D );
double error2 = polarDistanceGeo - D;
double absError2 = Math.Abs( error2 );
double errorRatio2 = absError2 / D;
if (errorRatio2 > t_maxPolarErrorRatio)
{
if (polarDistanceGeo > 0)
{
if (errorRatio2 > t_maxPolarErrorRatio * 1.1)
Helper.test();
t_maxPolarErrorRatio = errorRatio2;
}
else
Helper.dubious();
}
// Compare Pythagorean Theorem with Latitude Adjustment.
double pythagoreanDistanceGeo = ApproxDistanceGeo_PythagoreanCosLatitude( lon1Xdeg, lat1Ydeg, lon2Xdeg, lat2Ydeg, D );
double error3 = pythagoreanDistanceGeo - D;
double absError3 = Math.Abs( error3 );
double errorRatio3 = absError3 / D;
if (errorRatio3 > t_maxPythagoreanErrorRatio)
{
if (D < 2000)
{
if (errorRatio3 > t_maxPythagoreanErrorRatio * 1.05)
Helper.test();
t_maxPythagoreanErrorRatio = errorRatio3;
}
}
}
return D;
}
// As a fraction of the distance.
private static double t_maxHaversineErrorRatio, t_maxPolarErrorRatio, t_maxPythagoreanErrorRatio;
// Average of equatorial and polar radii (meters).
public const double EarthAvgRadius = 6371000;
public const double EarthAvgCircumference = EarthAvgRadius * 2 * PI;
// CAUTION: This is an average of great circles; won't be the actual distance of any longitude or latitude degree.
public const double EarthAvgMeterPerGreatCircleDegree = EarthAvgCircumference / 360;
// Haversine formula (assumes Earth is sphere).
// "deg" = degrees.
// Perhaps based on Haversine Formula in https://cs.nyu.edu/visual/home/proj/tiger/gisfaq.html
public static double HaversineApproxDistanceGeo(double lon1Xdeg, double lat1Ydeg, double lon2Xdeg, double lat2Ydeg)
{
double lon1 = degreesToRadians( lon1Xdeg );
double lat1 = degreesToRadians( lat1Ydeg );
double lon2 = degreesToRadians( lon2Xdeg );
double lat2 = degreesToRadians( lat2Ydeg );
double dlon = lon2 - lon1;
double dlat = lat2 - lat1;
double sinDLat2 = Sin( dlat / 2 );
double sinDLon2 = Sin( dlon / 2 );
double a = sinDLat2 * sinDLat2 + Cos( lat1 ) * Cos( lat2 ) * sinDLon2 * sinDLon2;
double c = 2 * Atan2( Sqrt( a ), Sqrt( 1 - a ) );
double d = EarthAvgRadius * c;
return d;
}
// From https://stackoverflow.com/a/19772119/199364
// Based on Polar Coordinate Flat Earth in https://cs.nyu.edu/visual/home/proj/tiger/gisfaq.html
public static double ApproxDistanceGeo_Polar( double lon1deg, double lat1deg, double lon2deg, double lat2deg, double D = 0 )
{
double approxUnitDistSq = ApproxUnitDistSq_Polar(lon1deg, lat1deg, lon2deg, lat2deg, D);
double c = Sqrt( approxUnitDistSq );
return EarthAvgRadius * c;
}
// Might be useful to avoid taking Sqrt, when comparing to some threshold.
// Threshold would have to be adjusted to match: Power(threshold / EarthAvgRadius, 2)
private static double ApproxUnitDistSq_Polar(double lon1deg, double lat1deg, double lon2deg, double lat2deg, double D = 0 )
{
const double HalfPi = PI / 2; //1.5707963267949;
double lon1 = degreesToRadians(lon1deg);
double lat1 = degreesToRadians(lat1deg);
double lon2 = degreesToRadians(lon2deg);
double lat2 = degreesToRadians(lat2deg);
double a = HalfPi - lat1;
double b = HalfPi - lat2;
double u = a * a + b * b;
double dlon21 = lon2 - lon1;
double cosDeltaLon = Cos( dlon21 );
double v = -2 * a * b * cosDeltaLon;
// TBD: Is "Abs" necessary? That is, is "u + v" ever negative?
// (I think not; "v" looks like a secondary term. Though might be round-off issue near zero when a~=b.)
double approxUnitDistSq = Abs(u + v);
//if (approxUnitDistSq.nearlyEquals(0, 1E-16))
// Helper.dubious();
//else if (D > 0)
//{
// double dba = b - a;
// double unitD = D / EarthAvgRadius;
// double unitDSq = unitD * unitD;
// if (approxUnitDistSq > 2 * unitDSq)
// Helper.dubious();
// else if (approxUnitDistSq * 2 < unitDSq)
// Helper.dubious();
//}
return approxUnitDistSq;
}
// Pythagorean Theorem with Latitude Adjustment - from Ewan Todd - https://stackoverflow.com/a/1664836/199364
// Refined by ToolmakerSteve - https://stackoverflow.com/a/53468745/199364
public static double ApproxDistanceGeo_PythagoreanCosLatitude( double lon1deg, double lat1deg, double lon2deg, double lat2deg, double D = 0 )
{
double approxDegreesSq = ApproxDegreesSq_PythagoreanCosLatitude( lon1deg, lat1deg, lon2deg, lat2deg );
// approximate degrees on the great circle between the points.
double d_degrees = Sqrt( approxDegreesSq );
return d_degrees * EarthAvgMeterPerGreatCircleDegree;
}
public static double ApproxDegreesSq_PythagoreanCosLatitude( double lon1deg, double lat1deg, double lon2deg, double lat2deg )
{
double avgLatDeg = average( lat1deg , lat2deg );
double avgLat = degreesToRadians( avgLatDeg );
double d_ew = (lon2deg - lon1deg) * Cos( avgLat );
double d_ns = (lat2deg - lat1deg);
double approxDegreesSq = d_ew * d_ew + d_ns * d_ns;
return approxDegreesSq;
}