Ошибка в мозаике растров разного размера с использованием пакета terra в R?

У меня есть папка, содержащая около 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, 25000, 1  (nrow, ncol, nlyr)
resolution  : 2, 2  (x, y)
extent      : -1800000, -1750000, -550000, -5e+05  (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_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, 1  (nrow, ncol, nlyr)
resolution  : 30, 30  (x, y)
extent      : -1800000, -1749990, -600010, -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      : memory 
name        : 35_23_1_1_2m_v3.0_reg_dem 
min value   :                 -15.62178 
max value   :                  233.6489 

[[2]]
class       : SpatRaster 
dimensions  : 1667, 1667, 1  (nrow, ncol, nlyr)
resolution  : 30, 30  (x, y)
extent      : -1800000, -1749990, -550010, -5e+05  (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        : 35_23_1_2_2m_v3.0_reg_dem 
min value   :                 -15.27713 
max value   :                  243.0772

r gis raster spatial terra
person apple    schedule 19.04.2021    source источник
comment
Один из способов сделать это - создать файл vrt со всеми плитками, а затем записать его как tif.   -  person Elia    schedule 20.04.2021

Ответы (1)


arrow_upward
3
arrow_downward

Я могу воспроизвести сообщение вот так

library(terra)
r1 <- rast(nrow=1667, ncol=1667, ext=c(-1800000, -1749990, -600010, -550000))
r2 <- rast(nrow=1667, ncol=1667, ext=c(-1800000, -1749990, -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 (в настоящее время это версия для разработки) теперь я получаю

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, ncol=1667, ext=c(-1800000, -1749990, -600010, -550000), crs=crs)
r2 <- rast(nrow=1667, ncol=1667, ext=c(-1800000, -1749990, -550010, -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, -1749990, -600010, -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, 1667, 1  (nrow, ncol, nlyr)
#resolution  : 30, 30  (x, y)
#extent      : -1800000, -1749990, -600020, -5e+05  (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      : test.vrt 
#name        :    test 
#min value   :       1 
#max value   : 2778889
person Robert Hijmans    schedule 20.04.2021
comment
Большое спасибо, спасибо! - person apple; 20.04.2021
comment
Я обновил свой ответ - person Robert Hijmans; 20.04.2021
comment
Спасибо! У меня были проблемы с загрузкой terra / 1.2-1 на HPC, поэтому мне придется попробовать маршрут vrt. Функция vrt совпадает с функцией gdalbuildvrt из пакета gdalUtils? В противном случае я не могу найти функцию vrt. - person apple; 21.04.2021
comment
vrt был представлен в версии 1.1-17 (на CRAN) и основан на той же базовой функциональности GDAL, что и gdalUtils; но детали реализации могут отличаться. - person Robert Hijmans; 21.04.2021