diff --git a/spatial/include/spatial/core/types.hpp b/spatial/include/spatial/core/types.hpp index 79133d70..e99d78f9 100644 --- a/spatial/include/spatial/core/types.hpp +++ b/spatial/include/spatial/core/types.hpp @@ -15,6 +15,9 @@ struct GeoTypes { static LogicalType GEOMETRY(); static LogicalType WKB_BLOB(); + static LogicalType RASTER(); + static LogicalType RASTER_COORD(); + static void Register(DatabaseInstance &db); static LogicalType CreateEnumType(const string &name, const vector &members); diff --git a/spatial/include/spatial/gdal/file_handler.hpp b/spatial/include/spatial/gdal/file_handler.hpp index cd3fd9bc..0253cac7 100644 --- a/spatial/include/spatial/gdal/file_handler.hpp +++ b/spatial/include/spatial/gdal/file_handler.hpp @@ -1,6 +1,7 @@ #pragma once #include "spatial/common.hpp" +#include "spatial/gdal/raster/raster_registry.hpp" #include "duckdb/main/database.hpp" namespace spatial { @@ -12,13 +13,17 @@ class DuckDBFileSystemHandler; class GDALClientContextState : public ClientContextState { string client_prefix; DuckDBFileSystemHandler *fs_handler; + std::map> registries; + std::mutex lock; public: explicit GDALClientContextState(ClientContext &context); ~GDALClientContextState() override; void QueryEnd() override; + void QueryEnd(ClientContext &context) override; string GetPrefix(const string &value) const; static GDALClientContextState &GetOrCreate(ClientContext &context); + RasterRegistry &GetRasterRegistry(ClientContext &context); }; } // namespace gdal diff --git a/spatial/include/spatial/gdal/functions.hpp b/spatial/include/spatial/gdal/functions.hpp index 0a4bb2a5..8751c9e6 100644 --- a/spatial/include/spatial/gdal/functions.hpp +++ b/spatial/include/spatial/gdal/functions.hpp @@ -34,6 +34,22 @@ struct GdalTableFunction : ArrowTableFunction { static void Register(DatabaseInstance &db); }; +struct GdalRasterTableFunction { +private: + static unique_ptr Bind(ClientContext &context, TableFunctionBindInput &input, + vector &return_types, vector &names); + + static void Execute(ClientContext &context, TableFunctionInput &input, DataChunk &output); + + static unique_ptr Cardinality(ClientContext &context, const FunctionData *data); + + static unique_ptr ReplacementScan(ClientContext &context, ReplacementScanInput &input, + optional_ptr data); + +public: + static void Register(DatabaseInstance &db); +}; + struct GdalDriversTableFunction { struct BindData : public TableFunctionData { @@ -61,6 +77,14 @@ struct GdalMetadataFunction { static void Register(DatabaseInstance &db); }; +struct GdalRasterMetadataFunction { + static void Register(DatabaseInstance &db); +}; + +struct GdalRasterCopyFunction { + static void Register(DatabaseInstance &db); +}; + } // namespace gdal } // namespace spatial \ No newline at end of file diff --git a/spatial/include/spatial/gdal/functions/aggregate.hpp b/spatial/include/spatial/gdal/functions/aggregate.hpp new file mode 100644 index 00000000..a349f0ec --- /dev/null +++ b/spatial/include/spatial/gdal/functions/aggregate.hpp @@ -0,0 +1,22 @@ +#pragma once +#include "spatial/common.hpp" + +namespace spatial { + +namespace gdal { + +struct GdalAggregateFunctions { +public: + static void Register(DatabaseInstance &db) { + RegisterStRasterMosaicAgg(db); + RegisterStRasterUnionAgg(db); + } + +private: + static void RegisterStRasterMosaicAgg(DatabaseInstance &db); + static void RegisterStRasterUnionAgg(DatabaseInstance &db); +}; + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/include/spatial/gdal/functions/cast.hpp b/spatial/include/spatial/gdal/functions/cast.hpp new file mode 100644 index 00000000..addbaa81 --- /dev/null +++ b/spatial/include/spatial/gdal/functions/cast.hpp @@ -0,0 +1,15 @@ +#pragma once +#include "spatial/common.hpp" + +namespace spatial { + +namespace gdal { + +struct GdalCastFunctions { +public: + static void Register(DatabaseInstance &db); +}; + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/include/spatial/gdal/functions/raster_agg.hpp b/spatial/include/spatial/gdal/functions/raster_agg.hpp new file mode 100644 index 00000000..1e9a1ade --- /dev/null +++ b/spatial/include/spatial/gdal/functions/raster_agg.hpp @@ -0,0 +1,114 @@ +#pragma once +#include "duckdb/function/aggregate_function.hpp" +#include "spatial/common.hpp" + +class GDALDataset; + +namespace spatial { + +namespace gdal { + +struct RasterAggState { + bool is_set; + std::vector datasets; +}; + +struct RasterAggUnaryOperation { + + template + static void Initialize(STATE &state) { + state.is_set = false; + } + + template + static void Combine(const STATE &source, STATE &target, AggregateInputData &) { + if (!source.is_set) { + return; + } + if (!target.is_set) { + target = source; + return; + } + target.datasets.insert(target.datasets.end(), source.datasets.begin(), source.datasets.end()); + } + + template + static void Operation(STATE &state, const INPUT_TYPE &input, AggregateUnaryInput &) { + GDALDataset *dataset = reinterpret_cast(input); + state.is_set = true; + state.datasets.emplace_back(dataset); + } + + template + static void ConstantOperation(STATE &state, const INPUT_TYPE &input, AggregateUnaryInput &agg, idx_t) { + Operation(state, input, agg); + } + + static bool IgnoreNull() { + return true; + } +}; + +struct RasterAggBinaryOperation { + + template + static void Initialize(STATE &state) { + state.is_set = false; + } + + template + static void Combine(const STATE &source, STATE &target, AggregateInputData &) { + if (!source.is_set) { + return; + } + if (!target.is_set) { + target = source; + return; + } + target.datasets.insert(target.datasets.end(), source.datasets.begin(), source.datasets.end()); + } + + template + static void Operation(STATE &state, const INPUT_TYPE &input, const OPTS_TYPE &opts, AggregateBinaryInput &) { + GDALDataset *dataset = reinterpret_cast(input); + state.is_set = true; + state.datasets.emplace_back(dataset); + } + + template + static void ConstantOperation(STATE &state, const INPUT_TYPE &input, const OPTS_TYPE &opts, + AggregateBinaryInput &agg, idx_t) { + Operation(state, input, opts, agg); + } + + static bool IgnoreNull() { + return true; + } +}; + +struct RasterAggBindData : public FunctionData { + //! The client context for the function call + ClientContext &context; + //! The list of options for the function + std::vector options; + + explicit RasterAggBindData(ClientContext &context, std::vector options) + : context(context), options(options) { + } + + unique_ptr Copy() const override { + return make_uniq(context, options); + } + + bool Equals(const FunctionData &other_p) const override { + auto &other = other_p.Cast(); + return options == other.options; + } +}; + +unique_ptr BindRasterAggOperation(ClientContext &context, AggregateFunction &function, + vector> &arguments); + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/include/spatial/gdal/functions/scalar.hpp b/spatial/include/spatial/gdal/functions/scalar.hpp new file mode 100644 index 00000000..60b79fd7 --- /dev/null +++ b/spatial/include/spatial/gdal/functions/scalar.hpp @@ -0,0 +1,84 @@ +#pragma once +#include "spatial/common.hpp" + +namespace spatial { + +namespace gdal { + +struct GdalScalarFunctions { +public: + static void Register(DatabaseInstance &db) { + RegisterStGetSRID(db); + RegisterStGetGeometry(db); + RegisterStGetMetadata(db); + RegisterStGetHasNoBand(db); + RegisterStBandPixelType(db); + RegisterStBandPixelTypeName(db); + RegisterStBandNoDataValue(db); + RegisterStBandColorInterp(db); + RegisterStBandColorInterpName(db); + RegisterStRasterToWorldCoord(db); + RegisterStWorldToRasterCoord(db); + RegisterStGetValue(db); + RegisterStRasterFromFile(db); + RegisterStRasterAsFile(db); + RegisterStRasterWarp(db); + RegisterStRasterClip(db); + } + +private: + // ST_SRID + static void RegisterStGetSRID(DatabaseInstance &db); + + // ST_GetGeometry + static void RegisterStGetGeometry(DatabaseInstance &db); + + // Raster Accessors: + // ST_Width, ST_Height, ST_NumBands, + // ST_UpperLeftX, ST_UpperLeftY, ST_ScaleX, ST_ScaleY, ST_SkewX, ST_SkewY, + // ST_PixelWidth, ST_PixelHeight + static void RegisterStGetMetadata(DatabaseInstance &db); + + // ST_HasNoBand + static void RegisterStGetHasNoBand(DatabaseInstance &db); + + // ST_GetBandPixelType + static void RegisterStBandPixelType(DatabaseInstance &db); + + // ST_GetBandPixelTypeName + static void RegisterStBandPixelTypeName(DatabaseInstance &db); + + // ST_GetBandNoDataValue + static void RegisterStBandNoDataValue(DatabaseInstance &db); + + // ST_GetBandColorInterp + static void RegisterStBandColorInterp(DatabaseInstance &db); + + // ST_GetBandColorInterpName + static void RegisterStBandColorInterpName(DatabaseInstance &db); + + // ST_RasterToWorldCoord[XY] + static void RegisterStRasterToWorldCoord(DatabaseInstance &db); + + // ST_WorldToRasterCoord[XY] + static void RegisterStWorldToRasterCoord(DatabaseInstance &db); + + // ST_Value + static void RegisterStGetValue(DatabaseInstance &db); + + // ST_RasterFromFile + static void RegisterStRasterFromFile(DatabaseInstance &db); + + // ST_RasterAsFile + static void RegisterStRasterAsFile(DatabaseInstance &db); + + // ST_RasterWarp + static void RegisterStRasterWarp(DatabaseInstance &db); + + // ST_RasterClip + static void RegisterStRasterClip(DatabaseInstance &db); +}; + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/include/spatial/gdal/raster/raster.hpp b/spatial/include/spatial/gdal/raster/raster.hpp new file mode 100644 index 00000000..0d1ce753 --- /dev/null +++ b/spatial/include/spatial/gdal/raster/raster.hpp @@ -0,0 +1,90 @@ +#pragma once +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/core/geometry/geometry.hpp" +#include "spatial/core/geometry/geometry_type.hpp" +#include "spatial/gdal/types.hpp" + +class GDALDataset; + +using namespace spatial::core; + +namespace spatial { + +namespace gdal { + +//! A simple wrapper of GDALDataset with useful methods to manage raster data. +//! Does not take ownership of the pointer. +class Raster { +public: + //! Constructor + Raster(GDALDataset *dataset); + + //! Returns the pointer to the dataset managed + GDALDataset *operator->() const noexcept { + return dataset_; + } + //! Returns the pointer to the dataset managed + GDALDataset *get() const noexcept { + return dataset_; + } + + //! Returns the raster width in pixels + int GetRasterXSize() const; + + //! Returns the raster height in pixels + int GetRasterYSize() const; + + //! Returns the number of raster bands + int GetRasterCount() const; + + //! Returns the spatial reference identifier of the raster + int32_t GetSrid() const; + + //! Gets the geometric transform matrix (double[6]) of the raster + bool GetGeoTransform(double *matrix) const; + + //! Gets the inverse geometric transform matrix (double[6]) of the raster + bool GetInvGeoTransform(double *inv_matrix) const; + + //! Returns the polygon representation of the extent of the raster + Geometry GetGeometry(ArenaAllocator &allocator) const; + + //! Returns the geometric X and Y (longitude and latitude) given a column and row + bool RasterToWorldCoord(PointXY &point, int32_t col, int32_t row) const; + + //! Returns the upper left corner as column and row given geometric X and Y + bool WorldToRasterCoord(RasterCoord &coord, double x, double y) const; + + //! Returns the value of a given band in a given col and row pixel + bool GetValue(double &value, int32_t band_num, int32_t col, int32_t row) const; + +public: + //! Returns the geometric X and Y (longitude and latitude) given a column and row + static bool RasterToWorldCoord(PointXY &point, double matrix[], int32_t col, int32_t row); + + //! Returns the upper left corner as column and row given geometric X and Y + static bool WorldToRasterCoord(RasterCoord &coord, double inv_matrix[], double x, double y); + + //! Builds a VRT from a list of Rasters + static GDALDataset *BuildVRT(const std::vector &datasets, + const std::vector &options = std::vector()); + + //! Performs mosaicing, reprojection and/or warping on a raster + static GDALDataset *Warp(GDALDataset *dataset, + const std::vector &options = std::vector()); + + //! Returns a raster that is clipped by the input geometry + static GDALDataset *Clip(GDALDataset *dataset, const geometry_t &geometry, + const std::vector &options = std::vector()); + + //! Get the last error message. + static string GetLastErrorMsg(); + +private: + GDALDataset *dataset_; +}; + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/include/spatial/gdal/raster/raster_factory.hpp b/spatial/include/spatial/gdal/raster/raster_factory.hpp new file mode 100644 index 00000000..a9fa3579 --- /dev/null +++ b/spatial/include/spatial/gdal/raster/raster_factory.hpp @@ -0,0 +1,33 @@ +#pragma once +#include "spatial/common.hpp" + +class GDALDataset; + +namespace spatial { + +namespace gdal { + +//! A loader of GDALDataset from several data source types. +//! Does not take ownership of the pointer. +class RasterFactory { +public: + //! Given a file path, returns a GDALDataset + static GDALDataset *FromFile(const std::string &file_path, + const std::vector &allowed_drivers = std::vector(), + const std::vector &open_options = std::vector(), + const std::vector &sibling_files = std::vector()); + + //! Writes a GDALDataset to a file path + static bool WriteFile(GDALDataset *dataset, const std::string &file_path, const std::string &driver_name = "COG", + const std::vector &write_options = std::vector()); + + //! Transforms a vector of strings as a vector of const char pointers. + static std::vector FromVectorOfStrings(const std::vector &input); + //! Transforms a map of params as a vector of const char pointers. + static std::vector FromNamedParameters(const named_parameter_map_t &input, + const std::string &keyname); +}; + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/include/spatial/gdal/raster/raster_registry.hpp b/spatial/include/spatial/gdal/raster/raster_registry.hpp new file mode 100644 index 00000000..37ba9829 --- /dev/null +++ b/spatial/include/spatial/gdal/raster/raster_registry.hpp @@ -0,0 +1,28 @@ +#pragma once +#include "spatial/common.hpp" + +#include "gdal_priv.h" + +namespace spatial { + +namespace gdal { + +//! A registry of Rasters (GDALDataset) where items are released. +//! This takes ownership of items registered. +class RasterRegistry { +public: + //! Constructor + RasterRegistry(); + //! Destructor + ~RasterRegistry(); + + //! Register a raster + void RegisterRaster(GDALDataset *dataset); + +private: + std::vector datasets_; +}; + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/include/spatial/gdal/raster/raster_value.hpp b/spatial/include/spatial/gdal/raster/raster_value.hpp new file mode 100644 index 00000000..95710dd1 --- /dev/null +++ b/spatial/include/spatial/gdal/raster/raster_value.hpp @@ -0,0 +1,24 @@ +#pragma once +#include "spatial/common.hpp" + +class GDALDataset; + +namespace spatial { + +namespace gdal { + +//! This Value object holds a Raster instance +class RasterValue : public Value { +public: + //! Returns the pointer to the dataset + GDALDataset *operator->() const; + //! Returns the pointer to the dataset + GDALDataset *get() const; + + //! Create a RASTER value + static Value CreateValue(GDALDataset *dataset); +}; + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/include/spatial/gdal/types.hpp b/spatial/include/spatial/gdal/types.hpp new file mode 100644 index 00000000..3ab41fb0 --- /dev/null +++ b/spatial/include/spatial/gdal/types.hpp @@ -0,0 +1,86 @@ +#pragma once +#include "spatial/common.hpp" + +namespace spatial { + +namespace gdal { + +//! Supported Pixel data types (GDALDataType). +typedef enum { + Unknown = 0, /**< Unknown or unspecified type */ + Byte = 1, /**< Eight bit unsigned integer */ + Int8 = 14, /**< 8-bit signed integer */ + UInt16 = 2, /**< Sixteen bit unsigned integer */ + Int16 = 3, /**< Sixteen bit signed integer */ + UInt32 = 4, /**< Thirty two bit unsigned integer */ + Int32 = 5, /**< Thirty two bit signed integer */ + UInt64 = 12, /**< 64 bit unsigned integer */ + Int64 = 13, /**< 64 bit signed integer */ + Float32 = 6, /**< Thirty two bit floating point */ + Float64 = 7, /**< Sixty four bit floating point */ + CInt16 = 8, /**< Complex Int16 */ + CInt32 = 9, /**< Complex Int32 */ + CFloat32 = 10, /**< Complex Float32 */ + CFloat64 = 11, /**< Complex Float64 */ + TypeCount = 15 /**< maximum type # + 1 */ +} PixelType; + +//! Returns the name of given PixelType +std::string GetPixelTypeName(const PixelType &pixel_type); + +//! Supported Types of color interpretation for raster bands (GDALColorInterp). +typedef enum { + Undefined = 0, /**< Undefined */ + GrayIndex = 1, /**< Greyscale */ + PaletteIndex = 2, /**< Paletted (see associated color table) */ + RedBand = 3, /**< Red band of RGBA image */ + GreenBand = 4, /**< Green band of RGBA image */ + BlueBand = 5, /**< Blue band of RGBA image */ + AlphaBand = 6, /**< Alpha (0=transparent, 255=opaque) */ + HueBand = 7, /**< Hue band of HLS image */ + SaturationBand = 8, /**< Saturation band of HLS image */ + LightnessBand = 9, /**< Lightness band of HLS image */ + CyanBand = 10, /**< Cyan band of CMYK image */ + MagentaBand = 11, /**< Magenta band of CMYK image */ + YellowBand = 12, /**< Yellow band of CMYK image */ + BlackBand = 13, /**< Black band of CMYK image */ + YCbCr_YBand = 14, /**< Y Luminance */ + YCbCr_CbBand = 15, /**< Cb Chroma */ + YCbCr_CrBand = 16 /**< Cr Chroma */ +} ColorInterp; + +//! Returns the name of given ColorInterp +std::string GetColorInterpName(const ColorInterp &color_interp); + +//! Supported Warp Resampling Algorithm (GDALResampleAlg). +typedef enum { + NearestNeighbour = 0, /**< Nearest neighbour (select on one input pixel) */ + Bilinear = 1, /**< Bilinear (2x2 kernel) */ + Cubic = 2, /**< Cubic Convolution Approximation (4x4 kernel) */ + CubicSpline = 3, /**< Cubic B-Spline Approximation (4x4 kernel) */ + Lanczos = 4, /**< Lanczos windowed sinc interpolation (6x6 kernel) */ + Average = 5, /**< Average (computes the weighted average of all non-NODATA contributing pixels) */ + Mode = 6, /**< Mode (selects the value which appears most often of all the sampled points) */ + Max = 8, /**< Max (selects maximum of all non-NODATA contributing pixels) */ + Min = 9, /**< Min (selects minimum of all non-NODATA contributing pixels) */ + Med = 10, /**< Med (selects median of all non-NODATA contributing pixels) */ + Q1 = 11, /**< Q1 (selects first quartile of all non-NODATA contributing pixels) */ + Q3 = 12, /**< Q3 (selects third quartile of all non-NODATA contributing pixels) */ + Sum = 13, /**< Sum (weighed sum of all non-NODATA contributing pixels) */ + RMS = 14, /**< RMS (weighted root mean square (quadratic mean) of all non-NODATA contributing pixels) */ +} ResampleAlg; + +//! Returns the name of given ResampleAlg +std::string GetResampleAlgName(const ResampleAlg &resample_alg); + +//! Position of a cell in a Raster (upper left corner as column and row) +struct RasterCoord { + int32_t col; + int32_t row; + explicit RasterCoord(int32_t col, int32_t row) : col(col), row(row) { + } +}; + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/core/types.cpp b/spatial/src/spatial/core/types.cpp index d33be0b9..665d02c4 100644 --- a/spatial/src/spatial/core/types.cpp +++ b/spatial/src/spatial/core/types.cpp @@ -64,6 +64,18 @@ LogicalType GeoTypes::WKB_BLOB() { return blob_type; } +LogicalType GeoTypes::RASTER() { + auto type = LogicalType(LogicalTypeId::POINTER); + type.SetAlias("RASTER"); + return type; +} + +LogicalType GeoTypes::RASTER_COORD() { + auto type = LogicalType::STRUCT({{"col", LogicalType::INTEGER}, {"row", LogicalType::INTEGER}}); + type.SetAlias("RASTER_COORD"); + return type; +} + LogicalType GeoTypes::CreateEnumType(const string &name, const vector &members) { auto varchar_vector = Vector(LogicalType::VARCHAR, members.size()); auto varchar_data = FlatVector::GetData(varchar_vector); @@ -101,6 +113,12 @@ void GeoTypes::Register(DatabaseInstance &db) { // WKB_BLOB ExtensionUtil::RegisterType(db, "WKB_BLOB", GeoTypes::WKB_BLOB()); + + // RASTER + ExtensionUtil::RegisterType(db, "RASTER", GeoTypes::RASTER()); + + // RASTER_COORD + ExtensionUtil::RegisterType(db, "RASTER_COORD", GeoTypes::RASTER_COORD()); } } // namespace core diff --git a/spatial/src/spatial/gdal/CMakeLists.txt b/spatial/src/spatial/gdal/CMakeLists.txt index 8dce5a4f..9b8ab307 100644 --- a/spatial/src/spatial/gdal/CMakeLists.txt +++ b/spatial/src/spatial/gdal/CMakeLists.txt @@ -1,7 +1,9 @@ add_subdirectory(functions) +add_subdirectory(raster) set(EXTENSION_SOURCES ${EXTENSION_SOURCES} ${CMAKE_CURRENT_SOURCE_DIR}/module.cpp ${CMAKE_CURRENT_SOURCE_DIR}/file_handler.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/types.cpp PARENT_SCOPE ) \ No newline at end of file diff --git a/spatial/src/spatial/gdal/file_handler.cpp b/spatial/src/spatial/gdal/file_handler.cpp index cb50144d..19f032c5 100644 --- a/spatial/src/spatial/gdal/file_handler.cpp +++ b/spatial/src/spatial/gdal/file_handler.cpp @@ -1,4 +1,5 @@ #include "spatial/gdal/file_handler.hpp" +#include "spatial/gdal/raster/raster_registry.hpp" #include "duckdb/common/mutex.hpp" #include "duckdb/main/client_context.hpp" @@ -376,6 +377,18 @@ void GDALClientContextState::QueryEnd() { }; +void GDALClientContextState::QueryEnd(ClientContext &context) { + std::lock_guard glock(lock); + auto keyid = CastPointerToValue(&context.GetExecutor()); + + auto entry = registries.find(keyid); + if (entry != registries.end()) { + std::unique_ptr ®istry = entry->second; + registry.reset(); + registries.erase(entry); + } +} + string GDALClientContextState::GetPrefix(const string &value) const { // If the user explicitly asked for a VSI prefix, we don't add our own if (StringUtil::StartsWith(value, "/vsi")) { @@ -384,6 +397,21 @@ string GDALClientContextState::GetPrefix(const string &value) const { return client_prefix + value; } +RasterRegistry &GDALClientContextState::GetRasterRegistry(ClientContext &context) { + std::lock_guard glock(lock); + auto keyid = CastPointerToValue(&context.GetExecutor()); + + auto entry = registries.find(keyid); + if (entry != registries.end()) { + std::unique_ptr ®istry = entry->second; + return *registry.get(); + } else { + std::unique_ptr entry = make_uniq(); + registries[keyid] = std::move(entry); + return *registries[keyid].get(); + } +} + GDALClientContextState &GDALClientContextState::GetOrCreate(ClientContext &context) { if (!context.registered_state["gdal"]) { context.registered_state["gdal"] = make_uniq(context); diff --git a/spatial/src/spatial/gdal/functions/CMakeLists.txt b/spatial/src/spatial/gdal/functions/CMakeLists.txt index 72a0eaab..2330e382 100644 --- a/spatial/src/spatial/gdal/functions/CMakeLists.txt +++ b/spatial/src/spatial/gdal/functions/CMakeLists.txt @@ -1,8 +1,14 @@ +add_subdirectory(aggregate) +add_subdirectory(scalar) set(EXTENSION_SOURCES ${EXTENSION_SOURCES} + ${CMAKE_CURRENT_SOURCE_DIR}/cast.cpp ${CMAKE_CURRENT_SOURCE_DIR}/st_drivers.cpp ${CMAKE_CURRENT_SOURCE_DIR}/st_read.cpp ${CMAKE_CURRENT_SOURCE_DIR}/st_read_meta.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/st_read_raster.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/st_read_raster_meta.cpp ${CMAKE_CURRENT_SOURCE_DIR}/st_write.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/st_write_raster.cpp PARENT_SCOPE ) \ No newline at end of file diff --git a/spatial/src/spatial/gdal/functions/aggregate/CMakeLists.txt b/spatial/src/spatial/gdal/functions/aggregate/CMakeLists.txt new file mode 100644 index 00000000..0b85406c --- /dev/null +++ b/spatial/src/spatial/gdal/functions/aggregate/CMakeLists.txt @@ -0,0 +1,7 @@ +set(EXTENSION_SOURCES + ${EXTENSION_SOURCES} + ${CMAKE_CURRENT_SOURCE_DIR}/raster_agg.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/st_mosaic_agg.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/st_union_agg.cpp + PARENT_SCOPE +) diff --git a/spatial/src/spatial/gdal/functions/aggregate/raster_agg.cpp b/spatial/src/spatial/gdal/functions/aggregate/raster_agg.cpp new file mode 100644 index 00000000..686e09e9 --- /dev/null +++ b/spatial/src/spatial/gdal/functions/aggregate/raster_agg.cpp @@ -0,0 +1,44 @@ +#include "spatial/gdal/functions/raster_agg.hpp" + +#include "gdal_priv.h" + +namespace spatial { + +namespace gdal { + +unique_ptr BindRasterAggOperation(ClientContext &context, AggregateFunction &function, + vector> &arguments) { + + std::vector options; + + // First arg is a raster, next ones are optional arguments + for (std::size_t i = 1; i < arguments.size(); i++) { + auto &arg = arguments[i]; + + if (arg->HasParameter()) { + throw ParameterNotResolvedException(); + } + if (!arg->IsFoldable()) { + throw BinderException("raster_agg: arguments must be constant"); + } + if (arg->alias == "options" || arg->return_type.id() == LogicalTypeId::LIST) { + + Value param_list = ExpressionExecutor::EvaluateScalar(context, *arg); + auto ¶ms = ListValue::GetChildren(param_list); + + for (std::size_t j = 0; j < params.size(); j++) { + const auto ¶m_value = params[j]; + const auto option = param_value.ToString(); + options.push_back(option); + } + } else { + throw BinderException(StringUtil::Format("raster_agg: Unknown argument '%s'", arg->alias.c_str())); + } + } + + return make_uniq(context, options); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/functions/aggregate/st_mosaic_agg.cpp b/spatial/src/spatial/gdal/functions/aggregate/st_mosaic_agg.cpp new file mode 100644 index 00000000..d8a58f98 --- /dev/null +++ b/spatial/src/spatial/gdal/functions/aggregate/st_mosaic_agg.cpp @@ -0,0 +1,123 @@ +#include "duckdb/parser/parsed_data/create_scalar_function_info.hpp" +#include "duckdb/parser/parsed_data/create_aggregate_function_info.hpp" + +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/gdal/file_handler.hpp" +#include "spatial/gdal/functions/aggregate.hpp" +#include "spatial/gdal/functions/raster_agg.hpp" +#include "spatial/gdal/raster/raster.hpp" + +#include "gdal_priv.h" + +using namespace spatial::core; + +namespace spatial { + +namespace gdal { + +//------------------------------------------------------------------------ +// ST_RasterMosaic_Agg +//------------------------------------------------------------------------ + +template +static void RasterMosaicFunction(STATE &state, T &target, AggregateFinalizeData &finalize_data) { + if (!state.is_set) { + finalize_data.ReturnNull(); + } else { + auto datasets = state.datasets; + auto &bind_data = finalize_data.input.bind_data->Cast(); + auto &context = bind_data.context; + auto &options = bind_data.options; + + GDALDataset *result = Raster::BuildVRT(datasets, options); + + if (result == nullptr) { + auto error = Raster::GetLastErrorMsg(); + throw IOException("Could not make mosaic: (" + error + ")"); + } + + auto &ctx_state = GDALClientContextState::GetOrCreate(context); + auto &raster_registry = ctx_state.GetRasterRegistry(context); + raster_registry.RegisterRaster(result); + + target = CastPointerToValue(result); + } +} + +struct MosaicAggUnaryOperation : RasterAggUnaryOperation { + + template + static void Finalize(STATE &state, T &target, AggregateFinalizeData &finalize_data) { + RasterMosaicFunction(state, target, finalize_data); + } +}; + +struct MosaicAggBinaryOperation : RasterAggBinaryOperation { + + template + static void Finalize(STATE &state, T &target, AggregateFinalizeData &finalize_data) { + RasterMosaicFunction(state, target, finalize_data); + } +}; + +//------------------------------------------------------------------------ +// Documentation +//------------------------------------------------------------------------ + +static constexpr const char *DOC_DESCRIPTION = R"( + Returns a mosaic of a set of raster tiles into a single raster. + + Tiles are considered as source rasters of a larger mosaic and the result dataset has as many bands as one of the input files. + + `options` is optional, an array of parameters like [GDALBuildVRT](https://gdal.org/programs/gdalbuildvrt.html). +)"; + +static constexpr const char *DOC_EXAMPLE = R"( + WITH __input AS ( + SELECT + 1 AS raster_id, + ST_RasterFromFile(file) AS raster + FROM + glob('./test/data/mosaic/*.tiff') + ), + SELECT + ST_RasterMosaic_Agg(raster, options => ['-r', 'bilinear']) AS r + FROM + __input + GROUP BY + raster_id + ; +)"; + +static constexpr DocTag DOC_TAGS[] = {{"ext", "spatial"}, {"category", "aggregation"}}; + +//------------------------------------------------------------------------ +// Register +//------------------------------------------------------------------------ + +void GdalAggregateFunctions::RegisterStRasterMosaicAgg(DatabaseInstance &db) { + + AggregateFunctionSet st_mosaic_agg("ST_RasterMosaic_Agg"); + + auto fun01 = AggregateFunction::UnaryAggregate( + GeoTypes::RASTER(), GeoTypes::RASTER()); + fun01.bind = BindRasterAggOperation; + + st_mosaic_agg.AddFunction(fun01); + + auto fun02 = AggregateFunction::BinaryAggregate( + GeoTypes::RASTER(), LogicalType::LIST(LogicalType::VARCHAR), GeoTypes::RASTER()); + fun02.bind = BindRasterAggOperation; + + st_mosaic_agg.AddFunction(fun02); + + ExtensionUtil::RegisterFunction(db, st_mosaic_agg); + + DocUtil::AddDocumentation(db, "ST_RasterMosaic_Agg", DOC_DESCRIPTION, DOC_EXAMPLE, DOC_TAGS); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/functions/aggregate/st_union_agg.cpp b/spatial/src/spatial/gdal/functions/aggregate/st_union_agg.cpp new file mode 100644 index 00000000..750ac132 --- /dev/null +++ b/spatial/src/spatial/gdal/functions/aggregate/st_union_agg.cpp @@ -0,0 +1,127 @@ +#include "duckdb/parser/parsed_data/create_scalar_function_info.hpp" +#include "duckdb/parser/parsed_data/create_aggregate_function_info.hpp" + +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/gdal/file_handler.hpp" +#include "spatial/gdal/functions/aggregate.hpp" +#include "spatial/gdal/functions/raster_agg.hpp" +#include "spatial/gdal/raster/raster.hpp" + +#include "gdal_priv.h" + +using namespace spatial::core; + +namespace spatial { + +namespace gdal { + +//------------------------------------------------------------------------ +// ST_RasterUnion_Agg +//------------------------------------------------------------------------ + +template +static void RasterUnionFunction(STATE &state, T &target, AggregateFinalizeData &finalize_data) { + if (!state.is_set) { + finalize_data.ReturnNull(); + } else { + auto datasets = state.datasets; + auto &bind_data = finalize_data.input.bind_data->Cast(); + auto &context = bind_data.context; + auto &options = bind_data.options; + + std::vector vrt_options; + vrt_options.push_back("-separate"); + vrt_options.insert(vrt_options.end(), options.begin(), options.end()); + + GDALDataset *result = Raster::BuildVRT(datasets, vrt_options); + + if (result == nullptr) { + auto error = Raster::GetLastErrorMsg(); + throw IOException("Could not make union: (" + error + ")"); + } + + auto &ctx_state = GDALClientContextState::GetOrCreate(context); + auto &raster_registry = ctx_state.GetRasterRegistry(context); + raster_registry.RegisterRaster(result); + + target = CastPointerToValue(result); + } +} + +struct UnionAggUnaryOperation : RasterAggUnaryOperation { + + template + static void Finalize(STATE &state, T &target, AggregateFinalizeData &finalize_data) { + RasterUnionFunction(state, target, finalize_data); + } +}; + +struct UnionAggBinaryOperation : RasterAggBinaryOperation { + + template + static void Finalize(STATE &state, T &target, AggregateFinalizeData &finalize_data) { + RasterUnionFunction(state, target, finalize_data); + } +}; + +//------------------------------------------------------------------------ +// Documentation +//------------------------------------------------------------------------ + +static constexpr const char *DOC_DESCRIPTION = R"( + Returns the union of a set of raster tiles into a single raster composed of at least one band. + + Each tiles goes into a separate band in the result dataset. + + `options` is optional, an array of parameters like [GDALBuildVRT](https://gdal.org/programs/gdalbuildvrt.html). +)"; + +static constexpr const char *DOC_EXAMPLE = R"( + WITH __input AS ( + SELECT + 1 AS raster_id, + ST_RasterFromFile(file) AS raster + FROM + glob('./test/data/bands/*.tiff') + ), + SELECT + ST_RasterUnion_Agg(raster, options => ['-resolution', 'highest']) AS r + FROM + __input + GROUP BY + raster_id + ; +)"; + +static constexpr DocTag DOC_TAGS[] = {{"ext", "spatial"}, {"category", "aggregation"}}; + +//------------------------------------------------------------------------ +// Register +//------------------------------------------------------------------------ + +void GdalAggregateFunctions::RegisterStRasterUnionAgg(DatabaseInstance &db) { + + AggregateFunctionSet st_union_agg("ST_RasterUnion_Agg"); + + auto fun01 = AggregateFunction::UnaryAggregate( + GeoTypes::RASTER(), GeoTypes::RASTER()); + fun01.bind = BindRasterAggOperation; + + st_union_agg.AddFunction(fun01); + + auto fun02 = + AggregateFunction::BinaryAggregate( + GeoTypes::RASTER(), LogicalType::LIST(LogicalType::VARCHAR), GeoTypes::RASTER()); + fun02.bind = BindRasterAggOperation; + + st_union_agg.AddFunction(fun02); + + ExtensionUtil::RegisterFunction(db, st_union_agg); + + DocUtil::AddDocumentation(db, "ST_RasterUnion_Agg", DOC_DESCRIPTION, DOC_EXAMPLE, DOC_TAGS); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/functions/cast.cpp b/spatial/src/spatial/gdal/functions/cast.cpp new file mode 100644 index 00000000..890df543 --- /dev/null +++ b/spatial/src/spatial/gdal/functions/cast.cpp @@ -0,0 +1,82 @@ +#include "duckdb/common/vector_operations/generic_executor.hpp" +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/core/functions/common.hpp" +#include "spatial/core/geometry/geometry.hpp" +#include "spatial/gdal/functions/cast.hpp" +#include "spatial/gdal/raster/raster.hpp" + +#include "gdal_priv.h" + +using namespace spatial::core; + +namespace spatial { + +namespace gdal { + +//------------------------------------------------------------------------------ +// RASTER -> VARCHAR +//------------------------------------------------------------------------------ + +static bool RasterToVarcharCast(Vector &source, Vector &result, idx_t count, CastParameters ¶meters) { + UnaryExecutor::Execute(source, result, count, + [&](uintptr_t &input) { return string_t("RASTER"); }); + return true; +} + +//------------------------------------------------------------------------------ +// RASTER -> GEOMETRY +//------------------------------------------------------------------------------ + +static bool RasterToGeometryCast(Vector &source, Vector &result, idx_t count, CastParameters ¶meters) { + auto &lstate = GeometryFunctionLocalState::ResetAndGet(parameters); + + UnaryExecutor::Execute(source, result, count, [&](uintptr_t &input) { + Raster raster(reinterpret_cast(input)); + return Geometry::Serialize(raster.GetGeometry(lstate.arena), result); + }); + return true; +} + +//------------------------------------------------------------------------------ +// RASTER_COORD -> VARCHAR +//------------------------------------------------------------------------------ + +static bool RasterCoordToVarcharCast(Vector &source, Vector &result, idx_t count, CastParameters ¶meters) { + + using COORD_TYPE = StructTypeBinary; + using VARCHAR_TYPE = PrimitiveType; + + GenericExecutor::ExecuteUnary(source, result, count, [&](COORD_TYPE &input) { + auto col = input.a_val; + auto row = input.b_val; + return StringVector::AddString(result, StringUtil::Format("COORD (%d, %d)", col, row)); + }); + return true; +} + +//------------------------------------------------------------------------------ +// Register +//------------------------------------------------------------------------------ + +void GdalCastFunctions::Register(DatabaseInstance &db) { + + ExtensionUtil::RegisterCastFunction(db, GeoTypes::RASTER(), LogicalType::VARCHAR, RasterToVarcharCast, 1); + + ExtensionUtil::RegisterCastFunction( + db, GeoTypes::RASTER(), GeoTypes::GEOMETRY(), + BoundCastInfo(RasterToGeometryCast, nullptr, GeometryFunctionLocalState::InitCast), 1); + + // POINTER -> RASTER is implicitly castable + ExtensionUtil::RegisterCastFunction(db, LogicalType::POINTER, GeoTypes::RASTER(), DefaultCasts::ReinterpretCast, 1); + + // RASTER -> POINTER is implicitly castable + ExtensionUtil::RegisterCastFunction(db, GeoTypes::RASTER(), LogicalType::POINTER, DefaultCasts::ReinterpretCast, 1); + + ExtensionUtil::RegisterCastFunction(db, GeoTypes::RASTER_COORD(), LogicalType::VARCHAR, RasterCoordToVarcharCast, + 1); +}; + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/functions/scalar/CMakeLists.txt b/spatial/src/spatial/gdal/functions/scalar/CMakeLists.txt new file mode 100644 index 00000000..43771e73 --- /dev/null +++ b/spatial/src/spatial/gdal/functions/scalar/CMakeLists.txt @@ -0,0 +1,18 @@ +set(EXTENSION_SOURCES + ${EXTENSION_SOURCES} + ${CMAKE_CURRENT_SOURCE_DIR}/st_band_color_interp.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/st_band_nodata_value.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/st_band_pixel_type.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/st_has_no_band.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/st_metadata.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/st_raster_as_file.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/st_raster_clip.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/st_raster_from_file.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/st_raster_geometry.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/st_raster_to_world_coord.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/st_raster_warp.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/st_srid.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/st_value.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/st_world_to_raster_coord.cpp + PARENT_SCOPE +) diff --git a/spatial/src/spatial/gdal/functions/scalar/st_band_color_interp.cpp b/spatial/src/spatial/gdal/functions/scalar/st_band_color_interp.cpp new file mode 100644 index 00000000..b981e87f --- /dev/null +++ b/spatial/src/spatial/gdal/functions/scalar/st_band_color_interp.cpp @@ -0,0 +1,150 @@ +#include "duckdb/parser/parsed_data/create_scalar_function_info.hpp" +#include "duckdb/common/vector_operations/generic_executor.hpp" +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/gdal/types.hpp" +#include "spatial/gdal/functions/scalar.hpp" + +#include "gdal_priv.h" + +using namespace spatial::core; + +namespace spatial { + +namespace gdal { + +//------------------------------------------------------------------------------ +// ST_GetBandColorInterp +//------------------------------------------------------------------------------ + +static void RasterGetColorInterpFunction(DataChunk &args, ExpressionState &state, Vector &result) { + D_ASSERT(args.data.size() == 2); + + BinaryExecutor::Execute( + args.data[0], args.data[1], result, args.size(), [&](uintptr_t input, int32_t band_num) { + GDALDataset *dataset = reinterpret_cast(input); + + if (band_num < 1) { + throw InvalidInputException("BandNum must be greater than 0"); + } + if (dataset->GetRasterCount() < band_num) { + throw InvalidInputException("Dataset only has %d RasterBands", dataset->GetRasterCount()); + } + GDALRasterBand *raster_band = dataset->GetRasterBand(band_num); + return (int32_t)raster_band->GetColorInterpretation(); + }); +} + +//------------------------------------------------------------------------------ +// ST_GetBandColorInterpName +//------------------------------------------------------------------------------ + +static void RasterGetColorInterpNameFunction(DataChunk &args, ExpressionState &state, Vector &result) { + D_ASSERT(args.data.size() == 2); + + BinaryExecutor::Execute( + args.data[0], args.data[1], result, args.size(), [&](uintptr_t input, int32_t band_num) { + GDALDataset *dataset = reinterpret_cast(input); + + if (band_num < 1) { + throw InvalidInputException("BandNum must be greater than 0"); + } + if (dataset->GetRasterCount() < band_num) { + throw InvalidInputException("Dataset only has %d RasterBands", dataset->GetRasterCount()); + } + GDALRasterBand *raster_band = dataset->GetRasterBand(band_num); + return GetColorInterpName((ColorInterp)raster_band->GetColorInterpretation()); + }); +} + +//------------------------------------------------------------------------ +// Documentation +//------------------------------------------------------------------------ + +static constexpr const char *DOC_DESCRIPTION_1 = R"( + Returns the color interpretation of a band in the raster. + + This is a code in the enumeration: + + + Undefined = 0: Undefined + + GrayIndex = 1: Greyscale + + PaletteIndex = 2: Paletted (see associated color table) + + RedBand = 3: Red band of RGBA image + + GreenBand = 4: Green band of RGBA image + + BlueBand = 5: Blue band of RGBA image + + AlphaBand = 6: Alpha (0=transparent, 255=opaque) + + HueBand = 7: Hue band of HLS image + + SaturationBand = 8: Saturation band of HLS image + + LightnessBand = 9: Lightness band of HLS image + + CyanBand = 10: Cyan band of CMYK image + + MagentaBand = 11: Magenta band of CMYK image + + YellowBand = 12: Yellow band of CMYK image + + BlackBand = 13: Black band of CMYK image + + YCbCr_YBand = 14: Y Luminance + + YCbCr_CbBand = 15: Cb Chroma + + YCbCr_CrBand = 16: Cr Chroma +)"; + +static constexpr const char *DOC_EXAMPLE_1 = R"( + SELECT ST_GetBandColorInterp(raster, 1) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_1[] = {{"ext", "spatial"}, {"category", "property"}}; + +static constexpr const char *DOC_DESCRIPTION_2 = R"( + Returns the color interpretation name of a band in the raster. + + This is a string in the enumeration: + + + Undefined: Undefined + + Greyscale: Greyscale + + Paletted: Paletted (see associated color table) + + Red: Red band of RGBA image + + Green: Green band of RGBA image + + Blue: Blue band of RGBA image + + Alpha: Alpha (0=transparent, 255=opaque) + + Hue: Hue band of HLS image + + Saturation: Saturation band of HLS image + + Lightness: Lightness band of HLS image + + Cyan: Cyan band of CMYK image + + Magenta: Magenta band of CMYK image + + Yellow: Yellow band of CMYK image + + Black: Black band of CMYK image + + YLuminance: Y Luminance + + CbChroma: Cb Chroma + + CrChroma: Cr Chroma +)"; + +static constexpr const char *DOC_EXAMPLE_2 = R"( + SELECT ST_GetBandColorInterpName(raster, 1) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_2[] = {{"ext", "spatial"}, {"category", "property"}}; + +//------------------------------------------------------------------------------ +// Register functions +//------------------------------------------------------------------------------ + +void GdalScalarFunctions::RegisterStBandColorInterp(DatabaseInstance &db) { + + ScalarFunctionSet set("ST_GetBandColorInterp"); + set.AddFunction( + ScalarFunction({GeoTypes::RASTER(), LogicalType::INTEGER}, LogicalType::INTEGER, RasterGetColorInterpFunction)); + ExtensionUtil::RegisterFunction(db, set); + + DocUtil::AddDocumentation(db, "ST_GetBandColorInterp", DOC_DESCRIPTION_1, DOC_EXAMPLE_1, DOC_TAGS_1); +} + +void GdalScalarFunctions::RegisterStBandColorInterpName(DatabaseInstance &db) { + + ScalarFunctionSet set("ST_GetBandColorInterpName"); + set.AddFunction(ScalarFunction({GeoTypes::RASTER(), LogicalType::INTEGER}, LogicalType::VARCHAR, + RasterGetColorInterpNameFunction)); + ExtensionUtil::RegisterFunction(db, set); + + DocUtil::AddDocumentation(db, "ST_GetBandColorInterpName", DOC_DESCRIPTION_2, DOC_EXAMPLE_2, DOC_TAGS_2); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/functions/scalar/st_band_nodata_value.cpp b/spatial/src/spatial/gdal/functions/scalar/st_band_nodata_value.cpp new file mode 100644 index 00000000..648cac98 --- /dev/null +++ b/spatial/src/spatial/gdal/functions/scalar/st_band_nodata_value.cpp @@ -0,0 +1,68 @@ +#include "duckdb/parser/parsed_data/create_scalar_function_info.hpp" +#include "duckdb/common/vector_operations/generic_executor.hpp" +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/gdal/types.hpp" +#include "spatial/gdal/functions/scalar.hpp" + +#include "gdal_priv.h" + +using namespace spatial::core; + +namespace spatial { + +namespace gdal { + +//------------------------------------------------------------------------------ +// ST_GetBandNoDataValue +//------------------------------------------------------------------------------ + +static void RasterGetBandNoDataFunction(DataChunk &args, ExpressionState &state, Vector &result) { + D_ASSERT(args.data.size() == 2); + + BinaryExecutor::Execute( + args.data[0], args.data[1], result, args.size(), [&](uintptr_t input, int32_t band_num) { + GDALDataset *dataset = reinterpret_cast(input); + + if (band_num < 1) { + throw InvalidInputException("BandNum must be greater than 0"); + } + if (dataset->GetRasterCount() < band_num) { + throw InvalidInputException("Dataset only has %d RasterBands", dataset->GetRasterCount()); + } + GDALRasterBand *raster_band = dataset->GetRasterBand(band_num); + return raster_band->GetNoDataValue(); + }); +} + +//------------------------------------------------------------------------ +// Documentation +//------------------------------------------------------------------------ + +static constexpr const char *DOC_DESCRIPTION = R"( + Returns the NODATA value of a band in the raster. +)"; + +static constexpr const char *DOC_EXAMPLE = R"( + SELECT ST_GetBandNoDataValue(raster, 1) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS[] = {{"ext", "spatial"}, {"category", "property"}}; + +//------------------------------------------------------------------------------ +// Register functions +//------------------------------------------------------------------------------ + +void GdalScalarFunctions::RegisterStBandNoDataValue(DatabaseInstance &db) { + + ScalarFunctionSet set("ST_GetBandNoDataValue"); + set.AddFunction( + ScalarFunction({GeoTypes::RASTER(), LogicalType::INTEGER}, LogicalType::DOUBLE, RasterGetBandNoDataFunction)); + ExtensionUtil::RegisterFunction(db, set); + + DocUtil::AddDocumentation(db, "ST_GetBandNoDataValue", DOC_DESCRIPTION, DOC_EXAMPLE, DOC_TAGS); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/functions/scalar/st_band_pixel_type.cpp b/spatial/src/spatial/gdal/functions/scalar/st_band_pixel_type.cpp new file mode 100644 index 00000000..fe008505 --- /dev/null +++ b/spatial/src/spatial/gdal/functions/scalar/st_band_pixel_type.cpp @@ -0,0 +1,146 @@ +#include "duckdb/parser/parsed_data/create_scalar_function_info.hpp" +#include "duckdb/common/vector_operations/generic_executor.hpp" +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/gdal/types.hpp" +#include "spatial/gdal/functions/scalar.hpp" + +#include "gdal_priv.h" + +using namespace spatial::core; + +namespace spatial { + +namespace gdal { + +//------------------------------------------------------------------------------ +// ST_GetBandPixelType +//------------------------------------------------------------------------------ + +static void RasterGetPixelTypeFunction(DataChunk &args, ExpressionState &state, Vector &result) { + D_ASSERT(args.data.size() == 2); + + BinaryExecutor::Execute( + args.data[0], args.data[1], result, args.size(), [&](uintptr_t input, int32_t band_num) { + GDALDataset *dataset = reinterpret_cast(input); + + if (band_num < 1) { + throw InvalidInputException("BandNum must be greater than 0"); + } + if (dataset->GetRasterCount() < band_num) { + throw InvalidInputException("Dataset only has %d RasterBands", dataset->GetRasterCount()); + } + GDALRasterBand *raster_band = dataset->GetRasterBand(band_num); + return (int32_t)raster_band->GetRasterDataType(); + }); +} + +//------------------------------------------------------------------------------ +// ST_GetBandPixelTypeName +//------------------------------------------------------------------------------ + +static void RasterGetPixelTypeNameFunction(DataChunk &args, ExpressionState &state, Vector &result) { + D_ASSERT(args.data.size() == 2); + + BinaryExecutor::Execute( + args.data[0], args.data[1], result, args.size(), [&](uintptr_t input, int32_t band_num) { + GDALDataset *dataset = reinterpret_cast(input); + + if (band_num < 1) { + throw InvalidInputException("BandNum must be greater than 0"); + } + if (dataset->GetRasterCount() < band_num) { + throw InvalidInputException("Dataset only has %d RasterBands", dataset->GetRasterCount()); + } + GDALRasterBand *raster_band = dataset->GetRasterBand(band_num); + return GetPixelTypeName((PixelType)raster_band->GetRasterDataType()); + }); +} + +//------------------------------------------------------------------------ +// Documentation +//------------------------------------------------------------------------ + +static constexpr const char *DOC_DESCRIPTION_1 = R"( + Returns the pixel type of a band in the raster. + + This is a code in the enumeration: + + + Unknown = 0: Unknown or unspecified type + + Byte = 1: Eight bit unsigned integer + + Int8 = 14: 8-bit signed integer + + UInt16 = 2: Sixteen bit unsigned integer + + Int16 = 3: Sixteen bit signed integer + + UInt32 = 4: Thirty two bit unsigned integer + + Int32 = 5: Thirty two bit signed integer + + UInt64 = 12: 64 bit unsigned integer + + Int64 = 13: 64 bit signed integer + + Float32 = 6: Thirty two bit floating point + + Float64 = 7: Sixty four bit floating point + + CInt16 = 8: Complex Int16 + + CInt32 = 9: Complex Int32 + + CFloat32 = 10: Complex Float32 + + CFloat64 = 11: Complex Float64 +)"; + +static constexpr const char *DOC_EXAMPLE_1 = R"( + SELECT ST_GetBandPixelType(raster, 1) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_1[] = {{"ext", "spatial"}, {"category", "property"}}; + +static constexpr const char *DOC_DESCRIPTION_2 = R"( + Returns the pixel type name of a band in the raster. + + This is a string in the enumeration: + + + Unknown: Unknown or unspecified type + + Byte: Eight bit unsigned integer + + Int8: 8-bit signed integer + + UInt16: Sixteen bit unsigned integer + + Int16: Sixteen bit signed integer + + UInt32: Thirty two bit unsigned integer + + Int32: Thirty two bit signed integer + + UInt64: 64 bit unsigned integer + + Int64: 64 bit signed integer + + Float32: Thirty two bit floating point + + Float64: Sixty four bit floating point + + CInt16: Complex Int16 + + CInt32: Complex Int32 + + CFloat32: Complex Float32 + + CFloat64: Complex Float64 +)"; + +static constexpr const char *DOC_EXAMPLE_2 = R"( + SELECT ST_GetBandPixelTypeName(raster, 1) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_2[] = {{"ext", "spatial"}, {"category", "property"}}; + +//------------------------------------------------------------------------------ +// Register functions +//------------------------------------------------------------------------------ + +void GdalScalarFunctions::RegisterStBandPixelType(DatabaseInstance &db) { + + ScalarFunctionSet set("ST_GetBandPixelType"); + set.AddFunction( + ScalarFunction({GeoTypes::RASTER(), LogicalType::INTEGER}, LogicalType::INTEGER, RasterGetPixelTypeFunction)); + ExtensionUtil::RegisterFunction(db, set); + + DocUtil::AddDocumentation(db, "ST_GetBandPixelType", DOC_DESCRIPTION_1, DOC_EXAMPLE_1, DOC_TAGS_1); +} + +void GdalScalarFunctions::RegisterStBandPixelTypeName(DatabaseInstance &db) { + + ScalarFunctionSet set("ST_GetBandPixelTypeName"); + set.AddFunction(ScalarFunction({GeoTypes::RASTER(), LogicalType::INTEGER}, LogicalType::VARCHAR, + RasterGetPixelTypeNameFunction)); + ExtensionUtil::RegisterFunction(db, set); + + DocUtil::AddDocumentation(db, "ST_GetBandPixelTypeName", DOC_DESCRIPTION_2, DOC_EXAMPLE_2, DOC_TAGS_2); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/functions/scalar/st_has_no_band.cpp b/spatial/src/spatial/gdal/functions/scalar/st_has_no_band.cpp new file mode 100644 index 00000000..2bb76eee --- /dev/null +++ b/spatial/src/spatial/gdal/functions/scalar/st_has_no_band.cpp @@ -0,0 +1,60 @@ +#include "duckdb/parser/parsed_data/create_scalar_function_info.hpp" +#include "duckdb/common/vector_operations/generic_executor.hpp" +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/gdal/functions/scalar.hpp" + +#include "gdal_priv.h" + +using namespace spatial::core; + +namespace spatial { + +namespace gdal { + +//------------------------------------------------------------------------------ +// ST_HasNoBand +//------------------------------------------------------------------------------ + +static void RasterHasNoBandFunction(DataChunk &args, ExpressionState &state, Vector &result) { + D_ASSERT(args.data.size() == 2); + + BinaryExecutor::Execute(args.data[0], args.data[1], result, args.size(), + [&](uintptr_t input, int32_t band_num) { + GDALDataset *dataset = reinterpret_cast(input); + return dataset->GetRasterCount() >= band_num; + }); +} + +//------------------------------------------------------------------------ +// Documentation +//------------------------------------------------------------------------ + +static constexpr const char *DOC_DESCRIPTION = R"( + Returns true if there is no band with given band number. + Band numbers start at 1 and band is assumed to be 1 if not specified. +)"; + +static constexpr const char *DOC_EXAMPLE = R"( + SELECT ST_HasNoBand(raster, 1) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS[] = {{"ext", "spatial"}, {"category", "property"}}; + +//------------------------------------------------------------------------------ +// Register functions +//------------------------------------------------------------------------------ + +void GdalScalarFunctions::RegisterStGetHasNoBand(DatabaseInstance &db) { + + ScalarFunctionSet set("ST_HasNoBand"); + set.AddFunction( + ScalarFunction({GeoTypes::RASTER(), LogicalType::INTEGER}, LogicalType::BOOLEAN, RasterHasNoBandFunction)); + ExtensionUtil::RegisterFunction(db, set); + + DocUtil::AddDocumentation(db, "ST_HasNoBand", DOC_DESCRIPTION, DOC_EXAMPLE, DOC_TAGS); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/functions/scalar/st_metadata.cpp b/spatial/src/spatial/gdal/functions/scalar/st_metadata.cpp new file mode 100644 index 00000000..98359858 --- /dev/null +++ b/spatial/src/spatial/gdal/functions/scalar/st_metadata.cpp @@ -0,0 +1,289 @@ +#include "duckdb/parser/parsed_data/create_scalar_function_info.hpp" +#include "duckdb/common/vector_operations/generic_executor.hpp" +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/gdal/functions/scalar.hpp" +#include "spatial/gdal/raster/raster.hpp" + +#include "gdal_priv.h" + +using namespace spatial::core; + +namespace spatial { + +namespace gdal { + +//------------------------------------------------------------------------------ +// Raster Accessors +//------------------------------------------------------------------------------ + +static void RasterGetWidthFunction(DataChunk &args, ExpressionState &state, Vector &result) { + + UnaryExecutor::Execute(args.data[0], result, args.size(), [&](uintptr_t input) { + GDALDataset *dataset = reinterpret_cast(input); + return dataset->GetRasterXSize(); + }); +} + +static void RasterGetHeightFunction(DataChunk &args, ExpressionState &state, Vector &result) { + + UnaryExecutor::Execute(args.data[0], result, args.size(), [&](uintptr_t input) { + GDALDataset *dataset = reinterpret_cast(input); + return dataset->GetRasterYSize(); + }); +} + +static void RasterGetNumBandsFunction(DataChunk &args, ExpressionState &state, Vector &result) { + + UnaryExecutor::Execute(args.data[0], result, args.size(), [&](uintptr_t input) { + GDALDataset *dataset = reinterpret_cast(input); + return dataset->GetRasterCount(); + }); +} + +static void RasterGetGeoTransformItemFunction(DataChunk &args, ExpressionState &state, Vector &result, + int32_t gt_index) { + + UnaryExecutor::Execute(args.data[0], result, args.size(), [&](uintptr_t input) { + Raster raster(reinterpret_cast(input)); + double gt[6] = {0}; + raster.GetGeoTransform(gt); + return gt[gt_index]; + }); +} + +static void RasterGetUpperLeftX(DataChunk &args, ExpressionState &state, Vector &result) { + RasterGetGeoTransformItemFunction(args, state, result, 0); +} +static void RasterGetUpperLeftY(DataChunk &args, ExpressionState &state, Vector &result) { + RasterGetGeoTransformItemFunction(args, state, result, 3); +} +static void RasterGetScaleX(DataChunk &args, ExpressionState &state, Vector &result) { + RasterGetGeoTransformItemFunction(args, state, result, 1); +} +static void RasterGetScaleY(DataChunk &args, ExpressionState &state, Vector &result) { + RasterGetGeoTransformItemFunction(args, state, result, 5); +} +static void RasterGetSkewX(DataChunk &args, ExpressionState &state, Vector &result) { + RasterGetGeoTransformItemFunction(args, state, result, 2); +} +static void RasterGetSkewY(DataChunk &args, ExpressionState &state, Vector &result) { + RasterGetGeoTransformItemFunction(args, state, result, 4); +} + +static void RasterGetPixelWidthFunction(DataChunk &args, ExpressionState &state, Vector &result) { + + UnaryExecutor::Execute(args.data[0], result, args.size(), [&](uintptr_t input) { + Raster raster(reinterpret_cast(input)); + double gt[6] = {0}; + raster.GetGeoTransform(gt); + return sqrt(gt[1] * gt[1] + gt[4] * gt[4]); + }); +} + +static void RasterGetPixelHeightFunction(DataChunk &args, ExpressionState &state, Vector &result) { + + UnaryExecutor::Execute(args.data[0], result, args.size(), [&](uintptr_t input) { + Raster raster(reinterpret_cast(input)); + double gt[6] = {0}; + raster.GetGeoTransform(gt); + return sqrt(gt[5] * gt[5] + gt[2] * gt[2]); + }); +} + +//------------------------------------------------------------------------ +// Documentation +//------------------------------------------------------------------------ + +static constexpr const char *DOC_DESCRIPTION_01 = R"( + Returns the width of the raster in pixels. +)"; + +static constexpr const char *DOC_EXAMPLE_01 = R"( + SELECT ST_Width(raster) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_01[] = {{"ext", "spatial"}, {"category", "property"}}; + +static constexpr const char *DOC_DESCRIPTION_02 = R"( + Returns the height of the raster in pixels. +)"; + +static constexpr const char *DOC_EXAMPLE_02 = R"( + SELECT ST_Height(raster) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_02[] = {{"ext", "spatial"}, {"category", "property"}}; + +static constexpr const char *DOC_DESCRIPTION_03 = R"( + Returns the number of bands in the raster. +)"; + +static constexpr const char *DOC_EXAMPLE_03 = R"( + SELECT ST_NumBands(raster) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_03[] = {{"ext", "spatial"}, {"category", "property"}}; + +static constexpr const char *DOC_DESCRIPTION_04 = R"( + Returns the upper left X coordinate of raster in projected spatial reference. +)"; + +static constexpr const char *DOC_EXAMPLE_04 = R"( + SELECT ST_UpperLeftX(raster) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_04[] = {{"ext", "spatial"}, {"category", "property"}}; + +static constexpr const char *DOC_DESCRIPTION_05 = R"( + Returns the upper left Y coordinate of raster in projected spatial reference. +)"; + +static constexpr const char *DOC_EXAMPLE_05 = R"( + SELECT ST_UpperLeftY(raster) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_05[] = {{"ext", "spatial"}, {"category", "property"}}; + +static constexpr const char *DOC_DESCRIPTION_06 = R"( + Returns the X component of the pixel width in units of coordinate reference system. + Refer to [World File](https://en.wikipedia.org/wiki/World_file) for more details. +)"; + +static constexpr const char *DOC_EXAMPLE_06 = R"( + SELECT ST_ScaleX(raster) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_06[] = {{"ext", "spatial"}, {"category", "property"}}; + +static constexpr const char *DOC_DESCRIPTION_07 = R"( + Returns the Y component of the pixel width in units of coordinate reference system. + Refer to [World File](https://en.wikipedia.org/wiki/World_file) for more details. +)"; + +static constexpr const char *DOC_EXAMPLE_07 = R"( + SELECT ST_ScaleY(raster) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_07[] = {{"ext", "spatial"}, {"category", "property"}}; + +static constexpr const char *DOC_DESCRIPTION_08 = R"( + Returns the georeference X skew (or rotation parameter). + Refer to [World File](https://en.wikipedia.org/wiki/World_file) for more details. +)"; + +static constexpr const char *DOC_EXAMPLE_08 = R"( + SELECT ST_SkewX(raster) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_08[] = {{"ext", "spatial"}, {"category", "property"}}; + +static constexpr const char *DOC_DESCRIPTION_09 = R"( + Returns the georeference Y skew (or rotation parameter). + Refer to [World File](https://en.wikipedia.org/wiki/World_file) for more details. +)"; + +static constexpr const char *DOC_EXAMPLE_09 = R"( + SELECT ST_SkewY(raster) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_09[] = {{"ext", "spatial"}, {"category", "property"}}; + +static constexpr const char *DOC_DESCRIPTION_10 = R"( + Returns the width of a pixel in geometric units of the spatial reference system. + In the common case where there is no skew, the pixel width is just the scale ratio between geometric coordinates and raster pixels. +)"; + +static constexpr const char *DOC_EXAMPLE_10 = R"( + SELECT ST_PixelWidth(raster) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_10[] = {{"ext", "spatial"}, {"category", "property"}}; + +static constexpr const char *DOC_DESCRIPTION_11 = R"( + Returns the height of a pixel in geometric units of the spatial reference system. + In the common case where there is no skew, the pixel height is just the scale ratio between geometric coordinates and raster pixels. +)"; + +static constexpr const char *DOC_EXAMPLE_11 = R"( + SELECT ST_PixelHeight(raster) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_11[] = {{"ext", "spatial"}, {"category", "property"}}; + +//------------------------------------------------------------------------------ +// Register functions +//------------------------------------------------------------------------------ + +void GdalScalarFunctions::RegisterStGetMetadata(DatabaseInstance &db) { + + ScalarFunctionSet set_01("ST_Width"); + set_01.AddFunction(ScalarFunction({GeoTypes::RASTER()}, LogicalType::INTEGER, RasterGetWidthFunction)); + ExtensionUtil::RegisterFunction(db, set_01); + + DocUtil::AddDocumentation(db, "ST_Width", DOC_DESCRIPTION_01, DOC_EXAMPLE_01, DOC_TAGS_01); + + ScalarFunctionSet set_02("ST_Height"); + set_02.AddFunction(ScalarFunction({GeoTypes::RASTER()}, LogicalType::INTEGER, RasterGetHeightFunction)); + ExtensionUtil::RegisterFunction(db, set_02); + + DocUtil::AddDocumentation(db, "ST_Height", DOC_DESCRIPTION_02, DOC_EXAMPLE_02, DOC_TAGS_02); + + ScalarFunctionSet set_03("ST_NumBands"); + set_03.AddFunction(ScalarFunction({GeoTypes::RASTER()}, LogicalType::INTEGER, RasterGetNumBandsFunction)); + ExtensionUtil::RegisterFunction(db, set_03); + + DocUtil::AddDocumentation(db, "ST_NumBands", DOC_DESCRIPTION_03, DOC_EXAMPLE_03, DOC_TAGS_03); + + ScalarFunctionSet set_04("ST_UpperLeftX"); + set_04.AddFunction(ScalarFunction({GeoTypes::RASTER()}, LogicalType::DOUBLE, RasterGetUpperLeftX)); + ExtensionUtil::RegisterFunction(db, set_04); + + DocUtil::AddDocumentation(db, "ST_UpperLeftX", DOC_DESCRIPTION_04, DOC_EXAMPLE_04, DOC_TAGS_04); + + ScalarFunctionSet set_05("ST_UpperLeftY"); + set_05.AddFunction(ScalarFunction({GeoTypes::RASTER()}, LogicalType::DOUBLE, RasterGetUpperLeftY)); + ExtensionUtil::RegisterFunction(db, set_05); + + DocUtil::AddDocumentation(db, "ST_UpperLeftY", DOC_DESCRIPTION_05, DOC_EXAMPLE_05, DOC_TAGS_05); + + ScalarFunctionSet set_06("ST_ScaleX"); + set_06.AddFunction(ScalarFunction({GeoTypes::RASTER()}, LogicalType::DOUBLE, RasterGetScaleX)); + ExtensionUtil::RegisterFunction(db, set_06); + + DocUtil::AddDocumentation(db, "ST_ScaleX", DOC_DESCRIPTION_06, DOC_EXAMPLE_06, DOC_TAGS_06); + + ScalarFunctionSet set_07("ST_ScaleY"); + set_07.AddFunction(ScalarFunction({GeoTypes::RASTER()}, LogicalType::DOUBLE, RasterGetScaleY)); + ExtensionUtil::RegisterFunction(db, set_07); + + DocUtil::AddDocumentation(db, "ST_ScaleY", DOC_DESCRIPTION_07, DOC_EXAMPLE_07, DOC_TAGS_07); + + ScalarFunctionSet set_08("ST_SkewX"); + set_08.AddFunction(ScalarFunction({GeoTypes::RASTER()}, LogicalType::DOUBLE, RasterGetSkewX)); + ExtensionUtil::RegisterFunction(db, set_08); + + DocUtil::AddDocumentation(db, "ST_SkewX", DOC_DESCRIPTION_08, DOC_EXAMPLE_08, DOC_TAGS_08); + + ScalarFunctionSet set_09("ST_SkewY"); + set_09.AddFunction(ScalarFunction({GeoTypes::RASTER()}, LogicalType::DOUBLE, RasterGetSkewY)); + ExtensionUtil::RegisterFunction(db, set_09); + + DocUtil::AddDocumentation(db, "ST_SkewY", DOC_DESCRIPTION_09, DOC_EXAMPLE_09, DOC_TAGS_09); + + ScalarFunctionSet set_10("ST_PixelWidth"); + set_10.AddFunction(ScalarFunction({GeoTypes::RASTER()}, LogicalType::DOUBLE, RasterGetPixelWidthFunction)); + ExtensionUtil::RegisterFunction(db, set_10); + + DocUtil::AddDocumentation(db, "ST_PixelWidth", DOC_DESCRIPTION_10, DOC_EXAMPLE_10, DOC_TAGS_10); + + ScalarFunctionSet set_11("ST_PixelHeight"); + set_11.AddFunction(ScalarFunction({GeoTypes::RASTER()}, LogicalType::DOUBLE, RasterGetPixelHeightFunction)); + ExtensionUtil::RegisterFunction(db, set_11); + + DocUtil::AddDocumentation(db, "ST_PixelHeight", DOC_DESCRIPTION_11, DOC_EXAMPLE_11, DOC_TAGS_11); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/functions/scalar/st_raster_as_file.cpp b/spatial/src/spatial/gdal/functions/scalar/st_raster_as_file.cpp new file mode 100644 index 00000000..fa86f28c --- /dev/null +++ b/spatial/src/spatial/gdal/functions/scalar/st_raster_as_file.cpp @@ -0,0 +1,148 @@ +#include "duckdb/parser/parsed_data/create_scalar_function_info.hpp" +#include "duckdb/common/vector_operations/generic_executor.hpp" +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/gdal/file_handler.hpp" +#include "spatial/gdal/functions/scalar.hpp" +#include "spatial/gdal/raster/raster.hpp" +#include "spatial/gdal/raster/raster_factory.hpp" +#include "spatial/gdal/raster/raster_registry.hpp" + +#include "gdal_priv.h" + +namespace spatial { + +namespace gdal { + +//------------------------------------------------------------------------------ +// ST_RasterAsFile +//------------------------------------------------------------------------------ + +static void RasterAsFileFunction_01(DataChunk &args, ExpressionState &state, Vector &result) { + auto &context = state.GetContext(); + + TernaryExecutor::Execute( + args.data[0], args.data[1], args.data[2], result, args.size(), + [&](uintptr_t input, string_t file_name, string_t driver_name) { + GDALDataset *dataset = reinterpret_cast(input); + + auto gdal_driver_name = driver_name.GetString(); + if (gdal_driver_name.empty()) { + throw InvalidInputException("Driver name must be specified"); + } + + auto raw_file_name = file_name.GetString(); + auto &client_ctx = GDALClientContextState::GetOrCreate(context); + auto prefixed_file_name = client_ctx.GetPrefix(raw_file_name); + + if (!RasterFactory::WriteFile(dataset, prefixed_file_name, gdal_driver_name)) { + auto error = Raster::GetLastErrorMsg(); + if (error.length()) { + throw IOException("Could not save file: " + raw_file_name + " (" + error + ")"); + } + return false; + } + return true; + }); +} + +static void RasterAsFileFunction_02(DataChunk &args, ExpressionState &state, Vector &result) { + auto &context = state.GetContext(); + + using POINTER_TYPE = PrimitiveType; + using STRING_TYPE = PrimitiveType; + using LIST_TYPE = PrimitiveType; + using BOOL_TYPE = PrimitiveType; + + auto &p1 = args.data[0]; + auto &p2 = args.data[1]; + auto &p3 = args.data[2]; + auto &p4 = args.data[3]; + auto &p4_entry = ListVector::GetEntry(p4); + + GenericExecutor::ExecuteQuaternary( + p1, p2, p3, p4, result, args.size(), [&](POINTER_TYPE p1, STRING_TYPE p2, STRING_TYPE p3, LIST_TYPE p4_offlen) { + auto input = p1.val; + auto file_name = p2.val; + auto driver_name = p3.val; + auto offlen = p4_offlen.val; + + GDALDataset *dataset = reinterpret_cast(input); + + auto gdal_driver_name = driver_name.GetString(); + if (gdal_driver_name.empty()) { + throw InvalidInputException("Driver name must be specified"); + } + + auto raw_file_name = file_name.GetString(); + auto &client_ctx = GDALClientContextState::GetOrCreate(context); + auto prefixed_file_name = client_ctx.GetPrefix(raw_file_name); + + auto options = std::vector(); + + for (idx_t i = offlen.offset; i < offlen.offset + offlen.length; i++) { + const auto &child_value = p4_entry.GetValue(i); + const auto option = child_value.ToString(); + options.emplace_back(option); + } + + if (!RasterFactory::WriteFile(dataset, prefixed_file_name, gdal_driver_name, options)) { + auto error = Raster::GetLastErrorMsg(); + if (error.length()) { + throw IOException("Could not save file: " + raw_file_name + " (" + error + ")"); + } + return false; + } + return true; + }); +} + +//------------------------------------------------------------------------ +// Documentation +//------------------------------------------------------------------------ + +static constexpr const char *DOC_DESCRIPTION = R"( + Writes a raster to a file path. + + `write_options` is optional, an array of parameters for the GDAL driver specified. +)"; + +static constexpr const char *DOC_EXAMPLE = R"( + WITH __input AS ( + SELECT + ST_RasterFromFile(file) AS raster + FROM + glob('./test/data/mosaic/SCL.tif-land-clip00.tiff') + ) + SELECT + ST_RasterAsFile(raster, './rasterasfile.tiff', 'Gtiff', ['COMPRESS=LZW']) AS result + FROM + __input + ; +)"; + +static constexpr DocTag DOC_TAGS[] = {{"ext", "spatial"}, {"category", "construction"}}; + +//------------------------------------------------------------------------------ +// Register functions +//------------------------------------------------------------------------------ + +void GdalScalarFunctions::RegisterStRasterAsFile(DatabaseInstance &db) { + + ScalarFunctionSet set("ST_RasterAsFile"); + + set.AddFunction(ScalarFunction({GeoTypes::RASTER(), LogicalType::VARCHAR, LogicalType::VARCHAR}, + LogicalType::BOOLEAN, RasterAsFileFunction_01)); + + set.AddFunction(ScalarFunction( + {GeoTypes::RASTER(), LogicalType::VARCHAR, LogicalType::VARCHAR, LogicalType::LIST(LogicalType::VARCHAR)}, + LogicalType::BOOLEAN, RasterAsFileFunction_02)); + + ExtensionUtil::RegisterFunction(db, set); + + DocUtil::AddDocumentation(db, "ST_RasterAsFile", DOC_DESCRIPTION, DOC_EXAMPLE, DOC_TAGS); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/functions/scalar/st_raster_clip.cpp b/spatial/src/spatial/gdal/functions/scalar/st_raster_clip.cpp new file mode 100644 index 00000000..e53f669f --- /dev/null +++ b/spatial/src/spatial/gdal/functions/scalar/st_raster_clip.cpp @@ -0,0 +1,159 @@ +#include "duckdb/parser/parsed_data/create_scalar_function_info.hpp" +#include "duckdb/common/vector_operations/generic_executor.hpp" +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/gdal/file_handler.hpp" +#include "spatial/core/functions/common.hpp" +#include "spatial/gdal/functions/scalar.hpp" +#include "spatial/gdal/raster/raster.hpp" + +#include "gdal_priv.h" + +using namespace spatial::core; + +namespace spatial { + +namespace gdal { + +//------------------------------------------------------------------------------ +// ST_RasterClip +//------------------------------------------------------------------------------ + +static void RasterClipFunction_01(DataChunk &args, ExpressionState &state, Vector &result) { + auto &context = state.GetContext(); + auto &ctx_state = GDALClientContextState::GetOrCreate(context); + + using POINTER_TYPE = PrimitiveType; + using GEOMETRY_TYPE = PrimitiveType; + + auto &p1 = args.data[0]; + auto &p2 = args.data[1]; + + GenericExecutor::ExecuteBinary( + p1, p2, result, args.size(), [&](POINTER_TYPE p1, GEOMETRY_TYPE p2) { + auto input = p1.val; + auto geometry = p2.val; + + GDALDataset *dataset = reinterpret_cast(input); + + if (dataset->GetRasterCount() == 0) { + throw InvalidInputException("Input Raster has no RasterBands"); + } + + GDALDataset *result = Raster::Clip(dataset, geometry); + + if (result == nullptr) { + auto error = Raster::GetLastErrorMsg(); + throw IOException("Could not clip raster (" + error + ")"); + } + + auto &raster_registry = ctx_state.GetRasterRegistry(context); + raster_registry.RegisterRaster(result); + + return CastPointerToValue(result); + }); +} + +static void RasterClipFunction_02(DataChunk &args, ExpressionState &state, Vector &result) { + auto &context = state.GetContext(); + auto &ctx_state = GDALClientContextState::GetOrCreate(context); + + using POINTER_TYPE = PrimitiveType; + using GEOMETRY_TYPE = PrimitiveType; + using LIST_TYPE = PrimitiveType; + + auto &p1 = args.data[0]; + auto &p2 = args.data[1]; + auto &p3 = args.data[2]; + auto &p3_entry = ListVector::GetEntry(p3); + + GenericExecutor::ExecuteTernary( + p1, p2, p3, result, args.size(), [&](POINTER_TYPE p1, GEOMETRY_TYPE p2, LIST_TYPE p3_offlen) { + auto input = p1.val; + auto geometry = p2.val; + auto offlen = p3_offlen.val; + + GDALDataset *dataset = reinterpret_cast(input); + + if (dataset->GetRasterCount() == 0) { + throw InvalidInputException("Input Raster has no RasterBands"); + } + + auto options = std::vector(); + + for (idx_t i = offlen.offset; i < offlen.offset + offlen.length; i++) { + const auto &child_value = p3_entry.GetValue(i); + const auto option = child_value.ToString(); + options.emplace_back(option); + } + + GDALDataset *result = Raster::Clip(dataset, geometry, options); + + if (result == nullptr) { + auto error = Raster::GetLastErrorMsg(); + throw IOException("Could not clip raster (" + error + ")"); + } + + auto &raster_registry = ctx_state.GetRasterRegistry(context); + raster_registry.RegisterRaster(result); + + return CastPointerToValue(result); + }); +} + +//------------------------------------------------------------------------ +// Documentation +//------------------------------------------------------------------------ + +static constexpr const char *DOC_DESCRIPTION = R"( + Returns a raster that is clipped by the input geometry. + + `options` is optional, an array of parameters like [GDALWarp](https://gdal.org/programs/gdalwarp.html). +)"; + +static constexpr const char *DOC_EXAMPLE = R"( + WITH __input AS ( + SELECT + ST_RasterFromFile(file) AS raster + FROM + glob('./test/data/mosaic/SCL.tif-land-clip00.tiff') + ), + __geometry AS ( + SELECT geom FROM ST_Read('./test/data/mosaic/CATAST_Pol_Township-PNA.gpkg') + ) + SELECT + ST_RasterClip(mosaic, + (SELECT geom FROM __geometry LIMIT 1), + options => + [ + '-r', 'bilinear', '-crop_to_cutline', '-wo', 'CUTLINE_ALL_TOUCHED=TRUE' + ] + ) AS clip + FROM + __input + ; +)"; + +static constexpr DocTag DOC_TAGS[] = {{"ext", "spatial"}, {"category", "construction"}}; + +//------------------------------------------------------------------------------ +// Register functions +//------------------------------------------------------------------------------ + +void GdalScalarFunctions::RegisterStRasterClip(DatabaseInstance &db) { + + ScalarFunctionSet set("ST_RasterClip"); + set.AddFunction( + ScalarFunction({GeoTypes::RASTER(), GeoTypes::GEOMETRY()}, GeoTypes::RASTER(), RasterClipFunction_01)); + + set.AddFunction(ScalarFunction({GeoTypes::RASTER(), GeoTypes::GEOMETRY(), LogicalType::LIST(LogicalType::VARCHAR)}, + GeoTypes::RASTER(), RasterClipFunction_02)); + + ExtensionUtil::RegisterFunction(db, set); + + DocUtil::AddDocumentation(db, "ST_RasterClip", DOC_DESCRIPTION, DOC_EXAMPLE, DOC_TAGS); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/functions/scalar/st_raster_from_file.cpp b/spatial/src/spatial/gdal/functions/scalar/st_raster_from_file.cpp new file mode 100644 index 00000000..4ea38bee --- /dev/null +++ b/spatial/src/spatial/gdal/functions/scalar/st_raster_from_file.cpp @@ -0,0 +1,78 @@ +#include "duckdb/parser/parsed_data/create_scalar_function_info.hpp" +#include "duckdb/common/vector_operations/generic_executor.hpp" +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/gdal/file_handler.hpp" +#include "spatial/gdal/functions/scalar.hpp" +#include "spatial/gdal/raster/raster.hpp" +#include "spatial/gdal/raster/raster_factory.hpp" +#include "spatial/gdal/raster/raster_registry.hpp" + +#include "gdal_priv.h" + +namespace spatial { + +namespace gdal { + +//------------------------------------------------------------------------------ +// ST_RasterFromFile +//------------------------------------------------------------------------------ + +static void RasterFromFileFunction(DataChunk &args, ExpressionState &state, Vector &result) { + auto &context = state.GetContext(); + + UnaryExecutor::Execute(args.data[0], result, args.size(), [&](string_t input) { + auto &ctx_state = GDALClientContextState::GetOrCreate(context); + + auto raw_file_name = input.GetString(); + auto prefixed_file_name = ctx_state.GetPrefix(raw_file_name); + + GDALDataset *dataset = RasterFactory::FromFile(prefixed_file_name); + if (dataset == nullptr) { + auto error = Raster::GetLastErrorMsg(); + throw IOException("Could not open file: " + raw_file_name + " (" + error + ")"); + } + + auto &raster_registry = ctx_state.GetRasterRegistry(context); + raster_registry.RegisterRaster(dataset); + + return CastPointerToValue(dataset); + }); +} + +//------------------------------------------------------------------------ +// Documentation +//------------------------------------------------------------------------ + +static constexpr const char *DOC_DESCRIPTION = R"( + Loads a raster from a file path. +)"; + +static constexpr const char *DOC_EXAMPLE = R"( + WITH __input AS ( + SELECT + ST_RasterFromFile(file) AS raster + FROM + glob('./test/data/mosaic/*.tiff') + ) + SELECT raster from __input; +)"; + +static constexpr DocTag DOC_TAGS[] = {{"ext", "spatial"}, {"category", "construction"}}; + +//------------------------------------------------------------------------------ +// Register functions +//------------------------------------------------------------------------------ + +void GdalScalarFunctions::RegisterStRasterFromFile(DatabaseInstance &db) { + + ScalarFunctionSet set("ST_RasterFromFile"); + set.AddFunction(ScalarFunction({LogicalType::VARCHAR}, GeoTypes::RASTER(), RasterFromFileFunction)); + ExtensionUtil::RegisterFunction(db, set); + + DocUtil::AddDocumentation(db, "ST_RasterFromFile", DOC_DESCRIPTION, DOC_EXAMPLE, DOC_TAGS); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/functions/scalar/st_raster_geometry.cpp b/spatial/src/spatial/gdal/functions/scalar/st_raster_geometry.cpp new file mode 100644 index 00000000..efccbdb4 --- /dev/null +++ b/spatial/src/spatial/gdal/functions/scalar/st_raster_geometry.cpp @@ -0,0 +1,62 @@ +#include "duckdb/parser/parsed_data/create_scalar_function_info.hpp" +#include "duckdb/common/vector_operations/generic_executor.hpp" +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/core/functions/common.hpp" +#include "spatial/core/geometry/geometry.hpp" +#include "spatial/gdal/functions/scalar.hpp" +#include "spatial/gdal/raster/raster.hpp" + +#include "gdal_priv.h" + +using namespace spatial::core; + +namespace spatial { + +namespace gdal { + +//------------------------------------------------------------------------------ +// ST_GetGeometry +//------------------------------------------------------------------------------ + +static void RasterGetGeometryFunction(DataChunk &args, ExpressionState &state, Vector &result) { + auto &lstate = GeometryFunctionLocalState::ResetAndGet(state); + + UnaryExecutor::Execute(args.data[0], result, args.size(), [&](uintptr_t input) { + Raster raster(reinterpret_cast(input)); + return Geometry::Serialize(raster.GetGeometry(lstate.arena), result); + }); +} + +//------------------------------------------------------------------------ +// Documentation +//------------------------------------------------------------------------ + +static constexpr const char *DOC_DESCRIPTION = R"( + Returns the polygon representation of the extent of the raster. +)"; + +static constexpr const char *DOC_EXAMPLE = R"( + SELECT ST_GetGeometry(raster) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS[] = {{"ext", "spatial"}, {"category", "property"}}; + +//------------------------------------------------------------------------------ +// Register functions +//------------------------------------------------------------------------------ + +void GdalScalarFunctions::RegisterStGetGeometry(DatabaseInstance &db) { + + ScalarFunctionSet set("ST_GetGeometry"); + set.AddFunction(ScalarFunction({GeoTypes::RASTER()}, GeoTypes::GEOMETRY(), RasterGetGeometryFunction, nullptr, + nullptr, nullptr, GeometryFunctionLocalState::Init)); + + ExtensionUtil::RegisterFunction(db, set); + + DocUtil::AddDocumentation(db, "ST_GetGeometry", DOC_DESCRIPTION, DOC_EXAMPLE, DOC_TAGS); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/functions/scalar/st_raster_to_world_coord.cpp b/spatial/src/spatial/gdal/functions/scalar/st_raster_to_world_coord.cpp new file mode 100644 index 00000000..28d28591 --- /dev/null +++ b/spatial/src/spatial/gdal/functions/scalar/st_raster_to_world_coord.cpp @@ -0,0 +1,143 @@ +#include "duckdb/parser/parsed_data/create_scalar_function_info.hpp" +#include "duckdb/common/vector_operations/generic_executor.hpp" +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/gdal/functions/scalar.hpp" +#include "spatial/gdal/raster/raster.hpp" + +#include "gdal_priv.h" + +using namespace spatial::core; + +namespace spatial { + +namespace gdal { + +//------------------------------------------------------------------------------ +// ST_RasterToWorldCoord[XY] +//------------------------------------------------------------------------------ + +static void RasterRasterToWorldCoordFunction(DataChunk &args, ExpressionState &state, Vector &result) { + D_ASSERT(args.data.size() == 3); + + using POINTER_TYPE = PrimitiveType; + using INT_TYPE = PrimitiveType; + using POINT_TYPE = StructTypeBinary; + + auto &p1 = args.data[0]; + auto &p2 = args.data[1]; + auto &p3 = args.data[2]; + + GenericExecutor::ExecuteTernary( + p1, p2, p3, result, args.size(), [&](POINTER_TYPE p1, INT_TYPE p2, INT_TYPE p3) { + auto input = p1.val; + auto col = p2.val; + auto row = p3.val; + Raster raster(reinterpret_cast(input)); + + PointXY coord(0, 0); + if (!raster.RasterToWorldCoord(coord, col, row)) { + throw InternalException("Could not compute geotransform matrix"); + } + return POINT_TYPE {coord.x, coord.y}; + }); +} + +static void RasterRasterToWorldCoordXFunction(DataChunk &args, ExpressionState &state, Vector &result) { + D_ASSERT(args.data.size() == 3); + + TernaryExecutor::Execute( + args.data[0], args.data[1], args.data[2], result, args.size(), [&](uintptr_t input, int32_t col, int32_t row) { + Raster raster(reinterpret_cast(input)); + + PointXY coord(0, 0); + if (!raster.RasterToWorldCoord(coord, col, row)) { + throw InternalException("Could not compute geotransform matrix"); + } + return coord.x; + }); +} + +static void RasterRasterToWorldCoordYFunction(DataChunk &args, ExpressionState &state, Vector &result) { + D_ASSERT(args.data.size() == 3); + + TernaryExecutor::Execute( + args.data[0], args.data[1], args.data[2], result, args.size(), [&](uintptr_t input, int32_t col, int32_t row) { + Raster raster(reinterpret_cast(input)); + + PointXY coord(0, 0); + if (!raster.RasterToWorldCoord(coord, col, row)) { + throw InternalException("Could not compute geotransform matrix"); + } + return coord.y; + }); +} + +//------------------------------------------------------------------------ +// Documentation +//------------------------------------------------------------------------ + +static constexpr const char *DOC_DESCRIPTION_1 = R"( + Returns the upper left corner as geometric X and Y (longitude and latitude) given a column and row. + Returned X and Y are in geometric units of the georeferenced raster. +)"; + +static constexpr const char *DOC_EXAMPLE_1 = R"( + SELECT ST_RasterToWorldCoord(raster) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_1[] = {{"ext", "spatial"}, {"category", "position"}}; + +static constexpr const char *DOC_DESCRIPTION_2 = R"( + Returns the upper left X coordinate of a raster column row in geometric units of the georeferenced raster. + Returned X is in geometric units of the georeferenced raster. +)"; + +static constexpr const char *DOC_EXAMPLE_2 = R"( + SELECT ST_RasterToWorldCoordX(raster) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_2[] = {{"ext", "spatial"}, {"category", "position"}}; + +static constexpr const char *DOC_DESCRIPTION_3 = R"( + Returns the upper left Y coordinate of a raster column row in geometric units of the georeferenced raster. + Returned Y is in geometric units of the georeferenced raster. +)"; + +static constexpr const char *DOC_EXAMPLE_3 = R"( + SELECT ST_RasterToWorldCoordY(raster) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_3[] = {{"ext", "spatial"}, {"category", "position"}}; + +//------------------------------------------------------------------------------ +// Register functions +//------------------------------------------------------------------------------ + +void GdalScalarFunctions::RegisterStRasterToWorldCoord(DatabaseInstance &db) { + + ScalarFunctionSet set01("ST_RasterToWorldCoord"); + set01.AddFunction(ScalarFunction({GeoTypes::RASTER(), LogicalType::INTEGER, LogicalType::INTEGER}, + GeoTypes::POINT_2D(), RasterRasterToWorldCoordFunction)); + ExtensionUtil::RegisterFunction(db, set01); + + DocUtil::AddDocumentation(db, "ST_RasterToWorldCoord", DOC_DESCRIPTION_1, DOC_EXAMPLE_1, DOC_TAGS_1); + + ScalarFunctionSet set02("ST_RasterToWorldCoordX"); + set02.AddFunction(ScalarFunction({GeoTypes::RASTER(), LogicalType::INTEGER, LogicalType::INTEGER}, + LogicalType::DOUBLE, RasterRasterToWorldCoordXFunction)); + ExtensionUtil::RegisterFunction(db, set02); + + DocUtil::AddDocumentation(db, "ST_RasterToWorldCoordX", DOC_DESCRIPTION_2, DOC_EXAMPLE_2, DOC_TAGS_2); + + ScalarFunctionSet set03("ST_RasterToWorldCoordY"); + set03.AddFunction(ScalarFunction({GeoTypes::RASTER(), LogicalType::INTEGER, LogicalType::INTEGER}, + LogicalType::DOUBLE, RasterRasterToWorldCoordYFunction)); + ExtensionUtil::RegisterFunction(db, set03); + + DocUtil::AddDocumentation(db, "ST_RasterToWorldCoordY", DOC_DESCRIPTION_3, DOC_EXAMPLE_3, DOC_TAGS_3); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/functions/scalar/st_raster_warp.cpp b/spatial/src/spatial/gdal/functions/scalar/st_raster_warp.cpp new file mode 100644 index 00000000..a70f4ca8 --- /dev/null +++ b/spatial/src/spatial/gdal/functions/scalar/st_raster_warp.cpp @@ -0,0 +1,110 @@ +#include "duckdb/parser/parsed_data/create_scalar_function_info.hpp" +#include "duckdb/common/vector_operations/generic_executor.hpp" +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/gdal/file_handler.hpp" +#include "spatial/core/functions/common.hpp" +#include "spatial/gdal/functions/scalar.hpp" +#include "spatial/gdal/raster/raster.hpp" + +#include "gdal_priv.h" + +using namespace spatial::core; + +namespace spatial { + +namespace gdal { + +//------------------------------------------------------------------------------ +// ST_RasterWarp +//------------------------------------------------------------------------------ + +static void RasterWarpFunction(DataChunk &args, ExpressionState &state, Vector &result) { + auto &context = state.GetContext(); + auto &ctx_state = GDALClientContextState::GetOrCreate(context); + + using POINTER_TYPE = PrimitiveType; + using LIST_TYPE = PrimitiveType; + + auto &p1 = args.data[0]; + auto &p2 = args.data[1]; + auto &p2_entry = ListVector::GetEntry(p2); + + GenericExecutor::ExecuteBinary( + p1, p2, result, args.size(), [&](POINTER_TYPE p1, LIST_TYPE p2_offlen) { + auto input = p1.val; + auto offlen = p2_offlen.val; + + GDALDataset *dataset = reinterpret_cast(input); + + if (dataset->GetRasterCount() == 0) { + throw InvalidInputException("Input Raster has no RasterBands"); + } + + auto options = std::vector(); + + for (idx_t i = offlen.offset; i < offlen.offset + offlen.length; i++) { + const auto &child_value = p2_entry.GetValue(i); + const auto option = child_value.ToString(); + options.emplace_back(option); + } + + GDALDataset *result = Raster::Warp(dataset, options); + + if (result == nullptr) { + auto error = Raster::GetLastErrorMsg(); + throw IOException("Could not warp raster (" + error + ")"); + } + + auto &raster_registry = ctx_state.GetRasterRegistry(context); + raster_registry.RegisterRaster(result); + + return CastPointerToValue(result); + }); +} + +//------------------------------------------------------------------------ +// Documentation +//------------------------------------------------------------------------ + +static constexpr const char *DOC_DESCRIPTION = R"( + Performs mosaicing, reprojection and/or warping on a raster. + + `options` is optional, an array of parameters like [GDALWarp](https://gdal.org/programs/gdalwarp.html). +)"; + +static constexpr const char *DOC_EXAMPLE = R"( + WITH __input AS ( + SELECT + raster + FROM + ST_ReadRaster('./test/data/mosaic/SCL.tif-land-clip00.tiff') + ), + SELECT + ST_RasterWarp(raster, options => ['-r', 'bilinear', '-tr', '40.0', '40.0']) AS warp + FROM + __input + ; +)"; + +static constexpr DocTag DOC_TAGS[] = {{"ext", "spatial"}, {"category", "construction"}}; + +//------------------------------------------------------------------------------ +// Register functions +//------------------------------------------------------------------------------ + +void GdalScalarFunctions::RegisterStRasterWarp(DatabaseInstance &db) { + + ScalarFunctionSet set("ST_RasterWarp"); + + set.AddFunction(ScalarFunction({GeoTypes::RASTER(), LogicalType::LIST(LogicalType::VARCHAR)}, GeoTypes::RASTER(), + RasterWarpFunction, nullptr, nullptr, nullptr, GeometryFunctionLocalState::Init)); + + ExtensionUtil::RegisterFunction(db, set); + + DocUtil::AddDocumentation(db, "ST_RasterWarp", DOC_DESCRIPTION, DOC_EXAMPLE, DOC_TAGS); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/functions/scalar/st_srid.cpp b/spatial/src/spatial/gdal/functions/scalar/st_srid.cpp new file mode 100644 index 00000000..efe52f46 --- /dev/null +++ b/spatial/src/spatial/gdal/functions/scalar/st_srid.cpp @@ -0,0 +1,57 @@ +#include "duckdb/parser/parsed_data/create_scalar_function_info.hpp" +#include "duckdb/common/vector_operations/generic_executor.hpp" +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/gdal/functions/scalar.hpp" +#include "spatial/gdal/raster/raster.hpp" + +using namespace spatial::core; + +namespace spatial { + +namespace gdal { + +//------------------------------------------------------------------------------ +// ST_SRID +//------------------------------------------------------------------------------ + +static void RasterGetSridFunction(DataChunk &args, ExpressionState &state, Vector &result) { + + UnaryExecutor::Execute(args.data[0], result, args.size(), [&](uintptr_t input) { + Raster raster(reinterpret_cast(input)); + return raster.GetSrid(); + }); +} + +//------------------------------------------------------------------------ +// Documentation +//------------------------------------------------------------------------ + +static constexpr const char *DOC_DESCRIPTION = R"( + Returns the spatial reference identifier (EPSG code) of the raster. + Refer to [EPSG](https://spatialreference.org/ref/epsg/) for more details. +)"; + +static constexpr const char *DOC_EXAMPLE = R"( + SELECT ST_SRID(raster) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS[] = {{"ext", "spatial"}, {"category", "property"}}; + +//------------------------------------------------------------------------------ +// Register functions +//------------------------------------------------------------------------------ + +void GdalScalarFunctions::RegisterStGetSRID(DatabaseInstance &db) { + + ScalarFunctionSet set("ST_SRID"); + set.AddFunction(ScalarFunction({GeoTypes::RASTER()}, LogicalType::INTEGER, RasterGetSridFunction)); + + ExtensionUtil::RegisterFunction(db, set); + + DocUtil::AddDocumentation(db, "ST_SRID", DOC_DESCRIPTION, DOC_EXAMPLE, DOC_TAGS); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/functions/scalar/st_value.cpp b/spatial/src/spatial/gdal/functions/scalar/st_value.cpp new file mode 100644 index 00000000..88c103e7 --- /dev/null +++ b/spatial/src/spatial/gdal/functions/scalar/st_value.cpp @@ -0,0 +1,94 @@ +#include "duckdb/parser/parsed_data/create_scalar_function_info.hpp" +#include "duckdb/common/vector_operations/generic_executor.hpp" +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/gdal/functions/scalar.hpp" +#include "spatial/gdal/raster/raster.hpp" + +#include "gdal_priv.h" + +namespace spatial { + +namespace gdal { + +//------------------------------------------------------------------------------ +// ST_Value +//------------------------------------------------------------------------------ + +static void RasterGetValueFunction(DataChunk &args, ExpressionState &state, Vector &result) { + D_ASSERT(args.data.size() == 4); + + using POINTER_TYPE = PrimitiveType; + using INT_TYPE = PrimitiveType; + using DOUBLE_TYPE = PrimitiveType; + + auto &p1 = args.data[0]; + auto &p2 = args.data[1]; + auto &p3 = args.data[2]; + auto &p4 = args.data[3]; + + GenericExecutor::ExecuteQuaternary( + p1, p2, p3, p4, result, args.size(), [&](POINTER_TYPE p1, INT_TYPE p2, INT_TYPE p3, INT_TYPE p4) { + auto input = p1.val; + auto band_num = p2.val; + auto col = p3.val; + auto row = p4.val; + + GDALDataset *dataset = reinterpret_cast(input); + auto cols = dataset->GetRasterXSize(); + auto rows = dataset->GetRasterYSize(); + + if (band_num < 1) { + throw InvalidInputException("BandNum must be greater than 0"); + } + if (dataset->GetRasterCount() < band_num) { + throw InvalidInputException("Dataset only has %d RasterBands", dataset->GetRasterCount()); + } + if (col < 0 || col >= cols || row < 0 || row >= rows) { + throw InvalidInputException( + "Attempting to get pixel value with out of range raster coordinates: (%d, %d)", col, row); + } + + Raster raster(dataset); + double value; + if (raster.GetValue(value, band_num, col, row)) { + return value; + } + throw InternalException("Failed attempting to get pixel value with raster coordinates: (%d, %d)", col, row); + }); +} + +//------------------------------------------------------------------------ +// Documentation +//------------------------------------------------------------------------ + +static constexpr const char *DOC_DESCRIPTION = R"( + Returns the value of a given band in a given column, row pixel. + Band numbers start at 1 and band is assumed to be 1 if not specified. +)"; + +static constexpr const char *DOC_EXAMPLE = R"( + SELECT ST_Value(raster, 1, 0, 0) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS[] = {{"ext", "spatial"}, {"category", "property"}}; + +//------------------------------------------------------------------------------ +// Register functions +//------------------------------------------------------------------------------ + +void GdalScalarFunctions::RegisterStGetValue(DatabaseInstance &db) { + + ScalarFunctionSet set("ST_Value"); + set.AddFunction( + ScalarFunction({GeoTypes::RASTER(), LogicalType::INTEGER, LogicalType::INTEGER, LogicalType::INTEGER}, + LogicalType::DOUBLE, RasterGetValueFunction)); + + ExtensionUtil::RegisterFunction(db, set); + + DocUtil::AddDocumentation(db, "ST_Value", DOC_DESCRIPTION, DOC_EXAMPLE, DOC_TAGS); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/functions/scalar/st_world_to_raster_coord.cpp b/spatial/src/spatial/gdal/functions/scalar/st_world_to_raster_coord.cpp new file mode 100644 index 00000000..dc06a2c5 --- /dev/null +++ b/spatial/src/spatial/gdal/functions/scalar/st_world_to_raster_coord.cpp @@ -0,0 +1,142 @@ +#include "duckdb/parser/parsed_data/create_scalar_function_info.hpp" +#include "duckdb/common/vector_operations/generic_executor.hpp" +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/gdal/types.hpp" +#include "spatial/gdal/functions/scalar.hpp" +#include "spatial/gdal/raster/raster.hpp" + +#include "gdal_priv.h" + +namespace spatial { + +namespace gdal { + +//------------------------------------------------------------------------------ +// ST_WorldToRasterCoord[XY] +//------------------------------------------------------------------------------ + +static void RasterWorldToRasterCoordFunction(DataChunk &args, ExpressionState &state, Vector &result) { + D_ASSERT(args.data.size() == 3); + + using POINTER_TYPE = PrimitiveType; + using DOUBLE_TYPE = PrimitiveType; + using COORD_TYPE = StructTypeBinary; + + auto &p1 = args.data[0]; + auto &p2 = args.data[1]; + auto &p3 = args.data[2]; + + GenericExecutor::ExecuteTernary( + p1, p2, p3, result, args.size(), [&](POINTER_TYPE p1, DOUBLE_TYPE p2, DOUBLE_TYPE p3) { + auto input = p1.val; + auto x = p2.val; + auto y = p3.val; + Raster raster(reinterpret_cast(input)); + + RasterCoord coord(0, 0); + if (!raster.WorldToRasterCoord(coord, x, y)) { + throw InternalException("Could not compute inverse geotransform matrix"); + } + return COORD_TYPE {coord.col, coord.row}; + }); +} + +static void RasterWorldToRasterCoordXFunction(DataChunk &args, ExpressionState &state, Vector &result) { + D_ASSERT(args.data.size() == 3); + + TernaryExecutor::Execute( + args.data[0], args.data[1], args.data[2], result, args.size(), [&](uintptr_t input, double_t x, double_t y) { + Raster raster(reinterpret_cast(input)); + + RasterCoord coord(0, 0); + if (!raster.WorldToRasterCoord(coord, x, y)) { + throw InternalException("Could not compute inverse geotransform matrix"); + } + return coord.col; + }); +} + +static void RasterWorldToRasterCoordYFunction(DataChunk &args, ExpressionState &state, Vector &result) { + D_ASSERT(args.data.size() == 3); + + TernaryExecutor::Execute( + args.data[0], args.data[1], args.data[2], result, args.size(), [&](uintptr_t input, double_t x, double_t y) { + Raster raster(reinterpret_cast(input)); + + RasterCoord coord(0, 0); + if (!raster.WorldToRasterCoord(coord, x, y)) { + throw InternalException("Could not compute inverse geotransform matrix"); + } + return coord.row; + }); +} + +//------------------------------------------------------------------------ +// Documentation +//------------------------------------------------------------------------ + +static constexpr const char *DOC_DESCRIPTION_1 = R"( + Returns the upper left corner as column and row given geometric X and Y (longitude and latitude). + Geometric X and Y must be expressed in the spatial reference coordinate system of the raster. +)"; + +static constexpr const char *DOC_EXAMPLE_1 = R"( + SELECT ST_WorldToRasterCoord(raster) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_1[] = {{"ext", "spatial"}, {"category", "position"}}; + +static constexpr const char *DOC_DESCRIPTION_2 = R"( + Returns the column in the raster given geometric X and Y (longitude and latitude). + Geometric X and Y must be expressed in the spatial reference coordinate system of the raster. +)"; + +static constexpr const char *DOC_EXAMPLE_2 = R"( + SELECT ST_WorldToRasterCoordX(raster) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_2[] = {{"ext", "spatial"}, {"category", "position"}}; + +static constexpr const char *DOC_DESCRIPTION_3 = R"( + Returns the row in the raster given geometric X and Y (longitude and latitude). + Geometric X and Y must be expressed in the spatial reference coordinate system of the raster. +)"; + +static constexpr const char *DOC_EXAMPLE_3 = R"( + SELECT ST_WorldToRasterCoordY(raster) FROM './test/data/mosaic/SCL.tif-land-clip00.tiff'; +)"; + +static constexpr DocTag DOC_TAGS_3[] = {{"ext", "spatial"}, {"category", "position"}}; + +//------------------------------------------------------------------------------ +// Register functions +//------------------------------------------------------------------------------ + +void GdalScalarFunctions::RegisterStWorldToRasterCoord(DatabaseInstance &db) { + + ScalarFunctionSet set01("ST_WorldToRasterCoord"); + set01.AddFunction(ScalarFunction({GeoTypes::RASTER(), LogicalType::DOUBLE, LogicalType::DOUBLE}, + GeoTypes::RASTER_COORD(), RasterWorldToRasterCoordFunction)); + ExtensionUtil::RegisterFunction(db, set01); + + DocUtil::AddDocumentation(db, "ST_WorldToRasterCoord", DOC_DESCRIPTION_1, DOC_EXAMPLE_1, DOC_TAGS_1); + + ScalarFunctionSet set02("ST_WorldToRasterCoordX"); + set02.AddFunction(ScalarFunction({GeoTypes::RASTER(), LogicalType::DOUBLE, LogicalType::DOUBLE}, + LogicalType::INTEGER, RasterWorldToRasterCoordXFunction)); + ExtensionUtil::RegisterFunction(db, set02); + + DocUtil::AddDocumentation(db, "ST_WorldToRasterCoordX", DOC_DESCRIPTION_2, DOC_EXAMPLE_2, DOC_TAGS_2); + + ScalarFunctionSet set03("ST_WorldToRasterCoordY"); + set03.AddFunction(ScalarFunction({GeoTypes::RASTER(), LogicalType::DOUBLE, LogicalType::DOUBLE}, + LogicalType::INTEGER, RasterWorldToRasterCoordYFunction)); + ExtensionUtil::RegisterFunction(db, set03); + + DocUtil::AddDocumentation(db, "ST_WorldToRasterCoordY", DOC_DESCRIPTION_3, DOC_EXAMPLE_3, DOC_TAGS_3); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/functions/st_drivers.cpp b/spatial/src/spatial/gdal/functions/st_drivers.cpp index 3c66be1e..1efe25ae 100644 --- a/spatial/src/spatial/gdal/functions/st_drivers.cpp +++ b/spatial/src/spatial/gdal/functions/st_drivers.cpp @@ -14,12 +14,14 @@ unique_ptr GdalDriversTableFunction::Bind(ClientContext &context, vector &return_types, vector &names) { return_types.emplace_back(LogicalType::VARCHAR); return_types.emplace_back(LogicalType::VARCHAR); + return_types.emplace_back(LogicalType::VARCHAR); return_types.emplace_back(LogicalType::BOOLEAN); return_types.emplace_back(LogicalType::BOOLEAN); return_types.emplace_back(LogicalType::BOOLEAN); return_types.emplace_back(LogicalType::VARCHAR); names.emplace_back("short_name"); names.emplace_back("long_name"); + names.emplace_back("type"); names.emplace_back("can_create"); names.emplace_back("can_copy"); names.emplace_back("can_open"); @@ -46,14 +48,23 @@ void GdalDriversTableFunction::Execute(ClientContext &context, TableFunctionInpu for (; state.current_idx < next_idx; state.current_idx++) { auto driver = GDALGetDriver((int)state.current_idx); - // Check if the driver is a vector driver - if (GDALGetMetadataItem(driver, GDAL_DCAP_VECTOR, nullptr) == nullptr) { + const char *vector_caps = GDALGetMetadataItem(driver, GDAL_DCAP_VECTOR, nullptr); + bool supports_vector = vector_caps && strcmp(vector_caps, "YES") == 0; + const char *raster_caps = GDALGetMetadataItem(driver, GDAL_DCAP_RASTER, nullptr); + bool supports_raster = raster_caps && strcmp(raster_caps, "YES") == 0; + + // Check if the driver supports vector or raster + if (!supports_vector && !supports_raster) { continue; } auto short_name = Value::CreateValue(GDALGetDriverShortName(driver)); auto long_name = Value::CreateValue(GDALGetDriverLongName(driver)); + std::string driver_type = supports_vector ? "vector" : ""; + if (supports_raster) + driver_type += (supports_vector) ? ",raster" : "raster"; + const char *create_flag = GDALGetMetadataItem(driver, GDAL_DCAP_CREATE, nullptr); auto create_value = Value::CreateValue(create_flag != nullptr); @@ -69,10 +80,11 @@ void GdalDriversTableFunction::Execute(ClientContext &context, TableFunctionInpu output.data[0].SetValue(count, short_name); output.data[1].SetValue(count, long_name); - output.data[2].SetValue(count, create_value); - output.data[3].SetValue(count, copy_value); - output.data[4].SetValue(count, open_value); - output.data[5].SetValue(count, help_topic_value); + output.data[2].SetValue(count, Value::CreateValue(driver_type)); + output.data[3].SetValue(count, create_value); + output.data[4].SetValue(count, copy_value); + output.data[5].SetValue(count, open_value); + output.data[6].SetValue(count, help_topic_value); count++; } output.SetCardinality(count); @@ -86,7 +98,7 @@ static constexpr DocTag DOC_TAGS[] = {{"ext", "spatial"}}; static constexpr const char *DOC_DESCRIPTION = R"( Returns the list of supported GDAL drivers and file formats - Note that far from all of these drivers have been tested properly, and some may require additional options to be passed to work as expected. If you run into any issues please first consult the [consult the GDAL docs](https://gdal.org/drivers/vector/index.html). + Note that far from all of these drivers have been tested properly, and some may require additional options to be passed to work as expected. If you run into any issues please first consult the GDAL docs for [vector drivers](https://gdal.org/drivers/vector/index.html) and [raster drivers](https://gdal.org/drivers/raster/index.html). )"; static constexpr const char *DOC_EXAMPLE = R"( diff --git a/spatial/src/spatial/gdal/functions/st_read_raster.cpp b/spatial/src/spatial/gdal/functions/st_read_raster.cpp new file mode 100644 index 00000000..3f41af2a --- /dev/null +++ b/spatial/src/spatial/gdal/functions/st_read_raster.cpp @@ -0,0 +1,203 @@ +#include "duckdb/parser/parsed_data/create_table_function_info.hpp" +#include "duckdb/parser/expression/constant_expression.hpp" +#include "duckdb/parser/expression/function_expression.hpp" +#include "duckdb/parser/tableref/table_function_ref.hpp" +#include "duckdb/function/function.hpp" +#include "duckdb/function/replacement_scan.hpp" + +#include "spatial/common.hpp" +#include "spatial/core/types.hpp" +#include "spatial/gdal/functions.hpp" +#include "spatial/gdal/file_handler.hpp" +#include "spatial/gdal/raster/raster.hpp" +#include "spatial/gdal/raster/raster_factory.hpp" +#include "spatial/gdal/raster/raster_registry.hpp" +#include "spatial/gdal/raster/raster_value.hpp" + +#include "gdal_priv.h" + +using namespace spatial::core; + +namespace spatial { + +namespace gdal { + +struct GdalRasterTableFunctionData : public TableFunctionData { + string file_name; + named_parameter_map_t parameters; + RasterRegistry raster_registry; + bool loaded; +}; + +//------------------------------------------------------------------------------ +// Bind +//------------------------------------------------------------------------------ + +unique_ptr GdalRasterTableFunction::Bind(ClientContext &context, TableFunctionBindInput &input, + vector &return_types, vector &names) { + return_types.emplace_back(LogicalType::VARCHAR); + return_types.emplace_back(GeoTypes::RASTER()); + names.emplace_back("path"); + names.emplace_back("raster"); + + auto raw_file_name = input.inputs[0].GetValue(); + auto parameters = input.named_parameters; + + auto result = make_uniq(); + result->file_name = raw_file_name; + result->parameters = parameters; + result->loaded = false; + return std::move(result); +} + +//------------------------------------------------------------------------------ +// Execute +//------------------------------------------------------------------------------ + +void GdalRasterTableFunction::Execute(ClientContext &context, TableFunctionInput &input, DataChunk &output) { + auto &bind_data = (GdalRasterTableFunctionData &)*input.bind_data; + + if (bind_data.loaded) { + output.SetCardinality(0); + return; + } + + auto &config = DBConfig::GetConfig(context); + if (!config.options.enable_external_access) { + throw PermissionException("Scanning GDAL files is disabled through configuration"); + } + + // First scan for "options" parameter + auto gdal_open_options = RasterFactory::FromNamedParameters(bind_data.parameters, "open_options"); + + auto gdal_allowed_drivers = RasterFactory::FromNamedParameters(bind_data.parameters, "allowed_drivers"); + + auto gdal_sibling_files = RasterFactory::FromNamedParameters(bind_data.parameters, "sibling_files"); + + // Now we can open the dataset + auto raw_file_name = bind_data.file_name; + auto &ctx_state = GDALClientContextState::GetOrCreate(context); + auto prefixed_file_name = ctx_state.GetPrefix(raw_file_name); + auto dataset = GDALDataset::Open(prefixed_file_name.c_str(), GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, + gdal_allowed_drivers.empty() ? nullptr : gdal_allowed_drivers.data(), + gdal_open_options.empty() ? nullptr : gdal_open_options.data(), + gdal_sibling_files.empty() ? nullptr : gdal_sibling_files.data()); + + if (dataset == nullptr) { + auto error = Raster::GetLastErrorMsg(); + throw IOException("Could not open file: " + raw_file_name + " (" + error + ")"); + } + + // Now we can bind the dataset + bind_data.raster_registry.RegisterRaster(dataset); + bind_data.loaded = true; + + // And fill the output + output.data[0].SetValue(0, Value::CreateValue(raw_file_name)); + output.data[1].SetValue(0, RasterValue::CreateValue(dataset)); + output.SetCardinality(1); +}; + +//------------------------------------------------------------------------------ +// Cardinality +//------------------------------------------------------------------------------ + +unique_ptr GdalRasterTableFunction::Cardinality(ClientContext &context, const FunctionData *data) { + auto result = make_uniq(); + result->has_estimated_cardinality = true; + result->estimated_cardinality = 1; + result->has_max_cardinality = true; + result->max_cardinality = 1; + return result; +} + +//------------------------------------------------------------------------------ +// ReplacementScan +//------------------------------------------------------------------------------ + +unique_ptr GdalRasterTableFunction::ReplacementScan(ClientContext &, ReplacementScanInput &input, + optional_ptr) { + + auto &table_name = input.table_name; + auto lower_name = StringUtil::Lower(table_name); + + // Check if the file name ends with some common raster file extensions + if (StringUtil::EndsWith(lower_name, ".img") || StringUtil::EndsWith(lower_name, ".tiff") || + StringUtil::EndsWith(lower_name, ".tif") || StringUtil::EndsWith(lower_name, ".vrt")) { + + auto table_function = make_uniq(); + vector> children; + children.push_back(make_uniq(Value(table_name))); + table_function->function = make_uniq("ST_ReadRaster", std::move(children)); + return std::move(table_function); + } + // else not something we can replace + return nullptr; +} + +//------------------------------------------------------------------------ +// Documentation +//------------------------------------------------------------------------ + +static constexpr const char *DOC_DESCRIPTION = R"( + The `ST_ReadRaster` table function is based on the [GDAL](https://gdal.org/index.html) translator library and enables reading spatial data from a variety of geospatial raster file formats as if they were DuckDB tables. + See `ST_Drivers` for a list of supported file formats and drivers. + + Except for the `path` parameter, all parameters are optional. + + | Parameter | Type | Description | + | --------- | -----| ----------- | + | `path` | VARCHAR | The path to the file to read. Mandatory | + | `open_options` | VARCHAR[] | A list of key-value pairs that are passed to the GDAL driver to control the opening of the file. | + | `allowed_drivers` | VARCHAR[] | A list of GDAL driver names that are allowed to be used to open the file. If empty, all drivers are allowed. | + | `sibling_files` | VARCHAR[] | A list of sibling files that are required to open the file. E.g., the ESRI Shapefile driver requires a .shx file to be present. Although most of the time these can be discovered automatically. | + + Note that GDAL is single-threaded, so this table function will not be able to make full use of parallelism. + + By using `ST_ReadRaster`, the spatial extension also provides “replacement scans” for common raster file formats, allowing you to query files of these formats as if they were tables directly. + In practice this is just syntax-sugar for calling ST_ReadRaster, so there is no difference in performance. If you want to pass additional options, you should use the ST_ReadRaster table function directly. + + The following formats are currently recognized by their file extension: + + | Format | Extension | + | ------ | --------- | + | GeoTiff COG | .tif, .tiff | + | Erdas Imagine | .img | + | GDAL Virtual | .vrt | +)"; + +static constexpr const char *DOC_EXAMPLE = R"( + SELECT * FROM ST_ReadRaster('some/file/path/filename.tiff'); + + SELECT * FROM './path/to/some/shapefile/dataset.tiff'; +)"; + +static constexpr DocTag DOC_TAGS[] = {{"ext", "spatial"}}; + +//------------------------------------------------------------------------------ +// Register +//------------------------------------------------------------------------------ + +void GdalRasterTableFunction::Register(DatabaseInstance &db) { + + TableFunctionSet set("ST_ReadRaster"); + + TableFunction func({LogicalType::VARCHAR}, Execute, Bind); + func.cardinality = GdalRasterTableFunction::Cardinality; + func.named_parameters["open_options"] = LogicalType::LIST(LogicalType::VARCHAR); + func.named_parameters["allowed_drivers"] = LogicalType::LIST(LogicalType::VARCHAR); + func.named_parameters["sibling_files"] = LogicalType::LIST(LogicalType::VARCHAR); + set.AddFunction(func); + + ExtensionUtil::RegisterFunction(db, set); + + DocUtil::AddDocumentation(db, "ST_ReadRaster", DOC_DESCRIPTION, DOC_EXAMPLE, DOC_TAGS); + + // Replacement scan + auto &config = DBConfig::GetConfig(db); + config.replacement_scans.emplace_back(GdalRasterTableFunction::ReplacementScan); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/functions/st_read_raster_meta.cpp b/spatial/src/spatial/gdal/functions/st_read_raster_meta.cpp new file mode 100644 index 00000000..aafd730e --- /dev/null +++ b/spatial/src/spatial/gdal/functions/st_read_raster_meta.cpp @@ -0,0 +1,177 @@ +#include "duckdb/parser/parsed_data/create_table_function_info.hpp" +#include "duckdb/parser/expression/constant_expression.hpp" +#include "duckdb/parser/expression/function_expression.hpp" +#include "duckdb/parser/tableref/table_function_ref.hpp" +#include "duckdb/common/multi_file_reader.hpp" + +#include "spatial/common.hpp" +#include "spatial/gdal/functions.hpp" +#include "spatial/gdal/file_handler.hpp" +#include "spatial/gdal/raster/raster.hpp" + +#include "gdal_priv.h" +#include + +namespace spatial { + +namespace gdal { + +//------------------------------------------------------------------------------ +// Bind +//------------------------------------------------------------------------------ + +struct GDALMetadataBindData : public TableFunctionData { + vector file_names; +}; + +static unique_ptr Bind(ClientContext &context, TableFunctionBindInput &input, + vector &return_types, vector &names) { + + auto file_name = input.inputs[0].GetValue(); + auto result = make_uniq(); + + auto multi_file_reader = MultiFileReader::Create(input.table_function); + result->file_names = + multi_file_reader->CreateFileList(context, input.inputs[0], FileGlobOptions::ALLOW_EMPTY)->GetAllFiles(); + + return_types.push_back(LogicalType::VARCHAR); + return_types.push_back(LogicalType::VARCHAR); + return_types.push_back(LogicalType::VARCHAR); + return_types.push_back(LogicalType::DOUBLE); + return_types.push_back(LogicalType::DOUBLE); + return_types.push_back(LogicalType::INTEGER); + return_types.push_back(LogicalType::INTEGER); + return_types.push_back(LogicalType::DOUBLE); + return_types.push_back(LogicalType::DOUBLE); + return_types.push_back(LogicalType::DOUBLE); + return_types.push_back(LogicalType::DOUBLE); + return_types.push_back(LogicalType::INTEGER); + return_types.push_back(LogicalType::INTEGER); + + names.emplace_back("file_name"); + names.emplace_back("driver_short_name"); + names.emplace_back("driver_long_name"); + names.emplace_back("upper_left_x"); + names.emplace_back("upper_left_y"); + names.emplace_back("width"); + names.emplace_back("height"); + names.emplace_back("scale_x"); + names.emplace_back("scale_y"); + names.emplace_back("skew_x"); + names.emplace_back("skew_y"); + names.emplace_back("srid"); + names.emplace_back("num_bands"); + + return std::move(result); +} + +//------------------------------------------------------------------------------ +// Init +//------------------------------------------------------------------------------ + +struct GDALMetadataState : public GlobalTableFunctionState { + idx_t current_file_idx = 0; +}; + +static unique_ptr Init(ClientContext &context, TableFunctionInitInput &input) { + auto result = make_uniq(); + return std::move(result); +} + +//------------------------------------------------------------------------------ +// Scan +//------------------------------------------------------------------------------ + +static void Scan(ClientContext &context, TableFunctionInput &input, DataChunk &output) { + auto &bind_data = input.bind_data->Cast(); + auto &state = input.global_state->Cast(); + + auto out_size = MinValue(STANDARD_VECTOR_SIZE, bind_data.file_names.size() - state.current_file_idx); + + for (idx_t out_idx = 0; out_idx < out_size; out_idx++, state.current_file_idx++) { + auto file_name = bind_data.file_names[state.current_file_idx]; + auto prefixed_file_name = GDALClientContextState::GetOrCreate(context).GetPrefix(file_name); + + GDALDatasetUniquePtr dataset; + try { + dataset = GDALDatasetUniquePtr( + GDALDataset::Open(prefixed_file_name.c_str(), GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR)); + } catch (...) { + // Just skip anything we cant open + out_idx--; + out_size--; + continue; + } + + Raster raster(dataset.get()); + double gt[6] = {0}; + raster.GetGeoTransform(gt); + + output.data[0].SetValue(out_idx, file_name); + output.data[1].SetValue(out_idx, dataset->GetDriver()->GetDescription()); + output.data[2].SetValue(out_idx, dataset->GetDriver()->GetMetadataItem(GDAL_DMD_LONGNAME)); + output.data[3].SetValue(out_idx, gt[0]); + output.data[4].SetValue(out_idx, gt[3]); + output.data[5].SetValue(out_idx, raster.GetRasterXSize()); + output.data[6].SetValue(out_idx, raster.GetRasterYSize()); + output.data[7].SetValue(out_idx, gt[1]); + output.data[8].SetValue(out_idx, gt[5]); + output.data[9].SetValue(out_idx, gt[2]); + output.data[10].SetValue(out_idx, gt[4]); + output.data[11].SetValue(out_idx, raster.GetSrid()); + output.data[12].SetValue(out_idx, raster.GetRasterCount()); + } + + output.SetCardinality(out_size); +} + +//------------------------------------------------------------------------ +// Documentation +//------------------------------------------------------------------------ + +static constexpr const char *DOC_DESCRIPTION = R"( + The `ST_Read_Meta` table function accompanies the [ST_ReadRaster](#st_readraster) table function, but instead of reading the contents of a file, this function scans the metadata instead. +)"; + +static constexpr const char *DOC_EXAMPLE = R"( + SELECT + driver_short_name, + driver_long_name, + upper_left_x, + upper_left_y, + width, + height, + scale_x, + scale_y, + skew_x, + skew_y, + srid, + num_bands + FROM + ST_ReadRaster_Meta('./test/data/mosaic/SCL.tif-land-clip00.tiff') + ; + + ┌───────────────────┬──────────────────┬──────────────┬──────────────┬───────┬────────┬─────────┬─────────┬────────┬────────┬───────┬───────────┐ + │ driver_short_name │ driver_long_name │ upper_left_x │ upper_left_y │ width │ height │ scale_x │ scale_y │ skew_x │ skew_y │ srid │ num_bands │ + │ varchar │ varchar │ double │ double │ int32 │ int32 │ double │ double │ double │ double │ int32 │ int32 │ + ├───────────────────┼──────────────────┼──────────────┼──────────────┼───────┼────────┼─────────┼─────────┼────────┼────────┼───────┼───────────┤ + │ GTiff │ GeoTIFF │ 541020.0 │ 4796640.0 │ 3438 │ 5322 │ 20.0 │ -20.0 │ 0.0 │ 0.0 │ 32630 │ 1 │ + └───────────────────┴──────────────────┴──────────────┴──────────────┴───────┴────────┴─────────┴─────────┴────────┴────────┴───────┴───────────┘ +)"; + +static constexpr DocTag DOC_TAGS[] = {{"ext", "spatial"}}; + +//------------------------------------------------------------------------------ +// Register +//------------------------------------------------------------------------------ + +void GdalRasterMetadataFunction::Register(DatabaseInstance &db) { + TableFunction func("ST_ReadRaster_Meta", {LogicalType::VARCHAR}, Scan, Bind, Init); + ExtensionUtil::RegisterFunction(db, MultiFileReader::CreateFunctionSet(func)); + + DocUtil::AddDocumentation(db, "ST_ReadRaster_Meta", DOC_DESCRIPTION, DOC_EXAMPLE, DOC_TAGS); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/functions/st_write_raster.cpp b/spatial/src/spatial/gdal/functions/st_write_raster.cpp new file mode 100644 index 00000000..b44de984 --- /dev/null +++ b/spatial/src/spatial/gdal/functions/st_write_raster.cpp @@ -0,0 +1,193 @@ +#include "duckdb/catalog/catalog.hpp" +#include "duckdb/common/types/value.hpp" +#include "duckdb/function/copy_function.hpp" +#include "duckdb/function/table_function.hpp" +#include "duckdb/main/client_context.hpp" +#include "duckdb/main/config.hpp" +#include "duckdb/parser/parsed_data/copy_info.hpp" +#include "duckdb/parser/parsed_data/create_copy_function_info.hpp" +#include "duckdb/parser/parsed_data/create_table_function_info.hpp" +#include "spatial/core/types.hpp" +#include "spatial/gdal/functions.hpp" +#include "spatial/gdal/file_handler.hpp" +#include "spatial/gdal/raster/raster.hpp" +#include "spatial/gdal/raster/raster_factory.hpp" + +#include "gdal_priv.h" + +namespace spatial { + +namespace gdal { + +struct BindRasterData : public TableFunctionData { + + string file_path; + vector field_sql_types; + vector field_names; + string driver_name; + vector creation_options; + + BindRasterData(string file_path, vector field_sql_types, vector field_names) + : file_path(std::move(file_path)), field_sql_types(std::move(field_sql_types)), + field_names(std::move(field_names)) { + } +}; + +struct LocalRasterState : public LocalFunctionData { + + explicit LocalRasterState(ClientContext &context) { + } +}; + +struct GlobalRasterState : public GlobalFunctionData { + + explicit GlobalRasterState(ClientContext &context) { + } +}; + +//===--------------------------------------------------------------------===// +// Bind +//===--------------------------------------------------------------------===// + +static unique_ptr Bind(ClientContext &context, CopyFunctionBindInput &input, const vector &names, + const vector &sql_types) { + + auto bind_data = make_uniq(input.info.file_path, sql_types, names); + + // check all the options in the copy info and set + for (auto &option : input.info.options) { + if (StringUtil::Upper(option.first) == "DRIVER") { + auto set = option.second.front(); + if (set.type().id() == LogicalTypeId::VARCHAR) { + bind_data->driver_name = set.GetValue(); + } else { + throw BinderException("Driver name must be a string"); + } + } else if (StringUtil::Upper(option.first) == "CREATION_OPTIONS") { + auto set = option.second; + for (auto &s : set) { + if (s.type().id() != LogicalTypeId::VARCHAR) { + throw BinderException("Creation options must be strings"); + } + bind_data->creation_options.push_back(s.GetValue()); + } + } else { + throw BinderException("Unknown option '%s'", option.first); + } + } + + if (bind_data->driver_name.empty()) { + throw BinderException("Driver name must be specified"); + } + + auto driver = GetGDALDriverManager()->GetDriverByName(bind_data->driver_name.c_str()); + if (!driver) { + throw BinderException("Unknown driver '%s'", bind_data->driver_name); + } + + // Try get the file extension from the driver + auto file_ext = driver->GetMetadataItem(GDAL_DMD_EXTENSION); + if (file_ext) { + input.file_extension = file_ext; + } else { + // Space separated list of file extensions + auto file_exts = driver->GetMetadataItem(GDAL_DMD_EXTENSIONS); + if (file_exts) { + auto exts = StringUtil::Split(file_exts, ' '); + if (!exts.empty()) { + input.file_extension = exts[0]; + } + } + } + + return std::move(bind_data); +} + +//===--------------------------------------------------------------------===// +// Init Local +//===--------------------------------------------------------------------===// + +static unique_ptr InitLocal(ExecutionContext &context, FunctionData &bind_data) { + + auto local_data = make_uniq(context.client); + return std::move(local_data); +} + +//===--------------------------------------------------------------------===// +// Init Global +//===--------------------------------------------------------------------===// + +static unique_ptr InitGlobal(ClientContext &context, FunctionData &bind_data, + const string &file_path) { + + auto global_data = make_uniq(context); + return std::move(global_data); +} + +//===--------------------------------------------------------------------===// +// Sink +//===--------------------------------------------------------------------===// + +static void Sink(ExecutionContext &context, FunctionData &bdata, GlobalFunctionData &gstate, LocalFunctionData &lstate, + DataChunk &input) { + auto &bind_data = bdata.Cast(); + + // Create the raster + input.Flatten(); + for (idx_t row_idx = 0; row_idx < input.size(); row_idx++) { + + for (idx_t col_idx = 0; col_idx < input.ColumnCount(); col_idx++) { + auto &type = bind_data.field_sql_types[col_idx]; + + if (type == core::GeoTypes::RASTER()) { + auto value = input.GetValue(col_idx, row_idx); + GDALDataset *dataset = reinterpret_cast(value.GetValueUnsafe()); + + auto raw_file_name = bind_data.file_path; + auto &client_ctx = GDALClientContextState::GetOrCreate(context.client); + auto prefixed_file_name = client_ctx.GetPrefix(raw_file_name); + auto driver_name = bind_data.driver_name; + auto creation_options = bind_data.creation_options; + + if (!RasterFactory::WriteFile(dataset, prefixed_file_name, driver_name, creation_options)) { + auto error = Raster::GetLastErrorMsg(); + throw IOException("Could not save file: " + raw_file_name + " (" + error + ")"); + } + break; + } + } + } +} + +//===--------------------------------------------------------------------===// +// Combine +//===--------------------------------------------------------------------===// + +static void Combine(ExecutionContext &context, FunctionData &bind_data, GlobalFunctionData &gstate, + LocalFunctionData &lstate) { +} + +//===--------------------------------------------------------------------===// +// Finalize +//===--------------------------------------------------------------------===// + +static void Finalize(ClientContext &context, FunctionData &bind_data, GlobalFunctionData &gstate) { +} + +void GdalRasterCopyFunction::Register(DatabaseInstance &db) { + // register the copy function + CopyFunction info("RASTER"); + info.copy_to_bind = Bind; + info.copy_to_initialize_local = InitLocal; + info.copy_to_initialize_global = InitGlobal; + info.copy_to_sink = Sink; + info.copy_to_combine = Combine; + info.copy_to_finalize = Finalize; + info.extension = "raster"; + + ExtensionUtil::RegisterFunction(db, info); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/module.cpp b/spatial/src/spatial/gdal/module.cpp index e0cee7ed..a8761e1c 100644 --- a/spatial/src/spatial/gdal/module.cpp +++ b/spatial/src/spatial/gdal/module.cpp @@ -1,5 +1,8 @@ #include "spatial/gdal/module.hpp" #include "spatial/gdal/functions.hpp" +#include "spatial/gdal/functions/aggregate.hpp" +#include "spatial/gdal/functions/cast.hpp" +#include "spatial/gdal/functions/scalar.hpp" #include "spatial/gdal/file_handler.hpp" #include "spatial/common.hpp" @@ -17,7 +20,7 @@ void GdalModule::Register(DatabaseInstance &db) { static std::once_flag loaded; std::call_once(loaded, [&]() { // Register all embedded drivers (dont go looking for plugins) - OGRRegisterAllInternal(); + GDALAllRegister(); // Set GDAL error handler @@ -62,9 +65,15 @@ void GdalModule::Register(DatabaseInstance &db) { // Register functions GdalTableFunction::Register(db); + GdalRasterTableFunction::Register(db); + GdalRasterMetadataFunction::Register(db); GdalDriversTableFunction::Register(db); GdalCopyFunction::Register(db); GdalMetadataFunction::Register(db); + GdalCastFunctions::Register(db); + GdalScalarFunctions::Register(db); + GdalAggregateFunctions::Register(db); + GdalRasterCopyFunction::Register(db); } } // namespace gdal diff --git a/spatial/src/spatial/gdal/raster/CMakeLists.txt b/spatial/src/spatial/gdal/raster/CMakeLists.txt new file mode 100644 index 00000000..74114cd6 --- /dev/null +++ b/spatial/src/spatial/gdal/raster/CMakeLists.txt @@ -0,0 +1,8 @@ +set(EXTENSION_SOURCES + ${EXTENSION_SOURCES} + ${CMAKE_CURRENT_SOURCE_DIR}/raster.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/raster_factory.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/raster_registry.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/raster_value.cpp + PARENT_SCOPE +) diff --git a/spatial/src/spatial/gdal/raster/raster.cpp b/spatial/src/spatial/gdal/raster/raster.cpp new file mode 100644 index 00000000..ee5a7553 --- /dev/null +++ b/spatial/src/spatial/gdal/raster/raster.cpp @@ -0,0 +1,348 @@ +#include "duckdb/common/types/uuid.hpp" +#include "spatial/core/types.hpp" +#include "spatial/core/geometry/wkb_writer.hpp" +#include "spatial/core/util/math.hpp" +#include "spatial/gdal/types.hpp" +#include "spatial/gdal/raster/raster.hpp" + +#include "gdal_priv.h" +#include "gdal_utils.h" +#include "gdalwarper.h" +#include /* for FLT_EPSILON */ + +using namespace spatial::core; + +namespace spatial { + +namespace gdal { + +Raster::Raster(GDALDataset *dataset) : dataset_(dataset) { +} + +int Raster::GetRasterXSize() const { + return dataset_->GetRasterXSize(); +} + +int Raster::GetRasterYSize() const { + return dataset_->GetRasterYSize(); +} + +int Raster::GetRasterCount() const { + return dataset_->GetRasterCount(); +} + +int32_t Raster::GetSrid() const { + + int32_t srid = 0; // SRID_UNKNOWN + + const char *proj_def = dataset_->GetProjectionRef(); + if (proj_def) { + OGRSpatialReference spatial_ref; + spatial_ref.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); + + if (spatial_ref.importFromWkt(proj_def) == OGRERR_NONE && spatial_ref.AutoIdentifyEPSG() == OGRERR_NONE) { + + const char *code = spatial_ref.GetAuthorityCode(nullptr); + if (code) { + srid = atoi(code); + } + } + } + return srid; +} + +bool Raster::GetGeoTransform(double *matrix) const { + + if (dataset_->GetGeoTransform(matrix) != CE_None) { + // Using default geotransform matrix (0, 1, 0, 0, 0, -1) + matrix[0] = 0; + matrix[1] = 1; + matrix[2] = 0; + matrix[3] = 0; + matrix[4] = 0; + matrix[5] = -1; + return false; + } + return true; +} + +bool Raster::GetInvGeoTransform(double *inv_matrix) const { + double gt[6] = {0}; + GetGeoTransform(gt); + + if (!GDALInvGeoTransform(gt, inv_matrix)) { + return false; + } + return true; +} + +static VertexXY rasterToWorldVertex(double matrix[], int32_t col, int32_t row) { + double xgeo = matrix[0] + matrix[1] * col + matrix[2] * row; + double ygeo = matrix[3] + matrix[4] * col + matrix[5] * row; + return VertexXY {xgeo, ygeo}; +} + +Geometry Raster::GetGeometry(ArenaAllocator &allocator) const { + auto cols = dataset_->GetRasterXSize(); + auto rows = dataset_->GetRasterYSize(); + + double gt[6] = {0}; + GetGeoTransform(gt); + + VertexXY vertex1 = rasterToWorldVertex(gt, 0, 0); + VertexXY vertex2 = rasterToWorldVertex(gt, cols, rows); + double minx = std::min(vertex1.x, vertex2.x); + double miny = std::min(vertex1.y, vertex2.y); + double maxx = std::max(vertex1.x, vertex2.x); + double maxy = std::max(vertex1.y, vertex2.y); + + return Polygon::CreateFromBox(allocator, minx, miny, maxx, maxy); +} + +bool Raster::RasterToWorldCoord(PointXY &point, int32_t col, int32_t row) const { + double gt[6] = {0}; + GetGeoTransform(gt); + return Raster::RasterToWorldCoord(point, gt, col, row); +} + +bool Raster::RasterToWorldCoord(PointXY &point, double matrix[], int32_t col, int32_t row) { + point.x = matrix[0] + matrix[1] * col + matrix[2] * row; + point.y = matrix[3] + matrix[4] * col + matrix[5] * row; + return true; +} + +bool Raster::WorldToRasterCoord(RasterCoord &coord, double x, double y) const { + double inv_gt[6] = {0}; + + if (GetInvGeoTransform(inv_gt)) { + return Raster::WorldToRasterCoord(coord, inv_gt, x, y); + } + return false; +} + +bool Raster::WorldToRasterCoord(RasterCoord &coord, double inv_matrix[], double x, double y) { + double xr = 0, yr = 0; + GDALApplyGeoTransform(inv_matrix, x, y, &xr, &yr); + + // Source: + // https://github.com/postgis/postgis/blob/stable-3.4/raster/rt_core/rt_raster.c#L808 + double rnd = 0; + +// Helper macro for symmetrical rounding +#define ROUND(x, y) (((x > 0.0) ? floor((x * pow(10, y) + 0.5)) : ceil((x * pow(10, y) - 0.5))) / pow(10, y)) +// Helper macro for consistent floating point equality checks +#define FLT_EQ(x, y) ((x == y) || (isnan(x) && isnan(y)) || (fabs(x - y) <= FLT_EPSILON)) + + rnd = ROUND(xr, 0); + if (FLT_EQ(rnd, xr)) + xr = rnd; + else + xr = floor(xr); + + rnd = ROUND(yr, 0); + if (FLT_EQ(rnd, yr)) + yr = rnd; + else + yr = floor(yr); + + coord.col = (int32_t)xr; + coord.row = (int32_t)yr; + return true; +} + +bool Raster::GetValue(double &value, int32_t band_num, int32_t col, int32_t row) const { + + GDALRasterBand *raster_band = dataset_->GetRasterBand(band_num); + double pixel_value = raster_band->GetNoDataValue(); + + if (raster_band->RasterIO(GF_Read, col, row, 1, 1, &pixel_value, 1, 1, GDT_Float64, 0, 0) == CE_None) { + value = pixel_value; + return true; + } + return false; +} + +GDALDataset *Raster::BuildVRT(const std::vector &datasets, const std::vector &options) { + + char **papszArgv = nullptr; + + for (auto it = options.begin(); it != options.end(); ++it) { + papszArgv = CSLAddString(papszArgv, (*it).c_str()); + } + + CPLErrorReset(); + + GDALBuildVRTOptions *psOptions = GDALBuildVRTOptionsNew(papszArgv, nullptr); + CSLDestroy(papszArgv); + + auto result = GDALDatasetUniquePtr(GDALDataset::FromHandle( + GDALBuildVRT(nullptr, datasets.size(), (GDALDatasetH *)&datasets[0], nullptr, psOptions, nullptr))); + + GDALBuildVRTOptionsFree(psOptions); + + if (result.get() != nullptr) { + result->FlushCache(); + } + return result.release(); +} + +GDALDataset *Raster::Warp(GDALDataset *dataset, const std::vector &options) { + + GDALDatasetH hDataset = GDALDataset::ToHandle(dataset); + + auto driver = GetGDALDriverManager()->GetDriverByName("MEM"); + if (!driver) { + throw InvalidInputException("Unknown driver 'MEM'"); + } + + char **papszArgv = nullptr; + papszArgv = CSLAddString(papszArgv, "-of"); + papszArgv = CSLAddString(papszArgv, "MEM"); + + for (auto it = options.begin(); it != options.end(); ++it) { + papszArgv = CSLAddString(papszArgv, (*it).c_str()); + } + + CPLErrorReset(); + + GDALWarpAppOptions *psOptions = GDALWarpAppOptionsNew(papszArgv, nullptr); + CSLDestroy(papszArgv); + + auto ds_name = UUID::ToString(UUID::GenerateRandomUUID()); + + auto result = GDALDatasetUniquePtr( + GDALDataset::FromHandle(GDALWarp(ds_name.c_str(), nullptr, 1, &hDataset, psOptions, nullptr))); + + GDALWarpAppOptionsFree(psOptions); + + if (result.get() != nullptr) { + result->FlushCache(); + } + return result.release(); +} + +//! Transformer of Geometries to pixel/line coordinates +class CutlineTransformer : public OGRCoordinateTransformation { +public: + void *hTransformArg = nullptr; + + explicit CutlineTransformer(void *hTransformArg) : hTransformArg(hTransformArg) { + } + virtual ~CutlineTransformer() { + GDALDestroyTransformer(hTransformArg); + } + + virtual const OGRSpatialReference *GetSourceCS() const override { + return nullptr; + } + virtual const OGRSpatialReference *GetTargetCS() const override { + return nullptr; + } + virtual OGRCoordinateTransformation *Clone() const override { + return nullptr; + } + virtual OGRCoordinateTransformation *GetInverse() const override { + return nullptr; + } + + virtual int Transform(int nCount, double *x, double *y, double *z, double * /* t */, int *pabSuccess) override { + return GDALGenImgProjTransform(hTransformArg, TRUE, nCount, x, y, z, pabSuccess); + } +}; + +GDALDataset *Raster::Clip(GDALDataset *dataset, const geometry_t &geometry, const std::vector &options) { + + GDALDatasetH hDataset = GDALDataset::ToHandle(dataset); + + auto driver = GetGDALDriverManager()->GetDriverByName("MEM"); + if (!driver) { + throw InvalidInputException("Unknown driver 'MEM'"); + } + + char **papszArgv = nullptr; + papszArgv = CSLAddString(papszArgv, "-of"); + papszArgv = CSLAddString(papszArgv, "MEM"); + + for (auto it = options.begin(); it != options.end(); ++it) { + papszArgv = CSLAddString(papszArgv, (*it).c_str()); + } + + // Add Bounds & Geometry in pixel/line coordinates to the options. + if (geometry.GetType() == GeometryType::POLYGON || geometry.GetType() == GeometryType::MULTIPOLYGON) { + + OGRGeometryUniquePtr ogr_geom; + + OGRSpatialReference srs; + srs.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER); + const char *proj_ref = dataset->GetProjectionRef(); + if (proj_ref) { + srs.importFromWkt(&proj_ref, nullptr); + } + + vector buffer; + WKBWriter::Write(geometry, buffer); + + OGRGeometry *ptr_geom = nullptr; + if (OGRGeometryFactory::createFromWkb(buffer.data(), &srs, &ptr_geom, buffer.size(), wkbVariantIso) != + OGRERR_NONE) { + + CSLDestroy(papszArgv); + throw InvalidInputException("Input Geometry could not imported"); + } else { + ogr_geom = OGRGeometryUniquePtr(ptr_geom); + } + + OGREnvelope envelope; + ogr_geom->getEnvelope(&envelope); + + CutlineTransformer transformer(GDALCreateGenImgProjTransformer2(hDataset, nullptr, nullptr)); + + if (ogr_geom->transform(&transformer) != OGRERR_NONE) { + CSLDestroy(papszArgv); + throw InvalidInputException("Transform of geometry to pixel/line coordinates failed"); + } + + char *pszWkt = nullptr; + if (ogr_geom->exportToWkt(&pszWkt) != OGRERR_NONE) { + CSLDestroy(papszArgv); + CPLFree(pszWkt); + throw InvalidInputException("Input Geometry could not loaded"); + } + std::string wkt_geom = pszWkt; + CPLFree(pszWkt); + + std::string wkt_option = "CUTLINE=" + wkt_geom; + papszArgv = CSLAddString(papszArgv, "-wo"); + papszArgv = CSLAddString(papszArgv, wkt_option.c_str()); + papszArgv = CSLAddString(papszArgv, "-te"); + papszArgv = CSLAddString(papszArgv, MathUtil::format_coord(envelope.MinX).c_str()); + papszArgv = CSLAddString(papszArgv, MathUtil::format_coord(envelope.MinY).c_str()); + papszArgv = CSLAddString(papszArgv, MathUtil::format_coord(envelope.MaxX).c_str()); + papszArgv = CSLAddString(papszArgv, MathUtil::format_coord(envelope.MaxY).c_str()); + } + + CPLErrorReset(); + + GDALWarpAppOptions *psOptions = GDALWarpAppOptionsNew(papszArgv, nullptr); + CSLDestroy(papszArgv); + + auto ds_name = UUID::ToString(UUID::GenerateRandomUUID()); + + auto result = GDALDatasetUniquePtr( + GDALDataset::FromHandle(GDALWarp(ds_name.c_str(), nullptr, 1, &hDataset, psOptions, nullptr))); + + GDALWarpAppOptionsFree(psOptions); + + if (result.get() != nullptr) { + result->FlushCache(); + } + return result.release(); +} + +string Raster::GetLastErrorMsg() { + return string(CPLGetLastErrorMsg()); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/raster/raster_factory.cpp b/spatial/src/spatial/gdal/raster/raster_factory.cpp new file mode 100644 index 00000000..6e599d46 --- /dev/null +++ b/spatial/src/spatial/gdal/raster/raster_factory.cpp @@ -0,0 +1,131 @@ +#include "spatial/common.hpp" +#include "spatial/gdal/raster/raster_factory.hpp" + +#include "gdal_priv.h" + +namespace spatial { + +namespace gdal { + +GDALDataset *RasterFactory::FromFile(const std::string &file_path, const std::vector &allowed_drivers, + const std::vector &open_options, + const std::vector &sibling_files) { + + auto gdal_allowed_drivers = RasterFactory::FromVectorOfStrings(allowed_drivers); + auto gdal_open_options = RasterFactory::FromVectorOfStrings(open_options); + auto gdal_sibling_files = RasterFactory::FromVectorOfStrings(sibling_files); + + GDALDataset *dataset = GDALDataset::Open(file_path.c_str(), GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR, + gdal_allowed_drivers.empty() ? nullptr : gdal_allowed_drivers.data(), + gdal_open_options.empty() ? nullptr : gdal_open_options.data(), + gdal_sibling_files.empty() ? nullptr : gdal_sibling_files.data()); + + return dataset; +} + +bool RasterFactory::WriteFile(GDALDataset *dataset, const std::string &file_path, const std::string &driver_name, + const std::vector &write_options) { + + auto driver = GetGDALDriverManager()->GetDriverByName(driver_name.c_str()); + if (!driver) { + throw InvalidInputException("Unknown driver '%s'", driver_name.c_str()); + } + + bool copy_available = CSLFetchBoolean(driver->GetMetadata(), GDAL_DCAP_CREATECOPY, FALSE); + auto gdal_write_options = RasterFactory::FromVectorOfStrings(write_options); + auto gdal_options = gdal_write_options.empty() ? nullptr : gdal_write_options.data(); + + GDALDatasetUniquePtr output; + CPLErrorReset(); + + if (copy_available) { + output = GDALDatasetUniquePtr(driver->CreateCopy(file_path.c_str(), dataset, FALSE, gdal_options, NULL, NULL)); + + if (output.get() == nullptr) { + return false; + } + } else { + int cols = dataset->GetRasterXSize(); + int rows = dataset->GetRasterYSize(); + int band_count = dataset->GetRasterCount(); + + if (band_count == 0) { + throw InvalidInputException("Input Raster has no RasterBands"); + } + + GDALRasterBand *raster_band = dataset->GetRasterBand(1); + GDALDataType data_type = raster_band->GetRasterDataType(); + int date_type_size = GDALGetDataTypeSize(data_type); + + output = + GDALDatasetUniquePtr(driver->Create(file_path.c_str(), cols, rows, band_count, data_type, gdal_options)); + + if (output.get() == nullptr) { + return false; + } + + double gt[6] = {0, 1, 0, 0, 0, -1}; + dataset->GetGeoTransform(gt); + output->SetGeoTransform(gt); + + output->SetProjection(dataset->GetProjectionRef()); + output->SetMetadata(dataset->GetMetadata()); + + void *pafScanline = CPLMalloc(date_type_size * cols * rows); + + for (int i = 1; i <= band_count; i++) { + GDALRasterBand *source_band = dataset->GetRasterBand(i); + GDALRasterBand *target_band = output->GetRasterBand(i); + + target_band->SetMetadata(source_band->GetMetadata()); + target_band->SetNoDataValue(source_band->GetNoDataValue()); + target_band->SetColorInterpretation(source_band->GetColorInterpretation()); + + if (source_band->RasterIO(GF_Read, 0, 0, cols, rows, pafScanline, cols, rows, data_type, 0, 0) != CE_None || + target_band->RasterIO(GF_Write, 0, 0, cols, rows, pafScanline, cols, rows, data_type, 0, 0) != + CE_None) { + CPLFree(pafScanline); + return false; + } + } + + CPLFree(pafScanline); + } + output->FlushCache(); + + return true; +} + +std::vector RasterFactory::FromVectorOfStrings(const std::vector &input) { + auto output = std::vector(); + + if (input.size()) { + output.reserve(input.size() + 1); + + for (auto it = input.begin(); it != input.end(); ++it) { + output.push_back((*it).c_str()); + } + output.push_back(nullptr); + } + return output; +} + +std::vector RasterFactory::FromNamedParameters(const named_parameter_map_t &input, + const std::string &keyname) { + auto output = std::vector(); + + auto input_param = input.find(keyname); + if (input_param != input.end()) { + output.reserve(input.size() + 1); + + for (auto ¶m : ListValue::GetChildren(input_param->second)) { + output.push_back(StringValue::Get(param).c_str()); + } + output.push_back(nullptr); + } + return output; +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/raster/raster_registry.cpp b/spatial/src/spatial/gdal/raster/raster_registry.cpp new file mode 100644 index 00000000..a5dea27d --- /dev/null +++ b/spatial/src/spatial/gdal/raster/raster_registry.cpp @@ -0,0 +1,31 @@ +#include "spatial/core/types.hpp" +#include "spatial/gdal/raster/raster.hpp" +#include "spatial/gdal/raster/raster_registry.hpp" + +#include "gdal_priv.h" + +using namespace spatial::core; + +namespace spatial { + +namespace gdal { + +RasterRegistry::RasterRegistry() { +} + +RasterRegistry::~RasterRegistry() { + + // Release items in reverse order, first children, then parent ones + for (auto it = datasets_.rbegin(); it != datasets_.rend(); ++it) { + auto datasetUniquePtr = std::move(*it); + datasetUniquePtr.reset(); + } +} + +void RasterRegistry::RegisterRaster(GDALDataset *dataset) { + datasets_.emplace_back(GDALDatasetUniquePtr(dataset)); +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/raster/raster_value.cpp b/spatial/src/spatial/gdal/raster/raster_value.cpp new file mode 100644 index 00000000..9695597a --- /dev/null +++ b/spatial/src/spatial/gdal/raster/raster_value.cpp @@ -0,0 +1,29 @@ +#include "spatial/core/types.hpp" +#include "spatial/gdal/raster/raster.hpp" +#include "spatial/gdal/raster/raster_value.hpp" + +using namespace spatial::core; + +namespace spatial { + +namespace gdal { + +Value RasterValue::CreateValue(GDALDataset *dataset) { + Value value = Value::POINTER(CastPointerToValue(dataset)); + value.Reinterpret(GeoTypes::RASTER()); + return value; +} + +GDALDataset *RasterValue::operator->() const { + GDALDataset *dataset = reinterpret_cast(GetValueUnsafe()); + return dataset; +} + +GDALDataset *RasterValue::get() const { + GDALDataset *dataset = reinterpret_cast(GetValueUnsafe()); + return dataset; +} + +} // namespace gdal + +} // namespace spatial diff --git a/spatial/src/spatial/gdal/types.cpp b/spatial/src/spatial/gdal/types.cpp new file mode 100644 index 00000000..15dab0c8 --- /dev/null +++ b/spatial/src/spatial/gdal/types.cpp @@ -0,0 +1,124 @@ +#include "spatial/gdal/types.hpp" + +#include "gdal_priv.h" + +namespace spatial { + +namespace gdal { + +std::string GetPixelTypeName(const PixelType &pixel_type) { + switch (pixel_type) { + case Byte: + return "Byte"; + case Int8: + return "Int8"; + case UInt16: + return "UInt16"; + case Int16: + return "Int16"; + case UInt32: + return "UInt32"; + case Int32: + return "Int32"; + case UInt64: + return "UInt64"; + case Int64: + return "Int64"; + case Float32: + return "Float32"; + case Float64: + return "Float64"; + case CInt16: + return "CInt16"; + case CInt32: + return "CInt32"; + case CFloat32: + return "CFloat32"; + case CFloat64: + return "CFloat64"; + case TypeCount: + case Unknown: + default: + return "Unknown"; + } +} + +std::string GetColorInterpName(const ColorInterp &color_interp) { + switch (color_interp) { + case Undefined: + return "Undefined"; + case GrayIndex: + return "Greyscale"; + case PaletteIndex: + return "Paletted"; + case RedBand: + return "Red"; + case GreenBand: + return "Green"; + case BlueBand: + return "Blue"; + case AlphaBand: + return "Alpha"; + case HueBand: + return "Hue"; + case SaturationBand: + return "Saturation"; + case LightnessBand: + return "Lightness"; + case CyanBand: + return "Cyan"; + case MagentaBand: + return "Magenta"; + case YellowBand: + return "Yellow"; + case BlackBand: + return "Black"; + case YCbCr_YBand: + return "YLuminance"; + case YCbCr_CbBand: + return "CbChroma"; + case YCbCr_CrBand: + return "CrChroma"; + default: + return "Unknown"; + } +} + +std::string GetResampleAlgName(const ResampleAlg &resample_alg) { + switch (resample_alg) { + case NearestNeighbour: + return "NearestNeighbour"; + case Bilinear: + return "Bilinear"; + case Cubic: + return "Cubic"; + case CubicSpline: + return "CubicSpline"; + case Lanczos: + return "Lanczos"; + case Average: + return "Average"; + case Mode: + return "Mode"; + case Max: + return "Maximum"; + case Min: + return "Minimun"; + case Med: + return "Median"; + case Q1: + return "Quartile1"; + case Q3: + return "Quartile3"; + case Sum: + return "Sum"; + case RMS: + return "RootMeanSquare"; + default: + return "Unknown"; + } +} + +} // namespace gdal + +} // namespace spatial diff --git a/test/data/mosaic/CATAST_Pol_Township-PNA.gpkg b/test/data/mosaic/CATAST_Pol_Township-PNA.gpkg new file mode 100644 index 00000000..edbe669e Binary files /dev/null and b/test/data/mosaic/CATAST_Pol_Township-PNA.gpkg differ diff --git a/test/data/mosaic/SCL.tif-land-clip00.tiff b/test/data/mosaic/SCL.tif-land-clip00.tiff new file mode 100644 index 00000000..c99014e9 Binary files /dev/null and b/test/data/mosaic/SCL.tif-land-clip00.tiff differ diff --git a/test/data/mosaic/SCL.tif-land-clip01.tiff b/test/data/mosaic/SCL.tif-land-clip01.tiff new file mode 100644 index 00000000..ec681c7d Binary files /dev/null and b/test/data/mosaic/SCL.tif-land-clip01.tiff differ diff --git a/test/data/mosaic/SCL.tif-land-clip10.tiff b/test/data/mosaic/SCL.tif-land-clip10.tiff new file mode 100644 index 00000000..1b34c780 Binary files /dev/null and b/test/data/mosaic/SCL.tif-land-clip10.tiff differ diff --git a/test/data/mosaic/SCL.tif-land-clip11.tiff b/test/data/mosaic/SCL.tif-land-clip11.tiff new file mode 100644 index 00000000..e5c6716f Binary files /dev/null and b/test/data/mosaic/SCL.tif-land-clip11.tiff differ diff --git a/test/sql/gdal/st_read_raster.test b/test/sql/gdal/st_read_raster.test new file mode 100644 index 00000000..a33ff02d --- /dev/null +++ b/test/sql/gdal/st_read_raster.test @@ -0,0 +1,31 @@ +require spatial + +query I +SELECT raster FROM ST_ReadRaster('__WORKING_DIRECTORY__/test/data/mosaic/SCL.tif-land-clip00.tiff'); +---- +RASTER + +query I +SELECT raster FROM '__WORKING_DIRECTORY__/test/data/mosaic/SCL.tif-land-clip00.tiff'; +---- +RASTER + +query IIIIIIIIIIII +SELECT + driver_short_name, + driver_long_name, + upper_left_x, + upper_left_y, + width, + height, + scale_x, + scale_y, + skew_x, + skew_y, + srid, + num_bands +FROM + ST_ReadRaster_Meta('__WORKING_DIRECTORY__/test/data/mosaic/SCL.tif-land-clip00.tiff') +; +---- +GTiff GeoTIFF 541020.0 4796640.0 3438 5322 20.0 -20.0 0.0 0.0 32630 1 diff --git a/test/sql/gdal/st_write_raster.test b/test/sql/gdal/st_write_raster.test new file mode 100644 index 00000000..f247ecea --- /dev/null +++ b/test/sql/gdal/st_write_raster.test @@ -0,0 +1,11 @@ +require spatial + +statement ok +COPY ( + SELECT * FROM '__WORKING_DIRECTORY__/test/data/mosaic/SCL.tif-land-clip00.tiff' +) +TO + '__TEST_DIR__/copytoraster.tiff' +WITH ( + FORMAT RASTER, DRIVER 'COG', CREATION_OPTIONS ('COMPRESS=LZW') +); diff --git a/test/sql/raster/st_getbandmetadata.test b/test/sql/raster/st_getbandmetadata.test new file mode 100644 index 00000000..717758b8 --- /dev/null +++ b/test/sql/raster/st_getbandmetadata.test @@ -0,0 +1,17 @@ +require spatial + +query IIIIIIII +SELECT + ST_HasNoBand(raster, 1) AS hb1, + ST_HasNoBand(raster, 2) AS hb2, + ST_HasNoBand(raster, 3) AS hb3, + ST_GetBandPixelType(raster, 1) AS pxt1, + ST_GetBandPixelTypeName(raster, 1) AS pxtn1, + ST_GetBandNoDataValue(raster, 1) AS ndv1, + ST_GetBandColorInterp(raster, 1) AS ci1, + ST_GetBandColorInterpName(raster, 1) AS cin1 +FROM + ST_ReadRaster('__WORKING_DIRECTORY__/test/data/mosaic/SCL.tif-land-clip00.tiff') +; +---- +1 0 0 3 Int16 -9999.0 1 Greyscale diff --git a/test/sql/raster/st_getgeometry.test b/test/sql/raster/st_getgeometry.test new file mode 100644 index 00000000..8b33618e --- /dev/null +++ b/test/sql/raster/st_getgeometry.test @@ -0,0 +1,11 @@ +require spatial + +query I +SELECT ST_GetGeometry(raster) FROM ST_ReadRaster('__WORKING_DIRECTORY__/test/data/mosaic/SCL.tif-land-clip00.tiff'); +---- +POLYGON ((541020 4690200, 541020 4796640, 609780 4796640, 609780 4690200, 541020 4690200)) + +query I +SELECT raster::geometry FROM ST_ReadRaster('__WORKING_DIRECTORY__/test/data/mosaic/SCL.tif-land-clip00.tiff'); +---- +POLYGON ((541020 4690200, 541020 4796640, 609780 4796640, 609780 4690200, 541020 4690200)) diff --git a/test/sql/raster/st_getmetadata.test b/test/sql/raster/st_getmetadata.test new file mode 100644 index 00000000..b6126c8a --- /dev/null +++ b/test/sql/raster/st_getmetadata.test @@ -0,0 +1,20 @@ +require spatial + +query IIIIIIIIIII +SELECT + ST_Width(raster) AS cols, + ST_Height(raster) AS rows, + ST_NumBands(raster) AS num_bands, + ST_UpperLeftX(raster) AS ulx, + ST_UpperLeftY(raster) AS uly, + ST_ScaleX(raster) AS scale_x, + ST_ScaleY(raster) AS scale_y, + ST_SkewX(raster) AS skew_x, + ST_SkewY(raster) AS skew_y, + ST_PixelWidth(raster) AS px_width, + ST_PixelHeight(raster) AS px_height +FROM + ST_ReadRaster('__WORKING_DIRECTORY__/test/data/mosaic/SCL.tif-land-clip00.tiff') +; +---- +3438 5322 1 541020.0 4796640.0 20.0 -20.0 0.0 0.0 20.0 20.0 diff --git a/test/sql/raster/st_getvalue.test b/test/sql/raster/st_getvalue.test new file mode 100644 index 00000000..da1bf077 --- /dev/null +++ b/test/sql/raster/st_getvalue.test @@ -0,0 +1,14 @@ +require spatial + +query IIIII +SELECT + ST_Value(raster, 1, (ST_Width(raster) / 2)::INT, (ST_Height(raster) / 2)::INT) AS valCC, + ST_Value(raster, 1, 0, 0) AS val00, + ST_Value(raster, 1, ST_Width(raster) - 1, 0) AS val10, + ST_Value(raster, 1, ST_Width(raster) - 1, ST_Height(raster) - 1) AS val11, + ST_Value(raster, 1, 0, ST_Height(raster) - 1) AS val01 +FROM + ST_ReadRaster('__WORKING_DIRECTORY__/test/data/mosaic/SCL.tif-land-clip00.tiff') +; +---- +1.0 -9999.0 -9999.0 15.0 -9999.0 diff --git a/test/sql/raster/st_rasterasfile.test b/test/sql/raster/st_rasterasfile.test new file mode 100644 index 00000000..dc5c6b37 --- /dev/null +++ b/test/sql/raster/st_rasterasfile.test @@ -0,0 +1,14 @@ +require spatial + +statement ok +WITH __input AS ( + SELECT + ST_RasterFromFile(file) AS raster + FROM + glob('__WORKING_DIRECTORY__/test/data/mosaic/SCL.tif-land-clip00.tiff') +) +SELECT + ST_RasterAsFile(raster, '__TEST_DIR__/rasterasfile.tiff', 'Gtiff') AS result +FROM + __input +; diff --git a/test/sql/raster/st_rasterclip.test b/test/sql/raster/st_rasterclip.test new file mode 100644 index 00000000..d0129872 --- /dev/null +++ b/test/sql/raster/st_rasterclip.test @@ -0,0 +1,40 @@ +require spatial + +query I +WITH __input AS ( + SELECT + 1 AS mosaic_id, + ST_RasterFromFile(file) AS raster + FROM + glob('__WORKING_DIRECTORY__/test/data/mosaic/*.tiff') +), +__mosaic AS ( + SELECT + ST_RasterMosaic_Agg(raster, options => ['-r', 'bilinear']) AS mosaic + FROM + __input + GROUP BY + mosaic_id +), +__geometry AS ( + SELECT geom FROM ST_Read('/__WORKING_DIRECTORY__/test/data/mosaic/CATAST_Pol_Township-PNA.gpkg') +), +__clip AS ( + SELECT + ST_RasterClip(mosaic, + (SELECT geom FROM __geometry LIMIT 1), + options => + [ + '-r', 'bilinear', '-crop_to_cutline', '-wo', 'CUTLINE_ALL_TOUCHED=TRUE' + ] + ) AS clip + FROM + __mosaic +) +SELECT + ST_Area(ST_GetGeometry(clip)) AS result +FROM + __clip +; +---- +44269454.49488351 diff --git a/test/sql/raster/st_rasterfromfile.test b/test/sql/raster/st_rasterfromfile.test new file mode 100644 index 00000000..6947ccf8 --- /dev/null +++ b/test/sql/raster/st_rasterfromfile.test @@ -0,0 +1,16 @@ +require spatial + +query I +WITH __input AS ( + SELECT + ST_RasterFromFile(file) AS raster + FROM + glob('__WORKING_DIRECTORY__/test/data/mosaic/SCL.tif-land-clip00.tiff') +) +SELECT + ST_GetGeometry(raster) AS geom +FROM + __input +; +---- +POLYGON ((541020 4690200, 541020 4796640, 609780 4796640, 609780 4690200, 541020 4690200)) diff --git a/test/sql/raster/st_rastermosaic.test b/test/sql/raster/st_rastermosaic.test new file mode 100644 index 00000000..2b7db98e --- /dev/null +++ b/test/sql/raster/st_rastermosaic.test @@ -0,0 +1,25 @@ +require spatial + +query I +WITH __input AS ( + SELECT + 1 AS mosaic_id, + ST_RasterFromFile(file) AS raster + FROM + glob('__WORKING_DIRECTORY__/test/data/mosaic/*.tiff') +), +__mosaic AS ( + SELECT + ST_RasterMosaic_Agg(raster, options => ['-r', 'bilinear']) AS mosaic + FROM + __input + GROUP BY + mosaic_id +) +SELECT + ST_GetGeometry(mosaic) AS result +FROM + __mosaic +; +---- +POLYGON ((541020 4640780, 541020 4796640, 685600 4796640, 685600 4640780, 541020 4640780)) diff --git a/test/sql/raster/st_rasterwarp.test b/test/sql/raster/st_rasterwarp.test new file mode 100644 index 00000000..5e8f0476 --- /dev/null +++ b/test/sql/raster/st_rasterwarp.test @@ -0,0 +1,23 @@ +require spatial + +query II +WITH __input AS ( + SELECT + raster + FROM + ST_ReadRaster('__WORKING_DIRECTORY__/test/data/mosaic/SCL.tif-land-clip00.tiff') +), +__warp AS ( + SELECT + ST_RasterWarp(raster, options => ['-r', 'bilinear', '-tr', '40.0', '40.0']) AS warp + FROM + __input +) +SELECT + ST_ScaleX(warp) AS scale_x, + ST_ScaleY(warp) AS scale_y +FROM + __warp +; +---- +40.0 -40.0 diff --git a/test/sql/raster/st_srid.test b/test/sql/raster/st_srid.test new file mode 100644 index 00000000..887a8962 --- /dev/null +++ b/test/sql/raster/st_srid.test @@ -0,0 +1,10 @@ +require spatial + +query I +SELECT + ST_Srid(raster) AS srid +FROM + ST_ReadRaster('__WORKING_DIRECTORY__/test/data/mosaic/SCL.tif-land-clip00.tiff') +; +---- +32630 diff --git a/test/sql/raster/st_transformcoords.test b/test/sql/raster/st_transformcoords.test new file mode 100644 index 00000000..d27efc9a --- /dev/null +++ b/test/sql/raster/st_transformcoords.test @@ -0,0 +1,15 @@ +require spatial + +query IIIIII +SELECT + ST_RasterToWorldCoord (raster, 0, 0) AS coord, + ST_RasterToWorldCoordX(raster, 0, 0) AS coord_x, + ST_RasterToWorldCoordY(raster, 0, 0) AS coord_y, + ST_WorldToRasterCoord (raster, ST_RasterToWorldCoordX(raster, 1, 2), ST_RasterToWorldCoordY(raster, 1, 2)) AS coord, + ST_WorldToRasterCoordX(raster, ST_RasterToWorldCoordX(raster, 1, 2), ST_RasterToWorldCoordY(raster, 1, 2)) AS col, + ST_WorldToRasterCoordY(raster, ST_RasterToWorldCoordX(raster, 1, 2), ST_RasterToWorldCoordY(raster, 1, 2)) AS row +FROM + ST_ReadRaster('__WORKING_DIRECTORY__/test/data/mosaic/SCL.tif-land-clip00.tiff') +; +---- +POINT (541020 4796640) 541020.0 4796640.0 COORD (1, 2) 1 2