[GIS] extract() function from r raster package returns different cells if “weights” parameter is TRUE or FALSE

extract-by-maskrraster

I am trying to extract cell values from a RasterLayer based on a SpatialPolygons object. I create the raster with

library(raster)

rasterValues <- matrix(rnorm(20),4,5)
r <- raster(rasterValues)
extent(r) <- extent(13.875,14.5,45.125,45.625)
projection(r) <- crs("+proj=longlat +datum=WGS84")

and then I create the SpatialPolygons with

x <- c(14.00, 14.27, 14.27, 14.00)
y <- c(45.25, 45.25, 45.5, 45.5)
spPoints <- cbind(x,y)
sPoly <- Polygon(spPoints)
sPolys <- Polygons(list(sPoly),1)
sPolyg <- SpatialPolygons(list(sPolys), proj4string = crs("+proj=longlat +datum=WGS84"))

As you can see from the picture, the polygon (colored in black) totally covers 4 cells and a small part of two cells of my raster.enter image description here

When I try to extract values with weights=TRUE, I get:

cellValue <- extract(r,sPolyg, cellnumbers=TRUE, weights=TRUE)
cellValue <- cellValue[[1]]
cellValue

      cell      value     weight
[1,]    7 -0.1658694 0.22727273
[2,]    8 -0.8996764 0.22727273
[3,]    9 -0.6951724 0.04545455
[4,]   12  0.5717841 0.22727273
[5,]   13  1.5263394 0.22727273
[6,]   14 -1.2750113 0.04545455

When I try to run the previous command with weights=FALSE, I get:

cellValue <- extract(r,sPolyg, cellnumbers=TRUE, weights=FALSE)
cellValue <- cellValue[[1]]
cellValue

cell      value
[1,]    7 -0.1658694
[2,]    8 -0.8996764
[3,]   12  0.5717841
[4,]   13  1.5263394

The two cells that are only partly covered by the SpatialPolygons object are not returned.

The question here is: why are some cells not returned, when we set the "weight" parameter to FALSE?

Best Answer

As described in ?extract,

A cell is covered if its center is inside the polygon (but see the weights option for considering partly covered cells; and argument small for getting values for small polygons anyway).

Therefore, if you run the following code using weights = FALSE (default), only values from the 4th and 6th cell are returned. On the other hand, values from the 7th and 9th cell are lacking since their respective centroids (small black crosses) are not covered by the single polygons. 

r <- raster(ncol = 36, nrow = 18)
r[] <- 1:ncell(r)
r <- aggregate(r, 8)

cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
polys <- spPolygons(cds1, cds2)

plot(r, col = rainbow(100))
plot(polys, add = TRUE)

## add centroids
pys <- rasterToPolygons(r)
cnt <- rgeos::gCentroid(pys, byid = TRUE)
plot(cnt, add = TRUE)

centroids

> extract(r, polys)
[[1]]

418.5 

[[2]]

154.5

You are explicitly required to set weights = TRUE in order to additionally include those partially covered cells in your analysis. 

> extract(r, polys, weights = TRUE)
[[1]]
     value    weight
[1,] 418.5 0.8548387
[2,] 426.5 0.1451613

[[2]]
     value weight
[1,] 154.5   0.35
[2,] 442.5   0.65
Related Question