3

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?

1

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
1

Short answer: because with weights=FALSE, only cells whose centroid falls into a polygon are accounted. With weights=TRUE, all cells that have at least a small part (but not necessarily the centroid) are reported.

Long answer: There are two aspects to take into account, weight and small, and hence multiple cases:

  • weight = FALSE (the case described by @fdetsch): return for each cell the polygon in which the centroid of the cell falls.

    • if small=FALSE: rule above. A raster cell will hence be attributed to at most one polygon. A polygon will be represented zero or more times.
    • if small=TRUE: add an exception for small polygons which do not contain any cell centroid, i.e. cells will be attributed possibly to zero, one or multiple polygons. A polygon will be represented one or more times.
  • weight=TRUE (your case): indicate the weight for all polygon found on the cell. Setting normalizeWeights=FALSE helps understand how much of each polygon a cell is covering. There will be as many weights as (parts of) polygons on the cell.

    • argument small: I don't think it plays a role in this case.

Illustration: To illustrate, let us add a very small polygon for which no cell's centroid fall into. And let us return the full data-frame (df=TRUE), and compare the four cases. enter image description here

### previous code
## small poly
x2 <- c(14.29, 14.31, 14.31, 14.29)
y2 <- c(45.48, 45.48, 45.49, 45.49)
sPolys2 <- Polygons(list(Polygon(cbind(x2,y2))),2)
sPolyg <- SpatialPolygons(list(sPolys, sPolys2), 
  proj4string =     crs("+proj=longlat +datum=WGS84"))

r_shp <- rasterToPolygons(r)

library(tmap)

tm_shape(r_shp)+
  tm_borders() +
  tm_shape(sPolyg) +
  tm_fill(col="blue") +
  tm_shape(r_shp)+
  tm_borders() +
  tm_text("layer")

Now compares the four cases:

cellValue_W_S <- raster::extract(r,sPolyg, cellnumbers=TRUE, df=TRUE, 
                             weights=TRUE, normalizeWeights=FALSE,
                             small=TRUE)
cellValue_W_nS <- raster::extract(r,sPolyg, cellnumbers=TRUE, df=TRUE, 
                              weights=TRUE, normalizeWeights=FALSE, 
                              small=FALSE)
cellValue_nW_S <- raster::extract(r,sPolyg, cellnumbers=TRUE, df=TRUE,
                              weights=FALSE,
                              small=TRUE)
cellValue_nW_nS <- raster::extract(r,sPolyg, cellnumbers=TRUE, df=TRUE,
                               weights=FALSE, 
                               small=FALSE)

Summarising for the cell 9, and showing the weights in each scenario:

   cell    ID `wgt=F, sml=F` `wgt=F, sml=T` `wgt=T, sml=F` `wgt=T, sml=T`
* <dbl> <dbl>          <dbl>          <dbl>          <dbl>          <dbl>
1     9     1              0              0           0.20           0.20
2     9     2              0              1           0.02           0.02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.