diff --git a/include/util/web_mercator.hpp b/include/util/web_mercator.hpp index 1bc4b1ca4..560d7884f 100644 --- a/include/util/web_mercator.hpp +++ b/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(); // ^ 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(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(clamped_latitude)); + return detail::RAD_TO_DEGREE * 0.5 * std::log((1 + f) / (1 - f)); } template 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(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(clamp(util::FloatLongitude{maxx})) * DEGREE_TO_PX; - maxy = latToY(clamp(util::FloatLatitude{maxy})) * DEGREE_TO_PX; + maxy = latToY(util::FloatLatitude{maxy}) * DEGREE_TO_PX; } } } diff --git a/unit_tests/util/coordinate_calculation.cpp b/unit_tests/util/coordinate_calculation.cpp index 7bd256b5b..103943337 100644 --- a/unit_tests/util/coordinate_calculation.cpp +++ b/unit_tests/util/coordinate_calculation.cpp @@ -118,8 +118,7 @@ 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::epsilon()}, FloatLatitude{0}}; @@ -127,14 +126,12 @@ BOOST_AUTO_TEST_CASE(compute_angle) 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::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::max()}, FloatLatitude{0})), + boost::numeric::positive_overflow); } // Regression test for bug captured in #1347 diff --git a/unit_tests/util/web_mercator.cpp b/unit_tests/util/web_mercator.cpp index 4a12c95cd..308a7c91b 100644 --- a/unit_tests/util/web_mercator.cpp +++ b/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()