R raster::extract function not working as expected

Question

I'm trying to do some polygon/raster work in R. The extract function isn't working as I'd expect. Here's a reproducible example.

library(sf)
library(raster)
library(maptools)

raster      <-     download.file('https://raw.githubusercontent.com/JimShady/london_osm_canyons/master/test.tif', 'downloaded_raster')
raster      <- raster('test.tif')

polygons    <- read_sf('https://raw.githubusercontent.com/JimShady/london_osm_canyons/master/test2.geojson')
polygons    <- st_transform(polygons, crs(raster)@projargs)
polygons$id <- 1
polygons    <- as(polygons, 'Spatial')

new_polygon <- elide(polygons, rotate=-50)
proj4string(new_polygon) <- crs(raster)@projargs

polygons    <- rbind(polygons, new_polygon)

plot(raster)
plot(polygons, add=T)

I want to get the number of cells that the each polygon intersects with, even if they are NA, and store them in the polygon as a column

polygons$cell_count     <- extract(raster, polygons, fun = function(x,...)length(x), na.rm=F)

But this returns the number 5 for polygon one, when clearly the polygon intersects more than 5 cells of the raster. What's going on with that?

I then want to get the weighted mean. if there are any NA cells intersected, then this should 'fail' and return NA. I think this works, but as the above is wrong, I'm not very confident of my method now.

polygons$mean           <- round(extract(raster, polygons, fun = mean, weights=T, na.rm=F),2)

I also want to know the number of cells that intersect that are less than one. I'm not sure how to do that. Anyone got a good idea?

polygons$almost_zero    <-

Thanks for any help.

Answer 1

From extract help:

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, you could use this to have the number of cells:

data <- extract(raster, polygons,  weights = T,  na.rm = F)[[1]]
nrow(data)

#[1] 18

Where data is:

      value      weight
 [1,]     8 0.003610108
 [2,]     6 0.025270758
 [3,]     0 0.010830325
 [4,]     7 0.104693141
 [5,]     0 0.066787004
 [6,]     0 0.088447653
 [7,]     0 0.001805054
 [8,]     0 0.104693141
 [9,]     0 0.010830325
[10,]     0 0.093862816
[11,]     0 0.030685921
[12,]     0 0.063176895
[13,]     0 0.057761733
[14,]     0 0.063176895
[15,]     7 0.003610108
[16,]     7 0.102888087
[17,]     0 0.093862816
[18,]     0 0.074007220

To have the number of cells that intersects and are less than 1:

sum(data[, 1] < 1)

#[1] 13

---

EDIT

In case of several polygons, you should first dissolve your polygons using rgeos::gUnaryUnion:

library(rgeos)

polygons$id <- 1
polygons <- gUnaryUnion(polygons, id = polygons$id)
# As suggested by @RobertHijmans, you can also use `raster::aggregate`
# polygons <- raster::aggregate(polygons, by = "id")

plot(raster)
plot(polygons, add = T)

Answer 2

Extract is behaving as it is meant to, by finding cells where the centroid lies within the polygon. One way to deal with partial coverage is to convert your polygon to a raster mask, using rasterize. Note that in the following I renamed your raster and polygon r1, and p1, respectively. Because its not good practice to use function names for objects:

r_mask = rasterize(p1, r1, getCover=TRUE)
r_mask[r_mask==0] = NA

Now we can use this mask to get the values you want:

r2 = mask(r1, r_mask)
cellStats(!is.na(r2), sum)
# [1] 18
cellStats(r2, mean)
# [1] 1.944444

To find how many cells are less than one, we can do

cellStats(r2<1, sum)
# [1] 13