R：使用“ rgdal”软件包裁剪GeoTiff栅格和“光栅” [英] R: Crop GeoTiff Raster using packages "rgdal" and "raster"

问题描述

我想使用提到的两个软件包 rgdal和 raster来裁剪GeoTiff栅格文件。一切工作正常，除了生成的输出tif的质量非常差且灰度而不是彩色之外。原始数据是来自瑞士联邦地形局的高质量栅格地图，可以下载示例文件此处。

I'd like to crop GeoTiff Raster Files using the two mentioned packages, "rgdal" and "raster". Everything works fine, except that the quality of the resulting output tif is very poor and in greyscale rather than colour. The original data are high quality raster maps from the swiss federal office of Topography, example files can be downloaded here.

这是我的代码：

## 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)

为了重现此内容例如，下载示例数据，并将其提取到文件夹 c：/ files /。奇怪的是，使用示例数据可以很好地裁剪图像，但是仍然是灰度。

In order to reproduce this example, download the sample data and extract it to the folder "c:/files/". Oddly enough, using the sample data the quality of the croped image is alright, but still greyscale.

我使用数据类型，格式选项玩耍，但是没有得到任何东西。有人可以指出解决方案吗？我应该提供输入数据的更多信息吗？

I played around using the options "datatype", "format", but didnt get anywhere with that. Can anybody point out a solution? Should I supply more information the the input data?

编辑：
Josh的示例与示例数据 2 。不幸的是，我的数据似乎比较旧，有些不同。您能否告诉我，如果您阅读以下GDALinfo，我会选择哪种选择：

Josh's example works superb with the sample data 2. Unfortunately, the data I have seems to be older and somewhat different. Can you tell me what option I choose if you read the following 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（）可能是最直接的解决方案：

If one simply wants to crop out a subset of an existing GeoTIFF and save the cropped part to a new *.tif file, using gdalUtils::gdal_translate() may be the most straightforward solution:

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)))

---

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

---

Looks like you need to change two details.

首先，您正在阅读的 *。tif 文件具有三个频段，因此应使用 stack（）。 （在其上使用 raster（）只能在单个波段（默认情况下为第一个波段）中读取，从而产生单色或灰度输出）。

First, the *.tif file you're reading in has three bands, so should be read in using stack(). (Using raster() on it will only read in a single band (the first one, by default) producing a monochromatic or 'greyscale' output).

第二（此处提到的原因） writeRaster（）默认情况下会将值写为实数（ Float64 我的机器）。为了明确地告诉您，您想使用字节，请给其参数 datatype ='INT1U'。

Second (for reasons mentioned here) writeRaster() will by default write out the values as real numbers (Float64 on my machine). To explicitly tell it you instead want to use bytes, give it the argument 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文件：

All of which outputs the following tiff file:

这篇关于R：使用“ rgdal”软件包裁剪GeoTiff栅格和“光栅”的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！