From e73209676d71b04e036cf15e2c68f23d9f3e6b40 Mon Sep 17 00:00:00 2001 From: Patrick Brosi Date: Fri, 21 Aug 2026 15:37:30 +0200 Subject: [PATCH 1/4] use the mean equatorial rad for the haversine, also add `vincenty()` (exact iterative method) and `andoyerLambert()` (refinement of haversine to correct for sphere approximation) and add `adaptiveMeterDist()` methods which try to select haversine, andoyerLambert and then vincenty (ordered here in computation cost) based on a given tolerance factor in meters (0.5 per default) AND ALSO use haversine for points which are nearly antipodal (all the spheriod methods fail catatrophically then) --- geo/Geo.h | 53 ++++++++++++++ geo/Geo.tpp | 203 +++++++++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 253 insertions(+), 3 deletions(-) diff --git a/geo/Geo.h b/geo/Geo.h index 9b5285c..51b89e9 100644 --- a/geo/Geo.h +++ b/geo/Geo.h @@ -11,8 +11,11 @@ #include #include +#include #include +#include #include +#include #include "util/Misc.h" #include "util/String.h" @@ -142,9 +145,20 @@ const static double AVERAGING_STEP = 20; const static double M_PER_DEG = 111319.4; +// earth radius at the equator const static double EQUATORIAL_RAD = 6378137.0; + +// mean earth radius, used for haversine +const static double MEAN_EARTH_RAD = 6371008.8; + const static double MIN_METERS_PER_LAT_RAD = 6356752.314245; +// WGS84 flattening f = 1 - b / a, used by the ellipsoidal distance measures +// (andoyerLambert, vincenty) +const static double FLATTENING = 1.0 / 298.257223563; +// max error of haversine against the true distance +const static double HAVERSINE_MAX_ERR = 0.00561437; + enum WKTType : uint8_t { NONE = 0, POINT = 1, @@ -1385,6 +1399,45 @@ double haversineWebMerc(T x1, T y1, T x2, T y2); template double haversineWebMerc(const Point& a, const Point& b); +template +double andoyerLambert(T lat1, T lon1, T lat2, T lon2); + +template +double andoyerLambert(const Point& a, const Point& b); + +template +double andoyerLambertWebMerc(T x1, T y1, T x2, T y2); + +template +double andoyerLambertWebMerc(const Point& a, const Point& b); + +template +double vincenty(T lat1, T lon1, T lat2, T lon2); + +template +double vincenty(const Point& a, const Point& b); + +template +double vincentyWebMerc(T x1, T y1, T x2, T y2); + +template +double vincentyWebMerc(const Point& a, const Point& b); + +// tol is the allows abs error in meters, 0.5 per default +template +double adaptiveMeterDist(T lat1, T lon1, T lat2, T lon2, double tol = 0.5); + +template +double adaptiveMeterDist(const Point& a, const Point& b, + double tol = 0.5); + +template +double adaptiveMeterDistWebMerc(T x1, T y1, T x2, T y2, double tol = 0.5); + +template +double adaptiveMeterDistWebMerc(const Point& a, const Point& b, + double tol = 0.5); + template Line densify(const Line& l, double d); diff --git a/geo/Geo.tpp b/geo/Geo.tpp index 1fff158..f796d14 100644 --- a/geo/Geo.tpp +++ b/geo/Geo.tpp @@ -6016,7 +6016,7 @@ double haversine(T lat1, T lon1, T lat2, T lon2) { const double sDLon = sin(dLon / 2); const double a = (sDLat * sDLat) + (sDLon * sDLon) * cos(lat1) * cos(lat2); - return EQUATORIAL_RAD * 2.0 * atan2(sqrt(a), sqrt(1.0 - a)); + return MEAN_EARTH_RAD * 2.0 * atan2(sqrt(a), sqrt(1.0 - a)); } // _____________________________________________________________________________ @@ -6042,7 +6042,7 @@ double haversineWebMerc(T x1, T y1, T x2, T y2) { const double cLat2 = 2.0 * t2 / q2; const double a = (sDLat * sDLat) + (sDLon * sDLon) * cLat1 * cLat2; - return EQUATORIAL_RAD * 2.0 * atan2(sqrt(a), sqrt(1.0 - a)); + return MEAN_EARTH_RAD * 2.0 * atan2(sqrt(a), sqrt(1.0 - a)); } // _____________________________________________________________________________ @@ -6051,6 +6051,203 @@ double haversineWebMerc(const Point& a, const Point& b) { return haversineWebMerc(a.getX(), a.getY(), b.getX(), b.getY()); } +// _____________________________________________________________________________ +template +double andoyerLambert(T lat1, T lon1, T lat2, T lon2) { + // see https://en.wikipedia.org/wiki/Geographical_distance#Andoyer-Lambert_formula_for_long_lines + double f1 = 1.0 - FLATTENING; + + double b1 = atan(f1 * tan(lat1 * RAD)); + double b2 = atan(f1 * tan(lat2 * RAD)); + + double cB1 = cos(b1); + double cB2 = cos(b2); + + double sDLat = sin((b2 - b1) / 2.0); + double sDLon = sin((lon2 - lon1) * RAD / 2.0); + + double a = (sDLat * sDLat) + (sDLon * sDLon) * cB1 * cB2; + double sig = 2.0 * atan2(sqrt(a), sqrt(1.0 - a)); + + if (sig == 0.0) return 0.0; + + double p = (b1 + b2) / 2.0; + double q = (b2 - b1) / 2.0; + + double sP = sin(p), cP = cos(p); + double sQ = sin(q), cQ = cos(q); + + double sSig = sin(sig); + double sSigH = sin(sig / 2.0), cSigH = cos(sig / 2.0); + + double x = (sig - sSig) * (sP * sP) * (cQ * cQ) / (cSigH * cSigH); + double y = (sig + sSig) * (cP * cP) * (sQ * sQ) / (sSigH * sSigH); + + return EQUATORIAL_RAD * (sig - (FLATTENING / 2.0) * (x + y)); +} + +// _____________________________________________________________________________ +template +double andoyerLambert(const Point& a, const Point& b) { + return andoyerLambert(a.getY(), a.getX(), b.getY(), b.getX()); +} + +// _____________________________________________________________________________ +template +double andoyerLambertWebMerc(T x1, T y1, T x2, T y2) { + const auto a = webMercToLatLng(x1, y1); + const auto b = webMercToLatLng(x2, y2); + return andoyerLambert(a.getY(), a.getX(), b.getY(), b.getX()); +} + +// _____________________________________________________________________________ +template +double andoyerLambertWebMerc(const Point& a, const Point& b) { + return andoyerLambertWebMerc(a.getX(), a.getY(), b.getX(), b.getY()); +} + +// _____________________________________________________________________________ +template +double vincenty(T lat1, T lon1, T lat2, T lon2) { + // see https://en.wikipedia.org/wiki/Vincenty's_formulae + size_t MAX_ITERS = 200; + double f1 = 1.0 - FLATTENING; + double b = EQUATORIAL_RAD * f1; + + double l = (lon2 - lon1) * RAD; + + double u1 = atan(f1 * tan(lat1 * RAD)); + double u2 = atan(f1 * tan(lat2 * RAD)); + + double sU1 = sin(u1), cU1 = cos(u1); + double sU2 = sin(u2), cU2 = cos(u2); + + double lambda = l; + double sSig = 0.0, cSig = 0.0, sig = 0.0, c2Alpha = 0.0, c2SigM = 0.0; + + bool converged = false; + + for (int i = 0; i < MAX_ITERS && !converged; ++i) { + double sLam = sin(lambda), cLam = cos(lambda); + + double t1 = cU2 * sLam; + double t2 = cU1 * sU2 - sU1 * cU2 * cLam; + + sSig = sqrt(t1 * t1 + t2 * t2); + + if (sSig == 0.0) return 0.0; + + cSig = sU1 * sU2 + cU1 * cU2 * cLam; + sig = atan2(sSig, cSig); + + double sAlpha = cU1 * cU2 * sLam / sSig; + c2Alpha = 1.0 - sAlpha * sAlpha; + + c2SigM = c2Alpha == 0.0 ? 0.0 : cSig - 2.0 * sU1 * sU2 / c2Alpha; + + double c = (FLATTENING / 16.0) * c2Alpha * + (4.0 + FLATTENING * (4.0 - 3.0 * c2Alpha)); + + double prev = lambda; + lambda = l + (1.0 - c) * FLATTENING * sAlpha * + (sig + + c * sSig * + (c2SigM + c * cSig * (-1.0 + 2.0 * c2SigM * c2SigM))); + + converged = fabs(lambda - prev) < 1e-12; + } + + if (!converged) return std::numeric_limits::quiet_NaN(); + + double uSq = c2Alpha * (EQUATORIAL_RAD * EQUATORIAL_RAD - b * b) / (b * b); + + double aa = 1.0 + (uSq / 16384.0) * + (4096.0 + uSq * (-768.0 + uSq * (320.0 - 175.0 * uSq))); + double bb = + (uSq / 1024.0) * (256.0 + uSq * (-128.0 + uSq * (74.0 - 47.0 * uSq))); + + double dSig = + bb * sSig * + (c2SigM + (bb / 4.0) * (cSig * (-1.0 + 2.0 * c2SigM * c2SigM) - + (bb / 6.0) * c2SigM * (-3.0 + 4.0 * sSig * sSig) * + (-3.0 + 4.0 * c2SigM * c2SigM))); + + return b * aa * (sig - dSig); +} + +// _____________________________________________________________________________ +template +double vincenty(const Point& a, const Point& b) { + return vincenty(a.getY(), a.getX(), b.getY(), b.getX()); +} + +// _____________________________________________________________________________ +template +double vincentyWebMerc(T x1, T y1, T x2, T y2) { + const auto a = webMercToLatLng(x1, y1); + const auto b = webMercToLatLng(x2, y2); + return vincenty(a.getY(), a.getX(), b.getY(), b.getX()); +} + +// _____________________________________________________________________________ +template +double vincentyWebMerc(const Point& a, const Point& b) { + return vincentyWebMerc(a.getX(), a.getY(), b.getX(), b.getY()); +} + +// _____________________________________________________________________________ +template +double adaptiveMeterDist(T lat1, T lon1, T lat2, T lon2, double tol) { + const double dHaversine = haversine(lat1, lon1, lat2, lon2); + + // the new methods below fail catastrophically if the two points are very + // close to being antipodes (on exactly the opposite side of the earth). The + // haversine distance doesnt fail in this case, so + // as soon as the haversine distance (exact to around 0.56%) is within 150km + // to half the earth cirumvernece, return the haversine + if ((M_PI * MEAN_EARTH_RAD) - dHaversine < 150000.0) return dHaversine; + + // dHaversine may have underestimated, so correct first by MAX_ERROR, then + // check whether the max error at that corrected distance ist still smaller + // than the tolerance if so, return haversine directly for speed + if (tol > HAVERSINE_MAX_ERR * ((HAVERSINE_MAX_ERR + 1.0) * dHaversine)) + return dHaversine; + + // andoyerLabert is correct to about a factor of 0.0000014 for dists within + // 100000 meters, pad this here a bit to be safe and return andoyerLambert + // directly if we are still within tolernace + if (dHaversine < 10000.0 * 1000.0 && tol > 0.000004 * dHaversine) + return andoyerLambert(lat1, lon1, lat2, lon2); + + // only now do iterative vincenty + const double v = vincenty(lat1, lon1, lat2, lon2); + + // if vincenty doesnt converge for some reason, return the haversine as + // fallback + return std::isnan(v) ? dHaversine : v; +} + +// _____________________________________________________________________________ +template +double adaptiveMeterDist(const Point& a, const Point& b, double tol) { + return adaptiveMeterDist(a.getY(), a.getX(), b.getY(), b.getX(), tol); +} + +// _____________________________________________________________________________ +template +double adaptiveMeterDistWebMerc(T x1, T y1, T x2, T y2, double tol) { + const auto a = webMercToLatLng(x1, y1); + const auto b = webMercToLatLng(x2, y2); + return adaptiveMeterDist(a.getY(), a.getX(), b.getY(), b.getX(), tol); +} + +// _____________________________________________________________________________ +template +double adaptiveMeterDistWebMerc(const Point& a, const Point& b, + double tol) { + return adaptiveMeterDistWebMerc(a.getX(), a.getY(), b.getX(), b.getY(), tol); +} + // _____________________________________________________________________________ template Line sparseify(const Line& l, double mind) { @@ -7232,7 +7429,7 @@ double webMercMaxEuclideanDist(const Box& boxA, const Box& boxB, auto scale = getMinMaxLocalScaleFactorsWebMerc(boxA, boxB, maxD); // use meters here directly, we are in web mercator world - return maxD / scale.first; + return (maxD * (1.0 + HAVERSINE_MAX_ERR)) / scale.first; } // _____________________________________________________________________________ From 6be905a8991c504d6994382c21cedb8e023255b0 Mon Sep 17 00:00:00 2001 From: Patrick Brosi Date: Fri, 21 Aug 2026 16:09:43 +0200 Subject: [PATCH 2/4] add tests for adaptiveMeterDistWebMerc --- tests/GeoTestDist.cpp | 123 ++++++++++++++++++++++++++++++++++++------ 1 file changed, 107 insertions(+), 16 deletions(-) diff --git a/tests/GeoTestDist.cpp b/tests/GeoTestDist.cpp index 96ffd07..1d8b4a4 100644 --- a/tests/GeoTestDist.cpp +++ b/tests/GeoTestDist.cpp @@ -378,13 +378,13 @@ static void testDistComplexGeoms(const LargeTestGeoms& g) { approx(6.5434)); TEST(util::geo::dist(g.germanyX, g.spainX), ==, approx(6.5434)); TEST(util::geo::webMercMeterDist(g.germanyMX, g.spainMX), ==, - approx(653276.57366)); + approx(652546.47237)); TEST(util::geo::withinDist(g.germany, g.spain, 10), ==, approx(6.5434)); TEST(util::geo::withinDist(g.germany, g.spain, 6.54341), ==, approx(6.5434)); TEST(util::geo::dist(g.germany, g.spain), ==, approx(6.5434)); TEST(util::geo::webMercMeterDist(g.germanyM, g.spainM), ==, - approx(653276.57366)); + approx(652546.47237)); TEST(util::geo::withinDist(g.germany, g.germany, 10), ==, approx(0)); TEST(util::geo::withinDist(g.germany, g.germany, 0), ==, approx(0)); @@ -431,14 +431,14 @@ static void testDistComplexGeoms(const LargeTestGeoms& g) { approx(7.00409)); TEST(util::geo::dist(g.spainX, g.flixbusX), ==, approx(7.00409)); TEST(util::geo::webMercMeterDist(g.spainMX, g.flixbusMX), ==, - approx(703461.25144)); + approx(702675.06380)); TEST(util::geo::withinDist(g.spain, g.flixbus, 10), ==, approx(7.00409)); TEST(util::geo::withinDist(g.spain, g.flixbus, 7.004091), ==, approx(7.00409)); TEST(util::geo::dist(g.spain, g.flixbus), ==, approx(7.00409)); TEST(util::geo::webMercMeterDist(g.spainM, g.flixbusM), ==, - approx(703461.25144)); + approx(702675.06380)); auto line = lineFromWKTProj("LINESTRING(7.8824970 48.0228303,7.8823288 48.0227874,7.8820604 48.0227417,7.8819946 48.0227305)", util::geo::projectToWebMerc); auto lineX = XSortedLine(line); @@ -463,10 +463,10 @@ static void testDistComplexGeoms(const LargeTestGeoms& g) { util::geo::webMercMeterDist(g.vaubanMX, lineX)); TEST(util::geo::webMercMeterDist(g.vaubanM, line), ==, - approx(6449.59555)); + approx(6442.38749)); TEST(util::geo::webMercMeterDist(line, g.vaubanM), ==, - approx(6449.59555)); + approx(6442.38749)); } // _____________________________________________________________________________ @@ -480,7 +480,7 @@ static void testDistHaversineNoPadding(const LargeTestGeoms& g) { double) -> double { return haversineWebMerc(a, b); }) * 10.0) / 10.0, - ==, approx(653276.6)); + ==, approx(652546.5)); TEST(std::round(util::geo::withinDist( g.germanyMX, g.germanyMX, 1000000, defaultPaddingFunc(), 1e8, @@ -496,7 +496,7 @@ static void testDistHaversineNoPadding(const LargeTestGeoms& g) { double) -> double { return haversineWebMerc(a, b); }) * 10.0) / 10.0, - ==, approx(1082465.4)); + ==, approx(1081255.6)); TEST(std::round( util::geo::withinDist( @@ -507,7 +507,7 @@ static void testDistHaversineNoPadding(const LargeTestGeoms& g) { double) -> double { return haversineWebMerc(a, b); }) * 10.0) / 10.0, - ==, approx(653276.6)); + ==, approx(652546.5)); TEST(std::round( util::geo::withinDist( @@ -518,7 +518,7 @@ static void testDistHaversineNoPadding(const LargeTestGeoms& g) { double) -> double { return haversineWebMerc(a, b); }) * 10.0) / 10.0, - ==, approx(653276.6)); + ==, approx(652546.5)); TEST(std::round( util::geo::withinDist( @@ -528,7 +528,7 @@ static void testDistHaversineNoPadding(const LargeTestGeoms& g) { double) -> double { return haversineWebMerc(a, b); }) * 10.0) / 10.0, - ==, approx(653276.6)); + ==, approx(652546.5)); } // _____________________________________________________________________________ @@ -544,7 +544,7 @@ static void testDistHaversineSmallPadding(const LargeTestGeoms& g) { double) -> double { return haversineWebMerc(a, b); }) * 10.0) / 10.0, - ==, approx(653276.6)); + ==, approx(652546.5)); TEST(std::round(util::geo::withinDist( g.germanyMX, g.flixbusMX, 1000000, [](double d, double, const Box&, @@ -565,7 +565,7 @@ static void testDistHaversineSmallPadding(const LargeTestGeoms& g) { double) -> double { return haversineWebMerc(a, b); }) * 10.0) / 10.0, - ==, approx(703461.3)); + ==, approx(702675.1)); } // _____________________________________________________________________________ @@ -574,13 +574,13 @@ static void testDistHaversineMeterDistPadding(const LargeTestGeoms& g) { TEST(std::round(util::geo::webMercMeterDist(g.germanyMX, g.spainMX) * 10.0) / 10.0, - ==, approx(653276.6)); + ==, approx(652546.5)); TEST(std::round(util::geo::webMercMeterDist(g.germanyMX, g.flixbusMX) * 10.0) / 10.0, ==, approx(0)); TEST(std::round(util::geo::webMercMeterDist(g.spainMX, g.flixbusMX) * 10.0) / 10.0, - ==, approx(703461.3)); + ==, approx(702675.1)); // takes too long // TEST(std::round(util::geo::meterDist(g.germanyMX, // XSortedCollection(Collection{g.spainM, g.saimaaM})) * 10.0) @@ -752,7 +752,7 @@ static void testDistOther() { TEST(util::geo::dist(point2, line2), ==, approx(9.01734)); TEST(util::geo::dist(line, point), ==, approx(5.5)); - TEST(util::geo::webMercMeterDist(point2M, line2M), ==, approx(999138.32522)); + TEST(util::geo::webMercMeterDist(point2M, line2M), ==, approx(999138.68916)); TEST(util::geo::webMercMeterDist(lineM, pointM), ==, approx(610368.37082)); // standard point/polygon @@ -1157,8 +1157,99 @@ static void testDistToSegmentExtreme() { } } +// _____________________________________________________________________________ +static void testDistAdaptiveMeterDist() { + auto p = [](double lng, double lat) { + return latLngToWebMerc(Point{lng, lat}); + }; + + // tolerance 0.1 m + // Freiburg, 10 m -> haversine + TEST(fabs(adaptiveMeterDistWebMerc(p(7.842, 47.998), + p(7.842066999, 47.998077887), 0.1) - + 10.0000) <= 0.1); + // Freiburg, 5 km -> andoyerLambert + TEST(fabs(adaptiveMeterDistWebMerc(p(7.842, 47.998), + p(7.899997428, 47.975501331), 0.1) - + 5000.0000) <= 0.1); + // Freiburg, 100 km -> vincenty + TEST(fabs(adaptiveMeterDistWebMerc(p(7.842, 47.998), + p(7.391038751, 47.151920583), 0.1) - + 99999.9999) <= 0.1); + + // tolerance 1.0 m + // equator, 100 m -> haversine + TEST(fabs(adaptiveMeterDistWebMerc(p(0.0, 0.0), p(0.0, 0.000904369), 1.0) - + 99.9999) <= 1.0); + // Helsinki, 50 km -> andoyerLambert + TEST(fabs(adaptiveMeterDistWebMerc(p(24.94, 60.17), + p(25.813024072, 60.283283669), 1.0) - + 50000.0000) <= 1.0); + // Helsinki, 1000 km -> vincenty + TEST(fabs(adaptiveMeterDistWebMerc(p(24.94, 60.17), + p(17.645254308, 52.150404569), 1.0) - + 1000000.0000) <= 1.0); + + // tolerance 10.0 m + // Cape Town, 1 km -> haversine + TEST(fabs(adaptiveMeterDistWebMerc(p(18.424, -33.925), + p(18.431646647, -33.918624893), 10.0) - + 1000.0000) <= 10.0); + // Cape Town, 500 km -> andoyerLambert + TEST(fabs(adaptiveMeterDistWebMerc(p(18.424, -33.925), + p(23.793744594, -34.590866672), 10.0) - + 500000.0000) <= 10.0); + // Cape Town, 5000 km -> vincenty + TEST(fabs(adaptiveMeterDistWebMerc(p(18.424, -33.925), + p(32.574341805, 9.183855187), 10.0) - + 5000000.0000) <= 10.0); + + // tolerance 100.0 m + // Tokyo, 10 km -> haversine + TEST(fabs(adaptiveMeterDistWebMerc(p(139.69, 35.69), + p(139.59426886, 35.735025795), 100.0) - + 10000.0000) <= 100.0); + // Tokyo, 5000 km -> andoyerLambert + TEST(fabs(adaptiveMeterDistWebMerc(p(139.69, 35.69), + p(101.657485041, 7.123964981), 100.0) - + 5000000.0000) <= 100.0); + // Tokyo, 12000 km -> vincenty + TEST(fabs(adaptiveMeterDistWebMerc(p(139.69, 35.69), + p(-97.87522015, 11.949948055), 100.0) - + 12000000.0000) <= 100.0); + + // tolerance 1000.0 m + // Buenos Aires, 100 km -> haversine + TEST(fabs(adaptiveMeterDistWebMerc(p(-58.38, -34.6), + p(-58.100741115, -33.728890554), 1000.0) - + 99999.9999) <= 1000.0); + // Buenos Aires, 5000 km -> andoyerLambert + TEST(fabs(adaptiveMeterDistWebMerc(p(-58.38, -34.6), + p(-72.504836992, 8.527465826), 1000.0) - + 5000000.0000) <= 1000.0); + // Buenos Aires, 12000 km -> vincenty + TEST(fabs(adaptiveMeterDistWebMerc(p(-58.38, -34.6), + p(35.958128097, 26.440931366), 1000.0) - + 12000000.0000) <= 1000.0); + + // tolerance 10000.0 m + // Reykjavik, 1000 km -> haversine + TEST(fabs(adaptiveMeterDistWebMerc(p(-21.942, 64.146), + p(-19.21998632, 55.273720353), 10000.0) - + 1000000.0000) <= 10000.0); + // Reykjavik, 5000 km -> andoyerLambert + TEST(fabs(adaptiveMeterDistWebMerc(p(-21.942, 64.146), + p(-36.852837223, 20.269614079), 10000.0) - + 5000000.0000) <= 10000.0); + // Reykjavik, 12000 km -> vincenty + TEST(fabs(adaptiveMeterDistWebMerc(p(-21.942, 64.146), + p(16.142559896, -39.804652894), 10000.0) - + 12000000.0000) <= 10000.0); +} + // _____________________________________________________________________________ void GeoTest::testDist() { + testDistAdaptiveMeterDist(); testDistCombinations(); testDistWithinRealisticMaxDist(); testDistWithinInfinityMaxDist(); From d0ab5069a06738c3bcc20794226346b7e169d905 Mon Sep 17 00:00:00 2001 From: Patrick Brosi Date: Fri, 21 Aug 2026 22:36:38 +0200 Subject: [PATCH 3/4] fix tests --- geo/Geo.tpp | 18 +++++++++++------- tests/GeoTestDist.cpp | 32 ++++++++++++++++---------------- 2 files changed, 27 insertions(+), 23 deletions(-) diff --git a/geo/Geo.tpp b/geo/Geo.tpp index f796d14..0f5e298 100644 --- a/geo/Geo.tpp +++ b/geo/Geo.tpp @@ -6054,7 +6054,8 @@ double haversineWebMerc(const Point& a, const Point& b) { // _____________________________________________________________________________ template double andoyerLambert(T lat1, T lon1, T lat2, T lon2) { - // see https://en.wikipedia.org/wiki/Geographical_distance#Andoyer-Lambert_formula_for_long_lines + // see + // https://en.wikipedia.org/wiki/Geographical_distance#Andoyer-Lambert_formula_for_long_lines double f1 = 1.0 - FLATTENING; double b1 = atan(f1 * tan(lat1 * RAD)); @@ -6605,7 +6606,7 @@ double webMercDistFactor(const G& a) { // euclidean distance on web mercator is in meters on equator, // and proportional to cos(lat) in both y directions double et = exp(a.getY() / EQUATORIAL_RAD); - return 2 * et / (et * et + 1); + return (MEAN_EARTH_RAD / EQUATORIAL_RAD) * 2 * et / (et * et + 1); } // _____________________________________________________________________________ @@ -7395,12 +7396,14 @@ std::pair getMinMaxLocalScaleFactors( -90.0 + util::geo::EPSILON, std::min(withinUp.getY() * 1.0, std::min(aUp.getY(), bUp.getY()))); - double a = cos(yRangeMin * util::geo::RAD); - double b = cos(yRangeMax * util::geo::RAD); + const double r = util::geo::MEAN_EARTH_RAD / util::geo::EQUATORIAL_RAD; - // if we crossed the equator, we encountered a scale factor of 1! + double a = r * cos(yRangeMin * util::geo::RAD); + double b = r * cos(yRangeMax * util::geo::RAD); + + // if we crossed the equator, we encountered the maximum scale factor! if (yRangeMin < 0 && yRangeMax > 0) { - return {std::min(a, b), std::max(1.0, std::max(a, b))}; + return {std::min(a, b), std::max(r, std::max(a, b))}; } return {std::min(a, b), std::max(a, b)}; @@ -7428,7 +7431,8 @@ double webMercMaxEuclideanDist(const Box& boxA, const Box& boxB, double maxD) { auto scale = getMinMaxLocalScaleFactorsWebMerc(boxA, boxB, maxD); - // use meters here directly, we are in web mercator world + // use meters here directly, we are in web mercator world, but acknowledge + // haversine error return (maxD * (1.0 + HAVERSINE_MAX_ERR)) / scale.first; } diff --git a/tests/GeoTestDist.cpp b/tests/GeoTestDist.cpp index 1d8b4a4..ed7795c 100644 --- a/tests/GeoTestDist.cpp +++ b/tests/GeoTestDist.cpp @@ -752,8 +752,8 @@ static void testDistOther() { TEST(util::geo::dist(point2, line2), ==, approx(9.01734)); TEST(util::geo::dist(line, point), ==, approx(5.5)); - TEST(util::geo::webMercMeterDist(point2M, line2M), ==, approx(999138.68916)); - TEST(util::geo::webMercMeterDist(lineM, pointM), ==, approx(610368.37082)); + TEST(util::geo::webMercMeterDist(point2M, line2M), ==, approx(998021.68916)); + TEST(util::geo::webMercMeterDist(lineM, pointM), ==, approx(609686.22369)); // standard point/polygon TEST(util::geo::withinDist(point, polyWithInner, 1), ==, approx(0.5)); @@ -763,9 +763,9 @@ static void testDistOther() { TEST(util::geo::dist(polyWithInner, point), ==, approx(0.5)); TEST(util::geo::webMercMeterDist(pointM, polyWithInnerM), ==, - approx(55488.16389)); + approx(55426.15037)); TEST(util::geo::webMercMeterDist(polyWithInnerM, pointM), ==, - approx(55488.16389)); + approx(55426.15037)); // standard line/polygon TEST(util::geo::withinDist(line2, polyWithInner, 100), ==, approx(0.25)); @@ -775,9 +775,9 @@ static void testDistOther() { TEST(util::geo::dist(polyWithInner, line2), ==, approx(0.25)); TEST(util::geo::webMercMeterDist(line2M, polyWithInnerM), ==, - approx(27744.08235)); + approx(27713.07559)); TEST(util::geo::webMercMeterDist(polyWithInnerM, line2M), ==, - approx(27744.08235)); + approx(27713.07559)); // standard line/line auto segLineA = lineFromWKT("LINESTRING(0 0, 1 0)"); @@ -893,9 +893,9 @@ static void testDistOther() { auto webMercLineA = lineFromWKT("LINESTRING(0 0, 10 0)"); auto webMercLineB = lineFromWKT("LINESTRING(20 0, 30 0)"); TEST(util::geo::webMercMeterDist(webMercLineA, webMercLineB), ==, - approx(10.0)); + approx(9.98882)); TEST(util::geo::webMercMeterDist(webMercLineB, webMercLineA), ==, - approx(10.0)); + approx(9.98882)); // haversine on web mercator coordinates must match haversine on the // corresponding lat/lng coordinates @@ -1020,7 +1020,7 @@ static void testDistLimitedPrecision() { util::geo::projectToWebMerc); TEST(util::geo::webMercMeterDist(germanyCoarseRaw, londonCoarseRaw), ==, - approx(426521.22769)); + approx(426044.54796)); TEST( util::geo::withinDist( @@ -1048,7 +1048,7 @@ static void testDistLimitedPrecision() { }), // NOTE: difference because of precision to only 10 cm because of coarse // projection - ==, approx(426521.18896)); + ==, approx(426044.50928)); } // _____________________________________________________________________________ @@ -1117,8 +1117,8 @@ static void testDistToSegmentExtreme() { latLngToWebMerc(Point{60.0, 67.410787}); double ref = haversineWebMerc(p, seg.front()); - TEST(ref, >, 2832051.0); - TEST(ref, <, 2832052.0); + TEST(ref, >, 2828886.0); + TEST(ref, <, 2828887.0); // FAILS with a significant error! // TEST(util::geo::webMercMeterDist(p, seg), ==, approx(ref)); @@ -1134,8 +1134,8 @@ static void testDistToSegmentExtreme() { double ref = haversineWebMerc(p, latLngToWebMerc(Point{0.278391, 76.962378})); - TEST(ref, >, 193127.0); - TEST(ref, <, 193128.0); + TEST(ref, >, 192911.0); + TEST(ref, <, 192912.0); // FAILS with a small error // TEST(util::geo::webMercMeterDist(p, seg), ==, approx(ref)); @@ -1149,8 +1149,8 @@ static void testDistToSegmentExtreme() { auto p = latLngToWebMerc(Point{150.0, 84.994120}); double ref = haversineWebMerc(p, seg.front()); - TEST(ref, >, 556393.0); - TEST(ref, <, 556395.0); + TEST(ref, >, 555772.0); + TEST(ref, <, 555773.0); // FAILS with a huge error, as expected // TEST(util::geo::webMercMeterDist(p, seg), ==, approx(ref)); From 5d7c2401dd7edfc75f01276240e76cd6afb9edb6 Mon Sep 17 00:00:00 2001 From: Patrick Brosi Date: Fri, 21 Aug 2026 22:40:39 +0200 Subject: [PATCH 4/4] fix typo comparison --- geo/Geo.tpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geo/Geo.tpp b/geo/Geo.tpp index 0f5e298..1abf339 100644 --- a/geo/Geo.tpp +++ b/geo/Geo.tpp @@ -6111,7 +6111,7 @@ double andoyerLambertWebMerc(const Point& a, const Point& b) { template double vincenty(T lat1, T lon1, T lat2, T lon2) { // see https://en.wikipedia.org/wiki/Vincenty's_formulae - size_t MAX_ITERS = 200; + int MAX_ITERS = 200; double f1 = 1.0 - FLATTENING; double b = EQUATORIAL_RAD * f1;