diff --git a/duckdb b/duckdb index 2ed9bf88..ff0f9595 160000 --- a/duckdb +++ b/duckdb @@ -1 +1 @@ -Subproject commit 2ed9bf887f61a0ac226ab8c8f1164601d985d607 +Subproject commit ff0f95954bdb4c7515481ac2b473261407bad18b diff --git a/src/sgl/sgl.cpp b/src/sgl/sgl.cpp index 27793af2..e614fec8 100644 --- a/src/sgl/sgl.cpp +++ b/src/sgl/sgl.cpp @@ -1,12 +1,15 @@ #include "sgl.hpp" #include +#include #include #include #include +#include #include #include #include // TODO: Remove this, it is only used for error messages in WKT reader +#include //====================================================================================================================== // Helpers @@ -745,6 +748,105 @@ void ops::visit_polygon_geometries(const geometry &geom, void *state, }); } +void ops::clip_point_to_box(const geometry &geom, const geometry &bbox, geometry &out_geom) { + SGL_ASSERT(geom.get_type() == geometry_type::POINT); + SGL_ASSERT(bbox.get_type() == geometry_type::POLYGON); + SGL_ASSERT(!geom.is_empty()); + + out_geom.set_type(geometry_type::POINT); + // TODO: determine if we want to keep Z. + out_geom.set_z(false); + out_geom.set_m(false); + + auto rect = sgl::extent_xy::smallest(); + if (sgl::ops::get_total_extent_xy(bbox, rect) == 0) { + out_geom.set_vertex_array(nullptr, 0); + return; + } + + const auto vertex = geom.get_vertex_xy(0); + if (rect.contains(vertex)) { + out_geom.set_vertex_array(geom.get_vertex_array(), 1); + } else { + out_geom.set_vertex_array(nullptr, 0); + } +} + +void ops::clip_multipoint_to_box(allocator &alloc, const geometry &geom, const geometry &bbox, geometry &out_geom) { + SGL_ASSERT(geom.get_type() == geometry_type::MULTI_POINT); + SGL_ASSERT(bbox.get_type() == geometry_type::POLYGON); + SGL_ASSERT(!geom.is_empty()); + + auto rect = sgl::extent_xy::smallest(); + if (sgl::ops::get_total_extent_xy(bbox, rect) == 0) { + out_geom.set_vertex_array(nullptr, 0); + return; + } + + // first pass + uint32_t num_points_to_return = 0; + visit_points(geom, [&rect, &num_points_to_return](const geometry &part) { + if (part.is_empty()) { + return; + } + const auto v_array = part.get_vertex_array(); + const auto v_width = part.get_vertex_width(); + + vertex_xyzm vertex = {0, 0, 0, 0}; + memcpy(&vertex, v_array, v_width); + if ((vertex.x < rect.min.x) || (vertex.x > rect.max.x) || (vertex.y < rect.min.y) || (vertex.y > rect.max.y)) { + return; + } + num_points_to_return += 1; + }); + + if (num_points_to_return == 0) { + out_geom.set_type(sgl::geometry_type::GEOMETRY_COLLECTION); + } else if (num_points_to_return == 1) { + out_geom.set_type(sgl::geometry_type::POINT); + // second pass + visit_points(geom, [&rect, &out_geom, &alloc](const geometry &part) { + if (part.is_empty()) { + return; + } + const auto v_array = part.get_vertex_array(); + const auto v_width = part.get_vertex_width(); + + vertex_xyzm vertex = {0, 0, 0, 0}; + memcpy(&vertex, v_array, v_width); + if ((vertex.x < rect.min.x) || (vertex.x > rect.max.x) || (vertex.y < rect.min.y) || (vertex.y > rect.max.y)) { + return; + } + const auto data_mem = static_cast(alloc.alloc(v_width)); + memcpy(data_mem, &vertex, v_width); + out_geom.set_vertex_array(data_mem, 1); + }); + } else { + out_geom.set_type(sgl::geometry_type::MULTI_POINT); + // second pass + visit_points(geom, [&rect, &out_geom, &alloc](const geometry &part) { + if (part.is_empty()) { + return; + } + const auto v_array = part.get_vertex_array(); + const auto v_width = part.get_vertex_width(); + + vertex_xyzm vertex = {0, 0, 0, 0}; + memcpy(&vertex, v_array, v_width); + if ((vertex.x < rect.min.x) || (vertex.x > rect.max.x) || (vertex.y < rect.min.y) || (vertex.y > rect.max.y)) { + return; + } + const auto data_mem = static_cast(alloc.alloc(v_width)); + memcpy(data_mem, &vertex, v_width); + auto point_mem = static_cast(alloc.alloc(sizeof(sgl::geometry))); + auto point_ptr = new (point_mem) sgl::geometry(geometry_type::POINT, part.has_z(), part.has_m()); + point_ptr->set_vertex_array(data_mem, 1); + out_geom.append_part(point_ptr); + }); + } +} + + } // namespace sgl //====================================================================================================================== diff --git a/src/sgl/sgl.hpp b/src/sgl/sgl.hpp index 29f2af62..469ab531 100644 --- a/src/sgl/sgl.hpp +++ b/src/sgl/sgl.hpp @@ -703,6 +703,9 @@ void locate_along(allocator &alloc, const geometry &linear_geom, double measure, void locate_between(allocator &alloc, const geometry &linear_geom, double measure_beg, double measure_end, double offset, geometry &out_geom); +void clip_point_to_box(const geometry &geom, const geometry &bbox, geometry &out_geom); +void clip_multipoint_to_box(allocator &alloc, const geometry &geom, const geometry &bbox, geometry &out_geom); + } // namespace ops // TODO: Move these diff --git a/src/spatial/modules/geos/geos_geometry.hpp b/src/spatial/modules/geos/geos_geometry.hpp index ea871a76..c8811850 100644 --- a/src/spatial/modules/geos/geos_geometry.hpp +++ b/src/spatial/modules/geos/geos_geometry.hpp @@ -91,6 +91,7 @@ class GeosGeometry { GeosGeometry get_coverage_invalid_edges(double tolerance) const; GeosGeometry get_coverage_simplified(double tolerance, bool preserve_boundary) const; GeosGeometry get_coverage_union() const; + GeosGeometry clip_by_rect(double x_min, double y_min, double x_max, double y_max) const; PreparedGeosGeometry get_prepared() const; @@ -480,6 +481,10 @@ inline PreparedGeosGeometry GeosGeometry::get_prepared() const { return PreparedGeosGeometry(handle, *this); } +inline GeosGeometry GeosGeometry::clip_by_rect(double x_min, double y_min, double x_max, double y_max) const { + return GeosGeometry(handle, GEOSClipByRect_r(handle, geom, x_min, y_min, x_max, y_max)); +} + //-- PreparedGeosGeometry --// inline bool PreparedGeosGeometry::contains(const GeosGeometry &other) const { diff --git a/src/spatial/modules/geos/geos_module.cpp b/src/spatial/modules/geos/geos_module.cpp index 0c0cb337..847bbe8c 100644 --- a/src/spatial/modules/geos/geos_module.cpp +++ b/src/spatial/modules/geos/geos_module.cpp @@ -2773,4 +2773,26 @@ void RegisterGEOSModule(ExtensionLoader &loader) { ST_CoverageSimplify_Agg::Register(loader); } +string_t GeosOperations::clip_to_rect(Vector &result, const string_t &blob, double x_min, double y_min, double x_max, double y_max) { + auto ctx = GEOS_init_r(); + + GEOSContext_setErrorMessageHandler_r(ctx, [](const char *message, void *) { throw InvalidInputException(message); }, nullptr); + + const auto geom = GeosSerde::Deserialize(ctx, blob.GetData(), blob.GetSize()); + if (geom == nullptr) { + GEOS_finish_r(ctx); + throw InvalidInputException("Could not deserialize geometry"); + } + auto geos = GeosGeometry(ctx, geom); + auto clipped = geos.clip_by_rect(x_min, y_min, x_max, y_max); + const auto raw = clipped.get_raw(); + const auto size = GeosSerde::GetRequiredSize(ctx, raw); + string_t serialized = StringVector::EmptyString(result, size); + const auto ptr = serialized.GetDataWriteable(); + GeosSerde::Serialize(ctx, raw, ptr, size); + serialized.Finalize(); + GEOS_finish_r(ctx); + return serialized; +} + } // namespace duckdb diff --git a/src/spatial/modules/geos/geos_module.hpp b/src/spatial/modules/geos/geos_module.hpp index fa6b69f0..d11d8c6e 100644 --- a/src/spatial/modules/geos/geos_module.hpp +++ b/src/spatial/modules/geos/geos_module.hpp @@ -1,9 +1,16 @@ #pragma once +#include "duckdb/common/types/string_type.hpp" + namespace duckdb { class ExtensionLoader; void RegisterGEOSModule(ExtensionLoader &loader); +class GeosOperations { +public: + static string_t clip_to_rect(Vector &result, const string_t &blob, double x_min, double y_min, double x_max, double y_max); +}; + } // namespace duckdb diff --git a/src/spatial/modules/main/spatial_functions_scalar.cpp b/src/spatial/modules/main/spatial_functions_scalar.cpp index 6ca172a4..d2715227 100644 --- a/src/spatial/modules/main/spatial_functions_scalar.cpp +++ b/src/spatial/modules/main/spatial_functions_scalar.cpp @@ -22,6 +22,10 @@ // Extra #include "yyjson.h" +#if SPATIAL_USE_GEOS +#include "spatial/modules/geos/geos_module.hpp" +#endif + namespace duckdb { namespace { @@ -1641,6 +1645,84 @@ struct ST_Centroid { } }; +//====================================================================================================================== +// ST_ClipByBox2D +//====================================================================================================================== + +struct ST_ClipByBox2D { + + //------------------------------------------------------------------------------------------------------------------ + // GEOMETRY + //------------------------------------------------------------------------------------------------------------------ + static void ExecuteGeometry(DataChunk &args, ExpressionState &state, Vector &result) { + auto &lstate = LocalState::ResetAndGet(state); + auto &alloc = lstate.GetAllocator(); + + BinaryExecutor::Execute( + args.data[0], args.data[1], result, args.size(), + [&](const string_t &geom_blob, const string_t &bbox_blob) { + sgl::geometry geom; + sgl::geometry bbox; + lstate.Deserialize(geom_blob, geom); + lstate.Deserialize(bbox_blob, bbox); + + if (geom.get_type() == sgl::geometry_type::POINT) { + sgl::geometry clipped; + sgl::ops::clip_point_to_box(geom, bbox, clipped); + return lstate.Serialize(result, clipped); + } else if (geom.get_type() == sgl::geometry_type::MULTI_POINT) { + sgl::geometry clipped; + sgl::ops::clip_multipoint_to_box(alloc, geom, bbox, clipped); + return lstate.Serialize(result, clipped); + } else { +#if SPATIAL_USE_GEOS + auto rect = sgl::extent_xy::smallest(); + if (sgl::ops::get_total_extent_xy(bbox, rect) == 0) { + sgl::geometry clipped; + clipped.set_vertex_array(nullptr, 0); + return lstate.Serialize(result, clipped); + } + return GeosOperations::clip_to_rect(result, geom_blob, rect.min.x, rect.min.y, rect.max.x, rect.max.y); +#else + throw InvalidInputException("ST_ClipByBox2D: input is not a POINT or MULTIPOINT, and GEOS fallback is not available"); +#endif + } + }); + } + + //------------------------------------------------------------------------------------------------------------------ + // Documentation + //------------------------------------------------------------------------------------------------------------------ + static constexpr auto DESCRIPTION = R"( + Clips a geometry by a rectangular bounding box (BBOX). + + This may produce a different kind of geometry. For example, a linestring can be converted into a multilinestring. + )"; + static constexpr auto EXAMPLE = ""; + + //------------------------------------------------------------------------------------------------------------------ + // Register + //------------------------------------------------------------------------------------------------------------------ + static void Register(ExtensionLoader &loader) { + FunctionBuilder::RegisterScalar(loader, "ST_ClipByBox2D", [](ScalarFunctionBuilder &func) { + func.AddVariant([](ScalarFunctionVariantBuilder &variant) { + variant.AddParameter("geom", GeoTypes::GEOMETRY()); + variant.AddParameter("box", GeoTypes::GEOMETRY()); + variant.SetReturnType(GeoTypes::GEOMETRY()); + + variant.SetFunction(ExecuteGeometry); + variant.SetInit(LocalState::Init); + }); + + func.SetDescription(DESCRIPTION); + func.SetExample(EXAMPLE); + + func.SetTag("ext", "spatial"); + func.SetTag("category", "construction"); + }); + } +}; + //====================================================================================================================== // ST_Collect //====================================================================================================================== @@ -9271,6 +9353,7 @@ void RegisterSpatialScalarFunctions(ExtensionLoader &loader) { ST_AsSVG::Register(loader); ST_Azimuth::Register(loader); ST_Centroid::Register(loader); + ST_ClipByBox2D::Register(loader); ST_Collect::Register(loader); ST_CollectionExtract::Register(loader); ST_Contains::Register(loader); diff --git a/test/sql/geometry/st_clipbyrect2d.test b/test/sql/geometry/st_clipbyrect2d.test new file mode 100644 index 00000000..81502dbb --- /dev/null +++ b/test/sql/geometry/st_clipbyrect2d.test @@ -0,0 +1,143 @@ +require spatial + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_MakePoint(0, -1), ST_MakeEnvelope(0, 0, 5, 5 ))); +---- +POINT EMPTY + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_MakePoint(-1, 0), ST_MakeEnvelope(0, 0, 5, 5 ))); +---- +POINT EMPTY + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_MakePoint(0, 0), ST_MakeEnvelope(0, 0, 5, 5 ))); +---- +POINT (0 0) + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_MakePoint(1, 3), ST_MakeEnvelope(0, 0, 5, 5 ))); +---- +POINT (1 3) + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_MakePoint(1, 6), ST_MakeEnvelope(0, 0, 5, 5 ))); +---- +POINT EMPTY + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_MakePoint(7, 3), ST_MakeEnvelope(0, 0, 5, 5 ))); +---- +POINT EMPTY + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_MakePoint(5, 4), ST_MakeEnvelope(0, 0, 5, 5 ))); +---- +POINT (5 4) + + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_GeomFromText('MULTIPOINT((1 2),(5 5),(10 5))'), ST_MakeEnvelope(0, 0, 4, 4))); +---- +POINT (1 2) + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_GeomFromText('MULTIPOINT((1 2),(5 5),(10 5))'), ST_MakeEnvelope(0, 0, 8, 8))); +---- +MULTIPOINT (1 2, 5 5) + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_GeomFromText('MULTIPOINT((1 2),(5 5),(10 5))'), ST_MakeEnvelope(2, 3, 4, 4))); +---- +GEOMETRYCOLLECTION EMPTY + + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_MakeLine(ST_Point(2, 3), ST_Point(4,5)), ST_MakeEnvelope(2, 2, 5, 5))); +---- +LINESTRING (2 3, 4 5) + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_GeomFromText('LINESTRING (2 3, 3 3, 3 2, 4 5)'), ST_MakeEnvelope(2, 2, 5, 5))); +---- +LINESTRING (2 3, 3 3, 3 2, 4 5) + +# This returns LINESTRING EMPTY on PostGIS +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_MakeLine(ST_Point(0, 0), ST_Point(1,1)), ST_MakeEnvelope(2, 2, 5, 5))); +---- +GEOMETRYCOLLECTION EMPTY + +# This returns LINESTRING EMPTY on PostGIS +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_MakeLine(ST_Point(0, 0), ST_Point(1, 3)), ST_MakeEnvelope(2, 2, 5, 5))); +---- +GEOMETRYCOLLECTION EMPTY + +# This returns LINESTRING EMPTY on PostGIS +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_MakeLine(ST_Point(0, 0), ST_Point(1, 7)), ST_MakeEnvelope(2, 2, 5, 5))); +---- +GEOMETRYCOLLECTION EMPTY + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_MakeLine(ST_Point(2.1, 2.1), ST_Point(2,4)), ST_MakeEnvelope(2, 2, 5, 5))); +---- +LINESTRING (2.1 2.1, 2 4) + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_GeomFromText('LINESTRING (0 0, 10 20)'), ST_MakeEnvelope(1, 2, 5, 5) ) ); +---- +LINESTRING (1 2, 2.5 5) + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_GeomFromText('LINESTRING (0 0, 20 20)'), ST_MakeBox2D(ST_Point(1, 2), ST_Point(5, 5) ) ) ); +---- +LINESTRING (2 2, 5 5) + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_GeomFromText('LINESTRING(0 0, 5 5, 10 5)'), ST_MakeEnvelope(1, 2, 8, 8))); +---- +LINESTRING (2 2, 5 5, 8 5) + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_GeomFromText('LINESTRING(0 0, 5 5, 10 5, 6 2)'), ST_MakeEnvelope(1, 2, 8, 8))); +---- +MULTILINESTRING ((2 2, 5 5, 8 5), (8 3.5, 6 2)) + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_MakeLine(ST_Point(0, 0), ST_Point(2,2)), ST_MakeEnvelope(2, 2, 5, 5))); +---- +GEOMETRYCOLLECTION EMPTY + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_GeomFromText('POLYGON((0 0, 5 5, 10 5, 5 0, 0 0))'), ST_MakeEnvelope(2, 2, 5, 5))); +---- +POLYGON ((2 2, 5 5, 5 2, 2 2)) + +# This returns POLYGON EMPTY on PostGIS +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_GeomFromText('POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))'), ST_MakeEnvelope(2, 2, 5, 5))); +---- +GEOMETRYCOLLECTION EMPTY + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_GeomFromText('POLYGON((0 0, 0 3, 4 3, 4 0, 0 0))'), ST_MakeEnvelope(2, 2, 5, 5))); +---- +POLYGON ((2 2, 2 3, 4 3, 4 2, 2 2)) + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_GeomFromText('POLYGON((0 0, 0 3, 4 3, 4 0, 0 0))'), ST_MakeEnvelope(0, 2, 2, 5))); +---- +POLYGON ((0 2, 0 3, 2 3, 2 2, 0 2)) + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_GeomFromText('MULTILINESTRING((0 0, 0 3, 4 3, 4 0, 0 0))'), ST_MakeEnvelope(2, 2, 5, 5))); +---- +LINESTRING (2 3, 4 3, 4 2) + +query I +SELECT ST_AsText(ST_ClipByBox2D(ST_GeomFromText('MULTILINESTRING((0 0, 0 3, 4 3, 4 0, 0 0),(3 3, 3 0))'), ST_MakeEnvelope(2, 2, 5, 5))); +---- +MULTILINESTRING ((2 3, 4 3, 4 2), (3 3, 3 2)) +