Fix possible division by zero by clamping latitude to 85.05°

Resolves #3530

Author: oxidase

include/util/web_mercator.hpp

@@ -20,7 +20,7 @@ const constexpr double EARTH_RADIUS_WGS84 = 6378137.0;
 // earth circumference devided by 2
 const constexpr double MAXEXTENT = EARTH_RADIUS_WGS84 * boost::math::constants::pi<double>();
 // ^ math functions are not constexpr since they have side-effects (setting errno) :(
-const constexpr double MAX_LATITUDE = 85.;
+const constexpr double EPSG3857_MAX_LATITUDE = 85.051128779806592378; // 90(4*atan(exp(pi))/pi-1)
 const constexpr double MAX_LONGITUDE = 180.0;
 }
 
@@ -29,6 +29,18 @@ const constexpr double DEGREE_TO_PX = detail::MAXEXTENT / 180.0;
 // This is the global default tile size for all Mapbox Vector Tiles
 const constexpr double TILE_SIZE = 256.0;
 
+inline FloatLatitude clamp(const FloatLatitude lat)
+{
+    return std::max(std::min(lat, FloatLatitude{detail::EPSG3857_MAX_LATITUDE}),
+                    FloatLatitude{-detail::EPSG3857_MAX_LATITUDE});
+}
+
+inline FloatLongitude clamp(const FloatLongitude lon)
+{
+    return std::max(std::min(lon, FloatLongitude{detail::MAX_LONGITUDE}),
+                    FloatLongitude{-detail::MAX_LONGITUDE});
+}
+
 inline FloatLatitude yToLat(const double y)
 {
     const auto clamped_y = std::max(-180., std::min(180., y));
@@ -41,10 +53,9 @@ inline FloatLatitude yToLat(const double y)
 inline double latToY(const FloatLatitude latitude)
 {
     // apparently this is the (faster) version of the canonical log(tan()) version
-    const double f = std::sin(detail::DEGREE_TO_RAD * static_cast<double>(latitude));
-    const double y = detail::RAD_TO_DEGREE * 0.5 * std::log((1 + f) / (1 - f));
-    const auto clamped_y = std::max(-180., std::min(180., y));
-    return clamped_y;
+    const auto clamped_latitude = clamp(latitude);
+    const double f = std::sin(detail::DEGREE_TO_RAD * static_cast<double>(clamped_latitude));
+    return detail::RAD_TO_DEGREE * 0.5 * std::log((1 + f) / (1 - f));
 }
 
 template <typename T> constexpr double horner(double, T an) { return an; }
@@ -91,18 +102,6 @@ inline double latToYapprox(const FloatLatitude latitude)
                   -3.23083224835967391884404730e-28);
 }
 
-inline FloatLatitude clamp(const FloatLatitude lat)
-{
-    return std::max(std::min(lat, FloatLatitude{detail::MAX_LATITUDE}),
-                    FloatLatitude{-detail::MAX_LATITUDE});
-}
-
-inline FloatLongitude clamp(const FloatLongitude lon)
-{
-    return std::max(std::min(lon, FloatLongitude{detail::MAX_LONGITUDE}),
-                    FloatLongitude{-detail::MAX_LONGITUDE});
-}
-
 inline void pixelToDegree(const double shift, double &x, double &y)
 {
     const double b = shift / 2.0;
@@ -166,9 +165,9 @@ inline void xyzToMercator(
     xyzToWGS84(x, y, z, minx, miny, maxx, maxy);
 
     minx = static_cast<double>(clamp(util::FloatLongitude{minx})) * DEGREE_TO_PX;
-    miny = latToY(clamp(util::FloatLatitude{miny})) * DEGREE_TO_PX;
+    miny = latToY(util::FloatLatitude{miny}) * DEGREE_TO_PX;
     maxx = static_cast<double>(clamp(util::FloatLongitude{maxx})) * DEGREE_TO_PX;
-    maxy = latToY(clamp(util::FloatLatitude{maxy})) * DEGREE_TO_PX;
+    maxy = latToY(util::FloatLatitude{maxy}) * DEGREE_TO_PX;
 }
 }
 }

unit_tests/util/coordinate_calculation.cpp

@@ -117,23 +117,20 @@ BOOST_AUTO_TEST_CASE(compute_angle)
     BOOST_CHECK_EQUAL(angle, 180);
 
     // Tiny changes below our calculation resolution
-    // This should be equivalent to having two points on the same
-    // spot.
+    // This should be equivalent to having two points on the same spot.
     first = Coordinate{FloatLongitude{0}, FloatLatitude{0}};
     middle = Coordinate{FloatLongitude{1}, FloatLatitude{0}};
     end = Coordinate{FloatLongitude{1 + std::numeric_limits<double>::epsilon()}, FloatLatitude{0}};
     angle = coordinate_calculation::computeAngle(first, middle, end);
     BOOST_CHECK_EQUAL(angle, 180);
 
     // Invalid values
-    /* TODO: Enable this when I figure out how to use BOOST_CHECK_THROW
-     *       and not have the whole test case fail...
-    first = Coordinate(FloatLongitude{0}, FloatLatitude{0});
-    middle = Coordinate(FloatLongitude{1}, FloatLatitude{0});
-    end = Coordinate(FloatLongitude(std::numeric_limits<double>::max()), FloatLatitude{0});
-    BOOST_CHECK_THROW( coordinate_calculation::computeAngle(first,middle,end),
-                       boost::numeric::positive_overflow);
-                       */
+    BOOST_CHECK_THROW(
+        coordinate_calculation::computeAngle(
+            Coordinate(FloatLongitude{0}, FloatLatitude{0}),
+            Coordinate(FloatLongitude{1}, FloatLatitude{0}),
+            Coordinate(FloatLongitude{std::numeric_limits<double>::max()}, FloatLatitude{0})),
+        boost::numeric::positive_overflow);
 }
 
 // Regression test for bug captured in #1347

unit_tests/util/web_mercator.cpp

@@ -73,12 +73,14 @@ BOOST_AUTO_TEST_CASE(xyz_to_mercator)
     double miny;
     double maxx;
     double maxy;
+
+    // http://tools.geofabrik.de/map/#13/85.0500/-175.5876&type=Geofabrik_Standard&grid=1
     web_mercator::xyzToMercator(100, 0, 13, minx, miny, maxx, maxy);
 
     BOOST_CHECK_CLOSE(minx, -19548311.361764118075, 0.0001);
-    BOOST_CHECK_CLOSE(miny, 19971868.8804085782, 0.0001);
+    BOOST_CHECK_CLOSE(miny, 20032616.372979045, 0.0001);
     BOOST_CHECK_CLOSE(maxx, -19543419.391953866929, 0.0001);
-    BOOST_CHECK_CLOSE(maxy, 19971868.880408578, 0.0001);
+    BOOST_CHECK_CLOSE(maxy, 20037508.342789277, 0.0001); // Mercator 6378137*pi, WGS 85.0511
 }
 
 BOOST_AUTO_TEST_SUITE_END()