如何解决了解 raster::extract 和 terra:extract
我在完全理解 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,1,NA,NA)
cds <- rbind(c(7.5,0),c(7.5,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,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,ymn= -7.31,ymx = 3.67,crs=crs(sps) )
ra2[] <- rep(20,24*31)
ra3 <- merge(ra,ra2)
rb <- raster(ncol=31,crs=crs(sps) )
rb[] <- rep(35,24*31)
rb2 <- raster(ncol=31,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,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,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
解决方法
更新答案
感谢您提出的扩展问题,以及您的坚持,并为这么长时间才给您回复表示歉意。这是 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,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
的{{1}}和raster
的{{1}}(默认)
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,na.rm=T)
# ID ra rb
#1 1 NaN NaN
#2 2 10 35
#3 3 10 35
(默认)和small=TRUE
用于raster
touches=TRUE
测试权重
terra
此问题已在 1.2.17 版中修复。您应该可以像这样在一小时内安装该版本
extract(r_raster,na.rm=T)
# ra rb
# [1,] 10 35
#[2,] 10 35
#[3,] 10 35
extract(r_terra,touches=TRUE)
# ID ra rb
#1 1 10 35
#2 2 10 35
#3 3 10 35
我将在未来几天进一步测试它;希望下周把它送到 CRAN。它曾经有效并且比我做得更快,但显然没有足够的测试用例。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。