diff --git a/include/mapnik/projection.hpp b/include/mapnik/projection.hpp index dd3183889..5a057b263 100644 --- a/include/mapnik/projection.hpp +++ b/include/mapnik/projection.hpp @@ -26,7 +26,7 @@ // mapnik #include #include - +#include #include MAPNIK_DISABLE_WARNING_PUSH #include @@ -79,7 +79,8 @@ class MAPNIK_DECL projection void inverse(double& x, double& y) const; std::string definition() const; std::string description() const; - void init_proj() const; + void init_proj(); + std::optional> area_of_use() const; private: void swap(projection& rhs); @@ -87,7 +88,8 @@ class MAPNIK_DECL projection private: std::string params_; bool defer_proj_init_; - mutable bool is_geographic_; + bool is_geographic_; + std::optional> area_of_use_; mutable PJ* proj_; mutable PJ_CONTEXT* proj_ctx_; }; diff --git a/src/proj_transform.cpp b/src/proj_transform.cpp index 6ae85fb5f..26ac44118 100644 --- a/src/proj_transform.cpp +++ b/src/proj_transform.cpp @@ -42,7 +42,6 @@ #include namespace mapnik { - namespace { // (local) // Returns points in clockwise order. This allows us to do anti-meridian checks. @@ -90,6 +89,36 @@ auto envelope_points(box2d const& env, std::size_t num_points) -> geometry::m return coords; } +std::vector> split_bounding_box(box2d const& bbox) +{ + std::vector> output; + if (bbox.minx() < 0 && bbox.maxx() > 0) + { + if (bbox.miny() < 0 && bbox.maxy() > 0) + { + output.emplace_back(bbox.minx(), bbox.miny(), 0, 0); + output.emplace_back(0, bbox.miny(), bbox.maxx(), 0); + output.emplace_back(0, 0, bbox.maxx(), bbox.maxy()); + output.emplace_back(bbox.minx(), 0, 0, bbox.maxy()); + } + else + { + output.emplace_back(bbox.minx(), bbox.miny(), 0, bbox.maxy()); + output.emplace_back(0, bbox.miny(), bbox.maxx(), bbox.maxy()); + } + } + else if (bbox.miny() < 0 && bbox.maxy() > 0) + { + output.emplace_back(bbox.minx(), bbox.miny(), bbox.maxx(), 0); + output.emplace_back(bbox.minx(), 0, bbox.maxx(), bbox.maxy()); + } + else + { + output.push_back(bbox); + } + return output; +} + } // namespace proj_transform::proj_transform(projection const& source, projection const& dest) @@ -398,7 +427,7 @@ bool proj_transform::backward(box2d& box) const } // More robust, but expensive, bbox transform -// in the face of proj4 out of bounds conditions. +// in the face of libproj out of bounds conditions. // Can result in 20 -> 10 r/s performance hit. // Alternative is to provide proper clipping box // in the target srs by setting map 'maximum-extent' @@ -413,27 +442,34 @@ bool proj_transform::backward(box2d& env, std::size_t points) const return backward(env); } - auto coords = envelope_points(env, points); // this is always clockwise - - for (auto& p : coords) - { - double z = 0; - if (!backward(p.x, p.y, z)) - return false; - } - box2d result; - boost::geometry::envelope(coords, result); - - if (is_source_longlat_ && !util::is_clockwise(coords)) + for (auto b : split_bounding_box(env)) { - // we've gone to a geographic CS, and our clockwise envelope has - // changed into an anticlockwise one. This means we've crossed the antimeridian, and - // need to expand the X direction to +/-180 to include all the data. Once we can deal - // with multiple bboxes in queries we can improve. - double miny = result.miny(); - result.expand_to_include(-180.0, miny); - result.expand_to_include(180.0, miny); + auto coords = envelope_points(b, points); // this is always clockwise + for (auto& p : coords) + { + double z = 0; + if (!backward(p.x, p.y, z)) + return false; + } + if (!result.valid()) + boost::geometry::envelope(coords, result); + else + { + box2d bb; + boost::geometry::envelope(coords, bb); + result.expand_to_include(bb); + } + if (is_source_longlat_ && !util::is_clockwise(coords)) + { + // we've gone to a geographic CS, and our clockwise envelope has + // changed into an anticlockwise one. This means we've crossed the antimeridian, and + // need to expand the X direction to +/-180 to include all the data. Once we can deal + // with multiple bboxes in queries we can improve. + double miny = result.miny(); + result.expand_to_include(-180.0, miny); + result.expand_to_include(180.0, miny); + } } env.re_center(result.center().x, result.center().y); @@ -453,32 +489,45 @@ bool proj_transform::forward(box2d& env, std::size_t points) const } auto coords = envelope_points(env, points); // this is always clockwise - + double z = 0; for (auto& p : coords) { - double z = 0; if (!forward(p.x, p.y, z)) return false; } box2d result; - boost::geometry::envelope(coords, result); - - if (is_dest_longlat_ && !util::is_clockwise(coords)) + for (auto b : split_bounding_box(env)) { - // we've gone to a geographic CS, and our clockwise envelope has - // changed into an anticlockwise one. This means we've crossed the antimeridian, and - // need to expand the X direction to +/-180 to include all the data. Once we can deal - // with multiple bboxes in queries we can improve. - double miny = result.miny(); - result.expand_to_include(-180.0, miny); - result.expand_to_include(180.0, miny); + auto coords = envelope_points(b, points); // this is always clockwise + for (auto& p : coords) + { + double z = 0; + if (!forward(p.x, p.y, z)) + return false; + } + if (!result.valid()) + boost::geometry::envelope(coords, result); + else + { + box2d bb; + boost::geometry::envelope(coords, bb); + result.expand_to_include(bb); + } + if (is_dest_longlat_ && !util::is_clockwise(coords)) + { + // we've gone to a geographic CS, and our clockwise envelope has + // changed into an anticlockwise one. This means we've crossed the antimeridian, and + // need to expand the X direction to +/-180 to include all the data. Once we can deal + // with multiple bboxes in queries we can improve. + double miny = result.miny(); + result.expand_to_include(-180.0, miny); + result.expand_to_include(180.0, miny); + } } - env.re_center(result.center().x, result.center().y); env.height(result.height()); env.width(result.width()); - return true; } diff --git a/src/projection.cpp b/src/projection.cpp index 29cffb4d8..6b434a175 100644 --- a/src/projection.cpp +++ b/src/projection.cpp @@ -94,7 +94,7 @@ bool projection::operator!=(const projection& other) const return !(*this == other); } -void projection::init_proj() const +void projection::init_proj() { #ifdef MAPNIK_USE_PROJ if (!proj_) @@ -117,6 +117,11 @@ void projection::init_proj() const } PJ_TYPE type = proj_get_type(proj_); is_geographic_ = (type == PJ_TYPE_GEOGRAPHIC_2D_CRS || type == PJ_TYPE_GEOGRAPHIC_3D_CRS) ? true : false; + double west_lon, south_lat, east_lon, north_lat; + if (proj_get_area_of_use(proj_ctx_, proj_, &west_lon, &south_lat, &east_lon, &north_lat, nullptr)) + { + area_of_use_ = box2d{west_lon, south_lat, east_lon, north_lat}; + } } #endif } @@ -131,6 +136,11 @@ bool projection::is_geographic() const return is_geographic_; } +std::optional> projection::area_of_use() const +{ + return area_of_use_; +} + std::optional projection::well_known() const { return is_well_known_srs(params_); diff --git a/test/unit/projection/proj_transform.cpp b/test/unit/projection/proj_transform.cpp index 7e3021ad1..d5ae1cbbf 100644 --- a/test/unit/projection/proj_transform.cpp +++ b/test/unit/projection/proj_transform.cpp @@ -16,6 +16,8 @@ TEST_CASE("projection transform") CHECK(proj_4326.is_geographic()); CHECK(!proj_3857.is_geographic()); + CHECK(*proj_4326.area_of_use() == mapnik::box2d(-180, -90, 180, 90)); + CHECK(*proj_3857.area_of_use() == mapnik::box2d(-180, -85.06, 180, 85.06)); mapnik::proj_transform prj_trans(proj_4326, proj_3857); double minx = -45.0; @@ -48,6 +50,7 @@ TEST_CASE("projection transform") { mapnik::projection proj_4269("epsg:4269"); mapnik::projection proj_3857("epsg:3857"); + mapnik::proj_transform prj_trans(proj_4269, proj_3857); mapnik::proj_transform prj_trans2(proj_3857, proj_4269); @@ -200,7 +203,42 @@ TEST_CASE("projection transform") } } } -#endif + SECTION("Test proj north/south poles bbox") + { + // WGS 84 / Arctic Polar Stereographic + + mapnik::projection prj_geog("epsg:4326"); + mapnik::projection prj_proj("epsg:3995"); + + mapnik::proj_transform prj_trans_fwd(prj_proj, prj_geog); + mapnik::proj_transform prj_trans_rev(prj_geog, prj_proj); + + // bounds + const mapnik::box2d bounds{-180.0, 60.0, 180.0, 90.0}; + { + // projected bounds + mapnik::box2d projected_bounds{-3299207.53, -3333134.03, 3299207.53, 3333134.03}; + CHECKED_IF(prj_trans_fwd.forward(projected_bounds, PROJ_ENVELOPE_POINTS)) + { + CHECK(projected_bounds.minx() == Approx(bounds.minx())); + CHECK(projected_bounds.miny() == Approx(48.65640)); // this is expected + CHECK(projected_bounds.maxx() == Approx(bounds.maxx())); + CHECK(projected_bounds.maxy() == Approx(bounds.maxy())); + } + } + + { + // check the same logic works for .backward() + mapnik::box2d projected_bounds{-3299207.53, -3333134.03, 3299207.53, 3333134.03}; + CHECKED_IF(prj_trans_rev.backward(projected_bounds, PROJ_ENVELOPE_POINTS)) + { + CHECK(projected_bounds.minx() == Approx(bounds.minx())); + CHECK(projected_bounds.miny() == Approx(48.65640)); // this is expected + CHECK(projected_bounds.maxx() == Approx(bounds.maxx())); + CHECK(projected_bounds.maxy() == Approx(bounds.maxy())); + } + } + } SECTION("proj_transform of coordinate arrays with stride > 1") { @@ -237,7 +275,6 @@ TEST_CASE("projection transform") } } -#ifdef MAPNIK_USE_PROJ SECTION("lonlat <-> New Zealand Transverse Mercator 2000") { mapnik::projection const proj_2193("epsg:2193");