extract_geom {gdalcubes}

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

cube

source data cube to extract values from

sf

object of class sf, see sf package

datetime

Date, POSIXt, or character vector containing per feature time information; length must be identical to the number of features in sf

time_column

name of the column in sf containing per feature time information

FUN

optional function to compute per feature summary statistics

merge

logical; return a combined data.frame with data cube values and labels, defaults to FALSE

drop_geom

logical; remove geometries from output, only used if merge is TRUE, defaults to FALSE

...

additional arguments passed to FUN

reduce_time

logical; if TRUE, time is ignored when FUN is applied

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