From ed7e22faf942895dad102eab6ffbb62a09f15e6f Mon Sep 17 00:00:00 2001 From: program-- Date: Fri, 28 Apr 2023 17:20:53 -0700 Subject: [PATCH] lots of refactoring -- streamlined wkb -> geojson types [no ci] --- include/geopackage/GeoPackage.hpp | 167 ++++++++++++++++++++++++++---- include/geopackage/wkb/reader.hpp | 135 +++++++++++++++--------- include/geopackage/wkb/wkt.hpp | 93 ----------------- test/geopackage/WKB_Test.cpp | 33 +++--- 4 files changed, 245 insertions(+), 183 deletions(-) delete mode 100644 include/geopackage/wkb/wkt.hpp diff --git a/include/geopackage/GeoPackage.hpp b/include/geopackage/GeoPackage.hpp index ece58c5c23..f226957d72 100644 --- a/include/geopackage/GeoPackage.hpp +++ b/include/geopackage/GeoPackage.hpp @@ -1,18 +1,28 @@ #ifndef NGEN_GEOPACKAGE_H #define NGEN_GEOPACKAGE_H +#include #include #include #include #include #include -#include -#include "SQLite.hpp" +#include +#include + #include "FeatureCollection.hpp" +#include "JSONGeometry.hpp" +#include "SQLite.hpp" +#include "features/CollectionFeature.hpp" +#include "wkb/reader.hpp" + +namespace bsrs = boost::geometry::srs; + namespace geopackage { +namespace { inline const geojson::FeatureType feature_type_map(const std::string& g) { if (g == "POINT") return geojson::FeatureType::Point; @@ -23,30 +33,66 @@ inline const geojson::FeatureType feature_type_map(const std::string& g) if (g == "MULTIPOLYGON") return geojson::FeatureType::MultiPolygon; return geojson::FeatureType::GeometryCollection; } +} // anonymous namespace inline geojson::geometry build_geometry(const sqlite_iter& row, const geojson::FeatureType geom_type, const std::string& geom_col) { const std::vector geometry_blob = row.get>(geom_col); + int index = 0; + if (geometry_blob[0] != 'G' && geometry_blob[1] != 'P') { + throw std::runtime_error("expected geopackage WKB, but found invalid format instead"); + } + index += 2; + + // skip version + index++; - geojson::geometry geometry; + // flags + const bool is_extended = geometry_blob[index] & 0x00100000; + const bool is_empty = geometry_blob[index] & 0x00010000; + const uint8_t indicator = (geometry_blob[index] & 0x00001110) >> 1; + const uint8_t endian = geometry_blob[index] & 0x00000001; + index++; - switch(geom_type) { - case geojson::FeatureType::Point: - break; - case geojson::FeatureType::LineString: - break; - case geojson::FeatureType::Polygon: - break; - case geojson::FeatureType::MultiPoint: - break; - case geojson::FeatureType::MultiLineString: - break; - case geojson::FeatureType::MultiPolygon: - break; - case geojson::FeatureType::GeometryCollection: - break; - default: - break; + // Read srs_id + uint32_t srs_id; + wkb::copy_from(geometry_blob, index, srs_id, endian); + + std::vector envelope; // may be unused + if (indicator > 0 & indicator < 5) { + // not an empty envelope + + envelope.resize(4); // only 4, not supporting Z or M dims + wkb::copy_from(geometry_blob, index, envelope[0], endian); + wkb::copy_from(geometry_blob, index, envelope[1], endian); + wkb::copy_from(geometry_blob, index, envelope[2], endian); + wkb::copy_from(geometry_blob, index, envelope[3], endian); + + // ensure `index` is at beginning of data + if (indicator == 2 || indicator == 3) { + index += 2 * sizeof(double); + } else if (indicator == 4) { + index += 4 * sizeof(double); + } + } + + const std::vector geometry_data(geometry_blob.begin() + index, geometry_blob.end()); + geojson::geometry geometry = wkb::read_wkb(geometry_data); + + if (srs_id != 4326) { + return geometry; + } else { + geojson::geometry projected; + + // project coordinates from whatever they are to 4326 + boost::geometry::srs::transformation<> tr{ + bsrs::epsg(srs_id), + bsrs::epsg(4326) + }; + + tr.forward(geometry, projected); + + return projected; } } @@ -101,9 +147,84 @@ inline geojson::Feature build_feature( const std::string& geom_col ) { - const auto type = feature_type_map(geom_type); - const geojson::PropertyMap properties = build_properties(row, geom_col); - const geojson::geometry geometry = build_geometry(row, type, geom_col); + const auto type = feature_type_map(geom_type); + const auto id = row.get("id"); + std::vector bounding_box = {0, 0, 0, 0}; // TODO + geojson::PropertyMap properties = std::move(build_properties(row, geom_col)); + geojson::geometry geometry = std::move(build_geometry(row, type, geom_col)); + + switch(type) { + case geojson::FeatureType::Point: + return std::make_shared(geojson::PointFeature( + boost::get(geometry), + id, + properties, + bounding_box, + std::vector(), + std::vector(), + {} + )); + case geojson::FeatureType::LineString: + return std::make_shared(geojson::LineStringFeature( + boost::get(geometry), + id, + properties, + bounding_box, + std::vector(), + std::vector(), + {} + )); + case geojson::FeatureType::Polygon: + return std::make_shared(geojson::PolygonFeature( + boost::get(geometry), + id, + properties, + bounding_box, + std::vector(), + std::vector(), + {} + )); + case geojson::FeatureType::MultiPoint: + return std::make_shared(geojson::MultiPointFeature( + boost::get(geometry), + id, + properties, + bounding_box, + std::vector(), + std::vector(), + {} + )); + case geojson::FeatureType::MultiLineString: + return std::make_shared(geojson::MultiLineStringFeature( + boost::get(geometry), + id, + properties, + bounding_box, + std::vector(), + std::vector(), + {} + )); + case geojson::FeatureType::MultiPolygon: + return std::make_shared(geojson::MultiPolygonFeature( + boost::get(geometry), + id, + properties, + bounding_box, + std::vector(), + std::vector(), + {} + )); + default: + return std::make_shared(geojson::CollectionFeature( + std::vector{geometry}, + id, + properties, + bounding_box, + std::vector(), + std::vector(), + {} + )); + } }; inline std::shared_ptr read(const std::string& gpkg_path, const std::string& layer = "", const std::vector& ids = {}) diff --git a/include/geopackage/wkb/reader.hpp b/include/geopackage/wkb/reader.hpp index 7dc51d13c7..cf773a425e 100644 --- a/include/geopackage/wkb/reader.hpp +++ b/include/geopackage/wkb/reader.hpp @@ -8,27 +8,14 @@ #include #include -#include + +#include "JSONGeometry.hpp" namespace geopackage { namespace wkb { using byte_t = uint8_t; using byte_vector = std::vector; -struct wkb_point { double x, y; }; -struct wkb_linestring { std::vector points; }; -struct wkb_polygon { std::vector rings; }; -struct wkb_multipoint { std::vector points; }; -struct wkb_multilinestring { std::vector lines; }; -struct wkb_multipolygon { std::vector polygons; }; -using wkb_geometry = boost::variant< - wkb_point, - wkb_linestring, - wkb_polygon, - wkb_multipoint, - wkb_multilinestring, - wkb_multipolygon ->; /** * @brief @@ -70,59 +57,106 @@ inline void copy_from(byte_vector src, T& index, S& dst, uint8_t order) template inline T read_wkb_internal(const byte_vector& buffer, int& index, uint8_t order); -namespace { -template -inline output_t read_wkb_compound_internal(const byte_vector& buffer, int& index, uint8_t order) -{ - uint32_t count; - copy_from(buffer, index, count, order); - std::vector children(count); - for (auto& child : children) { - child = read_wkb_internal(buffer, index, order); - } - return output_t{children}; -} -} // anonymous namespace - #define READ_WKB_INTERNAL_SIG(output_t) output_t read_wkb_internal(const byte_vector& buffer, int& index, uint8_t order) template<> -inline READ_WKB_INTERNAL_SIG(wkb_point) +inline READ_WKB_INTERNAL_SIG(geojson::coordinate_t) { double x, y; copy_from(buffer, index, x, order); copy_from(buffer, index, y, order); - return wkb_point{x, y}; + return geojson::coordinate_t{x, y}; }; template<> -inline READ_WKB_INTERNAL_SIG(wkb_linestring) +inline READ_WKB_INTERNAL_SIG(geojson::linestring_t) { - return read_wkb_compound_internal(buffer, index, order); + uint32_t count; + copy_from(buffer, index, count, order); + + geojson::linestring_t linestring; + linestring.resize(count); + + for (auto& child : linestring) { + child = read_wkb_internal(buffer, index, order); + } + + return linestring; } template<> -inline READ_WKB_INTERNAL_SIG(wkb_polygon) +inline READ_WKB_INTERNAL_SIG(geojson::polygon_t) { - return read_wkb_compound_internal(buffer, index, order); + uint32_t count; + copy_from(buffer, index, count, order); + + geojson::polygon_t polygon; + + if (count > 1) { + polygon.inners().resize(count - 1); + } + + auto outer = read_wkb_internal(buffer, index, order); + polygon.outer().reserve(outer.size()); + for (auto& p : outer) { + polygon.outer().push_back(p); + } + + for (uint32_t i = 1; i < count; i++) { + auto inner = read_wkb_internal(buffer, index, order); + polygon.inners().at(i).reserve(inner.size()); + for (auto& p : inner) { + polygon.inners().at(i).push_back(p); + } + } + + return polygon; } template<> -inline READ_WKB_INTERNAL_SIG(wkb_multipoint) +inline READ_WKB_INTERNAL_SIG(geojson::multipoint_t) { - return read_wkb_compound_internal(buffer, index, order); + uint32_t count; + copy_from(buffer, index, count, order); + + geojson::multipoint_t mp; + mp.resize(count); + for (auto& point : mp) { + point = read_wkb_internal(buffer, index, order); + } + + return mp; } template<> -inline READ_WKB_INTERNAL_SIG(wkb_multilinestring) +inline READ_WKB_INTERNAL_SIG(geojson::multilinestring_t) { - return read_wkb_compound_internal(buffer, index, order); + uint32_t count; + copy_from(buffer, index, count, order); + + geojson::multilinestring_t ml; + ml.resize(count); + for (auto& line : ml) { + auto l = read_wkb_internal(buffer, index, order); + + } + + return ml; } template<> -inline READ_WKB_INTERNAL_SIG(wkb_multipolygon) +inline READ_WKB_INTERNAL_SIG(geojson::multipolygon_t) { - return read_wkb_compound_internal(buffer, index, order); + uint32_t count; + copy_from(buffer, index, count, order); + + geojson::multipolygon_t mpl; + mpl.resize(count); + for (auto& polygon : mpl) { + polygon = read_wkb_internal(buffer, index, order); + } + + return mpl; } #undef READ_WKB_INTERNAL_SIG @@ -161,7 +195,7 @@ static inline T read_known_wkb(const byte_vector& buffer) * @param buffer vector of bytes * @return wkb::geometry geometry struct containing the parsed WKB values. */ -static inline wkb::wkb_geometry read_wkb(const byte_vector&buffer) +static inline geojson::geometry read_wkb(const byte_vector&buffer) { if (buffer.size() < 5) { throw std::runtime_error("buffer reached end before encountering WKB"); @@ -174,15 +208,14 @@ static inline wkb::wkb_geometry read_wkb(const byte_vector&buffer) uint32_t type; copy_from(buffer, index, type, order); - wkb_geometry g; + geojson::geometry g = geojson::coordinate_t{std::nan("0"), std::nan("0")}; switch(type) { - case 1: g = read_wkb_internal(buffer, index, order); break; - case 2: g = read_wkb_internal(buffer, index, order); break; - case 3: g = read_wkb_internal(buffer, index, order); break; - case 4: g = read_wkb_internal(buffer, index, order); break; - case 5: g = read_wkb_internal(buffer, index, order); break; - case 6: g = read_wkb_internal(buffer, index, order); break; - default: g = wkb_point{std::nan("0"), std::nan("0")}; break; + case 1: g = read_wkb_internal(buffer, index, order); break; + case 2: g = read_wkb_internal(buffer, index, order); break; + case 3: g = read_wkb_internal(buffer, index, order); break; + case 4: g = read_wkb_internal(buffer, index, order); break; + case 5: g = read_wkb_internal(buffer, index, order); break; + case 6: g = read_wkb_internal(buffer, index, order); break; } return g; @@ -191,4 +224,4 @@ static inline wkb::wkb_geometry read_wkb(const byte_vector&buffer) } // namespace wkb } // namespace geopackage -#endif // NGEN_GEOPACKAGE_WKB_READER_H \ No newline at end of file +#endif // NGEN_GEOPACKAGE_WKB_READER_H diff --git a/include/geopackage/wkb/wkt.hpp b/include/geopackage/wkb/wkt.hpp deleted file mode 100644 index 2f96da7bfc..0000000000 --- a/include/geopackage/wkb/wkt.hpp +++ /dev/null @@ -1,93 +0,0 @@ -#ifndef NGEN_GEOPACKAGE_WKB_VISITOR_WKT_H -#define NGEN_GEOPACKAGE_WKB_VISITOR_WKT_H - -#include "reader.hpp" - -namespace geopackage { -namespace wkb { - -class wkt_visitor : public boost::static_visitor -{ - public: - - /** - * @brief Get the WKT form from WKB structs - * - * @tparam T WKB geometry struct type - * @param g geometry object - * @return std::string @param{g} in WKT form - */ - template - std::string operator()(T& g) const - { - return std::string(this->wkt_type(g)) + " " + this->wkt_coords(g); - } - - private: - constexpr const char* wkt_type(const wkb_point&) const { return "POINT"; } - constexpr const char* wkt_type(const wkb_linestring&) const { return "LINESTRING"; } - constexpr const char* wkt_type(const wkb_polygon&) const { return "POLYGON"; } - constexpr const char* wkt_type(const wkb_multipoint&) const { return "MULTIPOINT"; } - constexpr const char* wkt_type(const wkb_multilinestring&) const { return "MULTILINESTRING"; } - constexpr const char* wkt_type(const wkb_multipolygon&) const { return "MULTIPOLYGON"; } - - std::string wkt_coords(const wkb_point& g) const - { - std::ostringstream out; - out.precision(3); - out << std::fixed << g.x << " " << g.y; - return std::move(out).str(); - } - - std::string wkt_coords(const wkb_linestring& g) const - { - return "(" + std::accumulate( - std::next(g.points.begin()), - g.points.end(), - wkt_coords(g.points[0]), - [this](const std::string& a, const wkb_point b) { return a + "," + this->wkt_coords(b); } - ) + ")"; - } - - std::string wkt_coords(const wkb_multipoint& g) const - { - return "(" + std::accumulate( - std::next(g.points.begin()), - g.points.end(), - wkt_coords(g.points[0]), - [this](const std::string& a, const wkb_point b) { return a + "," + this->wkt_coords(b); } - ) + ")"; - } - - std::string wkt_coords(const wkb_polygon& g) const - { - std::string output; - for (const auto& gg : g.rings) { - output += "(" + wkt_coords(gg) + ")"; - } - return output; - } - - std::string wkt_coords(const wkb_multilinestring& g) const - { - std::string output; - for (const auto& gg : g.lines) { - output += "(" + wkt_coords(gg) + ")"; - } - return output; - } - - std::string wkt_coords(const wkb_multipolygon& g) const - { - std::string output; - for (const auto& gg : g.polygons) { - output += "(" + wkt_coords(gg) + ")"; - } - return output; - } -}; - -} // namespace wkb -} // namespace geopackage - -#endif // NGEN_GEOPACKAGE_WKB_VISITOR_WKT_H \ No newline at end of file diff --git a/test/geopackage/WKB_Test.cpp b/test/geopackage/WKB_Test.cpp index 9d3c3d8131..0f495e4a1a 100644 --- a/test/geopackage/WKB_Test.cpp +++ b/test/geopackage/WKB_Test.cpp @@ -1,6 +1,8 @@ +#include "JSONGeometry.hpp" +#include #include + #include -#include using namespace geopackage; @@ -29,7 +31,7 @@ class WKB_Test : public ::testing::Test this->little_endian = { 0x01, - 0x00, 0x00, 0x00, 0x01, + 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x40, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x10, 0x40 }; @@ -51,10 +53,10 @@ class WKB_Test : public ::testing::Test TEST_F(WKB_Test, wkb_read_test) { - const wkb::wkb_geometry geom = wkb::read_wkb(this->wkb); + const geojson::geometry geom = wkb::read_wkb(this->wkb); EXPECT_EQ(geom.which() + 1, 3); // +1 since variant.which() is 0-based - const wkb::wkb_polygon& poly = boost::get(geom); + const geojson::polygon_t& poly = boost::get(geom); const std::vector> expected_coordinates = { {10.689, -25.092}, {34.595, -20.170}, @@ -64,20 +66,19 @@ TEST_F(WKB_Test, wkb_read_test) }; for (int i = 0; i < expected_coordinates.size(); i++) { - EXPECT_NEAR(poly.rings[0].points[i].x, expected_coordinates[i].first, 0.0001); - EXPECT_NEAR(poly.rings[0].points[i].y, expected_coordinates[i].second, 0.0001); + EXPECT_NEAR(poly.outer().at(i).get<0>(), expected_coordinates[i].first, 0.0001); + EXPECT_NEAR(poly.outer().at(i).get<1>(), expected_coordinates[i].second, 0.0001); } - - const wkb::wkt_visitor wkt_v; - const auto wkt = "POLYGON ((10.689 -25.092,34.595 -20.170,38.814 -35.639,13.502 -39.155,10.689 -25.092))"; - EXPECT_EQ(boost::apply_visitor(wkt_v, geom), wkt); - EXPECT_EQ(wkt_v(poly), wkt); } TEST_F(WKB_Test, wkb_endianness_test) { - const auto geom_big = wkb::read_known_wkb(this->big_endian); - const auto geom_little = wkb::read_known_wkb(this->little_endian); - EXPECT_NEAR(geom_big.x, geom_little.x, 0.000001); - EXPECT_NEAR(geom_big.y, geom_little.y, 0.000001); -} \ No newline at end of file + const auto geom_big = wkb::read_known_wkb(this->big_endian); + const auto geom_little = wkb::read_known_wkb(this->little_endian); + EXPECT_NEAR(geom_big.get<0>(), geom_little.get<0>(), 0.000001); + EXPECT_NEAR(geom_big.get<1>(), geom_little.get<1>(), 0.000001); + EXPECT_NEAR(geom_big.get<0>(), 2.0, 0.000001); + EXPECT_NEAR(geom_big.get<1>(), 4.0, 0.000001); + EXPECT_NEAR(geom_big.get<0>(), 2.0, 0.000001); + EXPECT_NEAR(geom_big.get<1>(), 4.0, 0.000001); +}