diff --git a/spatial/include/spatial/core/functions/scalar.hpp b/spatial/include/spatial/core/functions/scalar.hpp index 4220d204..99c8889a 100644 --- a/spatial/include/spatial/core/functions/scalar.hpp +++ b/spatial/include/spatial/core/functions/scalar.hpp @@ -30,7 +30,9 @@ struct CoreScalarFunctions { RegisterStIntersectsExtent(db); RegisterStIsEmpty(db); RegisterStLength(db); + RegisterStMakeEnvelope(db); RegisterStMakeLine(db); + RegisterStMakePolygon(db); RegisterStNGeometries(db); RegisterStNInteriorRings(db); RegisterStNPoints(db); @@ -114,9 +116,15 @@ struct CoreScalarFunctions { // ST_Length static void RegisterStLength(DatabaseInstance &db); + // ST_MakeEnvelope + static void RegisterStMakeEnvelope(DatabaseInstance &db); + // ST_MakeLine static void RegisterStMakeLine(DatabaseInstance &db); + // ST_MakePolygon + static void RegisterStMakePolygon(DatabaseInstance &db); + // ST_NGeometries static void RegisterStNGeometries(DatabaseInstance &db); diff --git a/spatial/src/spatial/core/functions/scalar/CMakeLists.txt b/spatial/src/spatial/core/functions/scalar/CMakeLists.txt index 5d96f43a..6cb7f1c3 100644 --- a/spatial/src/spatial/core/functions/scalar/CMakeLists.txt +++ b/spatial/src/spatial/core/functions/scalar/CMakeLists.txt @@ -21,7 +21,9 @@ set(EXTENSION_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/st_intersects.cpp ${CMAKE_CURRENT_SOURCE_DIR}/st_intersects_extent.cpp ${CMAKE_CURRENT_SOURCE_DIR}/st_length.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/st_makeenvelope.cpp ${CMAKE_CURRENT_SOURCE_DIR}/st_makeline.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/st_makepolygon.cpp ${CMAKE_CURRENT_SOURCE_DIR}/st_ngeometries.cpp ${CMAKE_CURRENT_SOURCE_DIR}/st_ninteriorrings.cpp ${CMAKE_CURRENT_SOURCE_DIR}/st_npoints.cpp diff --git a/spatial/src/spatial/core/functions/scalar/st_makeenvelope.cpp b/spatial/src/spatial/core/functions/scalar/st_makeenvelope.cpp new file mode 100644 index 00000000..10fa0f07 --- /dev/null +++ b/spatial/src/spatial/core/functions/scalar/st_makeenvelope.cpp @@ -0,0 +1,55 @@ +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/core/functions/scalar.hpp" +#include "spatial/core/functions/common.hpp" +#include "spatial/core/geometry/geometry.hpp" + +#include "duckdb/parser/parsed_data/create_scalar_function_info.hpp" +#include "duckdb/common/vector_operations/generic_executor.hpp" + +namespace spatial { + +namespace core { + +static void MakeEnvelopeFunction(DataChunk &args, ExpressionState &state, Vector &result) { + auto &lstate = GeometryFunctionLocalState::ResetAndGet(state); + auto count = args.size(); + + auto &min_x_vec = args.data[0]; + auto &min_y_vec = args.data[1]; + auto &max_x_vec = args.data[2]; + auto &max_y_vec = args.data[3]; + + using DOUBLE_TYPE = PrimitiveType; + using GEOMETRY_TYPE = PrimitiveType; + + GenericExecutor::ExecuteQuaternary( + min_x_vec, min_y_vec, max_x_vec, max_y_vec, result, count, + [&](DOUBLE_TYPE x_min, DOUBLE_TYPE y_min, DOUBLE_TYPE x_max, DOUBLE_TYPE y_max) { + uint32_t capacity = 5; + auto envelope_geom = lstate.factory.CreatePolygon(1, &capacity); + // Create the exterior ring in CCW order + envelope_geom.Shell().Add(Vertex(x_min.val, y_min.val)); + envelope_geom.Shell().Add(Vertex(x_min.val, y_max.val)); + envelope_geom.Shell().Add(Vertex(x_max.val, y_max.val)); + envelope_geom.Shell().Add(Vertex(x_max.val, y_min.val)); + envelope_geom.Shell().Add(Vertex(x_min.val, y_min.val)); + + return lstate.factory.Serialize(result, Geometry(envelope_geom)); + }); +} + +void CoreScalarFunctions::RegisterStMakeEnvelope(DatabaseInstance &db) { + + ScalarFunctionSet set("ST_MakeEnvelope"); + + set.AddFunction(ScalarFunction({LogicalType::DOUBLE, LogicalType::DOUBLE, LogicalType::DOUBLE, LogicalType::DOUBLE}, + GeoTypes::GEOMETRY(), MakeEnvelopeFunction, nullptr, nullptr, nullptr, + GeometryFunctionLocalState::Init)); + + ExtensionUtil::RegisterFunction(db, set); +} + +} // namespace core + +} // namespace spatial diff --git a/spatial/src/spatial/core/functions/scalar/st_makepolygon.cpp b/spatial/src/spatial/core/functions/scalar/st_makepolygon.cpp new file mode 100644 index 00000000..3c4ad112 --- /dev/null +++ b/spatial/src/spatial/core/functions/scalar/st_makepolygon.cpp @@ -0,0 +1,143 @@ +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/core/functions/scalar.hpp" +#include "spatial/core/functions/common.hpp" +#include "spatial/core/geometry/geometry.hpp" + +#include "duckdb/parser/parsed_data/create_scalar_function_info.hpp" +#include "duckdb/common/vector_operations/unary_executor.hpp" +#include "duckdb/common/vector_operations/binary_executor.hpp" + +namespace spatial { + +namespace core { + +static void MakePolygonFromRingsFunction(DataChunk &args, ExpressionState &state, Vector &result) { + auto &lstate = GeometryFunctionLocalState::ResetAndGet(state); + auto count = args.size(); + + auto &child_vec = ListVector::GetEntry(args.data[1]); + UnifiedVectorFormat format; + child_vec.ToUnifiedFormat(count, format); + + BinaryExecutor::Execute( + args.data[0], args.data[1], result, count, [&](string_t line_blob, list_entry_t &rings_list) { + // First, setup the shell + auto shell_geom = lstate.factory.Deserialize(line_blob); + if (shell_geom.Type() != GeometryType::LINESTRING) { + throw InvalidInputException("ST_MakePolygon only accepts LINESTRING geometries"); + } + + auto &shell = shell_geom.GetLineString(); + auto shell_vert_count = shell.Count(); + if (shell_vert_count < 4) { + throw InvalidInputException("ST_MakePolygon shell requires at least 4 vertices"); + } + + auto &shell_verts = shell.Vertices(); + if (shell_verts[0] != shell_verts[shell_vert_count - 1]) { + throw InvalidInputException( + "ST_MakePolygon shell must be closed (first and last vertex must be equal)"); + } + + // Validate and count the hole ring sizes + auto holes_offset = rings_list.offset; + auto holes_length = rings_list.length; + + vector rings_counts; + vector rings; + rings_counts.push_back(shell_vert_count); + rings.push_back(shell); + + for (idx_t hole_idx = 0; hole_idx < holes_length; hole_idx++) { + auto mapped_idx = format.sel->get_index(holes_offset + hole_idx); + if (!format.validity.RowIsValid(mapped_idx)) { + continue; + } + + auto geometry_blob = UnifiedVectorFormat::GetData(format)[mapped_idx]; + auto hole_geometry = lstate.factory.Deserialize(geometry_blob); + + if (hole_geometry.Type() != GeometryType::LINESTRING) { + throw InvalidInputException( + StringUtil::Format("ST_MakePolygon hole #%lu is not a LINESTRING geometry", hole_idx + 1)); + } + auto &hole = hole_geometry.GetLineString(); + auto hole_count = hole.Count(); + if (hole_count < 4) { + throw InvalidInputException( + StringUtil::Format("ST_MakePolygon hole #%lu requires at least 4 vertices", hole_idx + 1)); + } + + auto &ring_verts = hole.Vertices(); + if (ring_verts[0] != ring_verts[hole_count - 1]) { + throw InvalidInputException(StringUtil::Format( + "ST_MakePolygon hole #%lu must be closed (first and last vertex must be equal)", hole_idx + 1)); + } + + rings_counts.push_back(hole_count); + rings.push_back(hole); + } + + auto polygon = lstate.factory.CreatePolygon(rings_counts.size(), rings_counts.data()); + + for (auto ring_idx = 0; ring_idx < rings.size(); ring_idx++) { + auto &new_ring = rings[ring_idx]; + auto &poly_ring = polygon.Ring(ring_idx); + for (auto &v : new_ring.Vertices()) { + poly_ring.Add(v); + } + } + + return lstate.factory.Serialize(result, Geometry(polygon)); + }); +} + +static void MakePolygonFromShellFunction(DataChunk &args, ExpressionState &state, Vector &result) { + auto &lstate = GeometryFunctionLocalState::ResetAndGet(state); + auto count = args.size(); + + UnaryExecutor::Execute(args.data[0], result, count, [&](string_t &line_blob) { + auto line_geom = lstate.factory.Deserialize(line_blob); + + if (line_geom.Type() != GeometryType::LINESTRING) { + throw InvalidInputException("ST_MakePolygon only accepts LINESTRING geometries"); + } + + auto &line = line_geom.GetLineString(); + auto line_count = line.Count(); + if (line_count < 4) { + throw InvalidInputException("ST_MakePolygon shell requires at least 4 vertices"); + } + + auto &line_verts = line.Vertices(); + if (line_verts[0] != line_verts[line_count - 1]) { + throw InvalidInputException("ST_MakePolygon shell must be closed (first and last vertex must be equal)"); + } + + auto polygon = lstate.factory.CreatePolygon(1, &line_count); + for (auto &v : line_verts) { + polygon.Shell().Add(v); + } + + return lstate.factory.Serialize(result, Geometry(polygon)); + }); +} + +void CoreScalarFunctions::RegisterStMakePolygon(DatabaseInstance &db) { + + ScalarFunctionSet set("ST_MakePolygon"); + + set.AddFunction(ScalarFunction({GeoTypes::GEOMETRY(), LogicalType::LIST(GeoTypes::GEOMETRY())}, + GeoTypes::GEOMETRY(), MakePolygonFromRingsFunction, nullptr, nullptr, nullptr, + GeometryFunctionLocalState::Init)); + + set.AddFunction(ScalarFunction({GeoTypes::GEOMETRY()}, GeoTypes::GEOMETRY(), MakePolygonFromShellFunction, nullptr, + nullptr, nullptr, GeometryFunctionLocalState::Init)); + + ExtensionUtil::RegisterFunction(db, set); +} + +} // namespace core + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/file_handler.cpp b/spatial/src/spatial/gdal/file_handler.cpp index e7e9c13e..fa041ea6 100644 --- a/spatial/src/spatial/gdal/file_handler.cpp +++ b/spatial/src/spatial/gdal/file_handler.cpp @@ -40,21 +40,34 @@ static void *DuckDBOpen(void *, const char *file_name, const char *access) { if (len > 1 && access[1] == '+') { flags |= FileFlags::FILE_FLAGS_WRITE; } + if (len > 2 && access[2] == '+') { + // might be "rb+" + flags |= FileFlags::FILE_FLAGS_WRITE; + } } else if (access[0] == 'w') { flags = FileFlags::FILE_FLAGS_WRITE | FileFlags::FILE_FLAGS_FILE_CREATE_NEW; if (len > 1 && access[1] == '+') { flags |= FileFlags::FILE_FLAGS_READ; } + if (len > 2 && access[2] == '+') { + // might be "wb+" + flags |= FileFlags::FILE_FLAGS_READ; + } } else if (access[0] == 'a') { flags = FileFlags::FILE_FLAGS_APPEND; if (len > 1 && access[1] == '+') { flags |= FileFlags::FILE_FLAGS_READ; } + if (len > 2 && access[2] == '+') { + // might be "ab+" + flags |= FileFlags::FILE_FLAGS_READ; + } } else { throw InternalException("Unknown file access type"); } try { + string path(file_name); auto file = fs.OpenFile(file_name, flags); return file.release(); } catch (std::exception &ex) { diff --git a/spatial/src/spatial/gdal/functions/st_write.cpp b/spatial/src/spatial/gdal/functions/st_write.cpp index 02e44086..36a71daa 100644 --- a/spatial/src/spatial/gdal/functions/st_write.cpp +++ b/spatial/src/spatial/gdal/functions/st_write.cpp @@ -251,8 +251,11 @@ static unique_ptr InitGlobal(ClientContext &context, Functio if (!gdal_data.target_srs.empty()) { srs.SetFromUserInput(gdal_data.target_srs.c_str()); } + // Not all GDAL drivers check if the SRS is empty (cough cough GeoJSONSeq) + // so we have to pass nullptr if we want the default behavior. + OGRSpatialReference *srs_ptr = gdal_data.target_srs.empty() ? nullptr : &srs; - auto layer = dataset->CreateLayer(gdal_data.layer_name.c_str(), &srs, wkbUnknown, lco); + auto layer = dataset->CreateLayer(gdal_data.layer_name.c_str(), srs_ptr, wkbUnknown, lco); if (!layer) { throw IOException("Could not create layer"); } diff --git a/spatial/src/spatial/gdal/module.cpp b/spatial/src/spatial/gdal/module.cpp index 1c4e5967..aacd0d84 100644 --- a/spatial/src/spatial/gdal/module.cpp +++ b/spatial/src/spatial/gdal/module.cpp @@ -20,7 +20,12 @@ void GdalModule::Register(DatabaseInstance &db) { OGRRegisterAllInternal(); // Set GDAL error handler + CPLSetErrorHandler([](CPLErr e, int code, const char *msg) { + // DuckDB doesnt do warnings, so we only throw on errors + if (e != CE_Failure && e != CE_Fatal) { + return; + } switch (code) { case CPLE_NoWriteAccess: throw PermissionException("GDAL Error (%d): %s", code, msg); diff --git a/test/sql/gdal/st_write_basic.test b/test/sql/gdal/st_write_basic.test new file mode 100644 index 00000000..f1b3a91a --- /dev/null +++ b/test/sql/gdal/st_write_basic.test @@ -0,0 +1,15 @@ +require spatial + +statement ok +COPY (SELECT ST_GeomFromText('POINT (1 1)')) +TO '__TEST_DIR__/test_seq.json' +WITH (FORMAT GDAL, DRIVER 'GeoJSONSeq'); + +statement ok +COPY (SELECT ST_GeomFromText('POINT (1 1)')) +TO '__TEST_DIR__/test_seq.shp' +WITH (FORMAT GDAL, DRIVER 'ESRI ShapeFile'); + +# MVT is broken due to threading issues +#statement ok +#COPY (SELECT ST_GeomFromText('POINT (1 1)')) TO '__TEST_DIR__/test_mvt.mvt' WITH (FORMAT GDAL, DRIVER 'MVT'); \ No newline at end of file diff --git a/test/sql/geometry/st_makeenvelope.test b/test/sql/geometry/st_makeenvelope.test new file mode 100644 index 00000000..cbd87136 --- /dev/null +++ b/test/sql/geometry/st_makeenvelope.test @@ -0,0 +1,11 @@ +require spatial + +query I +SELECT ST_AsText(ST_MakeEnvelope(5, 5, 10, 10)); +---- +POLYGON ((5 5, 5 10, 10 10, 10 5, 5 5)) + +query I +SELECT ST_AsText(ST_MakeEnvelope(5, NULL, 10, 10)); +---- +NULL \ No newline at end of file diff --git a/test/sql/geometry/st_makepolygon.test b/test/sql/geometry/st_makepolygon.test new file mode 100644 index 00000000..047e5d93 --- /dev/null +++ b/test/sql/geometry/st_makepolygon.test @@ -0,0 +1,116 @@ +require spatial +# +## Single shell input tests (no holes) +#query I +#SELECT ST_AsText(ST_MakePolygon(ST_GeomFromText('LINESTRING(0 0, 0 10, 10 0, 0 0)'))); +#---- +#POLYGON ((0 0, 0 10, 10 0, 0 0)) +# +#query I +#SELECT ST_AsText(ST_MakePolygon(NULL)); +#---- +#NULL +# +#statement error +#SELECT st_makepolygon(ST_GeomFromText('POINT(0 0)')); +#---- +#ST_MakePolygon only accepts LINESTRING geometries +# +#statement error +#SELECT ST_MakePolygon(ST_GeomFromText('LINESTRING(0 0, 1 1)')); +#---- +#ST_MakePolygon shell requires at least 4 vertices +# +#statement error +#SELECT ST_MakePolygon(ST_GeomFromText('LINESTRING(0 0, 1 1, 2 2, 3 3)')); +#---- +#ST_MakePolygon shell must be closed (first and last vertex must be equal) + +# With holes tests +query I +SELECT ST_AsText( + ST_MakePolygon( + ST_GeomFromText('LINESTRING(0 0, 0 10, 10 0, 0 0)'), + [ + ST_GeomFromText('LINESTRING(1 1, 1 2, 2 2, 2 1, 1 1)'), + ST_GeomFromText('LINESTRING(3 3, 3 4, 4 4, 4 3, 3 3)') + ] + ) +); +---- +POLYGON ((0 0, 0 10, 10 0, 0 0), (1 1, 1 2, 2 2, 2 1, 1 1), (3 3, 3 4, 4 4, 4 3, 3 3)) + +# Normal null handling, null arg returns null +query I +SELECT ST_MakePolygon(ST_GeomFromText('LINESTRING(0 0, 0 10, 10 0, 0 0)'), NULL); +---- +NULL + +# Null rings are ignored +query I +SELECT ST_AsText( + ST_MakePolygon( + ST_GeomFromText('LINESTRING(0 0, 0 10, 10 0, 0 0)'), + [ + NULL, + ST_GeomFromText('LINESTRING(3 3, 3 4, 4 4, 4 3, 3 3)'), + NULL + ] + ) +); +---- +POLYGON ((0 0, 0 10, 10 0, 0 0), (3 3, 3 4, 4 4, 4 3, 3 3)) + +# Rings require at least 4 vertices +statement error +SELECT ST_AsText( + ST_MakePolygon( + ST_GeomFromText('LINESTRING(0 0, 0 10, 10 0, 0 0)'), + [ + ST_GeomFromText('LINESTRING(3 3, 3 4, 4 4)') + ] + ) +); +---- +ST_MakePolygon hole #1 requires at least 4 vertices + +# Rings must be closed +statement error +SELECT ST_AsText( + ST_MakePolygon( + ST_GeomFromText('LINESTRING(0 0, 0 10, 10 0, 0 0)'), + [ + ST_GeomFromText('LINESTRING(3 3, 3 4, 4 4, 5 5)') + ] + ) +); +---- +ST_MakePolygon hole #1 must be closed (first and last vertex must be equal) + +statement error +SELECT ST_AsText( + ST_MakePolygon( + ST_GeomFromText('LINESTRING(0 0, 0 10, 10 0, 0 0)'), + [ + ST_GeomFromText('LINESTRING(1 1, 1 2, 2 2, 2 1, 1 1)'), + ST_GeomFromText('LINESTRING(3 3, 3 4, 4 4, 4 3, 5 5)') + ] + ) +); +---- +ST_MakePolygon hole #2 must be closed (first and last vertex must be equal) + + +# Must be linestring geometries +statement error +SELECT ST_AsText( + ST_MakePolygon( + ST_GeomFromText('LINESTRING(0 0, 0 10, 10 0, 0 0)'), + [ + ST_Point(3, 3), + ST_GeomFromText('LINESTRING(1 1, 1 2, 2 2, 2 1, 1 1)') + ] + ) +); +---- +ST_MakePolygon hole #1 is not a LINESTRING geometry \ No newline at end of file