extract() data from raster with small polygons
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
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)
上一篇: 在Objective中以编程方式计算IRR（内部利率回报）和NPV
下一篇: 从具有小多边形的栅格提取（）数据