From 600cbc7b3bd221e9a0e51fd1634b063df52144e2 Mon Sep 17 00:00:00 2001 From: Marius Appel Date: Wed, 21 Feb 2024 22:10:09 +0100 Subject: [PATCH 1/4] add window_space [experimental] --- .Rbuildignore | 6 +- NAMESPACE | 2 +- R/RcppExports.R | 8 + R/window.R | 143 +++++++++--- man/window_space.Rd | 71 ++++++ man/window_time.Rd | 49 ++-- man/window_time.cube.Rd | 65 ------ src/Makevars.in | 1 + src/Makevars.ucrt | 1 + src/Makevars.win | 1 + src/RcppExports.cpp | 33 +++ src/gdalcubes.cpp | 36 +++ src/gdalcubes/src/cube.cpp | 95 ++++++++ src/gdalcubes/src/cube.h | 17 ++ src/gdalcubes/src/cube_factory.cpp | 22 ++ src/gdalcubes/src/gdalcubes.h | 1 + src/gdalcubes/src/window_space.cpp | 354 +++++++++++++++++++++++++++++ src/gdalcubes/src/window_space.h | 167 ++++++++++++++ 18 files changed, 956 insertions(+), 116 deletions(-) create mode 100644 man/window_space.Rd delete mode 100644 man/window_time.cube.Rd create mode 100644 src/gdalcubes/src/window_space.cpp create mode 100644 src/gdalcubes/src/window_space.h diff --git a/.Rbuildignore b/.Rbuildignore index 9830203..e44b4cb 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -39,6 +39,7 @@ ^src/gdalcubes/Dockerfile$ ^src/gdalcubes/docs\.cmake$ ^src/gdalcubes/formats$ +^src/gdalcubes/build$ ^src/gdalcubes/LICENSE$ ^src/gdalcubes/LICENSE_THIRDPARTY$ ^src/gdalcubes/README\.md$ @@ -75,6 +76,7 @@ ^\.pkgdown.R$ ^\.test.*$ ^\.*.db$ +^.*\.aux\.xml$ ^proj_conf_test$ ^proj_conf_test.cpp$ ^proj_conf_test.c$ @@ -85,4 +87,6 @@ ^gdalcubes.log$ ^.bash_history$ ^.vscode$ -^README.html$ \ No newline at end of file +^README.html$ +^inst/gdal$ +^inst/proj$ \ No newline at end of file diff --git a/NAMESPACE b/NAMESPACE index 86a7466..be1177c 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -18,7 +18,6 @@ S3method(reduce_space,array) S3method(reduce_space,cube) S3method(reduce_time,array) S3method(reduce_time,cube) -S3method(window_time,cube) export(add_collection_format) export(add_images) export(aggregate_space) @@ -73,6 +72,7 @@ export(srs) export(st_as_stars.cube) export(stac_image_collection) export(stack_cube) +export(window_space) export(window_time) export(write_chunk_from_array) export(write_ncdf) diff --git a/R/RcppExports.R b/R/RcppExports.R index db55471..8fb38b2 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -149,6 +149,14 @@ gc_create_window_time_cube_kernel <- function(pin, window, kernel) { .Call('_gdalcubes_gc_create_window_time_cube_kernel', PACKAGE = 'gdalcubes', pin, window, kernel) } +gc_create_window_space_cube_reduce <- function(pin, reducers, bands, win_size_y, win_size_x, keep_bands) { + .Call('_gdalcubes_gc_create_window_space_cube_reduce', PACKAGE = 'gdalcubes', pin, reducers, bands, win_size_y, win_size_x, keep_bands) +} + +gc_create_window_space_cube_kernel <- function(pin, kernel, win_size_y, win_size_x, keep_bands) { + .Call('_gdalcubes_gc_create_window_space_cube_kernel', PACKAGE = 'gdalcubes', pin, kernel, win_size_y, win_size_x, keep_bands) +} + gc_create_join_bands_cube <- function(pin_list, cube_names) { .Call('_gdalcubes_gc_create_join_bands_cube', PACKAGE = 'gdalcubes', pin_list, cube_names) } diff --git a/R/window.R b/R/window.R index cc9c1f5..7fca9f6 100644 --- a/R/window.R +++ b/R/window.R @@ -1,37 +1,3 @@ -#' Apply a moving window operation over time -#' -#' This generic function applies a reducer function over a moving window over the time dimension of a data cube, an R array, or other classes if implemented. -#' @param x object to be reduced -#' @param ... further arguments passed to specific implementations -#' @return value and type depend on the class of x -#' @seealso \code{\link{window_time.cube}} -#' @examples -#' # create image collection from example Landsat data only -#' # if not already done in other examples -#' if (!file.exists(file.path(tempdir(), "L8.db"))) { -#' L8_files <- list.files(system.file("L8NY18", package = "gdalcubes"), -#' ".TIF", recursive = TRUE, full.names = TRUE) -#' create_image_collection(L8_files, "L8_L1TP", file.path(tempdir(), "L8.db"), quiet = TRUE) -#' } -#' -#' L8.col = image_collection(file.path(tempdir(), "L8.db")) -#' v = cube_view(extent=list(left=388941.2, right=766552.4, -#' bottom=4345299, top=4744931, t0="2018-01", t1="2018-07"), -#' srs="EPSG:32618", nx = 400, dt="P1M") -#' L8.cube = raster_cube(L8.col, v) -#' L8.nir = select_bands(L8.cube, c("B05")) -#' window_time(L8.nir, window = c(2,2), "min(B05)") -#' window_time(L8.nir, kernel=c(-1,1), window=c(1,0)) -#' -#' \donttest{ -#' plot(window_time(L8.nir, kernel=c(-1,1), window=c(1,0)), key.pos=1) -#' } -#' @export -window_time <- function(x, ...) { - UseMethod("window_time") -} - - #' Apply a moving window function over the time dimension of a data cube #' #' Create a proxy data cube, which applies one ore more moving window functions to selected bands over pixel time series of a data cube. @@ -80,7 +46,7 @@ window_time <- function(x, ...) { #' #' #' @export -window_time.cube <- function(x, expr, ..., kernel, window) { +window_time <- function(x, expr, ..., kernel, window) { stopifnot(is.cube(x)) @@ -149,4 +115,111 @@ is.window_time_cube <- function(obj) { +#' Apply a moving window operation or convolution kernel over the spatial dimensions of a data cube +#' +#' Create a proxy data cube, which applies a convolution kernel or an aggregation functions on two-dimensional moving +#' windows sliding over spatial slices of a data cube. The function can either execute a predefined agggregation function or +#' apply a custom convolution kernel. Among others, use cases include image processing (edge detection, median filter noise reduction, etc.) and +#' enriching pixel values with local neighborhood properties (e.g. to use as predictor variables in ML models). +#' +#' @param x source data cube +#' @param kernel two dimensional kernel (matrix) applied as convolution (must have odd number of rows and columns) +#' @param expr either a single string, or a vector of strings, defining which reducers will be applied over which bands of the input cube +#' @param window integer vector with two elements defining the size of the window before and after a cell, the total size of the window is window[1] + 1 + window[2] +#' @param keep_bands logical; if FALSE (the default), original data cube bands will be dropped. +#' @param ... optional additional expressions (if expr is not a vector) +#' @return proxy data cube object +#' @note Implemented reducers will ignore any NAN values (as \code{na.rm = TRUE} does). +#' @note THIS FUNCTION IS STILL EXPERIMENTAL: The current implementation does not yet support any special edge handling. +#' @note THIS FUNCTION IS STILL EXPERIMENTAL: The current implementation will simply ignore NA values during the convolution, i.e. the sum will not be updated for NA pixels. +#' @examples +#' # create image collection from example Landsat data only +#' # if not already done in other examples +#' if (!file.exists(file.path(tempdir(), "L8.db"))) { +#' L8_files <- list.files(system.file("L8NY18", package = "gdalcubes"), +#' ".TIF", recursive = TRUE, full.names = TRUE) +#' create_image_collection(L8_files, "L8_L1TP", file.path(tempdir(), "L8.db"), quiet = TRUE) +#' } +#' +#' L8.col = image_collection(file.path(tempdir(), "L8.db")) +#' v = cube_view(extent=list(left=388941.2, right=766552.4, +#' bottom=4345299, top=4744931, t0="2018-04", t1="2018-06"), +#' srs="EPSG:32618", nx = 497, ny=526, dt="P1M") +#' L8.cube = raster_cube(L8.col, v, chunking = c(1,1000,1000)) +#' L8.cube = select_bands(L8.cube, c("B04", "B05")) +#' L8.cube.mean5x5 = window_space(L8.cube, kernel = matrix(1/25, 5, 5)) +#' L8.cube.mean5x5 +#' +#' L8.cube.med_sd = window_space(L8.cube, "median(B04)" ,"sd(B04)", "median(B05)", "sd(B05)", window = c(5,5), keep_bands = TRUE) +#' L8.cube.med_sd +#' +#' @note This function returns a proxy object, i.e., it will not start any computations besides deriving the shape of the result. +#' @details +#' The function either applies a kernel convolution (if the \code{kernel} argument is provided) or one or more built-in reducer function +#' over moving windows. In the former case, the kernel convolution will be applied over all bands of the input +#' cube, i.e., the output cube will have the same number of bands as the input cubes. +#' To apply one or more reducer functions, the window argument must be provided as a vector with two integer sizes in the order y, x. +#' Several string expressions can be provided to create multiple bands in the output cube. +#' +#' Notice that expressions have a very simple format: the reducer is followed by the name of a band in parentheses. You cannot add +#' more complex functions or arguments. Possible reducers currently include "min", "max", "sum", "prod", "count", "mean", "median", "var", and "sd". +#' +#' @export +window_space <- function(x, expr, ..., kernel, window, keep_bands = FALSE) { + stopifnot(is.cube(x)) + + if (!missing(kernel)) { + if (!missing(expr)) { + warning("argument expr will be ignored, applying kernel convolution") + } + if (length(list(...))> 0) { + warning("additional arguments will be ignored, applying kernel convolution") + } + if (!is.matrix(kernel)) { + stop("Kernel must be provided as a matrix") + } + x = gc_create_window_space_cube_kernel(x, as.double(kernel), as.integer(nrow(kernel)), as.integer(ncol(kernel)), keep_bands) + } + else { + stopifnot(is.character(expr)) + stopifnot(length(window) == 2) + stopifnot(window[1] %% 1 == 0) + stopifnot(window[2] %% 1 == 0) + + stopifnot(window[1] %% 2 == 1) + stopifnot(window[2] %% 2 == 1) + + if (length(list(...))> 0) { + stopifnot(all(sapply(list(...), is.character))) + expr = c(expr, unlist(list(...))) + } + + # parse expr to separate reducers and bands + reducers = gsub("\\(.*\\)", "", expr) + bands = gsub("[\\(\\)]", "", regmatches(expr, gregexpr("\\(.*?\\)", expr))) + stopifnot(length(reducers) == length(bands)) + x = gc_create_window_space_cube_reduce(x, reducers, bands, as.integer(window[1]), as.integer(window[2]), keep_bands) + } + class(x) <- c("window_space_cube", "cube", "xptr") + return(x) + +} + + + + +is.window_space_cube <- function(obj) { + if(!("window_space_cube" %in% class(obj))) { + return(FALSE) + } + if (gc_is_null(obj)) { + warning("GDAL data cube proxy object is invalid") + return(FALSE) + } + return(TRUE) +} + + + + diff --git a/man/window_space.Rd b/man/window_space.Rd new file mode 100644 index 0000000..f077d24 --- /dev/null +++ b/man/window_space.Rd @@ -0,0 +1,71 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/window.R +\name{window_space} +\alias{window_space} +\title{Apply a moving window operation or convolution kernel over the spatial dimensions of a data cube} +\usage{ +window_space(x, expr, ..., kernel, window, keep_bands = FALSE) +} +\arguments{ +\item{x}{source data cube} + +\item{expr}{either a single string, or a vector of strings, defining which reducers will be applied over which bands of the input cube} + +\item{...}{optional additional expressions (if expr is not a vector)} + +\item{kernel}{two dimensional kernel (matrix) applied as convolution (must have odd number of rows and columns)} + +\item{window}{integer vector with two elements defining the size of the window before and after a cell, the total size of the window is window[1] + 1 + window[2]} + +\item{keep_bands}{logical; if FALSE (the default), original data cube bands will be dropped.} +} +\value{ +proxy data cube object +} +\description{ +Create a proxy data cube, which applies a convolution kernel or an aggregation functions on two-dimensional moving +windows sliding over spatial slices of a data cube. The function can either execute a predefined agggregation function or +apply a custom convolution kernel. Among others, use cases include image processing (edge detection, median filter noise reduction, etc.) and +enriching pixel values with local neighborhood properties (e.g. to use as predictor variables in ML models). +} +\details{ +The function either applies a kernel convolution (if the \code{kernel} argument is provided) or one or more built-in reducer function +over moving windows. In the former case, the kernel convolution will be applied over all bands of the input +cube, i.e., the output cube will have the same number of bands as the input cubes. +To apply one or more reducer functions, the window argument must be provided as a vector with two integer sizes in the order y, x. +Several string expressions can be provided to create multiple bands in the output cube. + +Notice that expressions have a very simple format: the reducer is followed by the name of a band in parentheses. You cannot add +more complex functions or arguments. Possible reducers currently include "min", "max", "sum", "prod", "count", "mean", "median", "var", and "sd". +} +\note{ +Implemented reducers will ignore any NAN values (as \code{na.rm = TRUE} does). + +THIS FUNCTION IS STILL EXPERIMENTAL: The current implementation does not yet support any special edge handling. + +THIS FUNCTION IS STILL EXPERIMENTAL: The current implementation will simply ignore NA values during the convolution, i.e. the sum will not be updated for NA pixels. + +This function returns a proxy object, i.e., it will not start any computations besides deriving the shape of the result. +} +\examples{ +# create image collection from example Landsat data only +# if not already done in other examples +if (!file.exists(file.path(tempdir(), "L8.db"))) { + L8_files <- list.files(system.file("L8NY18", package = "gdalcubes"), + ".TIF", recursive = TRUE, full.names = TRUE) + create_image_collection(L8_files, "L8_L1TP", file.path(tempdir(), "L8.db"), quiet = TRUE) +} + +L8.col = image_collection(file.path(tempdir(), "L8.db")) +v = cube_view(extent=list(left=388941.2, right=766552.4, + bottom=4345299, top=4744931, t0="2018-04", t1="2018-06"), + srs="EPSG:32618", nx = 497, ny=526, dt="P1M") +L8.cube = raster_cube(L8.col, v, chunking = c(1,1000,1000)) +L8.cube = select_bands(L8.cube, c("B04", "B05")) +L8.cube.mean5x5 = window_space(L8.cube, kernel = matrix(1/25, 5, 5)) +L8.cube.mean5x5 + +L8.cube.med_sd = window_space(L8.cube, "median(B04)" ,"sd(B04)", "median(B05)", "sd(B05)", window = c(5,5), keep_bands = TRUE) +L8.cube.med_sd + +} diff --git a/man/window_time.Rd b/man/window_time.Rd index d513f23..03ce65c 100644 --- a/man/window_time.Rd +++ b/man/window_time.Rd @@ -2,20 +2,44 @@ % Please edit documentation in R/window.R \name{window_time} \alias{window_time} -\title{Apply a moving window operation over time} +\title{Apply a moving window function over the time dimension of a data cube} \usage{ -window_time(x, ...) +window_time(x, expr, ..., kernel, window) } \arguments{ -\item{x}{object to be reduced} +\item{x}{source data cube} -\item{...}{further arguments passed to specific implementations} +\item{expr}{either a single string, or a vector of strings defining which reducers wlil be applied over which bands of the input cube} + +\item{...}{optional additional expressions (if expr is not a vector)} + +\item{kernel}{numeric vector with elements of the kernel} + +\item{window}{integer vector with two elements defining the size of the window before and after a cell, the total size of the window is window[1] + 1 + window[2]} } \value{ -value and type depend on the class of x +proxy data cube object } \description{ -This generic function applies a reducer function over a moving window over the time dimension of a data cube, an R array, or other classes if implemented. +Create a proxy data cube, which applies one ore more moving window functions to selected bands over pixel time series of a data cube. +The fuction can either use a predefined agggregation function or apply a custom convolution kernel. +} +\details{ +The function either applies a kernel convolution (if the \code{kernel} argument is provided) or a general reducer function +over moving temporal windows. In the former case, the kernel convolution will be applied over all bands of the input +cube, i.e., the output cube will have the same number of bands as the input cubes. If a kernel is given and the \code{window} argument is missing, +the window will be symmetric to the center pixel with the size of the provided kernel. For general reducer functions, the window argument must be provided and +several expressions can be used to create multiple bands in the output cube. + +Notice that expressions have a very simple format: the reducer is followed by the name of a band in parantheses. You cannot add +more complex functions or arguments. + +Possible reducers currently are "min", "max", "sum", "prod", "count", "mean", "median". +} +\note{ +Implemented reducers will ignore any NAN values (as na.rm=TRUE does). + +This function returns a proxy object, i.e., it will not start any computations besides deriving the shape of the result. } \examples{ # create image collection from example Landsat data only @@ -32,13 +56,10 @@ v = cube_view(extent=list(left=388941.2, right=766552.4, srs="EPSG:32618", nx = 400, dt="P1M") L8.cube = raster_cube(L8.col, v) L8.nir = select_bands(L8.cube, c("B05")) -window_time(L8.nir, window = c(2,2), "min(B05)") -window_time(L8.nir, kernel=c(-1,1), window=c(1,0)) +L8.nir.min = window_time(L8.nir, window = c(2,2), "min(B05)") +L8.nir.min + +L8.nir.kernel = window_time(L8.nir, kernel=c(-1,1), window=c(1,0)) +L8.nir.kernel -\donttest{ -plot(window_time(L8.nir, kernel=c(-1,1), window=c(1,0)), key.pos=1) -} -} -\seealso{ -\code{\link{window_time.cube}} } diff --git a/man/window_time.cube.Rd b/man/window_time.cube.Rd deleted file mode 100644 index f52872c..0000000 --- a/man/window_time.cube.Rd +++ /dev/null @@ -1,65 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/window.R -\name{window_time.cube} -\alias{window_time.cube} -\title{Apply a moving window function over the time dimension of a data cube} -\usage{ -\method{window_time}{cube}(x, expr, ..., kernel, window) -} -\arguments{ -\item{x}{source data cube} - -\item{expr}{either a single string, or a vector of strings defining which reducers wlil be applied over which bands of the input cube} - -\item{...}{optional additional expressions (if expr is not a vector)} - -\item{kernel}{numeric vector with elements of the kernel} - -\item{window}{integer vector with two elements defining the size of the window before and after a cell, the total size of the window is window[1] + 1 + window[2]} -} -\value{ -proxy data cube object -} -\description{ -Create a proxy data cube, which applies one ore more moving window functions to selected bands over pixel time series of a data cube. -The fuction can either use a predefined agggregation function or apply a custom convolution kernel. -} -\details{ -The function either applies a kernel convolution (if the \code{kernel} argument is provided) or a general reducer function -over moving temporal windows. In the former case, the kernel convolution will be applied over all bands of the input -cube, i.e., the output cube will have the same number of bands as the input cubes. If a kernel is given and the \code{window} argument is missing, -the window will be symmetric to the center pixel with the size of the provided kernel. For general reducer functions, the window argument must be provided and -several expressions can be used to create multiple bands in the output cube. - -Notice that expressions have a very simple format: the reducer is followed by the name of a band in parantheses. You cannot add -more complex functions or arguments. - -Possible reducers currently are "min", "max", "sum", "prod", "count", "mean", "median". -} -\note{ -Implemented reducers will ignore any NAN values (as na.rm=TRUE does). - -This function returns a proxy object, i.e., it will not start any computations besides deriving the shape of the result. -} -\examples{ -# create image collection from example Landsat data only -# if not already done in other examples -if (!file.exists(file.path(tempdir(), "L8.db"))) { - L8_files <- list.files(system.file("L8NY18", package = "gdalcubes"), - ".TIF", recursive = TRUE, full.names = TRUE) - create_image_collection(L8_files, "L8_L1TP", file.path(tempdir(), "L8.db"), quiet = TRUE) -} - -L8.col = image_collection(file.path(tempdir(), "L8.db")) -v = cube_view(extent=list(left=388941.2, right=766552.4, - bottom=4345299, top=4744931, t0="2018-01", t1="2018-07"), - srs="EPSG:32618", nx = 400, dt="P1M") -L8.cube = raster_cube(L8.col, v) -L8.nir = select_bands(L8.cube, c("B05")) -L8.nir.min = window_time(L8.nir, window = c(2,2), "min(B05)") -L8.nir.min - -L8.nir.kernel = window_time(L8.nir, kernel=c(-1,1), window=c(1,0)) -L8.nir.kernel - -} diff --git a/src/Makevars.in b/src/Makevars.in index 570069e..befe49b 100644 --- a/src/Makevars.in +++ b/src/Makevars.in @@ -20,6 +20,7 @@ OBJECTS = gdalcubes/src/aggregate_time.o \ gdalcubes/src/select_time.o \ gdalcubes/src/reduce_time.o \ gdalcubes/src/reduce_space.o \ + gdalcubes/src/window_space.o \ gdalcubes/src/window_time.o \ gdalcubes/src/select_bands.o \ gdalcubes/src/slice_time.o \ diff --git a/src/Makevars.ucrt b/src/Makevars.ucrt index 727086a..94ec5b8 100644 --- a/src/Makevars.ucrt +++ b/src/Makevars.ucrt @@ -19,6 +19,7 @@ OBJECTS = gdalcubes/src/aggregate_time.o \ gdalcubes/src/select_time.o \ gdalcubes/src/reduce_time.o \ gdalcubes/src/reduce_space.o \ + gdalcubes/src/window_space.o \ gdalcubes/src/window_time.o \ gdalcubes/src/select_bands.o \ gdalcubes/src/slice_time.o \ diff --git a/src/Makevars.win b/src/Makevars.win index ef78589..cab8950 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -22,6 +22,7 @@ LIBGDALCUBES = gdalcubes/src/aggregate_time.o \ gdalcubes/src/select_time.o \ gdalcubes/src/reduce_time.o \ gdalcubes/src/reduce_space.o \ + gdalcubes/src/window_space.o \ gdalcubes/src/window_time.o \ gdalcubes/src/select_bands.o \ gdalcubes/src/slice_time.o \ diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index dcebf67..327e391 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -447,6 +447,37 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// gc_create_window_space_cube_reduce +SEXP gc_create_window_space_cube_reduce(SEXP pin, std::vector reducers, std::vector bands, int win_size_y, int win_size_x, bool keep_bands); +RcppExport SEXP _gdalcubes_gc_create_window_space_cube_reduce(SEXP pinSEXP, SEXP reducersSEXP, SEXP bandsSEXP, SEXP win_size_ySEXP, SEXP win_size_xSEXP, SEXP keep_bandsSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< SEXP >::type pin(pinSEXP); + Rcpp::traits::input_parameter< std::vector >::type reducers(reducersSEXP); + Rcpp::traits::input_parameter< std::vector >::type bands(bandsSEXP); + Rcpp::traits::input_parameter< int >::type win_size_y(win_size_ySEXP); + Rcpp::traits::input_parameter< int >::type win_size_x(win_size_xSEXP); + Rcpp::traits::input_parameter< bool >::type keep_bands(keep_bandsSEXP); + rcpp_result_gen = Rcpp::wrap(gc_create_window_space_cube_reduce(pin, reducers, bands, win_size_y, win_size_x, keep_bands)); + return rcpp_result_gen; +END_RCPP +} +// gc_create_window_space_cube_kernel +SEXP gc_create_window_space_cube_kernel(SEXP pin, std::vector kernel, int win_size_y, int win_size_x, bool keep_bands); +RcppExport SEXP _gdalcubes_gc_create_window_space_cube_kernel(SEXP pinSEXP, SEXP kernelSEXP, SEXP win_size_ySEXP, SEXP win_size_xSEXP, SEXP keep_bandsSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< SEXP >::type pin(pinSEXP); + Rcpp::traits::input_parameter< std::vector >::type kernel(kernelSEXP); + Rcpp::traits::input_parameter< int >::type win_size_y(win_size_ySEXP); + Rcpp::traits::input_parameter< int >::type win_size_x(win_size_xSEXP); + Rcpp::traits::input_parameter< bool >::type keep_bands(keep_bandsSEXP); + rcpp_result_gen = Rcpp::wrap(gc_create_window_space_cube_kernel(pin, kernel, win_size_y, win_size_x, keep_bands)); + return rcpp_result_gen; +END_RCPP +} // gc_create_join_bands_cube SEXP gc_create_join_bands_cube(Rcpp::List pin_list, std::vector cube_names); RcppExport SEXP _gdalcubes_gc_create_join_bands_cube(SEXP pin_listSEXP, SEXP cube_namesSEXP) { @@ -854,6 +885,8 @@ static const R_CallMethodDef CallEntries[] = { {"_gdalcubes_gc_create_reduce_space_cube", (DL_FUNC) &_gdalcubes_gc_create_reduce_space_cube, 4}, {"_gdalcubes_gc_create_window_time_cube_reduce", (DL_FUNC) &_gdalcubes_gc_create_window_time_cube_reduce, 4}, {"_gdalcubes_gc_create_window_time_cube_kernel", (DL_FUNC) &_gdalcubes_gc_create_window_time_cube_kernel, 3}, + {"_gdalcubes_gc_create_window_space_cube_reduce", (DL_FUNC) &_gdalcubes_gc_create_window_space_cube_reduce, 6}, + {"_gdalcubes_gc_create_window_space_cube_kernel", (DL_FUNC) &_gdalcubes_gc_create_window_space_cube_kernel, 5}, {"_gdalcubes_gc_create_join_bands_cube", (DL_FUNC) &_gdalcubes_gc_create_join_bands_cube, 2}, {"_gdalcubes_gc_create_select_bands_cube", (DL_FUNC) &_gdalcubes_gc_create_select_bands_cube, 2}, {"_gdalcubes_gc_create_select_time_cube", (DL_FUNC) &_gdalcubes_gc_create_select_time_cube, 2}, diff --git a/src/gdalcubes.cpp b/src/gdalcubes.cpp index 98bc924..9eceed6 100644 --- a/src/gdalcubes.cpp +++ b/src/gdalcubes.cpp @@ -1143,6 +1143,42 @@ SEXP gc_create_window_time_cube_kernel(SEXP pin, std::vector window, std::v } } +// [[Rcpp::export]] +SEXP gc_create_window_space_cube_reduce(SEXP pin, std::vector reducers, std::vector bands, int win_size_y, int win_size_x, bool keep_bands) { + try { + Rcpp::XPtr< std::shared_ptr > aa = Rcpp::as>>(pin); + + std::vector> reducer_bands; + for (uint16_t i=0; i* x = new std::shared_ptr(window_space_cube::create(*aa, reducer_bands, win_size_y, win_size_x, keep_bands)); + Rcpp::XPtr< std::shared_ptr > p(x, true) ; + + return p; + + } + catch (std::string s) { + Rcpp::stop(s); + } +} + +// [[Rcpp::export]] +SEXP gc_create_window_space_cube_kernel(SEXP pin, std::vector kernel, int win_size_y, int win_size_x, bool keep_bands) { + try { + Rcpp::XPtr< std::shared_ptr > aa = Rcpp::as>>(pin); + + std::shared_ptr* x = new std::shared_ptr(window_space_cube::create(*aa, kernel, win_size_y, win_size_x, keep_bands)); + Rcpp::XPtr< std::shared_ptr > p(x, true) ; + return p; + + } + catch (std::string s) { + Rcpp::stop(s); + } +} // [[Rcpp::export]] SEXP gc_create_join_bands_cube(Rcpp::List pin_list, std::vector cube_names) { diff --git a/src/gdalcubes/src/cube.cpp b/src/gdalcubes/src/cube.cpp index a3d9bca..fb8fd18 100644 --- a/src/gdalcubes/src/cube.cpp +++ b/src/gdalcubes/src/cube.cpp @@ -1799,6 +1799,101 @@ std::shared_ptr cube::to_double_array(std::shared_ptr cube::read_window(std::array lower, std::array upper) { + // This function allows negative coordinates + // This funtion guarantees that the output buffer has size prod(upper - lower + 1) + + if (upper[0] < lower[0] || upper[1] < lower[1] || upper[2] < lower[2]) { + GCBS_ERROR("ERROR in cube::read_window(): invalid window provided."); + throw std::string("ERROR in cube::read_window(): invalid window provided."); + } + + std::array chnk_min = {lower[0] / int32_t(_chunk_size[0]), lower[1] / int32_t(_chunk_size[1]), lower[2] / int32_t(_chunk_size[2])}; + std::array chnk_max = {upper[0] / int32_t(_chunk_size[0]), upper[1] / int32_t(_chunk_size[1]), upper[2] / int32_t(_chunk_size[2])}; + + std::shared_ptr out = std::make_shared(); + coords_nd size_tyx = {upper[0] - lower[0] + 1, upper[1] - lower[1] + 1, upper[2] - lower[2] + 1}; + + + coords_nd size_btyx = {uint32_t(_bands.count()), uint32_t(size_tyx[0]), uint32_t(size_tyx[1]), uint32_t(size_tyx[2])}; + out->size(size_btyx); + + // Fill buffers accordingly + out->buf(std::calloc(size_btyx[0] * size_btyx[1] * size_btyx[2] * size_btyx[3], sizeof(double))); + double* begin = (double*)out->buf(); + double* end = ((double*)out->buf()) + size_btyx[0] * size_btyx[1] * size_btyx[2] * size_btyx[3]; + std::fill(begin, end, NAN); + + // Areas outside the cube are always filled with NAN. Padding must be explicitly implemented in upstram code, if needed (e.g. window_space) + for (int32_t ict = std::max(chnk_min[0], 0); ict<= std::min(chnk_max[0], (int32_t)count_chunks_t()-1); ++ict) { + for (int32_t icy = std::max(chnk_min[1], 0); icy<= std::min(chnk_max[1], (int32_t)count_chunks_y()-1); ++icy) { + for (int32_t icx = std::max(chnk_min[2], 0); icx<= std::min(chnk_max[2], (int32_t)count_chunks_x()-1); ++icx) { + + // 1. Read chunk + chunkid_t chnk_id = chunk_id_from_coords({uint32_t(ict), uint32_t(icy), uint32_t(icx)}); // casting to unsigned needed? + std::shared_ptr cin = this->read_chunk(chnk_id); + if (cin->empty()) { + continue; + } + + + // 2. Iterate over chunk and copy values to correct position in the buffer + + // Offsets: Relative position between chunk pixels and target buffer pixels + int32_t offst_t = lower[0] - ict * _chunk_size[0]; + int32_t offst_y = lower[1] - icy * _chunk_size[1]; + int32_t offst_x = lower[2] - icx * _chunk_size[2]; + + int32_t start_t = std::max(offst_t, 0); + int32_t end_t = std::min(int32_t(cin->size()[1]) - 1, offst_t + (upper[0] - lower[0])); + int32_t start_y = std::max(offst_y, 0); + int32_t end_y = std::min(int32_t(cin->size()[2]) - 1, offst_y + (upper[1] - lower[1])); + int32_t start_x = std::max(offst_x, 0); + int32_t end_x = std::min(int32_t(cin->size()[3]) - 1, offst_x + (upper[2] - lower[2])); + + for (int32_t ib = 0; ib < int32_t(cin->size()[0]); ++ib) { + for (int32_t it = start_t; it<=end_t; ++it) { + for (int32_t iy = start_y; iy<=end_y; ++iy) { + for (int32_t ix = start_x; ix<=end_x; ++ix) { + ((double*)(out->buf()))[ + ib * size_tyx[0] * size_tyx[1] * size_tyx[2] + + (it - offst_t) * size_tyx[1] * size_tyx[2] + + (iy - offst_y) * size_tyx[2] + + (ix - offst_x)] = + ((double*)(cin->buf()))[ + ib * (int32_t)cin->size()[1] * (int32_t)cin->size()[2] * cin->size()[3] + + it * (int32_t)cin->size()[2] * (int32_t)cin->size()[3] + + iy * (int32_t)cin->size()[3] + + ix]; + } + } + } + } + } + } + } + + return out; + + +} + + + + + + + + + + + + + + void chunk_data::write_ncdf(std::string path, uint8_t compression_level, bool force) { if (filesystem::exists(path)) { GCBS_ERROR("File already exists"); diff --git a/src/gdalcubes/src/cube.h b/src/gdalcubes/src/cube.h index 86164df..a474cb5 100644 --- a/src/gdalcubes/src/cube.h +++ b/src/gdalcubes/src/cube.h @@ -760,6 +760,23 @@ class cube : public std::enable_shared_from_this { */ virtual std::shared_ptr read_chunk(chunkid_t id) = 0; + + /** + * @brief Read a window subset of a data cube to a buffer + * + * This funtion guarantees that the output buffer has size prod(upper - lower + 1). + * It can work with negative cube coordinates. + * + * Areas outside the cube are always filled with NAN. + * Padding must be explicitly implemented in upstram code, if needed (e.g. window_space) + * + * + * @param lower coordinates of the lower pixel (cube coordinates) + * @param upper coordinates of the upper point (cube coordinates) + * @return a smart pointer to chunk data + */ + std::shared_ptr read_window(std::array lower, std::array upper); + /** * @brief Write a data cube as a set of GeoTIFF files under a given directory * diff --git a/src/gdalcubes/src/cube_factory.cpp b/src/gdalcubes/src/cube_factory.cpp index 1734f75..c7ce051 100644 --- a/src/gdalcubes/src/cube_factory.cpp +++ b/src/gdalcubes/src/cube_factory.cpp @@ -53,6 +53,7 @@ #include "stream_apply_time.h" #include "stream_reduce_space.h" #include "stream_reduce_time.h" +#include "window_space.h" #include "window_time.h" namespace gdalcubes { @@ -143,6 +144,27 @@ void cube_factory::register_default() { j["win_size_l"].int_value(), j["win_size_r"].int_value()); } })); + cube_generators.insert(std::make_pair(json11::Json&)>>( + "window_space", [](json11::Json& j) { + if (!j["kernel"].is_null()) { + std::vector kernel; + for (uint16_t i = 0; i < j["kernel"].array_items().size(); ++i) { + kernel.push_back(j["kernel"][i].number_value()); + } + return window_space_cube::create(instance()->create_from_json(j["in_cube"]), kernel, + j["win_size_y"].int_value(), j["win_size_x"].int_value(), + j["keep_bands"].bool_value()); + } + else { + std::vector> band_reducers; + for (uint16_t i = 0; i < j["reducer_bands"].array_items().size(); ++i) { + band_reducers.push_back(std::make_pair(j["reducer_bands"][i][0].string_value(), j["reducer_bands"][i][1].string_value())); + } + return window_space_cube::create(instance()->create_from_json(j["in_cube"]), band_reducers, + j["win_size_y"].int_value(), j["win_size_x"].int_value(), + j["keep_bands"].bool_value()); + } + })); cube_generators.insert(std::make_pair(json11::Json&)>>( "select_bands", [](json11::Json& j) { std::vector bands; diff --git a/src/gdalcubes/src/gdalcubes.h b/src/gdalcubes/src/gdalcubes.h index 3670829..c221eea 100644 --- a/src/gdalcubes/src/gdalcubes.h +++ b/src/gdalcubes/src/gdalcubes.h @@ -59,6 +59,7 @@ #include "stream_reduce_time.h" #include "utils.h" #include "vector_queries.h" +#include "window_space.h" #include "window_time.h" #ifndef GDALCUBES_NO_SWARM diff --git a/src/gdalcubes/src/window_space.cpp b/src/gdalcubes/src/window_space.cpp new file mode 100644 index 0000000..ab32048 --- /dev/null +++ b/src/gdalcubes/src/window_space.cpp @@ -0,0 +1,354 @@ +/* + MIT License + + Copyright (c) 2024 Marius Appel + + Permission is hereby granted, free of charge, to any person obtaining a copy + of this software and associated documentation files (the "Software"), to deal + in the Software without restriction, including without limitation the rights + to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + copies of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be included in all + copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + SOFTWARE. +*/ +#include "window_space.h" + +namespace gdalcubes { + + + + + +struct window_reducer_singleband { + virtual ~window_reducer_singleband() {} + virtual void init() = 0; + virtual void update(double &v) = 0; + virtual double finalize() = 0; +}; + +struct window_reducer_mean : public window_reducer_singleband { + void init() override { + sum = 0.0; + count = 0; + } + void update(double &v) override { + sum += v; + count++; + } + double finalize() override { + return sum / double(count); + } + double sum; + uint32_t count; +}; + + +struct window_reducer_sum : public window_reducer_singleband { + void init() override { + sum = 0.0; + } + void update(double &v) override { + sum += v; + } + double finalize() override { + return sum; + } + double sum; +}; + +struct window_reducer_prod : public window_reducer_singleband { + void init() override { + prod = 0.0; + } + void update(double &v) override { + prod *= v; + } + double finalize() override { + return prod; + } + double prod; +}; + +struct window_reducer_max : public window_reducer_singleband { + void init() override { + max = NAN; + } + void update(double &v) override { + if (std::isnan(max) || v > max) { + max = v; + } + } + double finalize() override { + return max; + } + double max; +}; + +struct window_reducer_min : public window_reducer_singleband { + void init() override { + min = NAN; + } + void update(double &v) override { + if (std::isnan(min) || v < min) { + min = v; + } + } + double finalize() override { + return min; + } + double min; +}; + + +struct window_reducer_count : public window_reducer_singleband { + void init() override { + count = 0; + } + void update(double &v) override { + if (std::isfinite(v)) { + count++; + } + } + double finalize() override { + return count; + } + uint32_t count; +}; + +struct window_reducer_median : public window_reducer_singleband { + void init() override { + values.clear(); + } + void update(double &v) override { + if (std::isfinite(v)) { + values.push_back(v); + } + } + double finalize() override { + if (values.size() == 0) { + return NAN; + } + std::sort(values.begin(), values.end()); + if (values.size() % 2 == 1) { + return values[values.size() / 2]; + } else { + return (values[values.size() / 2] + values[values.size() / 2 - 1]) / ((double)2); + } + } + std::vector values; +}; + + + +struct window_reducer_var : public window_reducer_singleband { + void init() override { + mean = 0.0; + count = 0; + state = 0.0; + } + void update(double &v) override { + if (std::isfinite(v)) { + ++count; + double delta = v - mean; + mean += delta / count; + state += delta * (v - mean); + } + } + double finalize() override { + return count > 1 ? state / (count - 1) : NAN; + } + + double mean; + uint32_t count; + double state; +}; + + +struct window_reducer_sd : public window_reducer_var { + double finalize() override { + return count > 1 ? sqrt(state / (count - 1)) : NAN; + } +}; + + + + + +std::shared_ptr window_space_cube::read_chunk(chunkid_t id) { + GCBS_TRACE("window_space_cube::read_chunk(" + std::to_string(id) + ")"); + std::shared_ptr out = std::make_shared(); + if (id >= count_chunks()) + return out; // chunk is outside of the view, we don't need to read anything. + + coords_nd size_tyx = chunk_size(id); + + coords_nd size_btyx = {uint32_t(_bands.count()), size_tyx[0], size_tyx[1], size_tyx[2]}; + out->size(size_btyx); + + // Fill buffers accordingly + out->buf(std::calloc(size_btyx[0] * size_btyx[1] * size_btyx[2] * size_btyx[3], sizeof(double))); + double* begin = (double*)out->buf(); + double* end = ((double*)out->buf()) + size_btyx[0] * size_btyx[1] * size_btyx[2] * size_btyx[3]; + std::fill(begin, end, NAN); + + + // 1. Extract "subwindow" from input chunks --> function + // size: (winsize[0] + chunk_size[0] + winsize[2]) x (winsize[1] + chunk_size[1] + winsize[3]) + auto ccoords = chunk_coords_from_id(id); + std::array lower = {int32_t(ccoords[0] * _chunk_size[0]), + int32_t(ccoords[1] * _chunk_size[1] - ((_win_size_y - 1) / 2)), + int32_t(ccoords[2] * _chunk_size[2] - ((_win_size_x - 1) / 2))}; + std::array upper = {int32_t(ccoords[0] * _chunk_size[0] + size_tyx[0] - 1), + int32_t(ccoords[1] * _chunk_size[1] + size_tyx[1] - 1 + ((_win_size_y - 1) / 2)), + int32_t(ccoords[2] * _chunk_size[2] + size_tyx[2] - 1 + ((_win_size_x - 1) / 2))}; + + auto cwin = _in_cube->read_window(lower, upper); + + + + + if (!_kernel.empty()) { + // TODO: Apply padding + + // 2. Iterate over target buffer and apply kernel. + for (uint32_t ib = 0; ib < size_btyx[0]; ++ib) { + for (uint32_t it = 0; it < size_btyx[1]; ++it) { + for (uint32_t iy = 0; iy < size_btyx[2]; ++iy) { + for (uint32_t ix = 0; ix < size_btyx[3]; ++ix) { + + // apply kernel + double sum = NAN; // TODO: complete.cases only? + for (int16_t ky=0; ky < _win_size_y; ++ky) { + for (int16_t kx=0; kx < _win_size_x; ++kx) { + double v = ((double*)(cwin->buf()))[ + ib * cwin->size()[1] * cwin->size()[2] * cwin->size()[3] + + it * cwin->size()[2] * cwin->size()[3] + + (iy+ky) * cwin->size()[3] + + (ix + kx) + ]; + if (std::isnan(sum)) { + sum = _kernel[ky * _win_size_x + kx] * v; + } + else if (std::isfinite(v)) + sum += _kernel[ky * _win_size_x + kx] * v; + + } + } + + ((double*)(out->buf()))[ + ib * size_tyx[0] * size_tyx[1] * size_tyx[2] + + it * size_tyx[1] * size_tyx[2] + + iy * size_tyx[2] + + ix] = sum; + + } + } + } + } + } + else { + std::vector reducers; + for (uint16_t i = 0; i < _reducer_bands.size(); ++i) { + window_reducer_singleband *r = nullptr; + if (_reducer_bands[i].first == "min") { + r = new window_reducer_min(); + } else if (_reducer_bands[i].first == "max") { + r = new window_reducer_max(); + } else if (_reducer_bands[i].first == "mean") { + r = new window_reducer_mean(); + } else if (_reducer_bands[i].first == "median") { + r = new window_reducer_median(); + } else if (_reducer_bands[i].first == "sum") { + r = new window_reducer_sum(); + } else if (_reducer_bands[i].first == "count") { + r = new window_reducer_count(); + } else if (_reducer_bands[i].first == "prod") { + r = new window_reducer_prod(); + } else if (_reducer_bands[i].first == "var") { + r = new window_reducer_var(); + } else if (_reducer_bands[i].first == "sd") { + r = new window_reducer_sd(); + } else + throw std::string("ERROR in window_space_cube::read_chunk(): Unknown reducer given"); + + reducers.push_back(r); + } + + uint32_t ib=0; + if (_keep_bands) { + int32_t offst_y = (_win_size_y-1) / 2; + int32_t offst_x = (_win_size_x-1) / 2; + for (; ib < _in_cube->size_bands(); ++ib) { + for (uint32_t it = 0; it < size_btyx[1]; ++it) { + for (uint32_t iy = 0; iy < size_btyx[2]; ++iy) { + for (uint32_t ix = 0; ix < size_btyx[3]; ++ix) { + + double v = ((double*)(cwin->buf()))[ + ib * cwin->size()[1] * cwin->size()[2] * cwin->size()[3] + + it * cwin->size()[2] * cwin->size()[3] + + (iy + offst_y) * cwin->size()[3] + + (ix + offst_x) + ]; + + ((double*)(out->buf()))[ + ib * size_tyx[0] * size_tyx[1] * size_tyx[2] + + it * size_tyx[1] * size_tyx[2] + + iy * size_tyx[2] + + ix] = v; + } + } + } + } + } + int16_t offst_b = _keep_bands ? _in_cube->size_bands() : 0; + for (; ib < size_btyx[0]; ++ib) { + for (uint32_t it = 0; it < size_btyx[1]; ++it) { + for (uint32_t iy = 0; iy < size_btyx[2]; ++iy) { + for (uint32_t ix = 0; ix < size_btyx[3]; ++ix) { + // apply reducer + reducers[ib-offst_b]->init(); + for (int16_t ky=0; ky < _win_size_y; ++ky) { + for (int16_t kx=0; kx < _win_size_x; ++kx) { + double v = ((double*)(cwin->buf()))[ + _band_idx_in[ib-offst_b] * cwin->size()[1] * cwin->size()[2] * cwin->size()[3] + + it * cwin->size()[2] * cwin->size()[3] + + (iy+ky) * cwin->size()[3] + + (ix + kx) + ]; + reducers[ib-offst_b]->update(v); + } + } + + ((double*)(out->buf()))[ + ib * size_tyx[0] * size_tyx[1] * size_tyx[2] + + it * size_tyx[1] * size_tyx[2] + + iy * size_tyx[2] + + ix] = reducers[ib-offst_b]->finalize(); + + } + } + } + } + for (uint16_t i = 0; i < reducers.size(); ++i) { + if (reducers[i]) delete reducers[i]; + } + } + + // check if chunk is completely NAN and if yes, return empty chunk + if (out->all_nan()) { + out = std::make_shared(); + } + return out; +} + +} // namespace gdalcubes diff --git a/src/gdalcubes/src/window_space.h b/src/gdalcubes/src/window_space.h new file mode 100644 index 0000000..39041d2 --- /dev/null +++ b/src/gdalcubes/src/window_space.h @@ -0,0 +1,167 @@ +/* + MIT License + + Copyright (c) 2024 Marius Appel + + Permission is hereby granted, free of charge, to any person obtaining a copy + of this software and associated documentation files (the "Software"), to deal + in the Software without restriction, including without limitation the rights + to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + copies of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be included in all + copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + SOFTWARE. +*/ + +#ifndef WINDOW_SPACE_H +#define WINDOW_SPACE_H + +#include "cube.h" + +namespace gdalcubes { + +// TODO: add custom reducer function +// TODO: add boundary fill options (e.g. NA, wrap, repeat, mirror, ...) +/** + * @brief A data cube that applies a focal operation / kernel over sliding spatial windows + */ +class window_space_cube : public cube { + public: + + /** + * @brief Create a data cube that applies a reducer function on a given input data cube over time + * @note This static creation method should preferably be used instead of the constructors as + * the constructors will not set connections between cubes properly. + * @param in input data cube + * @param reducer reducer function + * @return a shared pointer to the created data cube instance + */ + static std::shared_ptr + create(std::shared_ptr in, std::vector> reducer_bands, + uint16_t win_size_y, uint16_t win_size_x, bool keep_bands) { + std::shared_ptr out = std::make_shared(in, reducer_bands, win_size_y, win_size_x, keep_bands); + in->add_child_cube(out); + out->add_parent_cube(in); + return out; + } + + /** + * @brief Create a data cube that applies a convolution kernel on a given input data cube over time + * @note This static creation method should preferably be used instead of the constructors as + * the constructors will not set connections between cubes properly. + * @param in input data cube + * @param reducer reducer function + * @return a shared pointer to the created data cube instance + */ + static std::shared_ptr + create(std::shared_ptr in, std::vector kernel, uint16_t win_size_y, uint16_t win_size_x, bool keep_bands) { + std::shared_ptr out = std::make_shared(in, kernel, win_size_y, win_size_x, keep_bands); + in->add_child_cube(out); + out->add_parent_cube(in); + return out; + } + + public: + + window_space_cube(std::shared_ptr in, std::vector> reducer_bands, + uint16_t win_size_y, uint16_t win_size_x, bool keep_bands) : cube(in->st_reference()->copy()), _in_cube(in), _reducer_bands(reducer_bands), _win_size_y(win_size_y), _win_size_x(win_size_x), _band_idx_in(), _kernel(), _keep_bands(keep_bands) { // it is important to duplicate st reference here, otherwise changes will affect input cube as well + _chunk_size[0] = _in_cube->chunk_size()[0]; + _chunk_size[1] = _in_cube->chunk_size()[1]; + _chunk_size[2] = _in_cube->chunk_size()[2]; + + if (win_size_y % 2 == 0 || win_size_x % 2 == 0) { + GCBS_ERROR("Window size must not be even"); + throw std::string( + "ERROR in window_space_cube::window_space_cube(): Window size must not be even"); + } + + if (_keep_bands) { + // TODO: check for band name conflicts here + for (uint16_t i = 0; i < _in_cube->size_bands(); ++i) { + _bands.add(_in_cube->bands().get(i)); + } + } + for (uint16_t i = 0; i < reducer_bands.size(); ++i) { + std::string reducerstr = reducer_bands[i].first; + std::string bandstr = reducer_bands[i].second; + + band b = in->bands().get(bandstr); + b.name = b.name + "_" + reducerstr; + _bands.add(b); + + _band_idx_in.push_back(in->bands().get_index(bandstr)); + } + } + + + window_space_cube(std::shared_ptr in, std::vector kernel, uint16_t win_size_y, uint16_t win_size_x, bool keep_bands) + : cube(in->st_reference()->copy()), _in_cube(in), _reducer_bands(), _win_size_y(win_size_y), _win_size_x(win_size_x), _band_idx_in(), _kernel(kernel), _keep_bands(keep_bands) { // it is important to duplicate st reference here, otherwise changes will affect input cube as well + _chunk_size[0] = _in_cube->chunk_size()[0]; + _chunk_size[1] = _in_cube->chunk_size()[1]; + _chunk_size[2] = _in_cube->chunk_size()[2]; + + if (win_size_y % 2 == 0 || win_size_x % 2 == 0) { + GCBS_ERROR("Window size must not be even"); + throw std::string( + "ERROR in window_space_cube::window_space_cube(): Window size must not be even"); + } + + if (kernel.size() != (win_size_y * win_size_x)) { + GCBS_ERROR("kernel size does not match window size"); + throw std::string( + "ERROR in window_space_cube::window_space_cube(): Kernel size does not match window size"); + } + // Important: keep_bands is ignored if a kernel is used + for (uint16_t i = 0; i < in->bands().count(); ++i) { + band b = in->bands().get(i); + _bands.add(b); + } + } + + public: + ~window_space_cube() {} + + std::shared_ptr read_chunk(chunkid_t id) override; + + json11::Json make_constructible_json() override { + json11::Json::object out; + out["cube_type"] = "window_space"; + if (!_kernel.empty()) { + out["kernel"] = _kernel; + } + else { + json11::Json::array rb; + for (uint16_t i = 0; i < _reducer_bands.size(); ++i) { + rb.push_back(json11::Json::array({_reducer_bands[i].first, _reducer_bands[i].second})); + } + out["reducer_bands"] = rb; + } + out["win_size_y"] = _win_size_y; + out["win_size_x"] = _win_size_x; + out["keep_bands"] = _keep_bands; + out["in_cube"] = _in_cube->make_constructible_json(); + return out; + } + + private: + std::shared_ptr _in_cube; + std::vector> _reducer_bands; + uint16_t _win_size_y; + uint16_t _win_size_x; + std::vector _band_idx_in; + std::vector _kernel; + bool _keep_bands; +}; + +} // namespace gdalcubes + +#endif // WINDOW_SPACE_H From 918142040c32e4f5f02b1444cc05cb6c27d9895a Mon Sep 17 00:00:00 2001 From: Marius Appel Date: Thu, 22 Feb 2024 17:31:35 +0100 Subject: [PATCH 2/4] add padding methods to window_time() [experimental] --- R/RcppExports.R | 8 +- R/window.R | 40 ++- man/window_space.Rd | 8 +- src/RcppExports.cpp | 20 +- src/gdalcubes.cpp | 8 +- src/gdalcubes/src/cube_factory.cpp | 4 +- src/gdalcubes/src/window_space.cpp | 380 ++++++++++++++++++++++++++++- src/gdalcubes/src/window_space.h | 70 +++++- 8 files changed, 497 insertions(+), 41 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 8fb38b2..ef94095 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -149,12 +149,12 @@ gc_create_window_time_cube_kernel <- function(pin, window, kernel) { .Call('_gdalcubes_gc_create_window_time_cube_kernel', PACKAGE = 'gdalcubes', pin, window, kernel) } -gc_create_window_space_cube_reduce <- function(pin, reducers, bands, win_size_y, win_size_x, keep_bands) { - .Call('_gdalcubes_gc_create_window_space_cube_reduce', PACKAGE = 'gdalcubes', pin, reducers, bands, win_size_y, win_size_x, keep_bands) +gc_create_window_space_cube_reduce <- function(pin, reducers, bands, win_size_y, win_size_x, keep_bands, pad_mode, pad_fill) { + .Call('_gdalcubes_gc_create_window_space_cube_reduce', PACKAGE = 'gdalcubes', pin, reducers, bands, win_size_y, win_size_x, keep_bands, pad_mode, pad_fill) } -gc_create_window_space_cube_kernel <- function(pin, kernel, win_size_y, win_size_x, keep_bands) { - .Call('_gdalcubes_gc_create_window_space_cube_kernel', PACKAGE = 'gdalcubes', pin, kernel, win_size_y, win_size_x, keep_bands) +gc_create_window_space_cube_kernel <- function(pin, kernel, win_size_y, win_size_x, keep_bands, pad_mode, pad_fill) { + .Call('_gdalcubes_gc_create_window_space_cube_kernel', PACKAGE = 'gdalcubes', pin, kernel, win_size_y, win_size_x, keep_bands, pad_mode, pad_fill) } gc_create_join_bands_cube <- function(pin_list, cube_names) { diff --git a/R/window.R b/R/window.R index 7fca9f6..ea57535 100644 --- a/R/window.R +++ b/R/window.R @@ -127,11 +127,11 @@ is.window_time_cube <- function(obj) { #' @param expr either a single string, or a vector of strings, defining which reducers will be applied over which bands of the input cube #' @param window integer vector with two elements defining the size of the window before and after a cell, the total size of the window is window[1] + 1 + window[2] #' @param keep_bands logical; if FALSE (the default), original data cube bands will be dropped. +#' @param pad Padding method applied to the borders; use NULL for no padding, a numeric a fill value, or one of "REPLICATE", "REFLECT", "REFLECT_PIXEL" #' @param ... optional additional expressions (if expr is not a vector) #' @return proxy data cube object #' @note Implemented reducers will ignore any NAN values (as \code{na.rm = TRUE} does). -#' @note THIS FUNCTION IS STILL EXPERIMENTAL: The current implementation does not yet support any special edge handling. -#' @note THIS FUNCTION IS STILL EXPERIMENTAL: The current implementation will simply ignore NA values during the convolution, i.e. the sum will not be updated for NA pixels. +#' @note THIS FUNCTION IS STILL EXPERIMENTAL. #' @examples #' # create image collection from example Landsat data only #' # if not already done in other examples @@ -165,9 +165,39 @@ is.window_time_cube <- function(obj) { #' more complex functions or arguments. Possible reducers currently include "min", "max", "sum", "prod", "count", "mean", "median", "var", and "sd". #' #' @export -window_space <- function(x, expr, ..., kernel, window, keep_bands = FALSE) { +window_space <- function(x, expr, ..., kernel, window, keep_bands = FALSE, pad = NA) { stopifnot(is.cube(x)) + + pad_fill = as.numeric(0) + pad_mode = "" + if (is.na(pad)) { + pad_mode = "" + } + else if (is.numeric(pad)) { + pad_fill = pad[1] + pad_mode = "CONSTANT" + } + else if (is.character(pad)) { + if (any(c("REPLICATE","replicate","edge") %in% pad)){ + pad_mode = "REPLICATE" + } + else if (any(c("REFLECT","reflect","symmetric") %in% pad)){ + pad_mode = "REFLECT" + } + else if (any(c("REFLECT_PIXEL","reflect_pixel") %in% pad)){ + pad_mode = "REFLECT_PIXEL" + } + else { + warning("Unknown padding method (argument pad) provided: falling back to default method (no padding)") + pad_mode = "" + } + } + else { + warning("Invalid padding method (argument pad) provided: falling back to default method (no padding)") + pad_mode = "" + } + if (!missing(kernel)) { if (!missing(expr)) { warning("argument expr will be ignored, applying kernel convolution") @@ -178,7 +208,7 @@ window_space <- function(x, expr, ..., kernel, window, keep_bands = FALSE) { if (!is.matrix(kernel)) { stop("Kernel must be provided as a matrix") } - x = gc_create_window_space_cube_kernel(x, as.double(kernel), as.integer(nrow(kernel)), as.integer(ncol(kernel)), keep_bands) + x = gc_create_window_space_cube_kernel(x, as.double(kernel), as.integer(nrow(kernel)), as.integer(ncol(kernel)), keep_bands, as.character(pad_mode), as.double(pad_fill)) } else { stopifnot(is.character(expr)) @@ -198,7 +228,7 @@ window_space <- function(x, expr, ..., kernel, window, keep_bands = FALSE) { reducers = gsub("\\(.*\\)", "", expr) bands = gsub("[\\(\\)]", "", regmatches(expr, gregexpr("\\(.*?\\)", expr))) stopifnot(length(reducers) == length(bands)) - x = gc_create_window_space_cube_reduce(x, reducers, bands, as.integer(window[1]), as.integer(window[2]), keep_bands) + x = gc_create_window_space_cube_reduce(x, reducers, bands, as.integer(window[1]), as.integer(window[2]), keep_bands, as.character(pad_mode), as.double(pad_fill)) } class(x) <- c("window_space_cube", "cube", "xptr") return(x) diff --git a/man/window_space.Rd b/man/window_space.Rd index f077d24..1a0f3e9 100644 --- a/man/window_space.Rd +++ b/man/window_space.Rd @@ -4,7 +4,7 @@ \alias{window_space} \title{Apply a moving window operation or convolution kernel over the spatial dimensions of a data cube} \usage{ -window_space(x, expr, ..., kernel, window, keep_bands = FALSE) +window_space(x, expr, ..., kernel, window, keep_bands = FALSE, pad = NA) } \arguments{ \item{x}{source data cube} @@ -18,6 +18,8 @@ window_space(x, expr, ..., kernel, window, keep_bands = FALSE) \item{window}{integer vector with two elements defining the size of the window before and after a cell, the total size of the window is window[1] + 1 + window[2]} \item{keep_bands}{logical; if FALSE (the default), original data cube bands will be dropped.} + +\item{pad}{Padding method applied to the borders; use NULL for no padding, a numeric a fill value, or one of "REPLICATE", "REFLECT", "REFLECT_PIXEL"} } \value{ proxy data cube object @@ -41,9 +43,7 @@ more complex functions or arguments. Possible reducers currently include "min", \note{ Implemented reducers will ignore any NAN values (as \code{na.rm = TRUE} does). -THIS FUNCTION IS STILL EXPERIMENTAL: The current implementation does not yet support any special edge handling. - -THIS FUNCTION IS STILL EXPERIMENTAL: The current implementation will simply ignore NA values during the convolution, i.e. the sum will not be updated for NA pixels. +THIS FUNCTION IS STILL EXPERIMENTAL. This function returns a proxy object, i.e., it will not start any computations besides deriving the shape of the result. } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 327e391..27119e5 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -448,8 +448,8 @@ BEGIN_RCPP END_RCPP } // gc_create_window_space_cube_reduce -SEXP gc_create_window_space_cube_reduce(SEXP pin, std::vector reducers, std::vector bands, int win_size_y, int win_size_x, bool keep_bands); -RcppExport SEXP _gdalcubes_gc_create_window_space_cube_reduce(SEXP pinSEXP, SEXP reducersSEXP, SEXP bandsSEXP, SEXP win_size_ySEXP, SEXP win_size_xSEXP, SEXP keep_bandsSEXP) { +SEXP gc_create_window_space_cube_reduce(SEXP pin, std::vector reducers, std::vector bands, int win_size_y, int win_size_x, bool keep_bands, std::string pad_mode, double pad_fill); +RcppExport SEXP _gdalcubes_gc_create_window_space_cube_reduce(SEXP pinSEXP, SEXP reducersSEXP, SEXP bandsSEXP, SEXP win_size_ySEXP, SEXP win_size_xSEXP, SEXP keep_bandsSEXP, SEXP pad_modeSEXP, SEXP pad_fillSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -459,13 +459,15 @@ BEGIN_RCPP Rcpp::traits::input_parameter< int >::type win_size_y(win_size_ySEXP); Rcpp::traits::input_parameter< int >::type win_size_x(win_size_xSEXP); Rcpp::traits::input_parameter< bool >::type keep_bands(keep_bandsSEXP); - rcpp_result_gen = Rcpp::wrap(gc_create_window_space_cube_reduce(pin, reducers, bands, win_size_y, win_size_x, keep_bands)); + Rcpp::traits::input_parameter< std::string >::type pad_mode(pad_modeSEXP); + Rcpp::traits::input_parameter< double >::type pad_fill(pad_fillSEXP); + rcpp_result_gen = Rcpp::wrap(gc_create_window_space_cube_reduce(pin, reducers, bands, win_size_y, win_size_x, keep_bands, pad_mode, pad_fill)); return rcpp_result_gen; END_RCPP } // gc_create_window_space_cube_kernel -SEXP gc_create_window_space_cube_kernel(SEXP pin, std::vector kernel, int win_size_y, int win_size_x, bool keep_bands); -RcppExport SEXP _gdalcubes_gc_create_window_space_cube_kernel(SEXP pinSEXP, SEXP kernelSEXP, SEXP win_size_ySEXP, SEXP win_size_xSEXP, SEXP keep_bandsSEXP) { +SEXP gc_create_window_space_cube_kernel(SEXP pin, std::vector kernel, int win_size_y, int win_size_x, bool keep_bands, std::string pad_mode, double pad_fill); +RcppExport SEXP _gdalcubes_gc_create_window_space_cube_kernel(SEXP pinSEXP, SEXP kernelSEXP, SEXP win_size_ySEXP, SEXP win_size_xSEXP, SEXP keep_bandsSEXP, SEXP pad_modeSEXP, SEXP pad_fillSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -474,7 +476,9 @@ BEGIN_RCPP Rcpp::traits::input_parameter< int >::type win_size_y(win_size_ySEXP); Rcpp::traits::input_parameter< int >::type win_size_x(win_size_xSEXP); Rcpp::traits::input_parameter< bool >::type keep_bands(keep_bandsSEXP); - rcpp_result_gen = Rcpp::wrap(gc_create_window_space_cube_kernel(pin, kernel, win_size_y, win_size_x, keep_bands)); + Rcpp::traits::input_parameter< std::string >::type pad_mode(pad_modeSEXP); + Rcpp::traits::input_parameter< double >::type pad_fill(pad_fillSEXP); + rcpp_result_gen = Rcpp::wrap(gc_create_window_space_cube_kernel(pin, kernel, win_size_y, win_size_x, keep_bands, pad_mode, pad_fill)); return rcpp_result_gen; END_RCPP } @@ -885,8 +889,8 @@ static const R_CallMethodDef CallEntries[] = { {"_gdalcubes_gc_create_reduce_space_cube", (DL_FUNC) &_gdalcubes_gc_create_reduce_space_cube, 4}, {"_gdalcubes_gc_create_window_time_cube_reduce", (DL_FUNC) &_gdalcubes_gc_create_window_time_cube_reduce, 4}, {"_gdalcubes_gc_create_window_time_cube_kernel", (DL_FUNC) &_gdalcubes_gc_create_window_time_cube_kernel, 3}, - {"_gdalcubes_gc_create_window_space_cube_reduce", (DL_FUNC) &_gdalcubes_gc_create_window_space_cube_reduce, 6}, - {"_gdalcubes_gc_create_window_space_cube_kernel", (DL_FUNC) &_gdalcubes_gc_create_window_space_cube_kernel, 5}, + {"_gdalcubes_gc_create_window_space_cube_reduce", (DL_FUNC) &_gdalcubes_gc_create_window_space_cube_reduce, 8}, + {"_gdalcubes_gc_create_window_space_cube_kernel", (DL_FUNC) &_gdalcubes_gc_create_window_space_cube_kernel, 7}, {"_gdalcubes_gc_create_join_bands_cube", (DL_FUNC) &_gdalcubes_gc_create_join_bands_cube, 2}, {"_gdalcubes_gc_create_select_bands_cube", (DL_FUNC) &_gdalcubes_gc_create_select_bands_cube, 2}, {"_gdalcubes_gc_create_select_time_cube", (DL_FUNC) &_gdalcubes_gc_create_select_time_cube, 2}, diff --git a/src/gdalcubes.cpp b/src/gdalcubes.cpp index 9eceed6..11c0f2c 100644 --- a/src/gdalcubes.cpp +++ b/src/gdalcubes.cpp @@ -1144,7 +1144,7 @@ SEXP gc_create_window_time_cube_kernel(SEXP pin, std::vector window, std::v } // [[Rcpp::export]] -SEXP gc_create_window_space_cube_reduce(SEXP pin, std::vector reducers, std::vector bands, int win_size_y, int win_size_x, bool keep_bands) { +SEXP gc_create_window_space_cube_reduce(SEXP pin, std::vector reducers, std::vector bands, int win_size_y, int win_size_x, bool keep_bands, std::string pad_mode, double pad_fill) { try { Rcpp::XPtr< std::shared_ptr > aa = Rcpp::as>>(pin); @@ -1154,7 +1154,7 @@ SEXP gc_create_window_space_cube_reduce(SEXP pin, std::vector reduc reducer_bands.push_back(std::make_pair(reducers[i], bands[i])); } - std::shared_ptr* x = new std::shared_ptr(window_space_cube::create(*aa, reducer_bands, win_size_y, win_size_x, keep_bands)); + std::shared_ptr* x = new std::shared_ptr(window_space_cube::create(*aa, reducer_bands, win_size_y, win_size_x, keep_bands, pad_mode, pad_fill)); Rcpp::XPtr< std::shared_ptr > p(x, true) ; return p; @@ -1166,11 +1166,11 @@ SEXP gc_create_window_space_cube_reduce(SEXP pin, std::vector reduc } // [[Rcpp::export]] -SEXP gc_create_window_space_cube_kernel(SEXP pin, std::vector kernel, int win_size_y, int win_size_x, bool keep_bands) { +SEXP gc_create_window_space_cube_kernel(SEXP pin, std::vector kernel, int win_size_y, int win_size_x, bool keep_bands, std::string pad_mode, double pad_fill) { try { Rcpp::XPtr< std::shared_ptr > aa = Rcpp::as>>(pin); - std::shared_ptr* x = new std::shared_ptr(window_space_cube::create(*aa, kernel, win_size_y, win_size_x, keep_bands)); + std::shared_ptr* x = new std::shared_ptr(window_space_cube::create(*aa, kernel, win_size_y, win_size_x, keep_bands, pad_mode, pad_fill)); Rcpp::XPtr< std::shared_ptr > p(x, true) ; return p; diff --git a/src/gdalcubes/src/cube_factory.cpp b/src/gdalcubes/src/cube_factory.cpp index c7ce051..a29368e 100644 --- a/src/gdalcubes/src/cube_factory.cpp +++ b/src/gdalcubes/src/cube_factory.cpp @@ -153,7 +153,7 @@ void cube_factory::register_default() { } return window_space_cube::create(instance()->create_from_json(j["in_cube"]), kernel, j["win_size_y"].int_value(), j["win_size_x"].int_value(), - j["keep_bands"].bool_value()); + j["keep_bands"].bool_value(), j["pad_str"].string_value(), j["pad_fill"].number_value()); } else { std::vector> band_reducers; @@ -162,7 +162,7 @@ void cube_factory::register_default() { } return window_space_cube::create(instance()->create_from_json(j["in_cube"]), band_reducers, j["win_size_y"].int_value(), j["win_size_x"].int_value(), - j["keep_bands"].bool_value()); + j["keep_bands"].bool_value(), j["pad_str"].string_value(), j["pad_fill"].number_value()); } })); cube_generators.insert(std::make_pair(json11::Json&)>>( diff --git a/src/gdalcubes/src/window_space.cpp b/src/gdalcubes/src/window_space.cpp index ab32048..385a897 100644 --- a/src/gdalcubes/src/window_space.cpp +++ b/src/gdalcubes/src/window_space.cpp @@ -183,7 +183,6 @@ struct window_reducer_sd : public window_reducer_var { - std::shared_ptr window_space_cube::read_chunk(chunkid_t id) { GCBS_TRACE("window_space_cube::read_chunk(" + std::to_string(id) + ")"); std::shared_ptr out = std::make_shared(); @@ -202,8 +201,7 @@ std::shared_ptr window_space_cube::read_chunk(chunkid_t id) { std::fill(begin, end, NAN); - // 1. Extract "subwindow" from input chunks --> function - // size: (winsize[0] + chunk_size[0] + winsize[2]) x (winsize[1] + chunk_size[1] + winsize[3]) + // 1. Extract "subwindow" from input chunks auto ccoords = chunk_coords_from_id(id); std::array lower = {int32_t(ccoords[0] * _chunk_size[0]), int32_t(ccoords[1] * _chunk_size[1] - ((_win_size_y - 1) / 2)), @@ -215,12 +213,373 @@ std::shared_ptr window_space_cube::read_chunk(chunkid_t id) { auto cwin = _in_cube->read_window(lower, upper); - - if (!_kernel.empty()) { - // TODO: Apply padding - // 2. Iterate over target buffer and apply kernel. + + // 2. Apply padding (if needed) + if (_pad.mode != padding::MODE::NONE) { + + double v; + if (_pad.mode == padding::MODE::CONSTANT) { + v = _pad.constant_value; + } + + + // lower y side + if (lower[1] < 0) { + int32_t wy = -lower[1]; + for (uint16_t ib=0; ib < cwin->size()[0]; ++ib) { + for (uint16_t it=0; it < cwin->size()[1]; ++it) { + for (int32_t ix=0; ix < int32_t(cwin->size()[3]); ++ix) { + for (int32_t off_y=1; off_y <= wy; ++off_y) { + int32_t tx = ix; + int32_t ty = wy - off_y; + int32_t sx = tx; + int32_t sy = ty; + if (_pad.mode == padding::MODE::REPLICATE) { + sx = tx; + sy = ty + off_y; + } + else if (_pad.mode == padding::MODE::REFLECT) { + sx = tx; + sy = ty + off_y + off_y - 1; + } + else if (_pad.mode == padding::MODE::REFLECT_PIXEL) { + sx = tx; + sy = ty + off_y + off_y; + } + if (_pad.mode != padding::MODE::CONSTANT) { + v = ((double*)(cwin->buf()))[ + ib * cwin->size()[1] * cwin->size()[2] * cwin->size()[3] + + it * cwin->size()[2] * cwin->size()[3] + + sy * cwin->size()[3] + + sx]; + } + ((double*)(cwin->buf()))[ + ib * cwin->size()[1] * cwin->size()[2] * cwin->size()[3] + + it * cwin->size()[2] * cwin->size()[3] + + ty * cwin->size()[3] + + tx ] = v; + } + } + } + } + } + + // upper y side + if (upper[1] >= int32_t(size_y())) { + int32_t wy = upper[1] - size_y() + 1; + for (uint16_t ib=0; ib < cwin->size()[0]; ++ib) { + for (uint16_t it=0; it < cwin->size()[1]; ++it) { + for (int32_t ix=0; ix < int32_t(cwin->size()[3]); ++ix) { + for (int32_t off_y=1; off_y <= wy; ++off_y) { + int32_t tx = ix; + int32_t ty = cwin->size()[2] - wy - 1 + off_y; + int32_t sx = tx; + int32_t sy = ty; + + if (_pad.mode == padding::MODE::REPLICATE) { + sx = tx; + sy = ty - off_y; + } + else if (_pad.mode == padding::MODE::REFLECT) { + sx = tx; + sy = ty - off_y - off_y + 1; + } + else if (_pad.mode == padding::MODE::REFLECT_PIXEL) { + sx = tx; + sy = ty - off_y - off_y; + } + if (_pad.mode != padding::MODE::CONSTANT) { + v = ((double*)(cwin->buf()))[ + ib * cwin->size()[1] * cwin->size()[2] * cwin->size()[3] + + it * cwin->size()[2] * cwin->size()[3] + + sy * cwin->size()[3] + + sx]; + } + ((double*)(cwin->buf()))[ + ib * cwin->size()[1] * cwin->size()[2] * cwin->size()[3] + + it * cwin->size()[2] * cwin->size()[3] + + ty * cwin->size()[3] + + tx ] = v; + } + } + } + } + } + + + + // lower x side + if (lower[2] < 0) { + int32_t wx = -lower[2]; + for (uint16_t ib=0; ib < cwin->size()[0]; ++ib) { + for (uint16_t it=0; it < cwin->size()[1]; ++it) { + for (int32_t iy=0; iy < int32_t(cwin->size()[2]); ++iy) { + for (int32_t off_x=1; off_x <= wx; ++off_x) { + + int32_t tx = wx - off_x; + int32_t ty = iy; + int32_t sx = tx; + int32_t sy = ty; + if (_pad.mode == padding::MODE::REPLICATE) { + sy = ty; + sx = tx + off_x; + } + else if (_pad.mode == padding::MODE::REFLECT) { + sx = tx + off_x + off_x - 1; + sy = ty; + } + else if (_pad.mode == padding::MODE::REFLECT_PIXEL) { + sx = tx + off_x + off_x; + sy = ty; + } + if (_pad.mode != padding::MODE::CONSTANT) { + v = ((double*)(cwin->buf()))[ + ib * cwin->size()[1] * cwin->size()[2] * cwin->size()[3] + + it * cwin->size()[2] * cwin->size()[3] + + sy * cwin->size()[3] + + sx]; + } + ((double*)(cwin->buf()))[ + ib * cwin->size()[1] * cwin->size()[2] * cwin->size()[3] + + it * cwin->size()[2] * cwin->size()[3] + + ty * cwin->size()[3] + + tx ] = v; + } + } + } + } + } + + // upper x side + if (upper[2] >= int32_t(size_x())) { + int32_t wx = upper[2] - size_x() + 1; + for (uint16_t ib=0; ib < cwin->size()[0]; ++ib) { + for (uint16_t it=0; it < cwin->size()[1]; ++it) { + for (int32_t iy=0; iy < int32_t(cwin->size()[2]); ++iy) { + for (int32_t off_x=1; off_x <= wx; ++off_x) { + int32_t tx = cwin->size()[3] - wx - 1 + off_x; + int32_t ty = iy; + int32_t sx = tx; + int32_t sy = ty; + + if (_pad.mode == padding::MODE::REPLICATE) { + sx = tx - off_x; + sy = ty; + } + else if (_pad.mode == padding::MODE::REFLECT) { + sx = tx - off_x - off_x + 1; + sy = ty; + } + else if (_pad.mode == padding::MODE::REFLECT_PIXEL) { + sx = tx - off_x - off_x; + sy = ty; + } + if (_pad.mode != padding::MODE::CONSTANT) { + v = ((double*)(cwin->buf()))[ + ib * cwin->size()[1] * cwin->size()[2] * cwin->size()[3] + + it * cwin->size()[2] * cwin->size()[3] + + sy * cwin->size()[3] + + sx]; + } + ((double*)(cwin->buf()))[ + ib * cwin->size()[1] * cwin->size()[2] * cwin->size()[3] + + it * cwin->size()[2] * cwin->size()[3] + + ty * cwin->size()[3] + + tx ] = v; + } + } + } + } + } + + // lower y, lower x side + if (lower[1] < 0 && lower[2] < 0) { + int32_t wx = -lower[2]; + int32_t wy = -lower[1]; + for (uint16_t ib=0; ib < cwin->size()[0]; ++ib) { + for (uint16_t it=0; it < cwin->size()[1]; ++it) { + for (int32_t off_y=1; off_y <= wy; ++off_y) { + for (int32_t off_x=1; off_x <= wx; ++off_x) { + int32_t tx = wx - off_x; + int32_t ty = wy - off_y; + int32_t sx = tx; + int32_t sy = ty; + if (_pad.mode == padding::MODE::REPLICATE) { + sx = tx + off_y; + sy = ty + off_y; + } + else if (_pad.mode == padding::MODE::REFLECT) { + sx = tx + off_x + off_x - 1; + sy = ty + off_y + off_y - 1; + } + else if (_pad.mode == padding::MODE::REFLECT_PIXEL) { + sx = tx + off_x + off_x; + sy = ty + off_y + off_y; + } + if (_pad.mode != padding::MODE::CONSTANT) { + v = ((double*)(cwin->buf()))[ + ib * cwin->size()[1] * cwin->size()[2] * cwin->size()[3] + + it * cwin->size()[2] * cwin->size()[3] + + sy * cwin->size()[3] + + sx]; + } + ((double*)(cwin->buf()))[ + ib * cwin->size()[1] * cwin->size()[2] * cwin->size()[3] + + it * cwin->size()[2] * cwin->size()[3] + + ty * cwin->size()[3] + + tx ] = v; + } + } + } + } + } + // upper y, lower x side + if (upper[1] >= int32_t(size_y()) && lower[2] < 0) { + int32_t wx = -lower[2]; + int32_t wy = upper[1] - size_y() + 1; + for (uint16_t ib=0; ib < cwin->size()[0]; ++ib) { + for (uint16_t it=0; it < cwin->size()[1]; ++it) { + for (int32_t off_y=1; off_y <= wy; ++off_y) { + for (int32_t off_x=1; off_x <= wx; ++off_x) { + int32_t tx = wx - off_x; + int32_t ty = cwin->size()[2] - wy - 1 + off_y; + int32_t sx = tx; + int32_t sy = ty; + if (_pad.mode == padding::MODE::REPLICATE) { + sx = tx + off_y; + sy = ty - off_y; + } + else if (_pad.mode == padding::MODE::REFLECT) { + sx = tx + off_x + off_x - 1; + sy = ty - off_y - off_y + 1; + } + else if (_pad.mode == padding::MODE::REFLECT_PIXEL) { + sx = tx + off_x + off_x; + sy = ty - off_y - off_y; + } + if (_pad.mode != padding::MODE::CONSTANT) { + v = ((double*)(cwin->buf()))[ + ib * cwin->size()[1] * cwin->size()[2] * cwin->size()[3] + + it * cwin->size()[2] * cwin->size()[3] + + sy * cwin->size()[3] + + sx]; + } + ((double*)(cwin->buf()))[ + ib * cwin->size()[1] * cwin->size()[2] * cwin->size()[3] + + it * cwin->size()[2] * cwin->size()[3] + + ty * cwin->size()[3] + + tx ] = v; + } + } + } + } + } + // upper y, upper x side + if (upper[1] >= int32_t(size_y()) && upper[2] >= int32_t(size_x())) { + int32_t wx = upper[2] - size_x() + 1; + int32_t wy = upper[1] - size_y() + 1; + for (uint16_t ib=0; ib < cwin->size()[0]; ++ib) { + for (uint16_t it=0; it < cwin->size()[1]; ++it) { + for (int32_t off_y=1; off_y <= wy; ++off_y) { + for (int32_t off_x=1; off_x <= wx; ++off_x) { + int32_t tx = cwin->size()[3] - wx - 1 + off_x; + int32_t ty = cwin->size()[2] - wy - 1 + off_y; + int32_t sx = tx; + int32_t sy = ty; + if (_pad.mode == padding::MODE::REPLICATE) { + sx = tx - off_x; + sy = ty - off_y; + } + else if (_pad.mode == padding::MODE::REFLECT) { + sx = tx - off_x - off_x + 1; + sy = ty - off_y - off_y + 1; + } + else if (_pad.mode == padding::MODE::REFLECT_PIXEL) { + sx = tx - off_x - off_x; + sy = ty - off_y - off_y; + } + if (_pad.mode != padding::MODE::CONSTANT) { + v = ((double*)(cwin->buf()))[ + ib * cwin->size()[1] * cwin->size()[2] * cwin->size()[3] + + it * cwin->size()[2] * cwin->size()[3] + + sy * cwin->size()[3] + + sx]; + } + ((double*)(cwin->buf()))[ + ib * cwin->size()[1] * cwin->size()[2] * cwin->size()[3] + + it * cwin->size()[2] * cwin->size()[3] + + ty * cwin->size()[3] + + tx ] = v; + } + } + } + } + } + // lower y, upper x side + if (upper[1] >= int32_t(size_y()) && upper[2] >= int32_t(size_x())) { + int32_t wx = upper[2] - size_x() + 1; + int32_t wy = -lower[1]; + for (uint16_t ib=0; ib < cwin->size()[0]; ++ib) { + for (uint16_t it=0; it < cwin->size()[1]; ++it) { + for (int32_t off_y=1; off_y <= wy; ++off_y) { + for (int32_t off_x=1; off_x <= wx; ++off_x) { + int32_t tx = cwin->size()[3] - wx - 1 + off_x; + int32_t ty = wy - off_y; + int32_t sx = tx; + int32_t sy = ty; + if (_pad.mode == padding::MODE::REPLICATE) { + sx = tx - off_x; + sy = ty + off_y; + } + else if (_pad.mode == padding::MODE::REFLECT) { + sx = tx - off_x - off_x + 1; + sy = ty + off_y + off_y - 1; + } + else if (_pad.mode == padding::MODE::REFLECT_PIXEL) { + sx = tx - off_x - off_x; + sy = ty + off_y + off_y; + } + if (_pad.mode != padding::MODE::CONSTANT) { + v = ((double*)(cwin->buf()))[ + ib * cwin->size()[1] * cwin->size()[2] * cwin->size()[3] + + it * cwin->size()[2] * cwin->size()[3] + + sy * cwin->size()[3] + + sx]; + } + ((double*)(cwin->buf()))[ + ib * cwin->size()[1] * cwin->size()[2] * cwin->size()[3] + + it * cwin->size()[2] * cwin->size()[3] + + ty * cwin->size()[3] + + tx ] = v; + } + } + } + } + } + } // end padding + + + + + + + + + + + + + + + + // 3. Apply kernel or aggregation + + // CASE 1: Convolution using a provided kernel + if (!_kernel.empty()) { + + // Iterate over target buffer and apply kernel. for (uint32_t ib = 0; ib < size_btyx[0]; ++ib) { for (uint32_t it = 0; it < size_btyx[1]; ++it) { for (uint32_t iy = 0; iy < size_btyx[2]; ++iy) { @@ -256,6 +615,12 @@ std::shared_ptr window_space_cube::read_chunk(chunkid_t id) { } } } + + + + + + // CASE 2: Aggregation of selected bands using built-in aggregation functions (see top of file) else { std::vector reducers; for (uint16_t i = 0; i < _reducer_bands.size(); ++i) { @@ -351,4 +716,5 @@ std::shared_ptr window_space_cube::read_chunk(chunkid_t id) { return out; } + } // namespace gdalcubes diff --git a/src/gdalcubes/src/window_space.h b/src/gdalcubes/src/window_space.h index 39041d2..9dd25c5 100644 --- a/src/gdalcubes/src/window_space.h +++ b/src/gdalcubes/src/window_space.h @@ -29,6 +29,11 @@ namespace gdalcubes { + + + + + // TODO: add custom reducer function // TODO: add boundary fill options (e.g. NA, wrap, repeat, mirror, ...) /** @@ -36,6 +41,14 @@ namespace gdalcubes { */ class window_space_cube : public cube { public: + + + struct padding { + enum MODE {NONE, CONSTANT, REPLICATE, REFLECT, REFLECT_PIXEL, WRAP}; // see https://processes.openeo.org/#apply_kernel + MODE mode; + double constant_value; + padding() : mode(NONE), constant_value(NAN) {} + }; /** * @brief Create a data cube that applies a reducer function on a given input data cube over time @@ -47,8 +60,8 @@ class window_space_cube : public cube { */ static std::shared_ptr create(std::shared_ptr in, std::vector> reducer_bands, - uint16_t win_size_y, uint16_t win_size_x, bool keep_bands) { - std::shared_ptr out = std::make_shared(in, reducer_bands, win_size_y, win_size_x, keep_bands); + uint16_t win_size_y, uint16_t win_size_x, bool keep_bands, std::string pad_str, double pad_fill=0.0) { + std::shared_ptr out = std::make_shared(in, reducer_bands, win_size_y, win_size_x, keep_bands, pad_str, pad_fill); in->add_child_cube(out); out->add_parent_cube(in); return out; @@ -63,8 +76,8 @@ class window_space_cube : public cube { * @return a shared pointer to the created data cube instance */ static std::shared_ptr - create(std::shared_ptr in, std::vector kernel, uint16_t win_size_y, uint16_t win_size_x, bool keep_bands) { - std::shared_ptr out = std::make_shared(in, kernel, win_size_y, win_size_x, keep_bands); + create(std::shared_ptr in, std::vector kernel, uint16_t win_size_y, uint16_t win_size_x, bool keep_bands, std::string pad_str, double pad_fill=0.0) { + std::shared_ptr out = std::make_shared(in, kernel, win_size_y, win_size_x, keep_bands, pad_str, pad_fill); in->add_child_cube(out); out->add_parent_cube(in); return out; @@ -73,7 +86,7 @@ class window_space_cube : public cube { public: window_space_cube(std::shared_ptr in, std::vector> reducer_bands, - uint16_t win_size_y, uint16_t win_size_x, bool keep_bands) : cube(in->st_reference()->copy()), _in_cube(in), _reducer_bands(reducer_bands), _win_size_y(win_size_y), _win_size_x(win_size_x), _band_idx_in(), _kernel(), _keep_bands(keep_bands) { // it is important to duplicate st reference here, otherwise changes will affect input cube as well + uint16_t win_size_y, uint16_t win_size_x, bool keep_bands, std::string pad_str, double pad_fill=0.0) : cube(in->st_reference()->copy()), _in_cube(in), _reducer_bands(reducer_bands), _win_size_y(win_size_y), _win_size_x(win_size_x), _band_idx_in(), _kernel(), _keep_bands(keep_bands), _pad_str(pad_str), _pad_fill(pad_fill), _pad() { // it is important to duplicate st reference here, otherwise changes will affect input cube as well _chunk_size[0] = _in_cube->chunk_size()[0]; _chunk_size[1] = _in_cube->chunk_size()[1]; _chunk_size[2] = _in_cube->chunk_size()[2]; @@ -100,11 +113,30 @@ class window_space_cube : public cube { _band_idx_in.push_back(in->bands().get_index(bandstr)); } + + if (pad_str == "CONSTANT") { + _pad.mode = padding::MODE::CONSTANT; + _pad.constant_value = pad_fill; + } + else if (pad_str == "REPLICATE") { + _pad.mode = padding::MODE::REPLICATE; + } + else if (pad_str == "REFLECT") { + _pad.mode = padding::MODE::REFLECT; + } + else if (pad_str == "REFLECT_PIXEL") { + _pad.mode = padding::MODE::REFLECT_PIXEL; + } + else { + if (!pad_str.empty()) + GCBS_WARN("Unknown padding method defined: falling back to default method (no padding)"); + _pad.mode = padding::MODE::NONE; + } } - window_space_cube(std::shared_ptr in, std::vector kernel, uint16_t win_size_y, uint16_t win_size_x, bool keep_bands) - : cube(in->st_reference()->copy()), _in_cube(in), _reducer_bands(), _win_size_y(win_size_y), _win_size_x(win_size_x), _band_idx_in(), _kernel(kernel), _keep_bands(keep_bands) { // it is important to duplicate st reference here, otherwise changes will affect input cube as well + window_space_cube(std::shared_ptr in, std::vector kernel, uint16_t win_size_y, uint16_t win_size_x, bool keep_bands, std::string pad_str, double pad_fill=0.0) + : cube(in->st_reference()->copy()), _in_cube(in), _reducer_bands(), _win_size_y(win_size_y), _win_size_x(win_size_x), _band_idx_in(), _kernel(kernel), _keep_bands(keep_bands), _pad_str(pad_str), _pad_fill(pad_fill), _pad() { // it is important to duplicate st reference here, otherwise changes will affect input cube as well _chunk_size[0] = _in_cube->chunk_size()[0]; _chunk_size[1] = _in_cube->chunk_size()[1]; _chunk_size[2] = _in_cube->chunk_size()[2]; @@ -125,6 +157,25 @@ class window_space_cube : public cube { band b = in->bands().get(i); _bands.add(b); } + + if (pad_str == "CONSTANT") { + _pad.mode = padding::MODE::CONSTANT; + _pad.constant_value = pad_fill; + } + else if (pad_str == "REPLICATE") { + _pad.mode = padding::MODE::REPLICATE; + } + else if (pad_str == "REFLECT") { + _pad.mode = padding::MODE::REFLECT; + } + else if (pad_str == "REFLECT_PIXEL") { + _pad.mode = padding::MODE::REFLECT_PIXEL; + } + else { + if (!pad_str.empty()) + GCBS_WARN("Unknown padding method defined: falling back to default method (no padding)"); + _pad.mode = padding::MODE::NONE; + } } public: @@ -147,6 +198,8 @@ class window_space_cube : public cube { } out["win_size_y"] = _win_size_y; out["win_size_x"] = _win_size_x; + out["pad_str"] = _pad_str; + out["pad_fill"] = _pad_fill; out["keep_bands"] = _keep_bands; out["in_cube"] = _in_cube->make_constructible_json(); return out; @@ -160,6 +213,9 @@ class window_space_cube : public cube { std::vector _band_idx_in; std::vector _kernel; bool _keep_bands; + std::string _pad_str; + double _pad_fill; + padding _pad; // TODO }; } // namespace gdalcubes From e1b19b45109c26ab18dca6a565f2b972adf11760 Mon Sep 17 00:00:00 2001 From: Marius Appel Date: Tue, 27 Feb 2024 23:29:43 +0100 Subject: [PATCH 3/4] add support for curvilinear grids in source images (e.g. S5P netCDFs) [experimental] --- src/gdalcubes/src/image_collection.cpp | 59 +++++- src/gdalcubes/src/warp.cpp | 261 ++++++++++++++++++++++++- src/gdalcubes/src/warp.h | 39 ++++ 3 files changed, 351 insertions(+), 8 deletions(-) diff --git a/src/gdalcubes/src/image_collection.cpp b/src/gdalcubes/src/image_collection.cpp index 1b4c896..ed9071a 100644 --- a/src/gdalcubes/src/image_collection.cpp +++ b/src/gdalcubes/src/image_collection.cpp @@ -852,12 +852,61 @@ void image_collection::add_with_collection_format(std::vector descr bbox.bottom = extent[1]; } - } else { - GDALClose((GDALDatasetH)dataset); - if (strict) throw std::string("ERROR in image_collection::add(): GDAL cannot derive spatial extent for '" + *it + "'."); - GCBS_WARN("Failed to derive spatial extent from " + *it); - continue; } + else { + char **slist = dataset->GetMetadata("GEOLOCATION"); + if (slist != NULL) { + std::string x_dataset = std::string(CSLFetchNameValue(slist, "X_DATASET")); + std::string y_dataset = std::string(CSLFetchNameValue(slist, "Y_DATASET")); + std::string geoloc_srs_str = std::string(CSLFetchNameValue(slist, "SRS")); + int16_t x_band = std::stoi(std::string(CSLFetchNameValue(slist, "X_BAND"))); + int16_t y_band = std::stoi(std::string(CSLFetchNameValue(slist, "Y_BAND"))); + + double adfMinMax[2]; + int bGotMin, bGotMax; + + GDALDataset* gd_x = (GDALDataset*)GDALOpen(x_dataset.c_str(), GA_ReadOnly); + if (!gd_x) { + GDALClose((GDALDatasetH)dataset); + if (strict) throw std::string("ERROR in image_collection::add(): GDAL cannot open '" + x_dataset + "'."); + GCBS_WARN("GDAL failed to open " + x_dataset); + continue; + } + GDALRasterBand* b_x = gd_x->GetRasterBand(x_band); + adfMinMax[0] = b_x->GetMinimum( &bGotMin ); + adfMinMax[1] = b_x->GetMaximum( &bGotMax ); + if(!(bGotMin && bGotMax)) + GDALComputeRasterMinMax((GDALRasterBandH)b_x, FALSE, adfMinMax); + bbox.left = adfMinMax[0]; + bbox.right = adfMinMax[1]; + GDALClose(gd_x); + + GDALDataset* gd_y = (GDALDataset*)GDALOpen(y_dataset.c_str(), GA_ReadOnly); + if (!gd_y) { + GDALClose((GDALDatasetH)dataset); + if (strict) throw std::string("ERROR in image_collection::add(): GDAL cannot open '" + y_dataset + "'."); + GCBS_WARN("GDAL failed to open " + y_dataset); + continue; + } + GDALRasterBand* b_y = gd_y->GetRasterBand(y_band); + adfMinMax[0] = b_y->GetMinimum( &bGotMin ); + adfMinMax[1] = b_y->GetMaximum( &bGotMax ); + if(!(bGotMin && bGotMax)) + GDALComputeRasterMinMax((GDALRasterBandH)b_y, FALSE, adfMinMax); + bbox.bottom = adfMinMax[0]; + bbox.top = adfMinMax[1]; + GDALClose(gd_y); + + bbox.transform(geoloc_srs_str, "EPSG:4326"); + } + else { // No extent??? + GDALClose((GDALDatasetH)dataset); + if (strict) throw std::string("ERROR in image_collection::add(): GDAL cannot derive spatial extent for '" + *it + "'."); + GCBS_WARN("Failed to derive spatial extent from " + *it); + continue; + } + + } } else { bbox.left = affine_in[0]; bbox.right = affine_in[0] + affine_in[1] * dataset->GetRasterXSize() + affine_in[2] * dataset->GetRasterYSize(); diff --git a/src/gdalcubes/src/warp.cpp b/src/gdalcubes/src/warp.cpp index 7511674..3316940 100644 --- a/src/gdalcubes/src/warp.cpp +++ b/src/gdalcubes/src/warp.cpp @@ -54,9 +54,25 @@ gdalwarp_client::gdal_transformation_cache::~gdal_transformation_cache() { } } + GDALDataset *gdalwarp_client::warp(GDALDataset *in, std::string s_srs, std::string t_srs, double te_left, double te_right, double te_top, double te_bottom, uint32_t ts_x, uint32_t ts_y, std::string resampling, std::vector srcnodata) { + double temp[6]; + if (in->GetGeoTransform(temp) == CE_None) { + return warp_simple(in, s_srs, t_srs, te_left, te_right, te_top, te_bottom, ts_x, ts_y, resampling, srcnodata); + } + else { // GCPs, RCPs, or GeolocationArrays + return warp_complex(in, s_srs, t_srs, te_left, te_right, te_top, te_bottom, ts_x, ts_y, resampling, srcnodata); + } +} + + + + +GDALDataset *gdalwarp_client::warp_simple(GDALDataset *in, std::string s_srs, std::string t_srs, double te_left, + double te_right, double te_top, double te_bottom, uint32_t ts_x, uint32_t ts_y, + std::string resampling, std::vector srcnodata) { char *wkt_out = NULL; OGRSpatialReference srs_out; @@ -87,6 +103,9 @@ GDALDataset *gdalwarp_client::warp(GDALDataset *in, std::string s_srs, std::stri if (s_srs.empty()) { // try reading from dataset s_srs = std::string(in->GetProjectionRef()); + if (s_srs.empty()) { + GCBS_DEBUG("Failed to read input crs."); + } } // if (std::string(in->GetProjectionRef()).empty()) { // if (in->SetProjection(s_srs.c_str()) != CE_None) { @@ -234,9 +253,7 @@ GDALDataset *gdalwarp_client::warp(GDALDataset *in, std::string s_srs, std::stri if (oOperation.Initialize(psWarpOptions) != CE_None) { GCBS_ERROR("Initialization of gdalwarp failed"); } - oOperation.ChunkAndWarpImage(0, 0, - GDALGetRasterXSize(out), - GDALGetRasterYSize(out)); + oOperation.ChunkAndWarpImage(0, 0, GDALGetRasterXSize(out), GDALGetRasterYSize(out)); destroy_transform((gdalwarp_client::gdalcubes_transform_info *)psWarpOptions->pTransformerArg); GDALDestroyWarpOptions(psWarpOptions); @@ -416,4 +433,242 @@ void gdalwarp_client::destroy_reprojection(gdalcubes_reprojection_info *reprojec } } +GDALDataset *gdalwarp_client::warp_complex(GDALDataset *in, std::string s_srs, std::string t_srs, double te_left, + double te_right, double te_top, double te_bottom, uint32_t ts_x, uint32_t ts_y, + std::string resampling, std::vector srcnodata) { + char *wkt_out = NULL; + + OGRSpatialReference srs_out; + srs_out.SetFromUserInput(t_srs.c_str()); + srs_out.exportToWkt(&wkt_out); + + double dst_geotransform[6]; + dst_geotransform[0] = te_left; + dst_geotransform[1] = (te_right - te_left) / double(ts_x); + dst_geotransform[2] = 0.0; + dst_geotransform[3] = te_top; + dst_geotransform[4] = 0.0; + dst_geotransform[5] = (te_bottom - te_top) / double(ts_y); + + GDALDriver *mem_driver = (GDALDriver *)GDALGetDriverByName("MEM"); + if (mem_driver == NULL) { + GCBS_ERROR("Cannot find GDAL MEM driver"); + throw std::string("Cannot find GDAL MEM driver"); + } + + GDALDataset *out = mem_driver->Create("", ts_x, ts_y, in->GetRasterCount(), GDT_Float64, NULL); + + for (uint16_t i=0; iGetRasterCount(); ++i) { + out->GetRasterBand(i+1)->Fill(NAN); // Avoid spurious output when warping fails. + } + + out->SetProjection(wkt_out); + out->SetGeoTransform(dst_geotransform); + + + // Setup warp options. + GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions(); + if (s_srs.empty()) { + // try reading from dataset + s_srs = std::string(in->GetProjectionRef()); + if (s_srs.empty()) { + GCBS_DEBUG("Failed to read input crs; assuming EPSG:4326"); + s_srs = "EPSG:4326"; + } + } + // if (std::string(in->GetProjectionRef()).empty()) { + // if (in->SetProjection(s_srs.c_str()) != CE_None) { + // GCBS_DEBUG("Failed to overwrite projection for dataset."); + // } + // } + // std::string s = (in->GetProjectionRef()); + // GCBS_DEBUG(s); + psWarpOptions->hSrcDS = in; + psWarpOptions->hDstDS = out; + psWarpOptions->pfnProgress = GDALDummyProgress; + + // see https://github.com/OSGeo/gdal/blob/64cf9b4e889c93e34177237665fe842186d1f581/alg/gdaltransformer.cpp#L1274C5-L1275C67 + GDALGenImgProjTransformInfo *psInfo = static_cast(GDALCreateGenImgProjTransformer(in, s_srs.c_str(), out, NULL, TRUE, 0.0, 1 )); + + // TODO: Do we need to check if t_srs and/or s_src are empty? + // Overwrite reprojection transform if needed, to benefit from cache + OGRSpatialReference srs_in; + srs_in.SetFromUserInput(s_srs.c_str()); + void *tmparg = psInfo->pReprojectArg; // must be freed later + if (!srs_in.IsSame(&srs_out)) { + psInfo->pReprojectArg = gdal_transformation_cache::instance()->get(s_srs, t_srs); + psInfo->pReproject = reproject; + } + + psWarpOptions->pTransformerArg = (void*)psInfo; + psWarpOptions->pfnTransformer = GDALGenImgProjTransform; + + // Derive best overview level to use + int n_ov = in->GetRasterBand(1)->GetOverviewCount(); + if (config::instance()->get_gdal_use_overviews() && n_ov > 0) { + double *x = (double *)std::malloc(sizeof(double) * 4); + double *y = (double *)std::malloc(sizeof(double) * 4); + int *succ = (int *)std::malloc(sizeof(int) * 4); + x[0] = 0; + x[1] = 0; + x[2] = ts_x; + x[3] = ts_x; + y[0] = ts_y; + y[1] = 0; + y[2] = ts_y; + y[3] = 0; + + // Transform corners + GDALGenImgProjTransform(psWarpOptions->pTransformerArg, 1, 4, x, y, NULL, succ); + + double minx = std::min(std::min(x[0], x[1]), std::min(x[2], x[3])); + double maxx = std::max(std::max(x[0], x[1]), std::max(x[2], x[3])); + double target_ratio = (maxx - minx) / double(ts_x); + + int16_t ilevel = 0; + while (ilevel < n_ov) { + double ov_ratio = double(in->GetRasterBand(1)->GetXSize()) / double(in->GetRasterBand(1)->GetOverview(ilevel)->GetXSize()); + if (ov_ratio > target_ratio) { + --ilevel; + break; + } + ++ilevel; + } + if (ilevel >= n_ov) { + ilevel = n_ov - 1; + } + if (ilevel >= 0) { + //GCBS_TRACE("Using overview level" + std::to_string(ilevel)); + char **oo = nullptr; + oo = CSLAddString(oo, ("OVERVIEW_LEVEL=" + std::to_string(ilevel)).c_str()); + + std::string descr = in->GetDescription(); + GDALClose(in); + in = (GDALDataset *)GDALOpenEx(descr.c_str(), GDAL_OF_RASTER | GDAL_OF_READONLY, NULL, oo, NULL); + if (in != NULL) { + destroy_transform((gdalwarp_client::gdalcubes_transform_info *)psWarpOptions->pTransformerArg); + psWarpOptions->pTransformerArg = create_transform(in, out, s_srs, t_srs); + psWarpOptions->hSrcDS = in; // TODO: close in_ov + } else { + GCBS_WARN("Failed to open GDAL overview dataset for '" + descr + "', using original full resolution image."); + } + CSLDestroy(oo); + } + + std::free(x); + std::free(y); + std::free(succ); + } + + psWarpOptions->eResampleAlg = GDALResampleAlg::GRA_NearestNeighbour; + if (resampling == "bilinear") { + psWarpOptions->eResampleAlg = GDALResampleAlg::GRA_Bilinear; + } else if (resampling == "cubic") { + psWarpOptions->eResampleAlg = GDALResampleAlg::GRA_Cubic; + } else if (resampling == "cubicspline") { + psWarpOptions->eResampleAlg = GDALResampleAlg::GRA_CubicSpline; + } else if (resampling == "lanczos") { + psWarpOptions->eResampleAlg = GDALResampleAlg::GRA_Lanczos; + } else if (resampling == "average") { + psWarpOptions->eResampleAlg = GDALResampleAlg::GRA_Average; + } else if (resampling == "mode") { + psWarpOptions->eResampleAlg = GDALResampleAlg::GRA_Mode; + } else if (resampling == "max") { + psWarpOptions->eResampleAlg = GDALResampleAlg::GRA_Max; + } else if (resampling == "min") { + psWarpOptions->eResampleAlg = GDALResampleAlg::GRA_Min; + } else if (resampling == "med") { + psWarpOptions->eResampleAlg = GDALResampleAlg::GRA_Med; + } else if (resampling == "q1") { + psWarpOptions->eResampleAlg = GDALResampleAlg::GRA_Q1; + } else if (resampling == "q3") { + psWarpOptions->eResampleAlg = GDALResampleAlg::GRA_Q3; + } + + psWarpOptions->nBandCount = in->GetRasterCount(); + psWarpOptions->panSrcBands = (int *)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount); + psWarpOptions->panDstBands = (int *)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount); + double *dst_nodata = (double *)CPLMalloc(sizeof(double) * psWarpOptions->nBandCount); + + // Due to issues on Windows with GDAL 2.2.3 if imag no data values are not set, + // we explicitly set then to 0 in the following. This might be removed in the future. + double *dst_nodata_img = (double *)CPLMalloc(sizeof(double) * psWarpOptions->nBandCount); + for (uint16_t i = 0; i < psWarpOptions->nBandCount; ++i) { + dst_nodata[i] = NAN; + dst_nodata_img[i] = 0.0; + psWarpOptions->panSrcBands[i] = i + 1; + psWarpOptions->panDstBands[i] = i + 1; + } + double *src_nodata = nullptr; + double *src_nodata_img = nullptr; + + if (!srcnodata.empty()) { + src_nodata = (double *)CPLMalloc(sizeof(double) * psWarpOptions->nBandCount); + src_nodata_img = (double *)CPLMalloc(sizeof(double) * psWarpOptions->nBandCount); + if (srcnodata.size() == 1) { + for (uint16_t i = 0; i < psWarpOptions->nBandCount; ++i) { + src_nodata[i] = srcnodata[0]; + src_nodata_img[i] = 0.0; + } + psWarpOptions->padfSrcNoDataReal = src_nodata; + psWarpOptions->padfSrcNoDataImag = src_nodata_img; + } else if (srcnodata.size() == (uint16_t)psWarpOptions->nBandCount) { + for (uint16_t i = 0; i < psWarpOptions->nBandCount; ++i) { + src_nodata[i] = srcnodata[i]; + src_nodata_img[i] = 0.0; + } + psWarpOptions->padfSrcNoDataReal = src_nodata; + psWarpOptions->padfSrcNoDataImag = src_nodata_img; + } else { + GCBS_DEBUG("Number of no data values does not match number of bands of source dataset, no data values will be ignored"); + std::free(src_nodata); + std::free(src_nodata_img); + } + } + psWarpOptions->padfDstNoDataReal = dst_nodata; + psWarpOptions->padfDstNoDataImag = dst_nodata_img; + + char **wo = nullptr; + wo = CSLAddString(wo, "INIT_DEST=nan"); + wo = CSLAddString(wo, ("NUM_THREADS=" + std::to_string(config::instance()->get_gdal_num_threads())).c_str()); + psWarpOptions->papszWarpOptions = wo; + + // Initialize and execute the warp operation. + GDALWarpOperation oOperation; + if (oOperation.Initialize(psWarpOptions) != CE_None) { + GCBS_ERROR("Initialization of gdalwarp failed"); + } + + if(oOperation.ChunkAndWarpImage(0, 0, + GDALGetRasterXSize(out), + GDALGetRasterYSize(out)) != CE_None) { + GCBS_DEBUG("gdalwarp returned error code #" + std::to_string(CPLGetLastErrorNo()) + ": " + std::string(CPLGetLastErrorMsg()) + + "[" + std::string(in->GetDescription()) + "]"); + + } + + static_cast(psWarpOptions->pTransformerArg)->pReprojectArg = tmparg; // safe release + GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg); + GDALDestroyWarpOptions(psWarpOptions); + + CPLFree(wkt_out); + + if (in) { + GDALClose(in); + } + return out; +} + + + + + + + } // namespace gdalcubes + + + + + + diff --git a/src/gdalcubes/src/warp.h b/src/gdalcubes/src/warp.h index 5404868..07cffad 100644 --- a/src/gdalcubes/src/warp.h +++ b/src/gdalcubes/src/warp.h @@ -96,6 +96,13 @@ class gdalwarp_client { */ static GDALDataset *warp(GDALDataset *in, std::string s_srs, std::string t_srs, double te_left, double te_right, double te_top, double te_bottom, uint32_t ts_x, uint32_t ts_y, std::string resampling, std::vector srcnodata); + + // See above, only for cases where source dataset has a GeoTransform (affine transformation) + static GDALDataset *warp_simple(GDALDataset *in, std::string s_srs, std::string t_srs, double te_left, double te_right, double te_top, double te_bottom, uint32_t ts_x, uint32_t ts_y, std::string resampling, std::vector srcnodata); + + // See above, but also working for source images with spatial reference by GCPs, RPCs, or GeoLocation arrays + static GDALDataset *warp_complex(GDALDataset *in, std::string s_srs, std::string t_srs, double te_left, double te_right, double te_top, double te_bottom, uint32_t ts_x, uint32_t ts_y, std::string resampling, std::vector srcnodata); + static gdalcubes_transform_info *create_transform(GDALDataset *in, GDALDataset *out, std::string srs_in_str, std::string srs_out_str); static void destroy_transform(gdalcubes_transform_info *transform); @@ -111,6 +118,38 @@ class gdalwarp_client { static int reproject(void *pTransformerArg, int bDstToSrc, int nPointCount, double *x, double *y, double *z = nullptr, int *panSuccess = nullptr); + + + private: + + // see https://github.com/OSGeo/gdal/blob/64cf9b4e889c93e34177237665fe842186d1f581/alg/gdaltransformer.cpp#L84 + typedef struct + { + + GDALTransformerInfo sTI; + + double adfSrcGeoTransform[6]; + double adfSrcInvGeoTransform[6]; + + void *pSrcTransformArg; + GDALTransformerFunc pSrcTransformer; + + void *pReprojectArg; + GDALTransformerFunc pReproject; + + double adfDstGeoTransform[6]; + double adfDstInvGeoTransform[6]; + + void *pDstTransformArg; + GDALTransformerFunc pDstTransformer; + + // Memorize the value of the CHECK_WITH_INVERT_PROJ at the time we + // instantiated the object, to be able to decide if + // GDALRefreshGenImgProjTransformer() must do something or not. + bool bCheckWithInvertPROJ; + + } GDALGenImgProjTransformInfo; + }; } // namespace gdalcubes From 4fab61618f0be7f4349c7636d0a009ef79ac87e7 Mon Sep 17 00:00:00 2001 From: Marius Appel Date: Thu, 29 Feb 2024 13:16:43 +0100 Subject: [PATCH 4/4] update NEWS.md; add window_space tests --- NEWS.md | 5 ++++ configure | 19 ++++++------ inst/tinytest/test_window_space.R | 49 +++++++++++++++++++++++++++++++ src/gdalcubes/src/utils.cpp | 10 +++++++ src/gdalcubes/src/utils.h | 1 + src/gdalcubes/src/warp.cpp | 12 ++++++-- 6 files changed, 84 insertions(+), 12 deletions(-) create mode 100644 inst/tinytest/test_window_space.R diff --git a/NEWS.md b/NEWS.md index ce63d3b..3ddfa17 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,6 +2,11 @@ * add `as.data.frame()` to easily convert data cubes to data frames * add `predict.cube()` to predict pixel values based on models (lm, glm, caret, tidymodels, and similar) +* add `window_space()` to apply (focal) moving window kernels or aggregation functions +* `extract()` now combines extracted values with input geometries and attributes (if `merge = TRUE`) +* add support for imagery with spatial reference from geolocation arrays (including curvilinear grids) +* `stac_image_collection()` now accepts STACItemCollection objects directly and should be more robust +* Windows build uses pkg-config if available # gdalcubes 0.6.4 (2023-04-14) diff --git a/configure b/configure index 6c728ea..d2ac17c 100755 --- a/configure +++ b/configure @@ -1,6 +1,6 @@ #! /bin/sh # Guess values for system-dependent variables and create Makefiles. -# Generated by GNU Autoconf 2.71 for gdalcubes 0.6.4. +# Generated by GNU Autoconf 2.71 for gdalcubes 0.6.9999. # # Report bugs to . # @@ -610,8 +610,8 @@ MAKEFLAGS= # Identity of this package. PACKAGE_NAME='gdalcubes' PACKAGE_TARNAME='gdalcubes' -PACKAGE_VERSION='0.6.4' -PACKAGE_STRING='gdalcubes 0.6.4' +PACKAGE_VERSION='0.6.9999' +PACKAGE_STRING='gdalcubes 0.6.9999' PACKAGE_BUGREPORT='marius.appel@hs-bochum.de' PACKAGE_URL='' @@ -1279,7 +1279,7 @@ if test "$ac_init_help" = "long"; then # Omit some internal or obsolete options to make the list less imposing. # This message is too long to be a string in the A/UX 3.1 sh. cat <<_ACEOF -\`configure' configures gdalcubes 0.6.4 to adapt to many kinds of systems. +\`configure' configures gdalcubes 0.6.9999 to adapt to many kinds of systems. Usage: $0 [OPTION]... [VAR=VALUE]... @@ -1341,7 +1341,7 @@ fi if test -n "$ac_init_help"; then case $ac_init_help in - short | recursive ) echo "Configuration of gdalcubes 0.6.4:";; + short | recursive ) echo "Configuration of gdalcubes 0.6.9999:";; esac cat <<\_ACEOF @@ -1439,7 +1439,7 @@ fi test -n "$ac_init_help" && exit $ac_status if $ac_init_version; then cat <<\_ACEOF -gdalcubes configure 0.6.4 +gdalcubes configure 0.6.9999 generated by GNU Autoconf 2.71 Copyright (C) 2021 Free Software Foundation, Inc. @@ -1595,7 +1595,7 @@ cat >config.log <<_ACEOF This file contains any messages produced by compilers while running configure, to aid debugging if configure makes a mistake. -It was created by gdalcubes $as_me 0.6.4, which was +It was created by gdalcubes $as_me 0.6.9999, which was generated by GNU Autoconf 2.71. Invocation command line was $ $0$ac_configure_args_raw @@ -4382,7 +4382,7 @@ cat >>$CONFIG_STATUS <<\_ACEOF || ac_write_fail=1 # report actual input values of CONFIG_FILES etc. instead of their # values after options handling. ac_log=" -This file was extended by gdalcubes $as_me 0.6.4, which was +This file was extended by gdalcubes $as_me 0.6.9999, which was generated by GNU Autoconf 2.71. Invocation command line was CONFIG_FILES = $CONFIG_FILES @@ -4437,7 +4437,7 @@ ac_cs_config_escaped=`printf "%s\n" "$ac_cs_config" | sed "s/^ //; s/'/'\\\\\\\\ cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1 ac_cs_config='$ac_cs_config_escaped' ac_cs_version="\\ -gdalcubes config.status 0.6.4 +gdalcubes config.status 0.6.9999 configured by $0, generated by GNU Autoconf 2.71, with options \\"\$ac_cs_config\\" @@ -4998,4 +4998,3 @@ printf "%s\n" "$as_me: WARNING: unrecognized options: $ac_unrecognized_opts" >&2 fi - diff --git a/inst/tinytest/test_window_space.R b/inst/tinytest/test_window_space.R new file mode 100644 index 0000000..0d59c4c --- /dev/null +++ b/inst/tinytest/test_window_space.R @@ -0,0 +1,49 @@ +library(gdalcubes) +v = cube_view(srs = "EPSG:4326", extent = list(left = 5, right = 15, bottom = 48, top = 58, + t0 = "2021-01-01", t1 = "2021-12-31"), dt = "P365D", + dx = 1, dy = 1) +v + +gdalcubes:::.raster_cube_dummy(v, 1, 1.0) |> + window_space("max(band1)","min(band1)", window = c(5,5)) |> + as_array() -> x + +expect_true(all(x == 1)) + + + +gdalcubes:::.raster_cube_dummy(v, 1, 1.0) |> + window_space("count(band1)", window = c(3,3)) |> + as_array() -> x +Xtrue = matrix(9, 10,10) +Xtrue[c(1,10),] = Xtrue[,c(1,10)] = 6 +Xtrue[1,1] = Xtrue[1,10] = Xtrue[10,1] = Xtrue[10,10] = 4 +expect_true(all(x[1,1,,] == Xtrue)) + + +K = matrix(1, 3,3) + +gdalcubes:::.raster_cube_dummy(v, 1, 1.0) |> + window_space(kernel = K, pad = 0) |> + as_array() -> x +expect_true(all(x[1,1,,] == Xtrue)) + +gdalcubes:::.raster_cube_dummy(v, 1, 1.0) |> + window_space(kernel = K, pad = "REFLECT") |> + as_array() -> x +expect_true(all(x == 9)) + + +gdalcubes:::.raster_cube_dummy(v, 1, 1.0) |> + window_space(kernel = K, pad = "REPLICATE") |> + as_array() -> x +expect_true(all(x == 9)) + + + + +gdalcubes:::.raster_cube_dummy(v, 1, 1.0, chunking = c(1,3,2)) |> + window_space(kernel = K, pad = "REPLICATE") |> + as_array() -> x +expect_true(all(x == 9)) + diff --git a/src/gdalcubes/src/utils.cpp b/src/gdalcubes/src/utils.cpp index 72b24d3..ba98f65 100644 --- a/src/gdalcubes/src/utils.cpp +++ b/src/gdalcubes/src/utils.cpp @@ -161,6 +161,16 @@ void utils::env::unset_all() { _vars.clear(); } + +std::string utils::env::get(std::string var_name, std::string default_value) { + std::string out = default_value; + char* r = getenv(var_name.c_str()); + if (r) { + out = r; + } + return out; +} + std::string utils::env::to_string() { std::string out; out = "{"; diff --git a/src/gdalcubes/src/utils.h b/src/gdalcubes/src/utils.h index bb79d30..7d5f69b 100644 --- a/src/gdalcubes/src/utils.h +++ b/src/gdalcubes/src/utils.h @@ -96,6 +96,7 @@ class utils { void set(std::map vars); void unset(std::set var_names); void unset_all(); + std::string get(std::string var_name, std::string default_value = ""); // Convert environment variable map to a JSON string std::string to_string(); diff --git a/src/gdalcubes/src/warp.cpp b/src/gdalcubes/src/warp.cpp index 3316940..07191d0 100644 --- a/src/gdalcubes/src/warp.cpp +++ b/src/gdalcubes/src/warp.cpp @@ -488,9 +488,17 @@ GDALDataset *gdalwarp_client::warp_complex(GDALDataset *in, std::string s_srs, s psWarpOptions->pfnProgress = GDALDummyProgress; // see https://github.com/OSGeo/gdal/blob/64cf9b4e889c93e34177237665fe842186d1f581/alg/gdaltransformer.cpp#L1274C5-L1275C67 - GDALGenImgProjTransformInfo *psInfo = static_cast(GDALCreateGenImgProjTransformer(in, s_srs.c_str(), out, NULL, TRUE, 0.0, 1 )); + CPLStringList trnsfrm_opts; // TODO: Do we need to check if t_srs and/or s_src are empty? + trnsfrm_opts.AddNameValue("SRC_SRS", s_srs.c_str()); + trnsfrm_opts.AddNameValue("GEOLOC_USE_TEMP_DATASETS", "NO"); + GDALGenImgProjTransformInfo *psInfo = static_cast(GDALCreateGenImgProjTransformer2(in, out, trnsfrm_opts.List())); + //GDALGenImgProjTransformInfo *psInfo = static_cast(GDALCreateGenImgProjTransformer(in, s_srs.c_str(), out, NULL, TRUE, 0.0, 1 )); + + // TODO DO we need to check if psInfo is NULL? + + // Overwrite reprojection transform if needed, to benefit from cache OGRSpatialReference srs_in; srs_in.SetFromUserInput(s_srs.c_str()); @@ -648,7 +656,7 @@ GDALDataset *gdalwarp_client::warp_complex(GDALDataset *in, std::string s_srs, s } static_cast(psWarpOptions->pTransformerArg)->pReprojectArg = tmparg; // safe release - GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg); + GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg); GDALDestroyWarpOptions(psWarpOptions); CPLFree(wkt_out);