R：使用包“rgdal”和“raster”裁剪 GeoTiff 栅格

我想使用提到的两个包“rgdal”和“raster”裁剪 GeoTiff 光栅文件。一切工作正常，除了生成的输出 tif 的质量非常差并且是灰度而不是彩色。原始数据是来自瑞士联邦地形局的高质量栅格地图，示例文件可以下载here http://www.swisstopo.admin.ch/internet/swisstopo/de/home/products/maps/national/digital/national/pk_25.html.

这是我的代码：

## install.packages("rgdal")
## install.packages("raster")
library("rgdal")
library("raster")

tobecroped <- raster("C:/files/krel_1129_2012_254dpi_LZW.tif")
ex  <- raster(xmn=648000, xmx=649000, ymn=224000, ymx=225000)
projection(ex) <- proj4string(tobecroped)
output <- "c:/files/output.tif"

crop(x = tobecroped, y = ex, filename = output)

为了重现此示例，请下载样本数据 http://www.swisstopo.admin.ch/internet/swisstopo/de/home/products/maps/national/digital/national/pk_25.parsys.60351.downloadList.20035.DownloadFile.tmp/krel11292012254dpilzw.zip并将其解压到文件夹“c:/files/”。奇怪的是，使用样本数据，裁剪后的图像质量还不错，但仍然是灰度的。

我尝试使用“数据类型”、“格式”选项，但没有得到任何结果。有人能指出解决方案吗？我应该提供输入数据的更多信息吗？

编辑： 乔希的示例对于示例数据效果非常好2 http://www.swisstopo.admin.ch/internet/swisstopo/de/home/products/maps/national/digital/national/pk_25.parsys.60351.downloadList.20035.DownloadFile.tmp/krel11292012254dpilzw.zip。不幸的是，我拥有的数据似乎更旧并且有些不同。如果您阅读以下 GDALinfo，您能告诉我我选择什么选项吗：

# packages same as above
OldInFile = "C:/files/krel1111.tif"
dataType(raster(OldInFile)
[1] "INT1U"

GDALinfo(OldInFile)

rows        4800 
columns     7000 
bands       1 
lower left origin.x        672500 
lower left origin.y        230000 
res.x       2.5 
res.y       2.5 
ysign       -1 
oblique.x   0 
oblique.y   0 
driver      GTiff 
projection  +proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333+k_0=1 +x_0=600000+y_0=200000 +ellps=bessel +units=m+no_defs 
file        C:/files/krel1111.tif 
apparent band summary:
  GDType hasNoDataValue NoDataValue blockSize1 blockSize2
1   Byte          FALSE           0          1       7000
apparent band statistics:
  Bmin Bmax Bmean Bsd
1    0  255    NA  NA
Metadata:
AREA_OR_POINT=Area 
TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch) 
TIFFTAG_XRESOLUTION=254 
TIFFTAG_YRESOLUTION=254 
Warning message:
statistics not supported by this driver

编辑（2015-03-10）：

如果只是想裁剪掉现有 GeoTIFF 的子集并将裁剪后的部分保存到新的*.tif文件，使用gdalUtils::gdal_translate()可能是最直接的解决方案：

library(raster)    # For extent(), xmin(), ymax(), et al.
library(gdalUtils) # For gdal_translate()

inFile <- "C:/files/krel_1129_2012_254dpi_LZW.tif"
outFile <- "subset.tif"
ex <- extent(c(686040.1, 689715.9, 238156.3, 241774.2))

gdal_translate(inFile, outFile,
               projwin=c(xmin(ex), ymax(ex), xmax(ex), ymin(ex)))

看来您需要更改两个细节。

首先，*.tif您正在读取的文件有三个波段，因此应该使用以下方式读取stack()。 （使用raster()它只会读取单个波段（默认情况下是第一个波段），产生单色或“灰度”输出。

Second (由于这里提到的原因 https://stat.ethz.ch/pipermail/r-sig-geo/2014-August/021624.html) writeRaster()默认情况下会将值写为实数（Float64在我的机器上）。要明确告诉它您想使用字节，请给它参数datatype='INT1U'.

library("rgdal")
library("raster")
inFile <- "C:/files/krel_1129_2012_254dpi_LZW.tif"
outFile <- "out.tif"

## Have a look at the format of your input file to:
## (1) Learn that it contains three bands (so should be read in as a RasterStack)
## (2) Contains values written as Bytes (so you should write output with datatype='INT1U')
GDALinfo(inFile)

## Read in as three separate layers (red, green, blue)
s <- stack(inFile)

## Crop the RasterStack to the desired extent
ex  <- raster(xmn=648000, xmx=649000, ymn=224000, ymx=225000)
projection(ex) <- proj4string(s)
s2 <- crop(s, ex)

## Write it out as a GTiff, using Bytes
writeRaster(s2, outFile, format="GTiff", datatype='INT1U', overwrite=TRUE)

所有这些都输出以下 tiff 文件：