Integrate PackedRTree (rbush)

.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

.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]

.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",

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"
+    }
+}

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)

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",

src/geobuf/geojson_cropping.hpp

@@ -20,6 +20,7 @@
 namespace cubao
 {
 using RowVectors = Eigen::Matrix<double, Eigen::Dynamic, 3, Eigen::RowMajor>;
+using RowVectorsNx2 = Eigen::Matrix<double, Eigen::Dynamic, 2, Eigen::RowMajor>;
 using BboxType = mapbox::geometry::box<double>;
 
 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<const RowVectorsNx2> &coords)
+{
+    using limits = std::numeric_limits<double>;
+    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];

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<const RowVectorsNx2> &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<mapbox::geojson::line_string>() ||
+                clipping_mode == "whole") {
+                // only check centroid
+                auto centroid = geometry_to_centroid(feature.geometry);
+                auto mask = point_in_polygon(
+                    Eigen::Map<const RowVectorsNx2>(&centroid.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<mapbox::geojson::line_string>();
+            auto polyline = Eigen::Map<const RowVectors>(&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<PolylineChunks::key_type> 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<FlatGeobuf::PackedRTree> rtree_;
+
+    template <typename G>
+    static FlatGeobuf::NodeItem envelope_2d(G const &geometry, uint64_t index)
+    {
+        // mapbox/geometry/envelope.hpp
+        using limits = std::numeric_limits<double>;
+        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<FlatGeobuf::NodeItem>{};
+        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