2

与R: Crop GeoTiff Raster using packages "rgdal" and "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 转换为 3 带 RGB rasterStack（如前面提到的帖子中的评论中所述），这显然可以保存颜色表。

但是，我不知道如何在保留颜色表的同时将单波段 *.tif 转换为 3 波段 RGB rasterStack。有谁知道如何进行这种转换，或者有没有人有更好的想法来解决这个问题？

4

1 回答 1

2

您可以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'))

在此处输入图像描述

于 2015-03-19T23:40:19.867 回答