Merge pull request #95 from appelmar/curvilinear

This pull request  mainly adds:
- support for images with spatial reference as geolocation arrays (including curvilinear grids)
- window_space operation

.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$
+^README.html$
+^inst/gdal$
+^inst/proj$

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)

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)

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, 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, 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) {
     .Call('_gdalcubes_gc_create_join_bands_cube', PACKAGE = 'gdalcubes', pin_list, cube_names)
 }

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,141 @@ 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 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.
+#' @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, 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")
+    }
+    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, as.character(pad_mode), as.double(pad_fill))
+  }
+  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, as.character(pad_mode), as.double(pad_fill))
+  }
+  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)
+}
+
+
+
+
 

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 <[email protected]>.
 #
@@ -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='[email protected]'
 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
 
 
-

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))
+