微信公众号搜"智元新知"关注
微信扫一扫可直接关注哦!

使用R中的terra包在不同程度上镶嵌栅格时出错?

如何解决使用R中的terra包在不同程度上镶嵌栅格时出错?

我有一个文件夹,其中包含大约 191 个 GeoTIFF 文件(每个文件都是更大区域的不同 DEM(高程)图块)。我想将所有图块合并到一个光栅文件中。我正在使用 terra 包并且能够成功加载每个栅格并将它们从 2 米分辨率聚合到 30 米分辨率。但是,当运行 mosaic 函数将它们全部合并时,我遇到了错误(请参阅下面的错误消息)。我已经能够在只有三个图块的较小子集上运行马赛克功能,但是当我放大到所有文件时,这就会成为一个问题。

通过调用栅格摘要(见下文),聚合确实稍微改变了范围 - 这可能是问题吗? resample 可能是一种选择,但每个单独的栅格都有不同的范围,我不确定如何实施此修复。

不确定样本数据集会有所帮助,因为我知道这些函数有效。我在高性能集群上运行此代码,因此运行小批量代码效率不高。

library(terra)

files <- as.list(list.files("./DEM_tiles",full.names = TRUE))

raster.list <- lapply(files,rast)
 
for(i in 1:length(raster.list)){
 raster.list[[i]] <- aggregate(raster.list[[i]],fact = 15)
 }

raster.mosaic <- do.call(mosaic,raster.list)

> Error: [mosaic] internal error: extents do not match ()
> Execution halted

下面是两个图块的示例:

### Before Aggregation
[[1]]
class       : SpatRaster 
dimensions  : 25000,25000,1  (nrow,ncol,nlyr)
resolution  : 2,2  (x,y)
extent      : -1800000,-1750000,-6e+05,-550000  (xmin,xmax,ymin,ymax)
coord. ref. : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=wgs84 +units=m +no_defs 
source      : 35_23_1_1_2m_v3.0_reg_dem.tif 
name        : 35_23_1_1_2m_v3.0_reg_dem 

[[2]]
class       : SpatRaster 
dimensions  : 25000,-550000,-5e+05  (xmin,ymax)
coord. ref. : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=wgs84 +units=m +no_defs 
source      : 35_23_1_2_2m_v3.0_reg_dem.tif 
name        : 35_23_1_2_2m_v3.0_reg_dem 


### After Aggregation
[[1]]
class       : SpatRaster 
dimensions  : 1667,1667,nlyr)
resolution  : 30,30  (x,-1749990,-600010,ymax)
coord. ref. : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=wgs84 +units=m +no_defs 
source      : memory 
name        : 35_23_1_1_2m_v3.0_reg_dem 
min value   :                 -15.62178 
max value   :                  233.6489 

[[2]]
class       : SpatRaster 
dimensions  : 1667,-550010,ymax)
coord. ref. : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=wgs84 +units=m +no_defs 
source      : memory 
name        : 35_23_1_2_2m_v3.0_reg_dem 
min value   :                 -15.27713 
max value   :                  243.0772 

解决方法

我可以复制这样的消息

library(terra)
r1 <- rast(nrow=1667,ncol=1667,ext=c(-1800000,-1749990,-600010,-550000))
r2 <- rast(nrow=1667,-550010,-5e+05))
values(r1) <- 1:ncell(r1)
values(r2) <- 1:ncell(r2)
m <- mosaic(r1,r2)
#Error: [mosaic] internal error: extents do not match ()

而且 merge 也失败了。这是一个我会调查的错误。它适用于 raster

library(raster)
x1 <- raster(r1)
x2 <- raster(r2)
m <- mosaic(x1,x2,fun=mean)
mm <- merge(x1,x2)

稍后:

错误是因为栅格未对齐。使用 terra 1.2-1(目前为 development version),我现在得到

library(terra)
#terra version 1.2.1
crs <- "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84 +units=m"
r1 <- rast(nrow=1667,-550000),crs=crs)
r2 <- rast(nrow=1667,-5e+05),crs=crs)
values(r1) <- 1:ncell(r1)
values(r2) <- 1:ncell(r2)
m <- mosaic(r1,r2)
#Warning message:
#[mosaic] rasters did not align and were resampled

m
#class       : SpatRaster 
#dimensions  : 3334,1667,1  (nrow,ncol,nlyr)
#resolution  : 30,30  (x,y)
#extent      : -1800000,-499990  (xmin,xmax,ymin,ymax)
#coord. ref. : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#source      : memory 
#name        :   lyr.1 
#min value   :       1 
#max value   : 2778889 

还有

@Elia 的建议是一个很好的解决方法:

r1 <- writeRaster(r1,"test1.tif",overwrite=TRUE)
r2 <- writeRaster(r2,"test2.tif",overwrite=TRUE)
v <- vrt(c("test1.tif","test2.tif"),"test.vrt",overwrite=TRUE)

v
#class       : SpatRaster 
#dimensions  : 3334,-600020,-5e+05  (xmin,ymax)
#coord. ref. : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#source      : test.vrt 
#name        :    test 
#min value   :       1 
#max value   : 2778889 

版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。

相关推荐


Selenium Web驱动程序和Java。元素在(x,y)点处不可单击。其他元素将获得点击?
Python-如何使用点“。” 访问字典成员?
Java 字符串是不可变的。到底是什么意思?
Java中的“ final”关键字如何工作?(我仍然可以修改对象。)
“loop:”在Java代码中。这是什么,为什么要编译?
java.lang.ClassNotFoundException:sun.jdbc.odbc.JdbcOdbcDriver发生异常。为什么?
这是用Java进行XML解析的最佳库。
Java的PriorityQueue的内置迭代器不会以任何特定顺序遍历数据结构。为什么?
如何在Java中聆听按键时移动图像。
Java“Program to an interface”。这是什么意思?