From 727e3a1ca75f0dd7eed0a8b405566e316f700ae3 Mon Sep 17 00:00:00 2001 From: Artem Pavlenko Date: Fri, 26 Jul 2024 14:28:41 +0100 Subject: [PATCH 1/4] Split bounding box along (0, 0) before reprojecting [WIP] [skip ci] (improves handling stereographic projections e.g epsg:3995) --- src/proj_transform.cpp | 119 +++++++++++++++++++++++++++++------------ 1 file changed, 84 insertions(+), 35 deletions(-) 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; } From 7358a66c2ffc288823bde5a5dc46395e245369fa Mon Sep 17 00:00:00 2001 From: Artem Pavlenko Date: Mon, 29 Jul 2024 11:21:54 +0100 Subject: [PATCH 2/4] mapnik::projection - add `area_of_use` method returning `std::optional>` (WGS84) [WIP] [skip ci] --- include/mapnik/projection.hpp | 9 +++++---- src/projection.cpp | 12 +++++++++++- 2 files changed, 16 insertions(+), 5 deletions(-) diff --git a/include/mapnik/projection.hpp b/include/mapnik/projection.hpp index dd3183889..032b1370b 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,15 +79,16 @@ 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); 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/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_); From 4f6a7a25715d6d50382b08643a8b81206e9207a4 Mon Sep 17 00:00:00 2001 From: Artem Pavlenko Date: Mon, 29 Jul 2024 11:23:39 +0100 Subject: [PATCH 3/4] Add basic forward/backward test for reprojection bounding box for `epsg:3995` (WGS 84 / Arctic Polar Stereographic) --- test/unit/projection/proj_transform.cpp | 42 +++++++++++++++++++++++-- 1 file changed, 40 insertions(+), 2 deletions(-) diff --git a/test/unit/projection/proj_transform.cpp b/test/unit/projection/proj_transform.cpp index 7e3021ad1..0b1829939 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,43 @@ 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 +276,6 @@ TEST_CASE("projection transform") } } -#ifdef MAPNIK_USE_PROJ SECTION("lonlat <-> New Zealand Transverse Mercator 2000") { mapnik::projection const proj_2193("epsg:2193"); From a982e6952577849e9f9a4c76f8ec8c18b912276d Mon Sep 17 00:00:00 2001 From: Artem Pavlenko Date: Mon, 29 Jul 2024 12:35:31 +0100 Subject: [PATCH 4/4] clang-format --- include/mapnik/projection.hpp | 1 + test/unit/projection/proj_transform.cpp | 9 ++++----- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/include/mapnik/projection.hpp b/include/mapnik/projection.hpp index 032b1370b..5a057b263 100644 --- a/include/mapnik/projection.hpp +++ b/include/mapnik/projection.hpp @@ -81,6 +81,7 @@ class MAPNIK_DECL projection std::string description() const; void init_proj(); std::optional> area_of_use() const; + private: void swap(projection& rhs); diff --git a/test/unit/projection/proj_transform.cpp b/test/unit/projection/proj_transform.cpp index 0b1829939..d5ae1cbbf 100644 --- a/test/unit/projection/proj_transform.cpp +++ b/test/unit/projection/proj_transform.cpp @@ -210,15 +210,14 @@ TEST_CASE("projection transform") 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}; + 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}; + // 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())); @@ -230,7 +229,7 @@ TEST_CASE("projection transform") { // check the same logic works for .backward() - mapnik::box2d projected_bounds{-3299207.53,-3333134.03, 3299207.53, 3333134.03}; + 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()));