extract() data from raster with small polygons - rounded weights too small

Question

Using R, I am trying to extract data from a raster layer using a polygon layer. The polygons are much smaller than the raster cells:

Now I call extract() from raster library:

a <- extract(raster, polygons, weights = TRUE, small = TRUE)
a
# ...
# [[1551]]
# value weight
#   209   0.03 # top left cell - more than 50% of the polygon area

There are two problems - the weight is the proportion of the cell area covered by the polygon, and the weights are rounded to 1/100. In my case, only the top left cell is present in the output (value 209) - the weight of 3 other cells was rounded to zero and they were excluded. However, the bottom left cell covers significant proportion of the polygon and should be included also!

I need a proper weighted mean. Can this be done somehow else using extract()? Or any other way?

PS: note aside: I think the weights in extract() are not designed very well - the weight should be the proportion of polygon area covered by the particular cell, not vice versa. Then, the weighted mean for the polygon would be also easier to compute (just multiply the two numbers in each row and sum up), and rounding to 1/100 wouldn't be a big problem.

Reproducible example - (download the files - simplified version, actual data are much bigger):

require(raster)
rast <- raster("my.tif")
poly <- readOGR(".", "socc_buff_Rx")
a <- extract(rast, poly, weights = TRUE, small = TRUE)
a

Related: Extract in R fails for small polygons and raster

Ledbetter

Can't you expand raster resolution to capture (or approach) the covered area proportion? – Blagoveshchensk

@PauloCardoso thanks - this was actually one of my 4 ideas how to solve problem: 1) persuade extract to sane behaviour, 2) convert raster cells to polygons, 3) read raster in R cell by cell, in fact implementing my own extract, 4) increase the raster resolution by further splitting the cells (your idea). So it depends which solution is in fact easiest and also computationally more efficient. And I don't know yet how to implement these solutions... – Ledbetter

A reproducible example would make this question useful. – Mudslinging

@Ledbetter I have the same question, have you found an answer? – Deme

I second the need for an answer to this question. Did @Mudslinging have any suggestion? – Blondie

I haven't found any solution to this (I'd post it if I did). My workaround was to enlarge the polygons (buffers). Real solution would probably be to send a change request to the authors of the extract function with a link to this page. Anyone is welcome to do that. – Ledbetter

Do you suggest any systematic way of doing it? Related either to the size of the polygon or the resolution of the polygons? I'm working on raster with resolution of 8x8 km. – Blondie

Answer 1

I think the easiest, albeit inelegant, solution is to disaggregate the RasterLayer first. I will have a look to see if I can change the extract function to do this automatically for very small (relative to the cells size) polygons.

library(raster)
r <- raster("my.tif")
pu <- shapefile("socc_buff_Rx.shp")
p <- spTransform(pu, crs(r))

extract(r, p, weights = TRUE, small = TRUE)
#[[1]]
# value weight
#   209   0.03

rr <- disaggregate(r, 10)
e <- extract(rr, p, weights = TRUE, small = TRUE)
lapply(e, function(x) { aggregate(x[,2,drop=F], list(value=x[,1]), sum ) } )

#[[1]]
#  value weight
#1   197   0.95
#2   209   3.44
#3   256   0.31
#4   293   0.04

plot(r, legend=F)
plot(p, add=T)
text(r)

Backflow