Title:	Earth Observation Data Cubes from Satellite Image Collections
Version:	0.7.1
Date:	2025-03-19
Description:	Processing collections of Earth observation images as on-demand multispectral, multitemporal raster data cubes. Users define cubes by spatiotemporal extent, resolution, and spatial reference system and let 'gdalcubes' automatically apply cropping, reprojection, and resampling using the 'Geospatial Data Abstraction Library' ('GDAL'). Implemented functions on data cubes include reduction over space and time, applying arithmetic expressions on pixel band values, moving window aggregates over time, filtering by space, time, bands, and predicates on pixel values, exporting data cubes as 'netCDF' or 'GeoTIFF' files, plotting, and extraction from spatial and or spatiotemporal features. All computational parts are implemented in C++, linking to the 'GDAL', 'netCDF', 'CURL', and 'SQLite' libraries. See Appel and Pebesma (2019) <doi:10.3390/data4030092> for further details.
Depends:	R (≥ 4.1.0)
Imports:	Rcpp, jsonlite, ncdf4
License:	MIT + file LICENSE
URL:	https://github.com/appelmar/gdalcubes
BugReports:	https://github.com/appelmar/gdalcubes/issues/
Encoding:	UTF-8
RoxygenNote:	7.3.2
LinkingTo:	Rcpp, BH
Suggests:	knitr, rmarkdown, stars, av, gifski, sf, tinytest, lubridate
VignetteBuilder:	knitr
Copyright:	file inst/COPYRIGHTS
NeedsCompilation:	yes
SystemRequirements:	GDAL (>= 2.0.1), PROJ (>= 4.8.0), netcdf, sqlite3
Packaged:	2025-03-19 07:42:02 UTC; marius
Author:	Marius Appel [aut, cre], Edzer Pebesma [ctb], Roger Bivand [ctb], Jeroen Ooms [ctb], Lewis Van Winkle [cph], Ole Christian Eidheim [cph], Howard Hinnant [cph], Adrian Colomitchi [cph], Florian Dang [cph], Paul Thompson [cph], Tomasz Kamiński [cph], Dropbox, Inc. [cph]
Maintainer:	Marius Appel <marius.appel@hs-bochum.de>
Repository:	CRAN
Date/Publication:	2025-03-20 01:10:02 UTC

gdalcubes: Earth Observation Data Cubes from Satellite Image Collections

Description

Processing collections of Earth observation images as on-demand multispectral, multitemporal raster data cubes. Users define cubes by spatiotemporal extent, resolution, and spatial reference system and let 'gdalcubes' automatically apply cropping, reprojection, and resampling using the 'Geospatial Data Abstraction Library' ('GDAL'). Implemented functions on data cubes include reduction over space and time, applying arithmetic expressions on pixel band values, moving window aggregates over time, filtering by space, time, bands, and predicates on pixel values, exporting data cubes as 'netCDF' or 'GeoTIFF' files, plotting, and extraction from spatial and or spatiotemporal features. All computational parts are implemented in C++, linking to the 'GDAL', 'netCDF', 'CURL', and 'SQLite' libraries. See Appel and Pebesma (2019) <doi:10.3390/data4030092> for further details.

Author(s)

Maintainer: Marius Appel marius.appel@hs-bochum.de (ORCID)

Other contributors:

• Edzer Pebesma (ORCID) [contributor]

• Roger Bivand [contributor]

• Jeroen Ooms (ORCID) [contributor]

• Lewis Van Winkle [copyright holder]

• Ole Christian Eidheim [copyright holder]

• Howard Hinnant [copyright holder]

• Adrian Colomitchi [copyright holder]

• Florian Dang [copyright holder]

• Paul Thompson [copyright holder]

• Tomasz Kamiński [copyright holder]

• Dropbox, Inc. [copyright holder]

See Also

Useful links:

• https://github.com/appelmar/gdalcubes

• Report bugs at https://github.com/appelmar/gdalcubes/issues/

Create a data cube proxy object copy

Description

Copy a data cube proxy object without copying any data

Usage

.copy_cube(cube)

Arguments

Details

This internal function copies the complete processing chain / graph of a data cube but does not copy any data It is used internally to avoid in-place modification for operations with potential side effects on source data cubes.

Value

copied data cube proxy object

Download and install an image collection format from a URL

Description

Download and install an image collection format from a URL

Usage

add_collection_format(url, name = NULL)

Arguments

Details

By default, the collection format name will be derived from the basename of the URL.

Examples

add_collection_format(
   "https://raw.githubusercontent.com/appelmar/gdalcubes_cpp/dev/formats/Sentinel1_IW_GRD.json")

Add images to an existing image collection

Description

This function adds provided files or GDAL dataset identifiers and to an existing image collection by extracting datetime, image identifiers, and band information according to the collection's format.

Usage

add_images(
  image_collection,
  files,
  unroll_archives = TRUE,
  out_file = "",
  quiet = FALSE
)

Arguments

Value

image collection proxy object, which can be used to create a data cube using raster_cube

Examples

L8_files <- list.files(system.file("L8NY18", package = "gdalcubes"),
                         ".TIF", recursive = TRUE, full.names = TRUE)
L8_col = create_image_collection(L8_files[1:12], "L8_L1TP") 
add_images(L8_col, L8_files[13:24])

Spatial aggregation of data cubes

Description

Create a proxy data cube, which applies an aggregation function to reduce the spatial resolution.

Usage

aggregate_space(cube, dx, dy, method = "mean", fact = NULL)

Arguments

Details

This function reduces the spatial resolution of a data cube by applying an aggregation function to smaller blocks of pixels.

The size of the cube may be expanded automatically in all directions if the original extent is not divisible by the new size of pixels.

Notice that if boundaries of the target cube do not align with the boundaries of the input cube (for example, if aggregating from 10m to 15m spatial resolution), pixels of the input cube will contribute to the output pixel that contains its center coordinate. If the center coordinate is exactly on a boundary, the input pixel will contribute to the right / bottom pixel of the output cube.

Note

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-12"),
              srs="EPSG:32618", dx = 500, dy=500, dt="P3M", aggregation = "median")
L8.cube = raster_cube(L8.col, v, mask=image_mask("BQA", bits=4, values=16))
L8.rgb = select_bands(L8.cube, c("B02", "B03", "B04"))
L8.5km = aggregate_space(L8.rgb, 5000,5000, "mean")
L8.5km


plot(L8.5km, rgb=3:1, zlim=c(5000,12000))

Aggregate data cube time series to lower temporal resolution

Description

Create a proxy data cube, which applies an aggregation function over pixel time series to lower temporal resolution.

Usage

aggregate_time(cube, dt, method = "mean", fact = NULL)

Arguments

Details

This function can be used to aggregate time series to lower resolution or to regularize a data cube with irregular (labeled) time axis. It is possible to change the unit of the temporal resolution (e.g. to create monthly composites from daily images). The size of the cube may be expanded automatically if the original temporal extent is not divisible by the new temporal size of pixels.

Note

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-12"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P3M", aggregation = "median")
L8.cube = raster_cube(L8.col, v, mask=image_mask("BQA", bits=4, values=16))
L8.rgb = select_bands(L8.cube, c("B02", "B03", "B04"))
L8.two_monthly = aggregate_time(L8.rgb, "P6M", "min")
L8.two_monthly


plot(L8.two_monthly, rgb=3:1, zlim=c(5000,12000))

Animate a data cube as an image time series

Description

This function can animate data cube time series as mp4 videos or animated GIFs. Depending on the desired output format, either the av or the gifski package is needed to create mp4 and GIF animations respectively.

Usage

animate(
  x,
  ...,
  fps = 1,
  loop = TRUE,
  width = 800,
  height = 800,
  save_as = tempfile(fileext = ".gif"),
  preview = interactive()
)

Arguments

Details

Animations can be created for single band data cubes or RGB plots of multi-band data cubes (by providing the argument rgb) only.

Value

character; path pointing to the the created file

See Also

plot.cube

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="P16D")

animate(select_bands(raster_cube(L8.col, v), c("B02", "B03", "B04")), rgb=3:1,
        zlim=c(0,20000), fps=1, loop=1)

animate(select_bands(raster_cube(L8.col, v), c("B05")), col=terrain.colors, key.pos=1)

Apply a function over (multi-band) pixels

Description

This generic function applies a function on pixels of a data cube, an R array, or other classes if implemented.

Usage

apply_pixel(x, ...)

Arguments

Value

return value and type depend on the class of x

See Also

apply_pixel.cube

apply_pixel.array

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

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.col = image_collection(file.path(tempdir(), "L8.db"))
apply_pixel(raster_cube(L8.col, v), "(B05-B04)/(B05+B04)", "NDVI") 



d <- c(4,16,128,128)
x <- array(rnorm(prod(d)), d)
y <- apply_pixel(x, function(v) {
  v[1] + v[2] + v[3] - v[4]
})

Apply a function over pixels in a four-dimensional (band, time, y, x) array

Description

Apply a function over pixels in a four-dimensional (band, time, y, x) array

Usage

## S3 method for class 'array'
apply_pixel(x, FUN, ...)

Arguments

Details

FUN is expected to produce a numeric vector (or scalar) where elements are interpreted as new bands in the result.

Note

This is a helper function that uses the same dimension ordering as gdalcubes. It can be used to simplify the application of R functions e.g. over time series in a data cube.

Examples

d <- c(4,16,32,32)
x <- array(rnorm(prod(d)), d)
y <- apply_pixel(x, function(v) {
  v[1] + v[2] + v[3] - v[4]
})
dim(y)

Apply arithmetic expressions over all pixels of a data cube

Description

Create a proxy data cube, which applies arithmetic expressions over all pixels of a data cube. Expressions may access band values by name.

Usage

## S3 method for class 'cube'
apply_pixel(
  x,
  expr,
  names = NULL,
  keep_bands = FALSE,
  ...,
  FUN,
  load_pkgs = FALSE,
  load_env = FALSE
)

Arguments

Details

The function can either apply simple arithmetic C expressions given as a character vector (expr argument), or apply a custom R reducer function if FUN is provided.

In the former case, gdalcubes uses the tinyexpr library to evaluate expressions in C / C++, you can look at the library documentation to see what kind of expressions you can execute. Pixel band values can be accessed by name. Predefined variables that can be used within the expression include integer pixel indexes (ix, iy, it), and pixel coordinates (left, right, top, bottom), t0, t1), where the last two values are provided seconds since epoch time.

FUN receives values of the bands from one pixel as a (named) vector and should return a numeric vector with identical length for all pixels. Elements of the result vectors will be interpreted as bands in the result data cube. Notice that by default, since FUN is executed in a separate R process, it cannot access any variables from outside and required packages must be loaded within FUN. To restore the current environment and automatically load packages, set load_env and/or load_pkgs to TRUE.

For more details and examples on how to write user-defined functions, please refer to the gdalcubes website at https://gdalcubes.github.io/source/concepts/udfs.html.

Value

a proxy data cube object

Note

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

# 1. Apply a C expression
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) 
L8.cube = select_bands(L8.cube, c("B04", "B05")) 
L8.ndvi = apply_pixel(L8.cube, "(B05-B04)/(B05+B04)", "NDVI") 
L8.ndvi


plot(L8.ndvi)


# 2. Apply a user defined R function
L8.ndvi.noisy = apply_pixel(L8.cube, names="NDVI_noisy", 
   FUN=function(x) {
       rnorm(1, 0, 0.1) + (x["B05"]-x["B04"])/(x["B05"]+x["B04"])
   })
L8.ndvi.noisy



plot(L8.ndvi.noisy)

 

Apply a function over (multi-band) pixel time series

Description

This generic function applies a function on pixel time series of a data cube, an R array, or other classes if implemented. The resulting object is expected to have the same spatial and temporal shape as the input, i.e., no reduction is performed.

Usage

apply_time(x, ...)

Arguments

Value

return value and type depend on the class of x

See Also

apply_time.cube

apply_time.array

Examples

# 1. input is data cube
# 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-06"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
L8.cube = raster_cube(L8.col, v) 
L8.cube = select_bands(L8.cube, c("B04", "B05")) 
L8.ndvi = apply_pixel(L8.cube, "(B05-B04)/(B05+B04)", "NDVI")

# Apply a user defined R function
apply_time(L8.ndvi, names="NDVI_residuals", 
   FUN=function(x) {
      y = x["NDVI",]
      if (sum(is.finite(y)) < 3) {
         return(rep(NA,ncol(x)))
      }
      t = 1:ncol(x)
      return(predict(lm(y ~ t)) -  x["NDVI",])})

# 2. input is array
d <- c(4,16,32,32)
x <- array(rnorm(prod(d)), d)
z <- apply_time(x, function(v) {
  y = matrix(NA, ncol=ncol(v), nrow=2)
  y[1,] = (v[1,] + v[2,]) / 2
  y[2,] = (v[3,] + v[4,]) / 2
  y
})
dim(z)

Apply a function over pixel time series in a four-dimensional (band, time, y, x) array

Description

Apply a function over pixel time series in a four-dimensional (band, time, y, x) array

Usage

## S3 method for class 'array'
apply_time(x, FUN, ...)

Arguments

Details

FUN is expected to produce a matrix (or vector if result has only one band) where rows are interpreted as new bands and columns represent time.

Note

This is a helper function that uses the same dimension ordering as gdalcubes. It can be used to simplify the application of R functions e.g. over time series in a data cube.

Examples

d <- c(4,16,32,32)
x <- array(rnorm(prod(d)), d)
z <- apply_time(x, function(v) {
  y = matrix(NA, ncol=ncol(v), nrow=2)
  y[1,] = (v[1,] + v[2,]) / 2
  y[2,] = (v[3,] + v[4,]) / 2
  y
})
dim(z)

Apply a user-defined R function over (multi-band) pixel time series

Description

Create a proxy data cube, which applies a user-defined R function over all pixel time series of a data cube. In contrast to reduce_time, the time dimension is not reduced, i.e., resulting time series must have identical length as the input data cube but may contain a different number of bands / variables. Example uses of this function may include time series decompositions, cumulative sums / products, smoothing, sophisticated NA filling, or similar.

Usage

## S3 method for class 'cube'
apply_time(
  x,
  names = NULL,
  keep_bands = FALSE,
  FUN,
  load_pkgs = FALSE,
  load_env = FALSE,
  ...
)

Arguments

Details

FUN receives a single (multi-band) pixel time series as a matrix with rows corresponding to bands and columns corresponding to time. In general, the function must return a matrix with the same number of columns. If the result contains only a single band, it may alternatively return a vector with length identical to the length of the input time series (number of columns of the input).

For more details and examples on how to write user-defined functions, please refer to the gdalcubes website at https://gdalcubes.github.io/source/concepts/udfs.html.

Value

a proxy data cube object

Note

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-06"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
L8.cube = raster_cube(L8.col, v) 
L8.cube = select_bands(L8.cube, c("B04", "B05")) 
L8.ndvi = apply_pixel(L8.cube, "(B05-B04)/(B05+B04)", "NDVI")

# Apply a user defined R function
L8.ndvi.resid = apply_time(L8.ndvi, names="NDVI_residuals", 
   FUN=function(x) {
      y = x["NDVI",]
      if (sum(is.finite(y)) < 3) {
         return(rep(NA,ncol(x)))
      }
      t = 1:ncol(x)
      return(predict(lm(y ~ t)) -  x["NDVI",])
   })
L8.ndvi.resid


plot(L8.ndvi.resid)

 

Convert a data cube to a data.frame

Description

Convert a data cube to a data.frame

Usage

## S3 method for class 'cube'
as.data.frame(x, ..., complete_only = FALSE)

Arguments

Value

A data.frame with bands / variables of the cube as columns and pixels as rows

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-05"),
              srs="EPSG:32618", nx = 100, ny=100, dt="P1M")
x = select_bands(raster_cube(L8.col, v), c("B04", "B05"))
df = as.data.frame(x, complete_only = TRUE)
head(df)

Convert a data cube to an in-memory R array

Description

Convert a data cube to an in-memory R array

Usage

as_array(x)

Arguments

Value

Four dimensional array with dimensions band, t, y, x

Note

Depending on the data cube size, this function may require substantial amounts of main memory, i.e. it makes sense for small data cubes only.

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-05"),
              srs="EPSG:32618", nx = 100, ny=100, dt="P1M")
x = as_array(select_bands(raster_cube(L8.col, v), c("B04", "B05")))
dim(x)
dimnames(x)

Query data cube properties

Description

gdalcubes internally uses a graph to serialize data cubes (including chained operations on cubes). This function derives a JSON representation, which can be used to save data cube objects without pixel data to disk.

Usage

as_json(obj, file = NULL)

Arguments

Value

If file is NULL, the function returns a JSON string representing a graph that can be used to recreate the same chain of gdalcubes operations even in a different R sessions.

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-04"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
cat(as_json(select_bands(raster_cube(L8.col, v), c("B04", "B05"))))

Query data cube properties

Description

Query data cube properties

Usage

bands(obj)

Arguments

Value

A data.frame with rows representing the bands and columns representing properties of a band (name, type, scale, offset, unit)

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")
bands(raster_cube(L8.col, v))

Apply an R function on chunks of a data cube

Description

Apply an R function on chunks of a data cube

Usage

chunk_apply(cube, f)

Arguments

Details

This function internally creates a gdalcubes stream data cube, which streams data of a chunk to a new R process. For reading data, the function typically calls x <- read_chunk_as_array() which then results in a 4 dimensional (band, time, y, x) array. Similarly write_chunk_from_array(x) will write a result array as a chunk in the resulting data cube. The chunk size of the input cube is important to control how the function will be exposed to the data cube. For example, if you want to apply an R function over complete pixel time series, you must define the chunk size argument in raster_cube to make sure that chunk contain the correct parts of the data.

Value

a proxy data cube object

Note

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-12"),
                          srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
L8.cube = raster_cube(L8.col, v)
L8.cube = select_bands(L8.cube, c("B04", "B05"))
f <- function() {
  x <- read_chunk_as_array()
  out <- reduce_time(x, function(x) {
    cor(x[1,], x[2,], use="na.or.complete", method = "kendall")
  }) 
  write_chunk_from_array(out)
}
L8.cor = chunk_apply(L8.cube, f)

List predefined image collection formats

Description

gdalcubes comes with some predefined collection formats e.g. to scan Sentinel 2 data. This function lists available formats including brief descriptions.

Usage

collection_formats(print = TRUE)

Arguments

Details

Image collection formats define how individual files / GDAL datasets relate to an image collection, i.e., which bands they contain, to which image they belong, and how to derive aquisition date/time. They are described as a set of regular expressions in a JSON file and used by gdalcubes to extract this information from the paths and/or filenames.

Value

data.frame with columns name and description where the former describes the unique identifier that can be used in create_image_collection and the latter gives a brief description of the format.

Examples

collection_formats()

Create an image collection from a set of GDAL datasets or files

Description

This function iterates over files or GDAL dataset identifiers and extracts datetime, image identifiers, and band information according to a given collection format.

Usage

create_image_collection(
  files,
  format = NULL,
  out_file = tempfile(fileext = ".sqlite"),
  date_time = NULL,
  band_names = NULL,
  use_subdatasets = FALSE,
  unroll_archives = TRUE,
  quiet = FALSE,
  one_band_per_file = NULL
)

Arguments

Details

An image collection is a simple index (a SQLite database) containing references to existing image files / GDAL dataset identifiers.

Collections can be created in two different ways: First, if a collection format is specified (argument format), date/time, bands, and metadata are automatically extracted from provided files. This is the most general approach but requires a collection format for the specific dataset.

Second, image collections can sometimes be created without collection format by manually specifying date/time of images (argument date_time) and names of bands (argument band_names). This is possible if either each image file contains all bands of the collection or only a single band. In the former case band_names simply contains the names of the bands or can be NULL to use default names. In the latter case (image files contain a single band only), the lengths of band_names and date_time must be identical. By default, the function assumes one band per file if length(band_names) == length(files). In the unlikely situation that this is not desired, it can be explicitly set using one_band_per_file.

Value

image collection proxy object, which can be used to create a data cube using raster_cube

Examples

# 1. create image collection using a collection format 
L8_files <- list.files(system.file("L8NY18", package = "gdalcubes"),
                       ".TIF", recursive = TRUE, full.names = TRUE)
x = create_image_collection(L8_files, "L8_L1TP")
x 


# 2. create image collection without format for a single band
L8_files_B4 <- list.files(system.file("L8NY18", package = "gdalcubes"),
                       "_B4.TIF", recursive = TRUE, full.names = TRUE)
d = as.Date(substr(basename(L8_files_B4), 18, 25), "%Y%m%d")
y = create_image_collection(L8_files_B4, date_time = d, band_names = "B4")
y


# 3. create image collection without format for all bands
d = as.Date(substr(basename(L8_files), 18, 25), "%Y%m%d")
fname = basename(tools::file_path_sans_ext(L8_files))
b = substr(fname, 27, nchar(fname))
z = create_image_collection(L8_files, date_time = d, band_names = b)
z

Crop data cube extent by space and/or time

Description

Create a proxy data cube, which crops a data cube by a spatial and/or temporal extent.

Usage

crop(cube, extent = NULL, iextent = NULL, snap = "near")

Arguments

Details

The new extent can be specified by spatial coordinates and datetime values (using the extent argument), or as zero-based integer indexes (using the iextent argument). In the former case, extent expects a list with numeric items left, right, top, bottom, t0, and t1, or a subset thereof. In the latter case, iextent is expected as a list with length-two integer vectors x, y, and t as items, defining the lower and upper cell indexes per dimension.

Notice that it is possible to crop only selected boundaries (e.g., only the right boundary) as missing boundaries in the extent or NA / NULL values in the iextent arguments are considered as "no change". It is, however, not possible to mix arguments extent and iextent.

If extent is given, the snap argument can be used to define what happens if the new boundary falls within a data cube cell.

Note

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-12"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P3M", aggregation = "median")
L8.cube = raster_cube(L8.col, v, mask=image_mask("BQA", bits=4, values=16))
L8.rgb = select_bands(L8.cube, c("B02", "B03", "B04"))

# crop by integer indexes
L8.cropped = crop(L8.rgb, iextent = list(x=c(0,400), y=c(0,400), t=c(1,1)))

# crop by spatiotemporal coordinates
L8.cropped = crop(L8.rgb, extent = list(left=388941.2, right=766552.4, 
   bottom=4345299, top=4744931, t0="2018-01", t1="2018-06"), snap = "in")
L8.cropped 

L8.cropped = crop(L8.rgb, extent = list(left=388941.2, right=766552.4, 
   bottom=4345299, top=4744931, t0="2018-01", t1="2018-06"), snap = "near")
L8.cropped 


plot(L8.cropped, rgb = 3:1, zlim=c(5000,10000))

Create or update a spatiotemporal data cube view

Description

Data cube views define the shape of a cube, i.e., the spatiotemporal extent, resolution, and spatial reference system (srs). They are used to access image collections as on-demand data cubes. The data cube will filter images based on the view's extent, read image data at the defined resolution, and warp / reproject images to the target srs automatically.

Usage

cube_view(
  view,
  extent,
  srs,
  nx,
  ny,
  nt,
  dx,
  dy,
  dt,
  aggregation,
  resampling,
  keep.asp = TRUE
)

Arguments

Details

The extent argument expects a simple list with elementes left, right, bottom, top, t0 (start date/time), t1 (end date/time) or an image collection object. In the latter case, the extent function is automatically called on the image collection object to get the full spatiotemporal extent of the collection. In the former case, datetimes are expressed as ISO8601 datetime strings.

The function can be used in two different ways. First, it can create data cube views from scratch by defining the extent, the spatial reference system, and for each dimension either the cell size (dx, dy, dt) or the total number of cells (nx, ny, nt). Second, the function can update an existing data cube view by overwriting specific fields. In this case, the extent or some elements of the extent may be missing.

In some cases, the extent of the view is automatically extended if the provided resolution would end within a pixel. For example, if the spatial extent covers an area of 1km x 1km and dx = dy = 300m, the extent would be enlarged to 1.2 km x 1.2km. The alignment will be reported to the user in a diagnostic message.

Value

A list with data cube view properties

Examples

 L8_files <- list.files(system.file("L8NY18", package = "gdalcubes"),
                        ".TIF", recursive = TRUE, full.names = TRUE)
 L8.col = create_image_collection(L8_files, "L8_L1TP")
 
 # 1. Create a new data cube view specification
 v = cube_view(extent=extent(L8.col,"EPSG:4326"), srs="EPSG:4326", dt="P1M", 
           nx=1000, ny=500, aggregation = "mean", resampling="bilinear")
 v

 # 2. overwrite parts of an existing data cube view
 vnew = cube_view(v, dt="P1M")

Query data cube properties

Description

Query data cube properties

Usage

## S3 method for class 'cube'
dim(x)

Arguments

Value

size of a data cube (number of cells) as integer vector in the order t, y, x

See Also

size

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")
dim(raster_cube(L8.col, v))

Query coordinate bounds for all dimensions of a data cube

Description

Dimension values give the coordinates bounds the spatial and temporal axes of a data cube.

Usage

dimension_bounds(obj, datetime_unit = NULL)

Arguments

Value

list with elements t,y,x, each a list with two elements, start and end

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")
dimension_bounds(raster_cube(L8.col, v))

Query coordinate values for all dimensions of a data cube

Description

Dimension values give the coordinates along the spatial and temporal axes of a data cube.

Usage

dimension_values(obj, datetime_unit = NULL)

Arguments

Value

list with elements t,y,x

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")
dimension_values(raster_cube(L8.col, v))

Query data cube properties

Description

Query data cube properties

Usage

dimensions(obj)

Arguments

Details

Elements of the returned list represent individual dimensions with properties such as dimension boundaries, names, and chunk size stored as inner lists

Value

Dimension information as a list

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")
dimensions(raster_cube(L8.col, v))

Derive the spatiotemporal extent of an image collection

Description

Derive the spatiotemporal extent of an image collection

Usage

extent(x, srs = "EPSG:4326")

Arguments

Value

a list with elements left, right, bottom, top, t0 (start date/time), and t1 (end date/time)

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"))
extent(L8.col,"EPSG:32618")
cube_view(extent=extent(L8.col,"EPSG:32618"),
          srs="EPSG:32618", nx = 497, ny=526, dt="P1M")

Extract values from a data cube by spatial or spatiotemporal features

Description

Extract pixel values of a data cube from a set of spatial or spatiotemporal features. Applications include the extraction of full time series at irregular points, extraction from spatiotemporal points, extraction of pixel values in polygons, and computing summary statistics over polygons.

Usage

extract_geom(
  cube,
  sf,
  datetime = NULL,
  time_column = NULL,
  FUN = NULL,
  merge = FALSE,
  drop_geom = FALSE,
  ...,
  reduce_time = FALSE
)

Arguments

Details

The geometry in sf can be of any simple feature type supported by GDAL, including POINTS, LINES, POLYGONS, MULTI*, and more. If no time information is provided in one of the arguments datetime or time_column, the full time series of pixels with regard to the features are returned.

Notice that feature identifiers in the FID column typically correspond to the row names / numbers of the provided sf object. This can be used to combine the output with the original geometries, e.g., using merge(). with gdalcubes > 0.6.4, this can be done automatically by setting merge=TRUE. In this case, the FID column is dropped from the result.

Pixels with missing values are automatically dropped from the result. It is hence not guaranteed that the result will contain rows for all input features.

Features are automatically reprojected if the coordinate reference system differs from the data cube.

Extracted values can be aggregated by features by providing a summary function. If reduce_time is FALSE (the default), the values are grouped by feature and time, i.e., the result will contain unique combinations of FID and time. To ignore time and produce a single value per feature, reduce_time can be set to TRUE.

Value

A data.frame with columns FID, time, and data cube bands / variables, see Details

Examples

# 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(srs="EPSG:32618", dy=1000, dx=1000, dt="P1M",
              aggregation = "median", resampling = "bilinear",
              extent=list(left=388941.2, right=766552.4,
                          bottom=4345299, top=4744931,
                          t0="2018-01-01", t1="2018-04-30"))
L8.cube = raster_cube(L8.col, v)
L8.cube = select_bands(L8.cube, c("B04", "B05"))
L8.ndvi = apply_pixel(L8.cube, "(B05-B04)/(B05+B04)", "NDVI")
L8.ndvi

if (gdalcubes_gdal_has_geos()) {
  if (requireNamespace("sf", quietly = TRUE)) {
  
    # create 50 random point locations
    x = runif(50, v$space$left, v$space$right)
    y = runif(50, v$space$bottom, v$space$top)
    t = sample(seq(as.Date("2018-01-01"),as.Date("2018-04-30"), by = 1),50, replace = TRUE)
    df = sf::st_as_sf(data.frame(x = x, y = y), coords = c("x", "y"), crs = v$space$srs)

    # 1. spatiotemporal points
    extract_geom(L8.ndvi, df, datetime = t)

    
    # 2. time series at spatial points
    extract_geom(L8.ndvi, df)
  
    # 3. summary statistics over polygons
    x = sf::st_read(system.file("nycd.gpkg", package = "gdalcubes"))
    zstats = extract_geom(L8.ndvi,x, FUN=median, reduce_time = TRUE, merge = TRUE)
    zstats
    plot(zstats["NDVI"])
    
  }
}

Fill NA data cube pixels by simple time series interpolation

Description

Create a proxy data cube, which fills NA pixels of a data cube by nearest neighbor or linear time series interpolation.

Usage

fill_time(cube, method = "near")

Arguments

Details

Please notice that completely empty (NA) time series will not be filled, i.e. the result cube might still contain NA values.

Value

a proxy data cube object

Note

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-12"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P3M", aggregation = "median")
L8.cube = raster_cube(L8.col, v, mask=image_mask("BQA", bits=4, values=16))
L8.rgb = select_bands(L8.cube, c("B02", "B03", "B04"))
L8.filled = fill_time(L8.rgb, "linear")
L8.filled


plot(L8.filled, rgb=3:1, zlim=c(5000,12000))

Filter data cube pixels by a polygon

Description

Create a proxy data cube, which filters pixels by a spatial (multi)polygon For all pixels whose center is within the polygon, the original

Usage

filter_geom(cube, geom, srs = NULL)

Arguments

Details

The resulting data cube will not be cropped but pixels outside of the polygon will be set to NAN.

If geom is provided as an sfc object with length > 1, geometries will be combined with sf::st_combine() before.

The geometry is automatically transformed to the data cube's spatial reference system if needed.

Value

a proxy data cube object

Note

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-06"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
L8.cube = raster_cube(L8.col, v) 
L8.cube = select_bands(L8.cube, c("B04", "B05")) 
L8.ndvi = apply_pixel(L8.cube, "(B05-B04)/(B05+B04)", "NDVI")
WKT = gsub(pattern='\\n',replacement="",x = 
  "Polygon ((-74.3541 40.9254, 
             -73.9813 41.2467, 
             -73.9997 41.4400, 
             -74.5362 41.1795, 
             -74.6286 40.9137, 
             -74.3541 40.9254))")
L8.ndvi.filtered = filter_geom(L8.ndvi, WKT, "EPSG:4326")
L8.ndvi.filtered

plot(L8.ndvi.filtered)

Filter data cube pixels by a user-defined predicate on band values

Description

Create a proxy data cube, which evaluates a predicate over all pixels of a data cube. For all pixels that fulfill the predicate, the original band values are returned. Other pixels are simply filled with NANs. The predicate may access band values by name.

Usage

filter_pixel(cube, pred)

Arguments

Details

gdalcubes uses and extends the tinyexpr library to evaluate expressions in C / C++, you can look at the library documentation to see what kind of expressions you can execute. Pixel band values can be accessed by name.

Value

a proxy data cube object

Note

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-06"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
L8.cube = raster_cube(L8.col, v) 
L8.cube = select_bands(L8.cube, c("B04", "B05")) 
L8.ndvi = apply_pixel(L8.cube, "(B05-B04)/(B05+B04)", "NDVI") 
L8.ndvi.filtered = filter_pixel(L8.ndvi, "NDVI > 0.5") 
L8.ndvi.filtered

plot(L8.ndvi.filtered)

Check if GDAL was built with GEOS

Description

Check if GDAL was built with GEOS

Usage

gdalcubes_gdal_has_geos()

Examples

gdalcubes_gdal_has_geos()

Get available GDAL drivers

Description

Get available GDAL drivers

Usage

gdalcubes_gdalformats()

Examples

gdalcubes_gdalformats()

Get the GDAL version used by gdalcubes

Description

Get the GDAL version used by gdalcubes

Usage

gdalcubes_gdalversion()

Examples

gdalcubes_gdalversion()

Set or read global options of the gdalcubes package

Description

Set global package options to change the default behavior of gdalcubes. These include how many parallel processes are used to process data cubes, how created netCDF files are compressed, and whether or not debug messages should be printed.

Usage

gdalcubes_options(
  ...,
  parallel,
  ncdf_compression_level,
  debug,
  cache,
  ncdf_write_bounds,
  use_overview_images,
  show_progress,
  default_chunksize,
  streaming_dir,
  log_file,
  threads
)

Arguments

Details

Data cubes can be processed in parallel where the number of chunks in a cube is distributed among parallel worker processes. The actual number of used workers can be lower if a data cube as less chunks. If parallel is TRUE, the number of available cores is used. Setting parallel = FALSE can be used to disable parallel processing. Notice that since version 0.6.0, separate processes are being used instead of parallel threads to avoid possible R session crashes due to some multithreading issues.

Caching has no effect on disk or memory consumption, it simply tries to reuse existing temporary files where possible. For example, changing only parameters to plot will void reprocessing the same data cube if cache is TRUE.

The streaming directory can be used to control the performance of user-defined functions, if disk IO is a bottleneck. Ideally, this can be set to a directory on a shared memory device.

Passing no arguments will return the current options as a list.

Examples

gdalcubes_options(parallel=4) # set the number 
gdalcubes_options() # print current options
gdalcubes_options(parallel=FALSE) # reset

Subsetting data cubes

Description

Subset data cube dimensions and bands / variables.

Usage

## S3 method for class 'cube'
x$name

## S3 method for class 'cube'
x[ib = TRUE, it = TRUE, iy = TRUE, ix = TRUE, ...]

Arguments

Details

The [] operator allows for flexible subsetting of data cubes by date, datetime, bounding box, spatial points, and band names. Depending on the arguments, it supports slicing (selecting one element of a dimension), cropping (selecting a subinterval of a dimension) and combinations thereof (e.g., selecting a spatial window and a temporal slice). Dimension subsets can be specified by integer indexes or coordinates / datetime values. Arguments are matched by type and order. For example, if the first argument is a length-two vector of type Date, the function will understand to subset the time dimension. Otherwise, arguments are treated in the order band, time, y, x.

Note

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-12"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P3M", aggregation = "median")
L8.cube = raster_cube(L8.col, v, mask=image_mask("BQA", bits=4, values=16))
L8.red = L8.cube$B04


plot(L8.red)


v = cube_view(extent=list(left=388941.2, right=766552.4,
                          bottom=4345299, top=4744931, t0="2018-01-01", t1="2018-12-31"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1D", aggregation = "median")
L8.cube = raster_cube(L8.col, v, mask=image_mask("BQA", bits=4, values=16))

L8.cube[c("B05","B04")] # select bands
L8.cube[as.Date(c("2018-01-10", "2018-01-20"))] # crop by time
L8.cube[as.Date("2018-01-10")] # slice by time
L8.cube["B05", "2018-01-10"] # select bands and slice by time
L8.cube["B05", c("2018-01-10","2018-01-17")] # select bands and crop by time
L8.cube[, c("2018-01-10","2018-01-17")] # crop by time

# crop by space (coordinates and integer indexes respectively)
L8.cube[list(left=388941.2 + 1e5, right=766552.4 - 1e5, bottom=4345299 + 1e5, top=4744931 - 1e5)]
L8.cube[,,c(1,100), c(1,100)] 

L8.cube[,c(1,2),,] # crop by time (integer indexes)

# subset by spatial point or bounding box
if (requireNamespace("sf", quietly = TRUE)) {
  s = sf::st_sfc(sf::st_point(c(500000, 4500000)), crs = "EPSG:32618")
  L8.cube[s]

  bbox =  sf::st_bbox(c(xmin = 388941.2 + 1e5, xmax = 766552.4 - 1e5,
                   ymax = 4744931 - 1e5, ymin = 4345299 + 1e5), crs = sf::st_crs(32618))
  L8.cube[bbox]
}

Set GDAL config options

Description

Set GDAL config options

Usage

gdalcubes_set_gdal_config(key, value)

Arguments

Details

Details and a list of possible options can be found at https://gdal.org/en/stable/user/configoptions.html.

Examples

gdalcubes_set_gdal_config("GDAL_NUM_THREADS", "ALL_CPUS")

Load an existing image collection from a file

Description

This function will load an image collection from an SQLite file. Image collection files index and reference existing imagery. To create a collection from files on disk, use create_image_collection.

Usage

image_collection(path)

Arguments

Value

an image collection proxy object, which can be used to create a data cube using raster_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"))
L8.col

Create a mask for images in a raster data cube

Description

Create an image mask based on a band and provided values to filter pixels of images read by raster_cube

Usage

image_mask(
  band,
  min = NULL,
  max = NULL,
  values = NULL,
  bits = NULL,
  invert = FALSE
)

Arguments

Details

Values of the selected mask band can be based on a range (by passing min and max) or on a set of values (by passing values). By default pixels with mask values contained in the range or in the values are masked out, i.e. set to NA. Setting invert = TRUE will invert the masking behavior. Passing values will override min and max.

Note

Notice that masks are applied per image while reading images as a raster cube. They can be useful to eliminate e.g. cloudy pixels before applying the temporal aggregation to merge multiple values for the same data cube pixel.

Examples

image_mask("SCL", values = c(3,8,9)) # Sentinel 2 L2A: mask cloud and cloud shadows
image_mask("BQA", bits=4, values=16) # Landsat 8: mask clouds
image_mask("B10", min = 8000, max=65000) 

Join bands of two identically shaped data cubes

Description

Create a proxy data cube, which joins the bands of two identically shaped data cubes. The resulting cube will have bands from both input cubes.

Usage

join_bands(cube_list, cube_names = NULL)

Arguments

Details

The number of provided cube_names must match the number of provided input cubes. If no cube_names are provided, bands of the output cube will adopt original names from the input cubes (without any prefix). If any two of the input bands have identical names, prefixes default prefixes ("X1", "X2", ...) will be used.

Value

proxy data cube object

Note

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-05"),
                          srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
L8.cube = raster_cube(L8.col, v)
L8.cube.b04 = select_bands(raster_cube(L8.col, v), c("B04"))
L8.cube.b05 = select_bands(raster_cube(L8.col, v), c("B05"))
join_bands(list(L8.cube.b04,L8.cube.b05))

plot(join_bands(list(L8.cube.b04,L8.cube.b05)))

Read a data cube from a json description file

Description

Read a data cube from a json description file

Usage

json_cube(json, path = NULL)

Arguments

Details

Data cubes can be stored as JSON description files. These files do not store any data but the recipe how a data cube is constructed, i.e., the chain (or graph) of processes involved.

Since data cube objects (as returned from raster_cube) cannot be saved with normal R methods, the combination of as_json and json_cube provides a cheap way to save virtual data cube objects across several R sessions, as in the examples.

Value

data cube proxy object

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-12"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
cube = raster_cube(L8.col, v) 

# save
fname = tempfile()
as_json(cube, fname)

# load
json_cube(path = fname)  

Query data cube properties

Description

Query data cube properties

Usage

memsize(obj, unit = "MiB")

Arguments

Value

Total data size of data cube values expressed in the given unit

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")
memsize(raster_cube(L8.col, v))

Query data cube properties

Description

Query data cube properties

Usage

## S3 method for class 'cube'
names(x)

Arguments

Value

Band names as character vector

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")
names(raster_cube(L8.col, v))

Query data cube properties

Description

Query data cube properties

Usage

nbands(obj)

Arguments

Value

Number of bands

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")
nbands(raster_cube(L8.col, v))

Read a data cube from an existing netCDF file

Description

Create a proxy data cube from a netCDF file that has been created using write_ncdf. This function does not read cubes from arbitrary netCDF files and can be used e.g., to load intermediate results and/or plotting existing netCDF cubes on disk without doing the data cube creation from image collections.

Usage

ncdf_cube(path, chunking = NULL, auto_unpack = TRUE)

Arguments

Value

a proxy data cube object

Note

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


ncfile = write_ncdf(select_bands(raster_cube(L8.col, v), c("B02", "B03", "B04")))
ncdf_cube(ncfile)
                           

Query data cube properties

Description

Query data cube properties

Usage

nt(obj)

Arguments

Value

Number of pixels in the time dimension

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")
nt(raster_cube(L8.col, v))

Query data cube properties

Description

Query data cube properties

Usage

nx(obj)

Arguments

Value

Number of pixels in the x dimension

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")
nx(raster_cube(L8.col, v))

Query data cube properties

Description

Query data cube properties

Usage

ny(obj)

Arguments

Value

Number of pixels in the y dimension

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")
ny(raster_cube(L8.col, v))

Helper function to define packed data exports by min / max values

Description

This function can be used to define packed exports in write_ncdf and write_tif. It will generate scale and offset values with maximum precision (unless simplify=TRUE).

Usage

pack_minmax(type = "int16", min, max, simplify = FALSE)

Arguments

Details

Nodata values will be mapped to the lowest value of the target data type.

Arguments min and max must have length 1 or length equal to the number of bands of the data cube to be exported. In the former case, the same values are used for all bands of the exported target cube, whereas the latter case allows to use different ranges for different bands.

Note

Using simplify=TRUE will round scale values to the next smaller power of 10.

Examples

ndvi_packing = pack_minmax(type="int16", min=-1, max=1)
ndvi_packing

Plot a gdalcubes data cube

Description

Plot a gdalcubes data cube

Usage

## S3 method for class 'cube'
plot(
  x,
  y,
  ...,
  nbreaks = 11,
  breaks = NULL,
  col = grey(1:(nbreaks - 1)/nbreaks),
  key.pos = NULL,
  bands = NULL,
  t = NULL,
  rgb = NULL,
  zlim = NULL,
  gamma = 1,
  periods.in.title = TRUE,
  join.timeseries = FALSE,
  axes = TRUE,
  ncol = NULL,
  nrow = NULL,
  downsample = TRUE,
  na.color = "#AAAAAA"
)

Arguments

Details

The style of the plot depends on provided parameters and on the shape of the cube, i.e., whether it is a pure time series and whether it contains multiple bands or not. Multi-band, multi-temporal images will be arranged with layout() such that bands are represented by columns and time is represented by rows. Time series plots can be combined to a single plot by setting join.timeseries = TRUE. The layout can be controlled with ncol and nrow, which define the number of rows and columns in the plot layout. Typically, only one of ncol and nrow is provided. For multi-band, multi-temporal plots, the actual number of rows or columns can be less if the input cube has less bands or time slices.

The downsample argument is used to speed-up plotting if the cube has much more pixels than the graphics device. If set to a scalar integer value > 1, the value is used to skip pixels in the spatial dimensions. For example, setting downsample = 4 means that every fourth pixel is used in the spatial dimensions. If TRUE (the default) downsample is derived automatically based on the sizes of the cube and the graphics device. If 1 or FALSE, no additional downsampling is performed. Notice that downsampling is only used for plotting. The size of the data cube (and hence the computation time to process the data cube) is not modified.

Note

If caching is enabled for the package (see gdalcubes_options), repeated calls of plot for the same data cube will not reevaluate the cube. Instead, the temporary result file will be reused, if possible.

Some parts of the function have been copied from the stars package (c) Edzer Pebesma

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")
              
plot(select_bands(raster_cube(L8.col, v), c("B02", "B03", "B04")), rgb=3:1)
              
L8.cube = select_bands(raster_cube(L8.col, v), c("B04", "B05")) 
L8.ndvi = apply_pixel(L8.cube, "(B05-B04)/(B05+B04)", "NDVI") 
plot(reduce_time(L8.ndvi, "median(NDVI)"), key.pos=1, zlim=c(0,1))

Model prediction

Description

Apply a trained model on all pixels of a data cube.

Usage

## S3 method for class 'cube'
predict(object, model, ..., output_names = c("pred"), keep_bands = FALSE)

Arguments

Details

The model-specific predict method will be automatically chosen based on the class of the provided model. It aims at supporting models from the packages tidymodels, caret, and simple models as from lm or glm.

For multiple output variables or output in form of lists or data.frames, output_names must be provided and match names of the columns / items of the result object returned from the underlying predict method. For example, predictions using tidymodels return a tibble (data.frame) with columns like .pred_class (classification case). This must be explicitly provided as output_names. Similarly, predict.lm and the like return lists if the standard error is requested by the user and output_names hence should be set to c("fit","se.fit").

For more complex cases or when predict expects something else than a data.frame, this function may not work at all.

Note

This function returns a proxy object, i.e., it will not immediately start any computations.

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

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="P3M")

L8.col = image_collection(file.path(tempdir(), "L8.db"))

x = sf::st_read(system.file("ny_samples.gpkg", package = "gdalcubes"))

raster_cube(L8.col, v) |>
  select_bands(c("B02","B03","B04","B05")) |>
  extract_geom(x) -> train

x$FID = rownames(x)  
train = merge(train, x, by = "FID")
train$iswater = as.factor(train$class == "water")

log_model <- glm(iswater ~ B02 + B03 + B04 + B05, data = train, family = "binomial")

raster_cube(L8.col, v) |>
  select_bands(c("B02","B03","B04","B05")) |>
  predict(model=log_model, type="response") |>
  plot(key.pos=1)

Print data cube information

Description

Prints information about the dimensions and bands of a data cube.

Usage

## S3 method for class 'cube'
print(x, ...)

Arguments

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-12"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
print(raster_cube(L8.col, v))

Print data cube view information

Description

Prints information about a data cube view, including its dimensions, spatial reference, aggregation method, and resampling method.

Usage

## S3 method for class 'cube_view'
print(x, ...)

Arguments

Examples

v = cube_view(extent=list(left=388941.2, right=766552.4, 
              bottom=4345299, top=4744931, t0="2018-01", t1="2018-12"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
print(v)

Print image collection information

Description

Prints information about images in an image collection.

Usage

## S3 method for class 'image_collection'
print(x, ..., n = 6)

Arguments

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"))
print(L8.col)

Query data cube properties

Description

Query data cube properties

Usage

proj4(obj)

Arguments

Value

The spatial reference system expressed as proj4 string

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")
proj4(raster_cube(L8.col, v))

Create a data cube from an image collection

Description

Create a proxy data cube, which loads data from a given image collection according to a data cube view

Usage

raster_cube(
  image_collection,
  view,
  mask = NULL,
  chunking = .pkgenv$default_chunksize,
  incomplete_ok = TRUE
)

Arguments

Details

The following steps will be performed when the data cube is requested to read data of a chunk:

1. Find images from the input collection that intersect with the spatiotemporal extent of the chunk 2. For all resulting images, apply gdalwarp to reproject, resize, and resample to an in-memory GDAL dataset 3. Read the resulting data to the chunk buffer and optionally apply a mask on the result 4. Update pixel-wise aggregator (as defined in the data cube view) to combine values of multiple images within the same data cube pixels

If chunking is provided as a function, it must accept exactly three arguments for the total size of the cube in t, y, and x axes (in this order).

Value

A proxy data cube object

Note

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-12"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
raster_cube(L8.col, v)
 
 # using a mask on the Landsat quality bit band to filter out clouds
 raster_cube(L8.col, v, mask=image_mask("BQA", bits=4, values=16))
 

Read chunk data of a data cube from stdin or a file

Description

This function can be used within function passed to chunk_apply in order to read a data cube chunk as a four-dimensional R array. It works only for R processes, which have been started from the gdalcubes C++ library. The resulting array has dimensions band, time, y, x (in this order).

Usage

read_chunk_as_array(with.dimnames = TRUE)

Arguments

Value

four-dimensional array

Note

Call this function ONLY from a function passed to chunk_apply.

This function only works in R sessions started from gdalcubes streaming.

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-12"),
                          srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
L8.cube = raster_cube(L8.col, v)
L8.cube = select_bands(L8.cube, c("B04", "B05"))
f <- function() {
  x <- read_chunk_as_array()
  out <- reduce_time(x, function(x) {
    cor(x[1,], x[2,], use="na.or.complete", method = "kendall")
  }) 
  write_chunk_from_array(out)
}
L8.cor = chunk_apply(L8.cube, f)
plot(L8.cor, zlim=c(0,1), key.pos=1)

Reduce multidimensional data over space

Description

This generic function applies a reducer function over a data cube, an R array, or other classes if implemented.

Usage

reduce_space(x, ...)

Arguments

Value

return value and type depend on the class of x

See Also

reduce_space.cube

reduce_space.array

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-12"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
reduce_space(raster_cube(L8.col, v) , "median(B02)")  


d <- c(4,16,32,32)
x <- array(rnorm(prod(d)), d)
y <- reduce_space(x, function(v) {
  apply(v, 1, mean)
})
 

Apply a function over space and bands in a four-dimensional (band, time, y, x) array and reduce spatial dimensions

Description

Apply a function over space and bands in a four-dimensional (band, time, y, x) array and reduce spatial dimensions

Usage

## S3 method for class 'array'
reduce_space(x, FUN, ...)

Arguments

Details

FUN is expected to produce a numeric vector (or scalar) where elements are interpreted as new bands in the result.

Note

This is a helper function that uses the same dimension ordering as gdalcubes streaming. It can be used to simplify the application of R functions e.g. over spatial slices in a data cube.

Examples

d <- c(4,16,32,32)
x <- array(rnorm(prod(d)), d)
# reduce individual bands over spatial slices 
y <- reduce_space(x, function(v) {
  apply(v, 1, mean)
})
dim(y)

Reduce a data cube over spatial (x,y or lat,lon) dimensions

Description

Create a proxy data cube, which applies one or more reducer functions to selected bands over spatial slices of a data cube

Usage

## S3 method for class 'cube'
reduce_space(
  x,
  expr,
  ...,
  FUN,
  names = NULL,
  load_pkgs = FALSE,
  load_env = FALSE
)

Arguments

Details

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".

For more details and examples on how to write user-defined functions, please refer to the gdalcubes website at https://gdalcubes.github.io/source/concepts/udfs.html.

Value

proxy data cube object

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-12"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
L8.cube = raster_cube(L8.col, v) 
L8.b02 = select_bands(L8.cube, c("B02"))
L8.b02.median = reduce_space(L8.b02, "median(B02)")  
L8.b02.median

plot(L8.b02.median)

Reduce multidimensional data over time

Description

This generic function applies a reducer function over a data cube, an R array, or other classes if implemented.

Usage

reduce_time(x, ...)

Arguments

Value

return value and type depend on the class of x

See Also

reduce_time.cube

reduce_time.array

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-06"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
reduce_time(raster_cube(L8.col, v) , "median(B02)", "median(B03)", "median(B04)")  
 
d <- c(4,16,32,32)
x <- array(rnorm(prod(d)), d)
y <- reduce_time(x, function(v) {
  apply(v, 1, mean)
})
 

Apply a function over time and bands in a four-dimensional (band, time, y, x) array and reduce time dimension

Description

Apply a function over time and bands in a four-dimensional (band, time, y, x) array and reduce time dimension

Usage

## S3 method for class 'array'
reduce_time(x, FUN, ...)

Arguments

Details

FUN is expected to produce a numeric vector (or scalar) where elements are interpreted as new bands in the result.

Note

This is a helper function that uses the same dimension ordering as gdalcubes streaming. It can be used to simplify the application of R functions e.g. over time series in a data cube.

Examples

d <- c(4,16,32,32)
x <- array(rnorm(prod(d)), d)
# reduce individual bands over pixel time series
y <- reduce_time(x, function(v) {
  apply(v, 1, mean)
})
dim(y)

Reduce a data cube over the time dimension

Description

Create a proxy data cube, which applies one or more reducer functions to selected bands over pixel time series of a data cube

Usage

## S3 method for class 'cube'
reduce_time(
  x,
  expr,
  ...,
  FUN,
  names = NULL,
  load_pkgs = FALSE,
  load_env = FALSE
)

Arguments

Details

The function can either apply a built-in reducer if expr is given, or apply a custom R reducer function if FUN is provided.

In the former case, 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", "var", "sd", "which_min", "which_max", "Q1" (1st quartile), and "Q3" (3rd quartile).

User-defined R reducer functions receive a two-dimensional array as input where rows correspond to the band and columns represent the time dimension. For example, one row is the time series of a specific band. FUN should always return a numeric vector with the same number of elements, which will be interpreted as bands in the result cube. Notice that it is recommended to specify the names of the output bands as a character vector. If names are missing, the number and names of output bands is tried to be derived automatically, which may fail in some cases.

For more details and examples on how to write user-defined functions, please refer to the gdalcubes website at https://gdalcubes.github.io/source/concepts/udfs.html.

Value

proxy data cube object

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-06"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
L8.cube = raster_cube(L8.col, v) 
L8.rgb = select_bands(L8.cube, c("B02", "B03", "B04"))
L8.rgb.median = reduce_time(L8.rgb, "median(B02)", "median(B03)", "median(B04)")  
L8.rgb.median


plot(L8.rgb.median, rgb=3:1)


# user defined reducer calculating interquartile ranges
L8.rgb.iqr = reduce_time(L8.rgb, names=c("iqr_R", "iqr_G","iqr_B"), FUN = function(x) {
    c(diff(quantile(x["B04",],c(0.25,0.75), na.rm=TRUE)),
      diff(quantile(x["B03",],c(0.25,0.75), na.rm=TRUE)),
      diff(quantile(x["B02",],c(0.25,0.75), na.rm=TRUE)))
})
L8.rgb.iqr

plot(L8.rgb.iqr, key.pos=1)

Rename bands of a data cube

Description

Create a proxy data cube, which renames specific bands of a data cube.

Usage

rename_bands(cube, ...)

Arguments

Details

The result data cube always contains the same number of bands. No subsetting is done if only names for some of the bands are provided. In this case, only provided bands are renamed whereas other bands keep their original name. Variable arguments must be named by the old band name and the new names must be provided as simple character values (see example).

Value

proxy data cube object

Note

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-07"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
L8.cube = raster_cube(L8.col, v) 
L8.rgb = select_bands(L8.cube, c("B02", "B03", "B04"))
L8.rgb
L8.rgb = rename_bands(L8.cube, B02 = "blue", B03 = "green", B04 = "red")
L8.rgb

Select bands of a data cube

Description

Create a proxy data cube, which selects specific bands of a data cube. The resulting cube will drop any other bands.

Usage

select_bands(cube, bands)

Arguments

Value

proxy data cube object

Note

This function returns a proxy object, i.e., it will not start any computations besides deriving the shape of the result.

For performance reasons, select_bands should always be called directly on a cube created with raster_cube and drop all unneded bands. This allows to reduce RasterIO and warp operations in GDAL.

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-07"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
L8.cube = raster_cube(L8.col, v) 
L8.rgb = select_bands(L8.cube, c("B02", "B03", "B04"))
L8.rgb

plot(L8.rgb, rgb=3:1)

Select time slices of a data cube

Description

Create a proxy data cube, which selects specific time slices of a data cube. The time dimension of the resulting cube will be irregular / labeled.

Usage

select_time(cube, t)

Arguments

Value

proxy data cube object

Note

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-07"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
L8.cube = raster_cube(L8.col, v) 
L8.rgb = select_bands(L8.cube, c("B02", "B03", "B04"))
L8.rgb = select_time(L8.rgb, c("2018-04", "2018-07"))
L8.rgb

plot(L8.rgb, rgb=3:1)

Query data cube properties

Description

Query data cube properties

Usage

size(obj)

Arguments

Value

size of a data cube (number of cells) as integer vector in the order t, y, x

See Also

dim.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-04", t1="2018-06"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
size(raster_cube(L8.col, v))

Extract a single time series (spatial slice) from a data cube

Description

Create a proxy data cube, which extracts a time series from a data cube defined by spatial coordinates or integer x and y indexes.

Usage

slice_space(cube, loc = NULL, i = NULL)

Arguments

Details

Either loc or i must be non-NULL. If both arguments are provided, integer indexes i are ignored.

Note

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-12"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P3M", aggregation = "median")
L8.cube = raster_cube(L8.col, v, mask=image_mask("BQA", bits=4, values=16))
L8.rgb = select_bands(L8.cube, c("B02", "B03", "B04"))
L8.ts = slice_space(L8.rgb, loc = c(5e05, 4400000))
L8.ts


plot(L8.ts, join.timeseries = TRUE)

Extract a single time slice from a data cube

Description

Create a proxy data cube, which extracts a time slice from a data cube defined by label (datetime string) or integer index.

Usage

slice_time(cube, datetime = NULL, it = NULL)

Arguments

Details

Either datetime or it must be non-NULL. If both arguments are provided, the integer index it is ignored.

Note

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-12"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M", aggregation = "median")
L8.cube = raster_cube(L8.col, v, mask=image_mask("BQA", bits=4, values=16))
L8.rgb = select_bands(L8.cube, c("B02", "B03", "B04"))
L8.slice = slice_time(L8.rgb, "2018-03")
L8.slice


plot(L8.slice, rgb=3:1, zlim=c(5000,12000))

Query data cube properties

Description

Query data cube properties

Usage

srs(obj)

Arguments

Value

The spatial reference system expressed as a string readable by GDAL

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")
srs(raster_cube(L8.col, v))

Coerce gdalcubes object into a stars object

Description

The function materializes a data cube as a temporary netCDF file and loads the file with the stars package.

Usage

st_as_stars.cube(.x, ...)

Arguments

Value

stars object

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-04"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
if(require("stars"))
  st_as_stars(select_bands(raster_cube(L8.col, v), c("B04", "B05")))

Create an image collection from a STAC feature collection

Description

This function creates an image collection from a STAC API collection response. It does not need to read any image data. Additionally, bands can be filtered and asset links can be transformed to make them readable for GDAL.

Usage

stac_image_collection(
  s,
  out_file = tempfile(fileext = ".db"),
  asset_names = NULL,
  asset_regex = NULL,
  url_fun = .default_url_fun,
  property_filter = NULL,
  skip_image_metadata = FALSE,
  srs = NULL,
  srs_overwrite = FALSE,
  duration = c("center", "start")
)

Arguments

Details

The property_filter argument can be used to filter images by metadata such as cloud coverage. The functions receives all properties of a STAC item (image) as input list and is expected to produce a single logical value, where an image will be ignored if the function returns FALSE.

Some STAC API endpoints may return items with duplicte IDs (image names), pointing to identical URLs. Such items are only added once during creation of the image collection.

Note

Currently, bbox results are expected to be WGS84 coordinates, even if bbox-crs is given in the STAC response.

Create a data cube from a set of images with the same spatial extent and spatial reference system

Description

Create a spatiotemporal data cube directly from images with identical spatial extent and spatial reference system, similar to a raster stack with an additional dimension supporting both, time and multiple bands / variables.

Usage

stack_cube(
  x,
  datetime_values,
  bands = NULL,
  band_names = NULL,
  chunking = c(1, 256, 256),
  dx = NULL,
  dy = NULL,
  incomplete_ok = TRUE
)

Arguments

Details

This function creates a four-dimensional (space, time, bands / variables) raster data cube from a set of provided files without the need to create an image collection before. This is possible if all images have the same spatial extent and spatial reference system and can be used for two different file organizations:

1. If all image files share the same bands / variables, the bands argument can be ignored (default NULL) can names of the bands can be specified using the band_names argument.

2. If image files represent different band / variable (e.g. individual files for red, green, and blue channels), the bands argument must be used to define the corresponding band / variable. Notice that in this case all files are expected to represent exactly one variable / band at one point in datetime. It is not possible to combine files with different numbers of variables / bands. If image files for different bands have different pixel sizes, the smallest size is used by default.

Notice that to avoid opening all image files in advance,no automatic check whether all images share the spatial extent and spatial reference system is performed.

Value

A proxy data cube object

Note

This function returns a proxy object, i.e., it will not start any computations besides deriving the shape of the result.

Examples

# toy example, repeating the same image as a daily time series
L8_file_nir <- 
system.file("L8NY18/LC08_L1TP_014032_20181122_20181129_01_T1/LC08_L1TP_014032_20181122_B5.TIF",
            package = "gdalcubes")
files = rep(L8_file_nir, 10)
datetime = as.Date("2018-11-22") + 1:10
stack_cube(files, datetime, band_names = "B05") 

# using a second band from different files
L8_file_red <- 
system.file("L8NY18/LC08_L1TP_014032_20181122_20181129_01_T1/LC08_L1TP_014032_20181122_B4.TIF",
            package = "gdalcubes")
files = rep(c(L8_file_nir, L8_file_red), each = 10)
datetime = rep(as.Date("2018-11-22") + 1:10, 2)
bands = rep(c("B5","B4"), each = 10)
stack_cube(files, datetime, bands = bands)         

Apply a moving window (focal) operation or a convolution kernel over spatial dimensions of a data cube.

Description

Create a proxy data cube, which applies a convolution kernel or aggregation functions over two-dimensional moving windows sliding over spatial slices of a data cube. The function can either execute one or more predefined aggregation functions or apply a custom convolution kernel. Among others, use cases include image processing (edge detection, noise reduction, etc.) and enriching pixel values with local neighborhood properties (e.g. to use as predictor variables in ML models).

Usage

window_space(x, expr, ..., kernel, window, keep_bands = FALSE, pad = NA)

Arguments

Details

The function either applies a kernel convolution (if the 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 aggregation functions over moving windows, 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, e.g, "mean(band1)". Possible reducers include "min", "max", "sum", "prod", "count", "mean", "median", "var", and "sd".

Padding methods "REPLICATE", "REFLECT", "REFLEX_PIXEL" are defined according to https://openeo.org/documentation/1.0/processes.html#apply_kernel.

Value

proxy data cube object

Note

Implemented reducers will ignore any NAN values (as na.rm = TRUE does).

Calling this function consecutively many times may result in long computation times depending on chunk and window sizes due to the need to read adjacent data cube chunks.

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


plot(L8.cube.mean5x5, key.pos=1)


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

plot(L8.cube.med_sd, key.pos=1)

Apply a moving window operation over the time dimension of a data cube

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 function can either apply a built-in aggregation function or apply a custom one-dimensional convolution kernel.

Usage

window_time(x, expr, ..., kernel, window)

Arguments

Details

The function either applies a kernel convolution (if the 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 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 include "min", "max", "sum", "prod", "count", "mean", and "median".

Value

proxy data cube object

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

Write chunk data of a cube to stdout or a file

Description

This function can be used within function passed to chunk_apply in order to pass four-dimensional R arrays as a data cube chunk to the gdalcubes C++ library. It works only for R processes, which have been started from the gdalcubes C++ library. The input array must have dimensions band, time, y, x (in this order).

Usage

write_chunk_from_array(v)

Arguments

Note

Call this function ONLY from a function passed to chunk_apply.

This function only works in R sessions started from gdalcubes streaming.

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-12"),
                          srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
L8.cube = raster_cube(L8.col, v)
L8.cube = select_bands(L8.cube, c("B04", "B05"))
f <- function() {
  x <- read_chunk_as_array()
  out <- reduce_time(x, function(x) {
    cor(x[1,], x[2,], use="na.or.complete", method = "kendall")
  }) 
  write_chunk_from_array(out)
}
L8.cor = chunk_apply(L8.cube, f)
plot(L8.cor, zlim=c(0,1), key.pos=1)

Export a data cube as netCDF file(s)

Description

This function will read chunks of a data cube and write them to a single (the default) or multitple (if chunked = TRUE) netCDF file(s). The resulting file(s) uses the enhanced netCDF-4 format, supporting chunking and compression.

Usage

write_ncdf(
  x,
  fname = tempfile(pattern = "gdalcubes", fileext = ".nc"),
  overwrite = FALSE,
  write_json_descr = FALSE,
  with_VRT = FALSE,
  pack = NULL,
  chunked = FALSE
)

Arguments

Details

The resulting netCDF file(s) contain three dimensions (t, y, x) and bands as variables.

If write_json_descr is TRUE, the function will write an addition file with the same name as the NetCDF file but ".json" suffix. This file includes a serialized description of the input data cube, including all chained data cube operations.

To reduce the size of created files, values can be packed by applying a scale factor and an offset value and using a smaller integer data type for storage (only supported if chunked = TRUE). The pack argument can be either NULL (the default), or a list with elements type, scale, offset, and nodata. type can be any of "uint8", "uint16" , "uint32", "int16", or "int32". scale, offset, and nodata must be numeric vectors with length one or length equal to the number of data cube bands (to use different values for different bands). The helper function pack_minmax can be used to derive offset and scale values with maximum precision from minimum and maximum data values on original scale.

If chunked = TRUE, names of the produced files will start with name (with removed extension), followed by an underscore and the internal integer chunk number.

Value

returns (invisibly) the path of the created netCDF file(s)

Note

Packing is currently ignored if chunked = TRUE

See Also

gdalcubes_options

pack_minmax

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-04"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
write_ncdf(select_bands(raster_cube(L8.col, v), c("B04", "B05")), fname=tempfile(fileext = ".nc"))

Export a data cube as a collection of GeoTIFF files

Description

This function will time slices of a data cube as GeoTIFF files in a given directory.

Usage

write_tif(
  x,
  dir = tempfile(pattern = ""),
  prefix = basename(tempfile(pattern = "cube_")),
  overviews = FALSE,
  COG = FALSE,
  rsmpl_overview = "nearest",
  creation_options = NULL,
  write_json_descr = FALSE,
  pack = NULL
)

Arguments

Details

If write_json_descr is TRUE, the function will write an additional file with name according to prefix (if not missing) or simply cube.json This file includes a serialized description of the input data cube, including all chained data cube operations.

Additional GDAL creation options for resulting GeoTIFF files must be passed as a named list of simple strings, where element names refer to the key. For example, creation_options = list("COMPRESS" = "DEFLATE", "ZLEVEL" = "5") would enable deflate compression at level 5.

To reduce the size of created files, values can be packed by applying a scale factor and an offset value and using a smaller integer data type for storage. The pack argument can be either NULL (the default), or a list with elements type, scale, offset, and nodata. type can be any of "uint8", "uint16" , "uint32", "int16", or "int32". scale, offset, and nodata must be numeric vectors with length one or length equal to the number of data cube bands (to use different values for different bands). The helper function pack_minmax can be used to derive offset and scale values with maximum precision from minimum and maximum data values on original scale.

If overviews=TRUE, the numbers of pixels are halved until the longer spatial dimensions counts less than 256 pixels. Setting COG=TRUE automatically sets overviews=TRUE.

Value

returns (invisibly) a vector of paths pointing to the created GeoTIFF files

See Also

pack_minmax

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-04"),
              srs="EPSG:32618", nx = 497, ny=526, dt="P1M")
write_tif(select_bands(raster_cube(L8.col, v), c("B04", "B05")), dir=tempdir())