我正在编写一个脚本,它将采用任意三个栅格,并尽可能地将它们裁剪掉。所有三个栅格都具有相同的分辨率和投影。但是,裁剪到最小程度不会改变三个栅格的范围。我尝试过setExtent,同样的事情发生了。如果有人能提出建议,我会非常感激。这是示例代码:
library(raster)
#Projection of all three rasters
newproj<- "+proj=utm +zone=4 +datum=WGS84 +units=m +no_defs +ellps=WGS84
+towgs84=0,0,0"
#Create three rasters with varying extents
raster1p<- raster(crs = newproj)
extent(raster1p)<- c(531247, 691847, 7856684, 7987884)
res(raster1p)<- c(100, 100)
values(raster1p)<- NA
raster2p<- raster(crs = newproj)
extent(raster2p)<- c(533550.8, 646550.8, 7881307, 7973807)
res(raster2p)<- c(100, 100)
values(raster2p)<- NA
raster3p<- raster(crs = newproj)
extent(raster3p)<- c(525739, 689839, 7857305, 7996505)
res(raster3p)<- c(100, 100)
values(raster3p)<- NA
#Find minimum extent
xmin1<- c(xmin(extent(raster1p)), xmin(extent(raster2p)), xmin(extent(raster3p)))
xmax1<- c(xmax(extent(raster1p)), xmax(extent(raster2p)), xmax(extent(raster3p)))
ymin1<- c(ymin(extent(raster1p)), ymin(extent(raster2p)), ymin(extent(raster3p)))
ymax1<- c(ymax(extent(raster1p)), ymax(extent(raster2p)), ymax(extent(raster3p)))
xmin_new<- min(xmin1)
xmax_new<- min(xmax1)
ymin_new<- min(ymin1)
ymax_new<- min(ymax1)
newextent=c(xmin_new, xmax_new, ymin_new, ymax_new)
#Crop rasters to minimum extent
crop(raster1p, newextent)
crop(raster2p, newextent)
crop(raster3p, newextent)
#Compare extents
extent_check<- c(extent(raster1p), extent(raster2p), extent(raster3p))
但是,当我查看extent_check以查看扩展区现在是否匹配时,我发现扩展区根本没有改变:
> extent_check
[[1]]
class : Extent
xmin : 531247
xmax : 691847
ymin : 7856684
ymax : 7987884
[[2]]
class : Extent
xmin : 533550.8
xmax : 646550.8
ymin : 7881307
ymax : 7973807
[[3]]
class : Extent
xmin : 525739
xmax : 689839
ymin : 7857305
ymax : 7996505
知道我可能做错了吗?感谢您的时间
我认为这不是一件坏事,而是一个误解(尽管你的代码中存在错误)。
示例数据
library(raster)
prj <- "+proj=utm +zone=4 +datum=WGS84"
r1 <- raster(res=100, ext=extent(c(531247, 691847, 7856684, 7987884)), crs=prj, vals=NA)
r2 <- raster(res=100, ext=extent(c(533550.8, 646550.8, 7881307, 7973807)), crs=prj, vals=NA)
r3 <- raster(res=100, ext=extent(c(525739, 689839, 7857305, 7996505)), crs=prj, vals=NA)
找到“最小范围”
e <- intersect(intersect(extent(r1), extent(r2)), extent(r3))
请注意,结果与您的结果不同,因为您使用
xmin_new <- min(xmin1)
和ymin_new <- min(ymin1)
它应该在哪里
xmin_new <- max(xmin1)
和ymin_new <- max(ymin1)
现在裁剪
r1e <- crop(r1, e)
r2e <- crop(r2, e)
r3e <- crop(r3, e)
检查生成的范围
t(sapply(c(r1e, r2e, r3e), function(i) as.vector(extent(i))))
# [,1] [,2] [,3] [,4]
#[1,] 533547.0 646547.0 7881284 7973784
#[2,] 533550.8 646550.8 7881307 7973807
#[3,] 533539.0 646539.0 7881305 7973805
它们并不完全相同,因为这是不可能的,因为栅格不对齐。他们的“起源”是不同的
t(sapply(c(r1e, r2e, r3e), origin))
# [,1] [,2]
#[1,] 47.0 -16
#[2,] -49.2 7
#[3,] 39.0 5
为了使它们对齐,你需要做这样的事情
r1e <- crop(r1, e)
r2e <- resample(r2, r1e)
r3e <- resample(r3, r1e)
t(sapply(c(r1e, r2e, r3e), function(i) as.vector(extent(i))))
# [,1] [,2] [,3] [,4]
#[1,] 533547 646547 7881284 7973784
#[2,] 533547 646547 7881284 7973784
#[3,] 533547 646547 7881284 7973784