Package 'gdalcubes'

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

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

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

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

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

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)

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

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)

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)

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)

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)

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)

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)

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

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)

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

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)

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

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

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

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

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

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

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

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

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

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

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

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)

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)

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.

Description

Check if GDAL was built with GEOS

Usage

gdalcubes_gdal_has_geos()

Examples

gdalcubes_gdal_has_geos()

Description

Get available GDAL drivers

Usage

gdalcubes_gdalformats()

Examples

gdalcubes_gdalformats()

Description

Get the GDAL version used by gdalcubes

Usage

gdalcubes_gdalversion()

Examples

gdalcubes_gdalversion()

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

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

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/user/configoptions.html.

Examples

gdalcubes_set_gdal_config("GDAL_NUM_THREADS", "ALL_CPUS")

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

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)

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

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)

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

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

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

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)

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

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

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

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

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

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)

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

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)