diff --git a/.gitignore b/.gitignore index c4d726c..55226fe 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,5 @@ tests/rapidjson.png.json __pycache__ benchmarks/*.json benchmarks/*.pbf* +# https://github.com/cubao/nano-fmm/blob/master/data/Makefile +data/suzhoubeizhan.json diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 432b377..b906381 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -45,6 +45,7 @@ repos: rev: 5.11.5 hooks: - id: isort + args: [--profile, black] # Upgrade older Python syntax - repo: https://github.com/asottile/pyupgrade @@ -61,7 +62,7 @@ repos: exclude: ^(docs|Makefile|benchmarks/Makefile) - repo: https://github.com/PyCQA/flake8 - rev: 3.9.2 + rev: 6.1.0 hooks: - id: flake8 additional_dependencies: [flake8-bugbear] diff --git a/.vscode/settings.json b/.vscode/settings.json index 0c4997d..5497d00 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -93,7 +93,14 @@ "regex": "cpp", "future": "cpp", "__bits": "cpp", - "bit": "cpp" + "bit": "cpp", + "charconv": "cpp", + "concepts": "cpp", + "memory_resource": "cpp", + "ranges": "cpp", + "shared_mutex": "cpp", + "span": "cpp", + "stop_token": "cpp" }, "workbench.colorCustomizations": { "activityBar.background": "#143328", diff --git a/data/suzhoubeizhan.pbf b/data/suzhoubeizhan.pbf new file mode 100644 index 0000000..b38185b Binary files /dev/null and b/data/suzhoubeizhan.pbf differ diff --git a/data/suzhoubeizhan_crossover.json b/data/suzhoubeizhan_crossover.json new file mode 100644 index 0000000..f0853b7 --- /dev/null +++ b/data/suzhoubeizhan_crossover.json @@ -0,0 +1,79 @@ +{ + "type": "Feature", + "properties": {}, + "geometry": { + "coordinates": [ + [ + [ + 120.65771383441779, + 31.406929844034934 + ], + [ + 120.65685254868885, + 31.405354623300482 + ], + [ + 120.65945178597582, + 31.403989410609327 + ], + [ + 120.66031307170232, + 31.402597923391937 + ], + [ + 120.66214330387294, + 31.401980936609206 + ], + [ + 120.66508090340801, + 31.403634977144364 + ], + [ + 120.66558844678389, + 31.4051183379097 + ], + [ + 120.66348137277282, + 31.406772323151472 + ], + [ + 120.66488096208019, + 31.407730571105034 + ], + [ + 120.66429651819169, + 31.408557544349165 + ], + [ + 120.66354289318042, + 31.409581405887053 + ], + [ + 120.66086675538332, + 31.409030097216544 + ], + [ + 120.66140505896414, + 31.407336771760555 + ], + [ + 120.66147901958857, + 31.407062414261205 + ], + [ + 120.66108904538675, + 31.406620542422047 + ], + [ + 120.65915262176173, + 31.40651150879127 + ], + [ + 120.65771383441779, + 31.406929844034934 + ] + ] + ], + "type": "Polygon" + } +} diff --git a/docs/about/release-notes.md b/docs/about/release-notes.md index 7d6e426..ba82d3e 100644 --- a/docs/about/release-notes.md +++ b/docs/about/release-notes.md @@ -10,6 +10,10 @@ To upgrade `pybind11-geobuf` to the latest version, use pip: pip install -U pybind11-geobuf ``` +## Version 0.1.7 (2023-11-11) + +* Integrate PackedRTree (rbush) + ## Version 0.1.6 (2023-07-02) * Crop geojson features by polygon (alpha release) diff --git a/headers b/headers index 64c2294..0c34001 160000 --- a/headers +++ b/headers @@ -1 +1 @@ -Subproject commit 64c2294d63e9e4d61f59e39fb31352cab7fc71c5 +Subproject commit 0c34001ac35ca9f89088ac176d1f1a739f69d124 diff --git a/setup.py b/setup.py index a4a39e9..a9c3d58 100644 --- a/setup.py +++ b/setup.py @@ -127,7 +127,7 @@ def build_extension(self, ext): # logic and declaration, and simpler if you include description/version in a file. setup( name="pybind11_geobuf", - version="0.1.6", + version="0.1.7", author="tzx", author_email="dvorak4tzx@gmail.com", url="https://geobuf-cpp.readthedocs.io", diff --git a/src/geobuf/geojson_cropping.hpp b/src/geobuf/geojson_cropping.hpp index 97524e8..34d0f81 100644 --- a/src/geobuf/geojson_cropping.hpp +++ b/src/geobuf/geojson_cropping.hpp @@ -20,6 +20,7 @@ namespace cubao { using RowVectors = Eigen::Matrix; +using RowVectorsNx2 = Eigen::Matrix; using BboxType = mapbox::geometry::box; inline BboxType row_vectors_to_bbox(const RowVectors &coords) @@ -30,7 +31,7 @@ inline BboxType row_vectors_to_bbox(const RowVectors &coords) auto bbox = BboxType({max_t, max_t, max_t}, {min_t, min_t, min_t}); auto &min = bbox.min; auto &max = bbox.max; - for (int i = 1, N = coords.rows(); i < N; ++i) { + for (int i = 0, N = coords.rows(); i < N; ++i) { double x = coords(i, 0); double y = coords(i, 1); double z = coords(i, 2); @@ -50,6 +51,30 @@ inline BboxType row_vectors_to_bbox(const RowVectors &coords) return bbox; } +inline BboxType +row_vectors_to_bbox(const Eigen::Ref &coords) +{ + using limits = std::numeric_limits; + double min_t = limits::has_infinity ? -limits::infinity() : limits::min(); + double max_t = limits::has_infinity ? limits::infinity() : limits::max(); + auto bbox = BboxType({max_t, max_t, 0.0}, {min_t, min_t, 0.0}); + auto &min = bbox.min; + auto &max = bbox.max; + for (int i = 0, N = coords.rows(); i < N; ++i) { + double x = coords(i, 0); + double y = coords(i, 1); + if (min.x > x) + min.x = x; + if (min.y > y) + min.y = y; + if (max.x < x) + max.x = x; + if (max.y < y) + max.y = y; + } + return bbox; +} + inline RowVectors bbox2row_vectors(const BboxType &bbox) { auto coords = RowVectors(5, 3); @@ -176,7 +201,7 @@ inline int geojson_cropping(const mapbox::geojson::feature &feature, double len2 = std::get<5>(k2) - std::get<2>(k2); return len1 < len2; }); - keys = {*itr}; + keys = {*itr}; // pick longest } for (auto &key : keys) { auto &coords = segs[key]; diff --git a/src/geobuf/planet.hpp b/src/geobuf/planet.hpp new file mode 100644 index 0000000..277b264 --- /dev/null +++ b/src/geobuf/planet.hpp @@ -0,0 +1,186 @@ +#ifndef CUBAO_PLANET_HPP +#define CUBAO_PLANET_HPP + +// https://github.com/microsoft/vscode-cpptools/issues/9692 +#if __INTELLISENSE__ +#undef __ARM_NEON +#undef __ARM_NEON__ +#endif + +#include "geojson_cropping.hpp" +#include "packedrtree.hpp" + +namespace cubao +{ +struct Planet +{ + using FeatureCollection = mapbox::geojson::feature_collection; + Planet() = default; + Planet(const FeatureCollection &features) { this->features(features); } + + const FeatureCollection &features() const { return features_; } + Planet &features(const FeatureCollection &features) + { + features_ = features; + rtree_.reset(); + return *this; + } + + void build() const { this->rtree(); } + + // TODO, query by style expression + Eigen::VectorXi query(const Eigen::Vector2d &min, + const Eigen::Vector2d &max) const + { + auto &tree = this->rtree(); + auto hits = tree.search(min[0], min[1], max[0], max[1]); + Eigen::VectorXi index(hits.size()); + for (int i = 0, N = hits.size(); i < N; ++i) { + index[i] = hits[i].offset; + } + return index; + } + FeatureCollection copy(const Eigen::VectorXi &index) const + { + auto fc = FeatureCollection(); + fc.reserve(index.size()); + for (int i = 0, N = index.size(); i < N; ++i) { + fc.push_back(features_[index[i]]); + } + return fc; + } + + FeatureCollection crop(const Eigen::Ref &polygon, + const std::string &clipping_mode = "longest", + bool strip_properties = false, + bool is_wgs84 = true) const + { + auto bbox = row_vectors_to_bbox(polygon); + auto hits = + this->query({bbox.min.x, bbox.min.y}, {bbox.max.x, bbox.max.y}); + auto fc = FeatureCollection(); + fc.reserve(hits.size()); + for (auto idx : hits) { + auto &feature = features_[idx]; + if (!feature.geometry.is() || + clipping_mode == "whole") { + // only check centroid + auto centroid = geometry_to_centroid(feature.geometry); + auto mask = point_in_polygon( + Eigen::Map(¢roid.x, 1, 2), + polygon); + if (mask[0]) { + fc.emplace_back( + feature.geometry, + strip_properties + ? mapbox::feature::property_map{{"index", idx}} + : feature.properties, + feature.id); + } + continue; + } + auto &line_string = + feature.geometry.get(); + auto polyline = Eigen::Map(&line_string[0].x, + line_string.size(), 3); + auto segs = polyline_in_polygon(polyline, polygon, is_wgs84); + if (segs.empty()) { + continue; + } + if (clipping_mode == "first") { + auto &coords = segs.begin()->second; + auto geom = mapbox::geojson::line_string(); + geom.resize(coords.rows()); + as_row_vectors(geom) = coords; + fc.emplace_back( + geom, + strip_properties + ? mapbox::feature::property_map{{"index", idx}} + : feature.properties, + feature.id); + continue; + } + // longest or all + std::vector keys; + keys.reserve(segs.size()); + for (auto &pair : segs) { + keys.push_back(pair.first); + } + if (clipping_mode == "longest") { // else assume all + auto itr = std::max_element( + keys.begin(), keys.end(), + [](const auto &k1, const auto &k2) { + double len1 = std::get<5>(k1) - std::get<2>(k1); + double len2 = std::get<5>(k2) - std::get<2>(k2); + return len1 < len2; + }); + keys = {*itr}; // pick longest + } + for (auto &key : keys) { + auto &coords = segs[key]; + auto geom = mapbox::geojson::line_string(); + geom.resize(coords.rows()); + as_row_vectors(geom) = coords; + fc.emplace_back( + geom, + strip_properties + ? mapbox::feature::property_map{{"index", idx}} + : feature.properties, + feature.id); + } + } + return fc; + } + + private: + FeatureCollection features_; + + mutable std::optional rtree_; + + template + static FlatGeobuf::NodeItem envelope_2d(G const &geometry, uint64_t index) + { + // mapbox/geometry/envelope.hpp + using limits = std::numeric_limits; + constexpr double min_t = + limits::has_infinity ? -limits::infinity() : limits::min(); + constexpr double max_t = + limits::has_infinity ? limits::infinity() : limits::max(); + double min_x = max_t; + double min_y = max_t; + double max_x = min_t; + double max_y = min_t; + mapbox::geometry::for_each_point( + geometry, [&](mapbox::geojson::point const &point) { + if (min_x > point.x) + min_x = point.x; + if (min_y > point.y) + min_y = point.y; + if (max_x < point.x) + max_x = point.x; + if (max_y < point.y) + max_y = point.y; + }); + return {min_x, min_y, max_x, max_y, index}; + } + + FlatGeobuf::PackedRTree &rtree() const + { + if (rtree_) { + return *rtree_; + } + auto nodes = std::vector{}; + nodes.reserve(features_.size()); + uint64_t index{0}; + for (auto &feature : features_) { + nodes.emplace_back(envelope_2d(feature.geometry, index++)); + } + auto extent = calcExtent(nodes); + hilbertSort(nodes, extent); + rtree_ = FlatGeobuf::PackedRTree(nodes, extent); + return *rtree_; + } +}; +} // namespace cubao + +#endif diff --git a/src/main.cpp b/src/main.cpp index c26c4b4..e148ae6 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -11,6 +11,7 @@ #include #include "geobuf/geobuf.hpp" +#include "geobuf/planet.hpp" #include "geobuf/pybind11_helpers.hpp" #include @@ -256,6 +257,25 @@ PYBIND11_MODULE(_pybind11_geobuf, m) auto geojson = m.def_submodule("geojson"); cubao::bind_geojson(geojson); + using Planet = cubao::Planet; + py::class_(m, "Planet", py::module_local()) + .def(py::init<>()) + .def(py::init()) + .def("features", py::overload_cast<>(&Planet::features, py::const_), + rvp::reference_internal) + .def("features", + py::overload_cast( + &Planet::features)) + .def("build", &Planet::build) + .def("query", &Planet::query, "min"_a, "max"_a) + .def("copy", &Planet::copy) + .def("crop", &Planet::crop, "polygon"_a, py::kw_only(), + "clipping_mode"_a = "longest", // + "strip_properties"_a = false, // + "is_wgs84"_a = true) + // + ; + cubao::bind_rapidjson(m); auto tf = m.def_submodule("tf"); diff --git a/tests/test_geobuf.py b/tests/test_geobuf.py index 8990ddd..a231406 100644 --- a/tests/test_geobuf.py +++ b/tests/test_geobuf.py @@ -15,6 +15,7 @@ from pybind11_geobuf import ( # noqa Decoder, Encoder, + Planet, geojson, normalize_json, pbf_decode, @@ -1400,7 +1401,7 @@ def test_geojson_feature(): assert props["list"].is_array() for x in list(props["list"].as_array()): assert isinstance(x, geojson.value) - assert type(x) == geojson.value + assert type(x) is geojson.value with pytest.raises(RuntimeError) as excinfo: props["list"].as_object() assert "in get()" in repr(excinfo) @@ -1412,7 +1413,7 @@ def test_geojson_feature(): for k, v in props["dict"].as_object().items(): assert isinstance(k, str) assert isinstance(v, geojson.value) - assert type(x) == geojson.value + assert type(x) is geojson.value with pytest.raises(RuntimeError) as excinfo: props["dict"].as_array() assert "in get()" in repr(excinfo) @@ -1755,6 +1756,27 @@ def test_rapidjson_normalize_non_geojson(): ) +def test_query(): + path = f"{__pwd}/../data/suzhoubeizhan.pbf" + fc = geojson.FeatureCollection().load(path) + planet = Planet(fc) + hits = planet.query([120.64094, 31.41515], [120.64137, 31.41534]) + assert len(hits) == 4 + + path = f"{__pwd}/../data/suzhoubeizhan_crossover.json" + polygon = geojson.Feature().load(path).to_numpy() + + cropped1 = planet.crop(polygon[:, :2]) + cropped2 = planet.crop(polygon[:, :2], strip_properties=True) + # cropped1.dump('cropped1.json', indent=True) + # cropped2.dump('cropped2.json', indent=True) + assert len(cropped1) == len(cropped2) == 54 + assert len(list(cropped1[0].properties().keys())) == 6 + assert list(cropped2[0].properties().keys()) == ["index"] + assert cropped2[0].properties()["index"]() == 438 + assert fc[438] == cropped1[0] + + if __name__ == "__main__": np.set_printoptions(suppress=True) pwd = os.path.abspath(os.path.dirname(__file__))