| extract_geom {gdalcubes} | R Documentation |
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.
extract_geom( cube, sf, datetime = NULL, time_column = NULL, FUN = NULL, ..., reduce_time = FALSE )
cube |
source data cube to extract values from |
sf |
object of class |
datetime |
Date, POSIXt, or character vector containing per feature time information; length must be identical to the number of features in |
time_column |
name of the column in |
FUN |
optional function to compute per feature summary statistics |
... |
additional arguments passed to |
reduce_time |
logical; if TRUE, time is ignored when |
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.
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.
A data.frame with columns FID, time, and data cube bands / variables
# 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"))
}
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)) {
x = runif(20, v$space$left, v$space$right)
y = runif(20, v$space$bottom, v$space$top)
t = sample(seq(as.Date("2018-01-01"),as.Date("2018-04-30"), by = 1),20, replace = TRUE)
df = sf::st_as_sf(data.frame(x = x, y = y), coords = c("x", "y"), crs = v$space$srs)
# spatiotemporal points
extract_geom(L8.ndvi, df, datetime = t)
# time series at spatial points
extract_geom(L8.ndvi, df)
# 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)
zstats
# combine with original sf object
x$FID = rownames(x)
x = merge(x, zstats, by = "FID")
x
# plot(x["NDVI"])
}
}