R: Преобразование одного полосового слоя rasterLayer с таблицей цветов в 3-полосный RGB rasterStack

Question

Аналогично вопросу, поднятому в R: Crop GeoTiff Raster с использованием пакетов "rgdal" и "raster". Я пытаюсь обрезать карту из Федерального управления топографии Швейцарии с пакетами "rgdal" и "raster", сохраняя при этом исходную таблицу цветов., Для одного *.tif файла с полосами обрезанное изображение теряет информацию таблицы цветов и, таким образом, не отображается должным образом (получающееся изображение является почти черным).

Входной файл можно скачать здесь и должен быть распакован в папку "C:/files". Это код:

## install.packages("rgdal")
## install.packages("raster")
library("rgdal")
library("raster")
input <- "C:/files/PK25_KOMB_20L_2004_1_1.tif"
output <- "C:/files/cropped.tif"
r <- raster(input)
ex  <- extent(c(600500, 601500, 196500, 197500))
cropped <- crop(r, ex)
writeRaster(cropped, output, format="GTiff", datatype='INT1U', overwrite=TRUE)

Решение, представленное в вышеупомянутом посте, работало только для 3-полосного *.tif, но не для 1-полосного *.tif (такого как файл примера).

Решение, которое должно работать, заключается в преобразовании одного полосового rasterLayer, который включает в себя таблицу цветов, в трехполосный RGB rasterStack (как обрисовано в комментарии в ранее упомянутом посте), который, очевидно, сохраняет таблицу цветов.

Тем не менее, я не знаю, как сделать преобразование одной полосы *.tif в 3-полосную RGB rasterStack при сохранении таблицы цветов. Кто-нибудь знает, как это преобразование может быть сделано, или у кого-то есть лучшая идея, чтобы решить проблему?

Answer 1

Ты можешь использовать gdalUtils::gdalwarp за это:

library(raster)
library(gdalUtils)

Загрузка данных:

download.file(file.path('http://www.swisstopo.admin.ch/internet/swisstopo/de',
                        'home/products/maps/national/digital/national',
                        'pk_25.parsys.89625.downloadList.82162.DownloadFile.tmp',
                        'pk25komblzw.zip'), f <- tempfile())
unzip(f, exdir=tempdir())

Обрезка с gdalwarp :

cropped <- gdalwarp(
  file.path(tempdir(), 'PK25_KOMB_20L_2004_1_1.tif'), 
  'cropped.tif', te=c(600500, 196500, 601500, 197500), output_Raster = TRUE)

Обратите внимание, что степень должна быть задана как c(xmin, ymin, xmax, ymax), который отличается от порядка, используемого для raster::extent,

Подтверждая, что это сработало:

plot(raster('cropped.tif'))