Extract and Redial Functions in R-Bitmap Package: Region-Weighted Values
I made this post in another forum , but since I really need an answer, I am posting it here again.
I am working in R and I want to calculate the value of a polygon obtained from intersecting raster cells. The value should take into account the weights at each intersecting cell. When I try to run the "extract" function with a raster and polygon sample, I get different weights, which I manually calculate, resulting in a different final value.
Here's my sample code:
require(raster)
r <- raster(nrow=2, ncol=2, xmn=-180, xmx=60, ymn=-30, ymx=90)
r[] <- c(1,2,4,5)
s <- raster(xmn=-120, xmx=-40, ymn=20, ymx=60, nrow=1, ncol=1)
s.pl <- as(s, 'SpatialPolygons')
w <- raster::extract(r, s.pl, method="simple",weights=T, normalizeWeights=F)
mean.value <- raster::extract(r, s.pl, method="simple",weights=T, fun=mean)
The value I get is 2.14, but according to the actual cell weights it should be 2. More specifically, for each part of the polygon that intersects with another cell, the data is:
Area Value
1800 1
600 2
600 4
200 5
So the end value of the polygon based on the above should be 2.
Could this be due to the projection being in lat / lon? But even when I assign the projection in meters, I have the same result. How can I get the value 2 I'm interested in? I also tried with the "resample" function, but I have other results as well.
My ultimate goal is to create a new raster at different resolutions and extents from the original and assign values based on the cell weights of the original raster that intersect with the cells in the new raster. But it seems that neither the rescues nor the extraction functions give the expected result.
source to share
Suppose we have a raster A
and two SpatialPolygon objects [B, C]
that are not rectangular (in this case, they will be hexagons). For demonstration purposes, the center of the hexagon B
is defined as the center of our raster A
(see left graph below). The hexagon is C
shifted to the right along the horizontal axis.
require(raster)
require(scales)
A <- raster(nrow=2, ncol=2, xmn=-180, xmx=180, ymn=-180, ymx=180)
A[] <- c(1,2,4,5)
A.pl <- as(A, 'SpatialPolygons')
B <- SpatialPolygons(list(Polygons(list(Polygon(cbind(c(0, 100, 100, 0, -100, -100, 0),
c(100, 50, -50, -100, -50, 50, 100)))), 'B')))
C <- SpatialPolygons(list(Polygons(list(Polygon(cbind(c(40, 140, 140, 40, -60, -60, 40),
c(100, 50, -50, -100, -50, 50, 100)))), 'C')))
Object B
Since the hexagon B
is in the center, the weights should be 0.25. You can see from the graph that the area of the hexagons is 30,000 (imagine a square that fits a hexagon (40,000) and subtracts 2 rectangles (-10000), each of which consists of 2 of 4 corners that you have to cut off), so each intersection area has size 7500 and7500/30000 = 0.25
# get intersections
intsct.B <- raster::intersect(B, A.pl)
intsct.C <- raster::intersect(C, A.pl)
### B
area.B <- B@polygons[[1]]@area
weights <- unlist(lapply(intsct.B@polygons, function(x) {
slot(x, 'area')/area.B
}))
weights
> [1] 0.25 0.25 0.25 0.25
Now we get the cell value that each intersection polygon falls into and calculate the average.
vals <- unlist(lapply(intsct.B@polygons, function(x) {
extract(A, data.frame(t(slot(x, 'labpt'))))
}))
sum(weights * vals)
> [1] 3
As you would expect, the mean c(1, 2, 4, 5)
is 3
.
Object C
Now let's do the same with the object. C
### C
area.C <- C@polygons[[1]]@area
weights <- unlist(lapply(intsct.C@polygons, function(x) {
slot(x, 'area')/area.C
}))
weights
> [1] 0.13 0.37 0.13 0.37
vals <- unlist(lapply(intsct.C@polygons, function(x) {
extract(A, data.frame(t(slot(x, 'labpt'))))
}))
sum(weights * vals)
> [1] 3.24
Again, as you would expect, the average is larger (since the weights for cells with values 2 and 5 are higher). Also, since we only shifted the hexagon along one axis, it makes sense that the two weights appear twice.
Raster with many cells
The following graph shows the intersections B
(left side) and C
(rhs) with a raster 4x4
with values c(1:8, 10:17)
. For B
there are 12 crossings and for C
8. Note again that the mean for B
is 9 due to symmetry.
This should work for any object SpatialPolygons
. Remember to use the same CRS for the objects you drop into intersect
.
source to share
Here's what I was able to do based on the answers in this post.
require(raster)
require(rgeos)
r <- raster(nrow=2, ncol=2, xmn=-180, xmx=60, ymn=-30, ymx=90)
r[] <- c(1,2,4,5)
r <- stack(r, r*2, r^2)
s <- raster(xmn=-120, xmx=-40, ymn=20, ymx=60, nrow=1, ncol=1)
s.pl <- as(s, 'SpatialPolygons')
r.s <- as(r, 'SpatialPolygonsDataFrame')
pi1 <- gIntersection(r.s, s.pl, byid = T)
areas1 <- data.frame(area=sapply(pi1@polygons, FUN=function(x) {slot(x, 'area')}))
row.names(areas1) <- sapply(pi1@polygons, FUN=function(x) {slot(x, 'ID')})
areas1$Pol.old <- as.numeric(vapply(strsplit(rownames(areas1), " "), `[`, 1, FUN.VALUE=character(1)))
areas1$pol.new <- as.numeric(vapply(strsplit(rownames(areas1), " "), `[`, 2, FUN.VALUE=character(1)))
f <- r.s@data
seqs <- match(areas1$Pol.old, rownames(f))
ar <- cbind(areas1, f[seqs,])
ar[,-(1:3)] <- ar[,-(1:3)]*ar$area
f <- aggregate.data.frame(ar, by=list(ar$pol.new), FUN=sum)
f[,-(1:4)] <- f[,-(1:4)]/f$area
ar.v <- as.matrix(f[, -c(1:4)])
s2 <- stack(s)
s1 <- setValues(s2, ar.v)
If anyone can suggest a more convenient and / or faster code, please let me know, because I don't really like my approach.
source to share