diff --git a/Makefile b/Makefile index 1356dad..7f48cd5 100644 --- a/Makefile +++ b/Makefile @@ -85,7 +85,7 @@ docs_serve: mkdocs serve -a 0.0.0.0:8088 DOCKER_TAG_WINDOWS ?= ghcr.io/cubao/build-env-windows-x64:v0.0.1 -DOCKER_TAG_LINUX ?= ghcr.io/cubao/build-env-manylinux2014-x64:v0.0.1 +DOCKER_TAG_LINUX ?= ghcr.io/cubao/build-env-manylinux2014-x64:v0.0.4 DOCKER_TAG_MACOS ?= ghcr.io/cubao/build-env-macos-arm64:v0.0.1 test_in_win: @@ -95,6 +95,15 @@ test_in_mac: test_in_linux: docker run --rm -w `pwd` -v `pwd`:`pwd` -v `pwd`/build/linux:`pwd`/build -it $(DOCKER_TAG_LINUX) bash +DEV_CONTAINER_NAME ?= $(USER)_$(subst /,_,$(PROJECT_NAME)____$(PROJECT_SOURCE_DIR)) +DEV_CONTAINER_IMAG ?= $(DOCKER_TAG_LINUX) +test_in_dev_container: + docker ps | grep $(DEV_CONTAINER_NAME) \ + && docker exec -it $(DEV_CONTAINER_NAME) bash \ + || docker run --rm --name $(DEV_CONTAINER_NAME) \ + --network host --security-opt seccomp=unconfined \ + -v `pwd`:`pwd` -w `pwd` -it $(DEV_CONTAINER_IMAG) bash + PYTHON ?= python3 python_install: $(PYTHON) setup.py install --force @@ -155,7 +164,7 @@ python_build_py310: PYTHON=python conda run --no-capture-output -n py310 make python_build python_build_all: python_build_py36 python_build_py37 python_build_py38 python_build_py39 python_build_py310 python_build_all_in_linux: - docker run --rm -w `pwd` -v `pwd`:`pwd` -v `pwd`/build/win:`pwd`/build -it $(DOCKER_TAG_LINUX) make python_build_all + docker run --rm -w `pwd` -v `pwd`:`pwd` -v `pwd`/build/linux:`pwd`/build -it $(DOCKER_TAG_LINUX) make python_build_all make repair_wheels && rm -rf dist/*.whl && mv wheelhouse/*.whl dist && rm -rf wheelhouse python_build_all_in_macos: python_build_py38 python_build_py39 python_build_py310 python_build_all_in_windows: python_build_all diff --git a/docs/about/release-notes.md b/docs/about/release-notes.md index b7dad5f..7d6e426 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.6 (2023-07-02) + +* Crop geojson features by polygon (alpha release) + ## Version 0.1.5 (2023-06-02) * Add `round_non_geojson` to `normalize_json` diff --git a/headers b/headers index 131c1ca..64c2294 160000 --- a/headers +++ b/headers @@ -1 +1 @@ -Subproject commit 131c1ca724002532479e966ee84670610b932456 +Subproject commit 64c2294d63e9e4d61f59e39fb31352cab7fc71c5 diff --git a/mkdocs.yml b/mkdocs.yml index 92f4742..cbad8ff 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -4,7 +4,7 @@ site_description: A compact binary encoding for geographic data, with python bin site_author: district10 repo_url: https://github.com/cubao/pybind11-geobuf -edit_uri: blob/master/docs/ +edit_uri: blob/dev/docs/ theme: readthedocs nav: diff --git a/pybind11_geobuf/crop.py b/pybind11_geobuf/crop.py new file mode 100755 index 0000000..3b33f65 --- /dev/null +++ b/pybind11_geobuf/crop.py @@ -0,0 +1,135 @@ +import os +from typing import List, Optional, Tuple, Union # noqa + +import numpy as np +from _pybind11_geobuf import geojson, tf +from loguru import logger + + +def bbox2polygon(bbox: np.ndarray): + lon0, lat0, lon1, lat1 = bbox + return np.array( + [ + [lon0, lat0, 0.0], + [lon1, lat0, 0.0], + [lon1, lat1, 0.0], + [lon0, lat1, 0.0], + [lon0, lat0, 0.0], + ] + ) + + +def crop_by_feature_id( + input_path: str, + output_path: str, + *, + feature_id: str, + buffer: Union[float, Tuple[float, float]] = 100.0, + clipping_mode: str = "longest", + max_z_offset: float = None, +) -> bool: + if not feature_id: + logger.info( + f"invalid feature id: {feature_id} (type: {type(feature_id)})" + ) # noqa + return False + g = geojson.GeoJSON().load(input_path) + if not g.is_feature_collection(): + logger.warning(f"{input_path} is not valid GeoJSON FeatureCollection") + return False + + fc = g.as_feature_collection() + bbox = None + height = None + for f in fc: + props = f.properties() + if "id" not in props: + continue + fid = props["id"]() + if fid == feature_id: + bbox = f.bbox() + if max_z_offset is not None: + height = f.bbox(with_z=True)[2::3].mean() + if bbox is None: + logger.error(f"not any feature matched by id: {feature_id}") + return False + + dlon, dlat = 1.0 / tf.cheap_ruler_k(bbox[1::2].mean())[:2] + if isinstance(buffer, (int, float, np.generic)): + dlon *= buffer + dlat *= buffer + else: + dlon *= buffer[0] + dlat *= buffer[1] + bbox += [-dlon, -dlat, dlon, dlat] + logger.info(f"bbox: {bbox}") + + polygon = bbox2polygon(bbox) + if height is not None: + polygon[:, 2] = height + logger.info(f"polygon:\n{polygon}") + + logger.info(f"writing to {output_path} ...") + os.makedirs(os.path.dirname(os.path.abspath(output_path)), exist_ok=True) + cropped = g.crop( + polygon, + clipping_mode=clipping_mode, + max_z_offset=max_z_offset, + ) + return cropped.to_rapidjson().sort_keys().dump(output_path, indent=True) + + +def crop_by_grid( + input_path: str, + output_dir: str, + *, + anchor_lla: Union[str, List[float]] = None, + grid_size: Union[float, Tuple[float, float]] = 1000.0, +): + os.makedirs(os.path.abspath(output_dir), exist_ok=True) + + +def crop_by_center( + input_path: str, + output_dir: str, + *, + anchor_lla: Union[str, List[float]] = None, + size: Union[float, Tuple[float, float]] = 1000.0, +): + os.makedirs(os.path.abspath(output_dir), exist_ok=True) + + +def crop_by_bbox( + input_path: str, + output_path: str, + *, + bbox: Union[str, List[float]], + z_center: float = None, + z_max_offset: float = None, +): + logger.info(f"wrote to {output_path}") + + +def crop_by_polygon( + input_path: str, + output_path: str, + *, + polygon: Union[str, np.ndarray], + z_max_offset: float = None, +): + pass + + +if __name__ == "__main__": + import fire + + fire.core.Display = lambda lines, out: print(*lines, file=out) + fire.Fire( + { + "by_feature_id": crop_by_feature_id, + "by_grid": crop_by_grid, + "by_center": crop_by_center, + "by_bbox": crop_by_bbox, + "by_polygon": crop_by_polygon, + } + ) diff --git a/setup.py b/setup.py index 0a1cc3d..a4a39e9 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.5", + version="0.1.6", 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 new file mode 100644 index 0000000..97524e8 --- /dev/null +++ b/src/geobuf/geojson_cropping.hpp @@ -0,0 +1,232 @@ +#ifndef CUBAO_GEOJSON_CROPPING_HPP +#define CUBAO_GEOJSON_CROPPING_HPP + +// https://github.com/microsoft/vscode-cpptools/issues/9692 +#if __INTELLISENSE__ +#undef __ARM_NEON +#undef __ARM_NEON__ +#endif + +#include "cubao/point_in_polygon.hpp" +#include "cubao/polyline_in_polygon.hpp" +#include "geojson_helpers.hpp" +#include +#include +#include +#include +#include +#include + +namespace cubao +{ +using RowVectors = Eigen::Matrix; +using BboxType = mapbox::geometry::box; + +inline BboxType row_vectors_to_bbox(const RowVectors &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, 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) { + double x = coords(i, 0); + double y = coords(i, 1); + double z = coords(i, 2); + if (min.x > x) + min.x = x; + if (min.y > y) + min.y = y; + if (min.z > z) + min.z = z; + if (max.x < x) + max.x = x; + if (max.y < y) + max.y = y; + if (max.z < z) + max.z = z; + } + return bbox; +} + +inline RowVectors bbox2row_vectors(const BboxType &bbox) +{ + auto coords = RowVectors(5, 3); + coords << bbox.min.x, bbox.min.y, 0.0, // + bbox.max.x, bbox.min.y, 0.0, // + bbox.max.x, bbox.max.y, 0.0, // + bbox.min.x, bbox.max.y, 0.0, // + bbox.min.x, bbox.min.y, 0.0; + return coords; +} + +inline RowVectors bbox2row_vectors(const Eigen::Vector4d &bbox) +{ + return bbox2row_vectors(BboxType({bbox[0], bbox[1]}, {bbox[2], bbox[3]})); +} + +inline mapbox::geojson::point +geometry_to_centroid(const mapbox::geojson::geometry &geom) +{ + auto centroid = mapbox::geojson::point(0, 0, 0); + int N = 0; + mapbox::geometry::for_each_point(geom, [&](auto &point) { + centroid.x += point.x; + centroid.y += point.y; + centroid.z += point.z; + ++N; + }); + centroid.x /= N; + centroid.y /= N; + centroid.z /= N; + return centroid; +} + +inline bool bbox_overlap(const BboxType &bbox1, const BboxType &bbox2, + bool check_z = false) +{ + if (check_z && (bbox1.max.z < bbox2.min.z || bbox2.min.z > bbox2.max.z)) { + return false; + } + if (bbox1.max.x < bbox2.min.x || bbox1.min.x > bbox2.max.x) { + return false; + } + if (bbox1.max.y < bbox2.min.y || bbox1.min.y > bbox2.max.y) { + return false; + } + return true; +} + +/* +there is a difference between +- "cropping" (more like in a bigger picture) and +- "clipping" (more like focused onto something specific, like in tunnel +vision). + +clipping_mode + - longest + - first + - all + - whole +*/ + +inline int geojson_cropping(const mapbox::geojson::feature &feature, + mapbox::geojson::feature_collection &output, + const RowVectors &polygon, + std::optional bbox, + const std::string &clipping_mode = "longest", + const std::optional max_z_offset = {}) +{ + if (!bbox) { + bbox = row_vectors_to_bbox(polygon); + } + if (max_z_offset) { + bbox->min.z -= *max_z_offset; + bbox->max.z += *max_z_offset; + } + + if (!bbox_overlap(mapbox::geometry::envelope(feature.geometry), // + *bbox, // + (bool)max_z_offset) // check_z? + ) { + return 0; + } + 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.leftCols<2>()); + if (mask[0]) { + output.push_back(feature); + return 1; + } else { + return 0; + } + } + 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.leftCols<2>()); + if (segs.empty()) { + return 0; + } + if (clipping_mode == "first") { + auto &coords = segs.begin()->second; + auto geom = mapbox::geojson::line_string(); + geom.resize(coords.rows()); + as_row_vectors(geom) = coords; + auto f = feature; + f.geometry = geom; + output.push_back(std::move(f)); + return 1; + } + // 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}; + } + for (auto &key : keys) { + auto &coords = segs[key]; + auto geom = mapbox::geojson::line_string(); + geom.resize(coords.rows()); + as_row_vectors(geom) = coords; + auto f = feature; + f.geometry = geom; + output.push_back(std::move(f)); + } + return keys.size(); +} + +inline mapbox::geojson::feature_collection +geojson_cropping(const mapbox::geojson::feature_collection &features, + const RowVectors &polygon, + const std::string &clipping_mode = "longest", + const std::optional max_z_offset = {}) +{ + auto bbox = row_vectors_to_bbox(polygon); + bbox.min.z = bbox.max.z = (bbox.min.z + bbox.max.z) / 2.0; + mapbox::geojson::feature_collection output; + for (auto &f : features) { + geojson_cropping(f, output, polygon, bbox, clipping_mode, max_z_offset); + } + return output; +} + +inline mapbox::geojson::feature_collection +geojson_cropping(const mapbox::geojson::geojson &geojson, + const RowVectors &polygon, + const std::string &clipping_mode = "longest", + const std::optional max_z_offset = {}) +{ + return geojson.match( + [&](const mapbox::geojson::geometry &g) { + return geojson_cropping( + mapbox::geojson::feature_collection{ + mapbox::geojson::feature_collection{g}}, + polygon, clipping_mode, max_z_offset); + }, + [&](const mapbox::geojson::feature &f) { + return geojson_cropping(mapbox::geojson::feature_collection{f}, + polygon, clipping_mode, max_z_offset); + }, + [&](const mapbox::geojson::feature_collection &fc) { + return geojson_cropping(fc, polygon, clipping_mode, max_z_offset); + }); +} + +} // namespace cubao + +#endif diff --git a/src/main.cpp b/src/main.cpp index 5c08a02..c26c4b4 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -26,6 +26,7 @@ namespace cubao { void bind_geojson(py::module &m); void bind_rapidjson(py::module &m); +void bind_crs_transform(py::module &m); } // namespace cubao PYBIND11_MODULE(_pybind11_geobuf, m) @@ -257,9 +258,19 @@ PYBIND11_MODULE(_pybind11_geobuf, m) cubao::bind_rapidjson(m); + auto tf = m.def_submodule("tf"); + cubao::bind_crs_transform(tf); + #ifdef VERSION_INFO m.attr("__version__") = MACRO_STRINGIFY(VERSION_INFO); #else m.attr("__version__") = "dev"; #endif } + +#define CUBAO_STATIC_LIBRARY +#ifndef CUBAO_ARGV_DEFAULT_NONE +#define CUBAO_ARGV_DEFAULT_NONE(argv) py::arg_v(#argv, std::nullopt, "None") +#endif + +#include "cubao/pybind11_crs_transform.hpp" diff --git a/src/pybind11_geojson.cpp b/src/pybind11_geojson.cpp index 02ed2db..09a3681 100644 --- a/src/pybind11_geojson.cpp +++ b/src/pybind11_geojson.cpp @@ -6,6 +6,7 @@ #include #include +#include "geobuf/geojson_cropping.hpp" #include "geobuf/geojson_helpers.hpp" #include "geobuf/pybind11_helpers.hpp" #include "geobuf/rapidjson_helpers.hpp" @@ -36,6 +37,24 @@ using rvp = py::return_value_policy; using PropertyMap = mapbox::geojson::value::object_type; +template Eigen::VectorXd geom2bbox(const T &t, bool with_z) +{ + auto bbox = mapbox::geometry::envelope(t); + if (with_z) { + return (Eigen::VectorXd(6) << // + bbox.min.x, + bbox.min.y, bbox.min.z, // + bbox.max.x, bbox.max.y, bbox.max.z) + .finished(); + } else { + return (Eigen::VectorXd(4) << // + bbox.min.x, + bbox.min.y, // + bbox.max.x, bbox.max.y) + .finished(); + } +} + inline bool endswith(const std::string &text, const std::string &suffix) { if (text.length() < suffix.length()) { @@ -146,6 +165,21 @@ void bind_geojson(py::module &geojson) "only_xy"_a = false, // "round_z"_a = std::nullopt) // + .def( + "crop", + [](mapbox::geojson::geojson &self, const RowVectors &polygon, + const std::string &clipping_mode, + std::optional max_z_offset) + -> mapbox::geojson::feature_collection { + return cubao::geojson_cropping(self, // + polygon, // + clipping_mode, // + max_z_offset); + }, + "polygon"_a, py::kw_only(), // + "clipping_mode"_a = "longest", // + "max_z_offset"_a = std::nullopt) + // .def( "load", [](mapbox::geojson::geojson &self, @@ -549,6 +583,9 @@ void bind_geojson(py::module &geojson) "sort_keys"_a = false, // "precision"_a = 8, // "only_xy"_a = false) + .def("bbox", [](const mapbox::geojson::geometry &self, bool with_z) -> Eigen::VectorXd { + return geom2bbox(self, with_z); + }, py::kw_only(), "with_z"_a = false) .def(py::self == py::self) // .def(py::self != py::self) // @@ -650,6 +687,11 @@ void bind_geojson(py::module &geojson) mapbox::geojson::geometry{self}, allocator); return json; }) + .def( + "bbox", + [](const mapbox::geojson::point &self, bool with_z) + -> Eigen::VectorXd { return geom2bbox(self, with_z); }, + py::kw_only(), "with_z"_a = false) .def(py::self == py::self) // .def(py::self != py::self) // .def( @@ -793,12 +835,18 @@ void bind_geojson(py::module &geojson) return self; \ }, \ rvp::reference_internal) \ - .def("to_rapidjson", [](const mapbox::geojson::geom_type &self) { \ - RapidjsonAllocator allocator; \ - auto json = mapbox::geojson::convert( \ - mapbox::geojson::geometry{self}, allocator); \ - return json; \ - }) + .def("to_rapidjson", \ + [](const mapbox::geojson::geom_type &self) { \ + RapidjsonAllocator allocator; \ + auto json = mapbox::geojson::convert( \ + mapbox::geojson::geometry{self}, allocator); \ + return json; \ + }) \ + .def( \ + "bbox", \ + [](const mapbox::geojson::geom_type &self, bool with_z) \ + -> Eigen::VectorXd { return geom2bbox(self, with_z); }, \ + py::kw_only(), "with_z"_a = false) py::class_>(geojson, "MultiPoint", @@ -956,7 +1004,12 @@ void bind_geojson(py::module &geojson) return self; \ }, \ py::kw_only(), "lon"_a = 8, "lat"_a = 8, "alt"_a = 3, \ - rvp::reference_internal) GEOMETRY_DEDUPLICATE_XYZ(geom_type) + rvp::reference_internal) GEOMETRY_DEDUPLICATE_XYZ(geom_type) \ + .def( \ + "bbox", \ + [](const mapbox::geojson::geom_type &self, bool with_z) \ + -> Eigen::VectorXd { return geom2bbox(self, with_z); }, \ + py::kw_only(), "with_z"_a = false) py::class_( @@ -1151,6 +1204,11 @@ void bind_geojson(py::module &geojson) mapbox::geojson::geometry{self}, allocator); return json; }) + .def( + "bbox", + [](const mapbox::geojson::multi_polygon &self, bool with_z) + -> Eigen::VectorXd { return geom2bbox(self, with_z); }, + py::kw_only(), "with_z"_a = false) .def(py::self == py::self) // .def(py::self != py::self) // // @@ -1720,6 +1778,11 @@ void bind_geojson(py::module &geojson) return self; }, rvp::reference_internal) + .def( + "bbox", + [](const mapbox::geojson::feature &self, bool with_z) + -> Eigen::VectorXd { return geom2bbox(self.geometry, with_z); }, + py::kw_only(), "with_z"_a = false) .def(py::self == py::self) // .def(py::self != py::self) // //