了解 raster::extract 和 terra:extract

Posted 2023-03-23

【中文标题】了解 raster::extract 和 terra:extract

【英文标题】：Understanding raster::extract and terra:extract

【发布时间】：2021-08-12 09:37:34

【问题描述】：

我在完全理解 terra:extract 时遇到问题。我希望提取管理 GADM 多边形的平均栅格值。我的栅格每个国家/地区都有一个值。我希望特定国家/地区内的每个行政多边形具有相同的值，并且包含某些国家边界的某些多边形被分配区域加权平均值。不幸的是，我当前的脚本并非如此。 raster::extract 似乎给出了合理的结果，但不是 terra:extract （请参阅下面的示例代码 - 提供具有不同值的输出）。有人可以根据我下面的代码解释我为什么吗？非常感谢。

## libraries
library(terra)
library(raster)

#===============================================    
## sample example - provides results as expected (1.333, that is (2*0.5+1*1)/1.5)

# sample raster and SpatialPolygons
r <- raster(ncol=2, nrow=3, xmn= 0, ymn= 0, xmx = 30,ymx = 30)
r[] <- c(2, 2, 2, 1, NA, NA)
cds <- rbind(c(7.5,0), c(7.5,20), c(30, 20),c(30,10))
library(sp)
p = Polygon(cds)
ps = Polygons(list(p),1)
sps = SpatialPolygons(list(ps))
plot(r)
plot(sps, add=T)

# test raster package
test1 <- raster::extract(r , sps, fun=mean, na.rm=T,  weights=TRUE) 
test1 # I get 1.333333 which is what I would expect

# test terra package
sps.spatv <- vect(sps)
r.spatR <-  rast(r) #conversion to SpatRaster class

test2 <- terra::extract(r.spatR, sps.spatv, fun=mean, na.rm=T,  weights=TRUE, exact=TRUE, touches=TRUE) 
test2 # I get 1.333333 which is what I would expect  

#===============================================    
## sample code that leads to different results between raster and terra packages - I wish to understand why such difference.
# sample SpatialPolygonsDataFrame 
ETH <- getData("GADM", country = 'ETH', level = 2)
SOM <- getData("GADM", country = 'SOM', level = 2)
sps <- bind(ETH, SOM)

# sample raster stack
ra <- raster(ncol=31, nrow=24, xmn= 33.3, ymn=  3.67, xmx = 47.5, ymx = 14.65, crs=crs(sps) )
ra[] <- rep(10, 24*31)
ra2 <- raster(ncol=31, nrow=24, xmn= 33.3, ymn= -7.31 , xmx = 47.5, ymx = 3.67, crs=crs(sps) )
ra2[] <- rep(20, 24*31)
ra3 <- merge(ra, ra2)

rb <- raster(ncol=31, nrow=24, xmn= 33.3, ymn=  3.67, xmx = 47.5, ymx = 14.65, crs=crs(sps) )
rb[] <- rep(35, 24*31)
rb2 <- raster(ncol=31, nrow=24, xmn= 33.3, ymn= -7.31 , xmx = 47.5, ymx = 3.67, crs=crs(sps) )
rb2[] <- rep(45, 24*31)
rb3 <- merge(rb, rb2)

stack.r <- stack(ra3, rb3)
names(stack.r) <- c("ra3", "rb3")

plot(stack.r[[1]])
plot(sps, add=T)

# raster::extract
rastR <- raster::extract(stack.r, sps, fun=mean, na.rm=T,  weights=TRUE)

# > head(rastR)
# [,1] [,2]
# [1,]   10   35
# [2,]   10   35
# [3,]   10   35
# [4,]   10   35
# [5,]   10   35
# [6,]   10   35

rastR2 <- rastR %>%
  cbind(sps@data["GID_2"]) # add ID

# terra::extract
sps.spatv <- vect(sps)
stack.r.spatR <-  rast(stack.r) 
rastT <- terra::extract(stack.r.spatR, sps.spatv, fun=mean, na.rm=T,  exact=TRUE)
# > head(rastT)
# ID ra3 rb3
# [1,]  1  10  10
# [2,]  2  10  10
# [3,]  3  10  10
# [4,]  4  10  10
# [5,]  5  10  10
# [6,]  6  10  10
rastT2 <- rastT %>%
  cbind(sps@data["GID_2"]) # add ID

【参考方案1】：

更新答案

感谢您提出的扩展问题和坚持，并很抱歉花了这么长时间才回复您。这是terra 中的一个错误，我没有立即发现。加权平均结果出现乱码（矩阵未按正确顺序填充）。现已修复：

您的简化示例数据

library(raster)
library(terra)
#terra version 1.2.17

sp <- getData("GADM", country = 'ETH', level = 2)[1:3,]
sv <- vect(sp)

ra <- raster(ncols=31, nrows=24, xmn= 33.3, ymn=  3.67, xmx = 47.5, ymx = 14.65, crs=crs(sp), vals=rep(10, 24*31))
rb <- raster(ncols=31, nrows=24, xmn= 33.3, ymn=  3.67, xmx = 47.5, ymx = 14.65, crs=crs(sv), vals=rep(35, 24*31))

r_raster <- stack(ra, rb)
names(r_raster) <- c("ra", "rb")
r_terra <-  rast(r_raster) 

测试无权重和small=FALSE 用于raster 和touches=FALSE 用于terra（默认）

extract(r_raster, sp, fun=mean, na.rm=T, small=FALSE)
#     [,1] [,2]
#[1,]   NA   NA
#[2,]   10   35
#[3,]   10   35

extract(r_terra, sv, fun=mean, na.rm=T)
#  ID  ra  rb
#1  1 NaN NaN
#2  2  10  35
#3  3  10  35

测试无权重和small=TRUE 用于raster（默认）和touches=TRUE 用于terra

extract(r_raster, sp, fun=mean, na.rm=T)
#      ra rb
# [1,] 10 35
#[2,] 10 35
#[3,] 10 35

extract(r_terra, sv, fun=mean, na.rm=T, touches=TRUE)
#  ID ra rb
#1  1 10 35
#2  2 10 35
#3  3 10 35
 

权重测试

extract(r_raster, sp, fun=mean, na.rm=T,  weights=TRUE)
#     ra rb
#[1,] 10 35
#[2,] 10 35
#[3,] 10 35

extract(r_terra, sv, fun=mean, na.rm=T,  weights=TRUE)
#     ID ra rb
#[1,]  1 10 35
#[2,]  2 10 35
#[3,]  3 10 35

这已在版本 1.2.17 中得到修复。你应该可以像这样在一小时内安装那个版本

install.packages('terra', repos='https://rspatial.r-universe.dev')

我将在接下来的几天里进一步测试它；希望下周能把它送到克兰。它曾经工作过，而且比我做得更快，但显然没有足够的测试用例。

【讨论】：

非常感谢，我在我的问题中添加了更多细节。你的例子对我来说很有意义。你的结果是我所期望的。然后我真的想知道为什么我的两个输出小标题 Nadmin2R 来自 raster extract 和 Nadmin2T 来自 terra extract 提供了如此不同的结果。 Nadmin2R 似乎是正确的，但 Nadmin2T 包含我无法解释的值......我可能误解了一些东西

你能简化你的代码，只显示埃塞俄比亚几个地区的两个提取物的输出吗？所有其他代码都很难看到发生了什么。或者更好的是，以eth <- getData("GADM", country="ETH", level=1")

之类的开头创建一个可重复的示例

感谢您的建议。我已经修改了我的问题以包含一个可重复的示例。

亲爱的罗伯特，在这个阶段，我的理解是要么我没有在光栅堆栈上正确使用 terra::extract 函数，要么 terra::extract 不能在光栅堆栈上使用。如果你再次看到这篇文章，我很想听听你的想法。谢谢

感谢您的坚持和道歉，我花了很长时间才明白您的报告内容