R:将带有颜色表的单个带状 rasterLayer 转换为 3 带状 RGB rasterStack

Question

类似于R: Crop GeoTiff Raster using packages "rgdal" and "raster"中提出的问题我正在尝试使用“rgdal”和“raster”包裁剪瑞士联邦地形局的 map ，同时保留原始颜色表。对于单个带状 *.tif 文件，裁剪后的图像会丢失颜色表信息，因此无法正确显示(生成的图像几乎是黑色的)。

输入文件可以下载here并应解压到文件夹“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)

前面提到的解决方案post仅适用于 3 频段 *.tif，不适用于 1 频段 *.tif(例如示例文件)。

一个可行的解决方案是将包含颜色表的单带rasterLayer转换为3带RGB rasterStack(如前面提到的post中的评论中所述)，这显然保留了颜色表。

但是，我不知道如何将单波段 *.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'))