| zonal_statistics {gdalcubes} | R Documentation |
This function will overlay spatial polygons with a data cube and compute summary statistics of pixel values within the polygons over all time slices.
zonal_statistics( x, geom, expr, out_path = tempfile(fileext = ".gpkg"), overwrite = FALSE, ogr_layer = NULL, as_stars = FALSE )
x |
source data cube |
geom |
Either an sf object, or a path to an OGR dataset (Shapefile, GeoPackage, or similar) with input (multi)polygon geometries |
expr |
character vector of summary statistics expressions, describing pairs of aggregation functions and data cube bands (e.g. "mean(band1)") |
out_path |
path to where resulting GeoPackage will be written to |
overwrite |
logical; overwrite |
ogr_layer |
If the input OGR dataset has multiple layers, a layer can be chosen by name |
as_stars |
logical; if TRUE, the created gpkg file will be loaded as a stars vector data cube |
The function creates a single GeoPackage output file containing:
A single layer "geom" containing the geometries (and feature identifiers) only.
Attribute tables (layers without geometry) for each time slice of the data cube containing summary statistics as columns. Corresponding layer names start with "attr_", followed by date and time.
Virtual spatial views for each time slice, joining the geometries and attribute tables. Corresponding layer names start with "map_", followed by date and time.
You will most-likely want to use the spatial view layers directly e.g. with the sf package.
Available summary statistics currently include "min", "max", "mean", "median", "count", "sum", "prod", "var", and "sd".
character length-one vector containing the path to the resulting GeoPackage file (see Details) or a stars object (if as_stars is TRUE)
Currently, the spatial reference systems of the data cube and the features must be identical.
This function requires GDAL with built-in GEOS support, which can checked with gdalcubes_gdal_has_geos)
# 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=300, dx=300, 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
# toy example: overlay NDVI data with NYC districts
if (gdalcubes_gdal_has_geos()) {
x = zonal_statistics(L8.ndvi, system.file("nycd.gpkg", package = "gdalcubes"),
expr = "median(NDVI)")
x
}