Extract Zonal Statistics

library(mapme.vegetation)
library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE
library(stringr)

This vignette showcases different ways to extract zonal and pixelwise statistics from a time series processed with the mapme.vegetation package. We assume that we have processed S2 data to a dense time-series following vignette xx. Let’s read in the files and our AOI.

ndvi_files = s2_files = list.files(system.file("extdata/filled", package = "mapme.vegetation"), ".tif", full.names = T)
bands = "NDVI"
times = str_sub(basename(ndvi_files), -14, -5)

aoi = st_read(system.file("extdata", "exp_region_wgs.gpkg", package = "mapme.vegetation"), quiet = T)
aoi_utm = st_transform(aoi, st_crs(32638))
aoi_utm$id = 1
bbox = st_bbox(aoi_utm)

As an example we will extract zonal statistics for the complete AOI. This function relies on gdalcubes which is why we need information of the spatio-temporal extent of the input files. These are the same as in the other vignettes. With this function we get information averaged over all pixels within a polygon for each time step. This information is useful for statistical analysis and might also be suitable for object based classification or regression problems which require information through the time axis. Specific metrics can be calculated. Thanks to the gdalcubes package, any R function that summarizes the data per zone can be supplied to a named list. It is important to declare a column as the primary identification column with the argument idcol. Note that the returned object will be in a time-long format. That also means that the geometry column will be repeated. Because of this, the resulting object might be quite large. You can specify a filepath with the outpath argument to write the object to a GeoPackage (.gpkg).

zs_aoi = extract_zonalstats(aoi = aoi_utm, 
                            idcol = "id", 
                            files = ndvi_files, 
                            times = times, 
                            bands = bands, 
                            bbox = bbox, 
                            zonalfuns = c(sum = sum, mean = mean, median = median, ncells = length), 
                            after = "2017-05-01", 
                            before = "2017-07-31", 
                            srs = "EPSG:32638", 
                            dx = 200, 
                            dy = 200, 
                            dt = "P14D", 
                            aggregation = "max", 
                            resampling = "near")
head(zs_aoi)
#> Simple feature collection with 6 features and 4 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 348194.9 ymin: 1056042 xmax: 362482.2 ymax: 1066928
#> Projected CRS: WGS 84 / UTM zone 38N
#>         time     NDVI stat id                           geom
#> 1 2017-05-01 440.8433  sum  1 POLYGON ((348194.9 1056042,...
#> 2 2017-05-15 555.8793  sum  1 POLYGON ((348194.9 1056042,...
#> 3 2017-05-29 570.3013  sum  1 POLYGON ((348194.9 1056042,...
#> 4 2017-06-12 569.3919  sum  1 POLYGON ((348194.9 1056042,...
#> 5 2017-06-26 578.9204  sum  1 POLYGON ((348194.9 1056042,...
#> 6 2017-07-10 600.0071  sum  1 POLYGON ((348194.9 1056042,...

On other occasions pixel-wise information for polygons might be of interest, e.g. when the goal of your analysis is the pixel-wise classification of a landscape. For this kind of extraction we rely on the terra package. For the resulting data.frame object to receive the right column names, however, it is important to hand over the band names and the associated time-stamp for each file.

ps_aoi = extract_pixels(files = ndvi_files,
                        bands = bands,
                        times = times, 
                        aoi = aoi_utm)
head(ps_aoi)
#>   ID NDVI-2017-05-01 NDVI-2017-05-15 NDVI-2017-05-29 NDVI-2017-06-12
#> 1  1      0.10098172       0.1437648       0.1508505       0.1471979
#> 2  1      0.09259786       0.1310685       0.1391816       0.1403345
#> 3  1      0.08804486       0.1242767       0.1325152       0.1372536
#> 4  1      0.08691528       0.1218589       0.1303333       0.1359665
#> 5  1      0.08776132       0.1230615       0.1315388       0.1365081
#> 6  1      0.09056787       0.1301583       0.1377396       0.1404345
#>   NDVI-2017-06-26 NDVI-2017-07-10 NDVI-2017-07-24 id
#> 1       0.1449453       0.1495825       0.1599105  1
#> 2       0.1381166       0.1425987       0.1544822  1
#> 3       0.1353973       0.1394072       0.1510371  1
#> 4       0.1338123       0.1368362       0.1471831  1
#> 5       0.1342536       0.1372983       0.1475986  1
#> 6       0.1394512       0.1461845       0.1615133  1

The resulting object is not an sf object. The spatially explicit information is list. However, the column ID indicates to to which feature in the aoi object the pixel information belongs in the order of rows of the input object. The data.frame object has a wide format. Each row represents a single pixel. The columns identify the band name as well as the time-stamp.

Last modified: 2023-02-10