Skip to content

Commit

Permalink
Integrate PackedRTree (rbush) (#24)
Browse files Browse the repository at this point in the history
* update headers

* update headers

* add planet

* add tree

* crop features

* query

* test query

* update to 0.1.7

* fix index

* bump flake: PyCQA/flake8#1885

* fix lint

---------

Co-authored-by: TANG ZHIXIONG <zhixiong.tang@momenta.ai>
  • Loading branch information
district10 and zhixiong-tang authored Nov 10, 2023
1 parent 9205a44 commit dd7cfee
Show file tree
Hide file tree
Showing 12 changed files with 354 additions and 8 deletions.
2 changes: 2 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -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
3 changes: 2 additions & 1 deletion .pre-commit-config.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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]
Expand Down
9 changes: 8 additions & 1 deletion .vscode/settings.json
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down
Binary file added data/suzhoubeizhan.pbf
Binary file not shown.
79 changes: 79 additions & 0 deletions data/suzhoubeizhan_crossover.json
Original file line number Diff line number Diff line change
@@ -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"
}
}
4 changes: 4 additions & 0 deletions docs/about/release-notes.md
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion headers
Submodule headers updated 774 files
2 changes: 1 addition & 1 deletion setup.py
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down
29 changes: 27 additions & 2 deletions src/geobuf/geojson_cropping.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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);
Expand All @@ -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);
Expand Down Expand Up @@ -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];
Expand Down
186 changes: 186 additions & 0 deletions src/geobuf/planet.hpp
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit dd7cfee

Please sign in to comment.