Weight raster cells by overlapping polygons to avoid over-counting when aggregating by polygons
overlap.weight(raster, polygons, count = FALSE, warn = TRUE)a RasterLayer object.
a SpatialPolygons, SpatialPolygonsDataFrame, or simple feature
collection with at least two features. The function will still work with only
one polygon, but values will be unchanged, and the result will be equivalent
to mask(raster, polygons).
a logical indicating whether to return a raster with the count of polygons intersecting each cell, or a raster with original values weighted by 1/number of intersecting polygons.
include warnings? Most common is that the returned raster will be
an intersection of the raster and the polygons. Default TRUE.
a RasterLayer object.
This function takes a raster and a set of polygons as arguments.
It counts the number of polygons that intersect each raster cell. It can
return either a raster with the count of the number of intersecting polygons
as cell values or the original raster with cell values weighted by 1 / the
number of intersecting polygons (the default behavior). Cells that do not
intersect any polygons will receive a value of NA. If the extent of the
polygons is less than the extent of the raster, then the function will warn
that it is cropping the raster to the polygons' extent.
library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.4.0, PROJ 8.1.1; sf_use_s2() is TRUE
library(raster)
#> Loading required package: sp
polys_t <- st_sfc(list(st_polygon(list(rbind(c(2,2), c(2,6),
c(6,6), c(6,2),
c(2, 2)))),
st_polygon(list(rbind(c(8,8), c(4,8),
c(4,4), c(8,4),
c(8,8))))),
crs = st_crs('OGC:CRS84'))
raster_t <- raster(nrows = 10, ncols = 10, xmn = 0,
xmx = 10, ymn = 0, ymx = 10,
vals = 1:100,
crs = CRS(st_crs(polys_t)$proj4string))
overlap.weight(raster_t, polys_t)
#> Warning: Raster objects have different extents. Result for their intersection is returned
#> class : RasterLayer
#> dimensions : 6, 6, 36 (nrow, ncol, ncell)
#> resolution : 1, 1 (x, y)
#> extent : 2, 8, 2, 8 (xmin, xmax, ymin, ymax)
#> crs : +proj=longlat +datum=WGS84 +no_defs
#> source : memory
#> names : layer
#> values : 22.5, Inf (min, max)
#>