extract_geom

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

Argument	Description
cube	source data cube to extract values from
sf	object of class sf, see [sf package](sf::st_as_sf)
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()](base::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

A data cube proxy object

Dimensions:
         low       high count pixel_size chunk_size
t 2018-01-01 2018-04-30     4        P1M          1
y    4345115    4745115   400       1000        320
x   388746.8   766746.8   378       1000        320

Bands:
  name offset scale nodata unit
1 NDVI      0     1    NaN     

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

Reading layer `nycd_1' from data source 
  `/home/marius/R/x86_64-pc-linux-gnu-library/4.3/gdalcubes/nycd.gpkg' 
  using driver `GPKG'
Simple feature collection with 71 features and 0 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 563069.9 ymin: 4483098 xmax: 609761.1 ymax: 4529895
Projected CRS: WGS 84 / UTM zone 18N