4

上記の 2 つのパッケージ "rgdal" と "raster" を使用して、GeoTiff ラスター ファイルをトリミングしたいと思います。結果として得られる出力 tif の品質が非常に悪く、カラーではなくグレースケールであることを除いて、すべてが正常に機能します。元のデータは、スイス連邦地形局の高品質ラスター マップです。サンプル ファイルは、ここからダウンロードできます。

これは私のコードです：

## 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/」フォルダに展開します。奇妙なことに、サンプル データを使用すると、トリミングされた画像の品質は問題ありませんが、グレースケールのままです。

オプション「datatype」、「format」を使用して遊んでいましたが、それでどこにも行きませんでした。誰でも解決策を指摘できますか？入力データにさらに情報を提供する必要がありますか?

編集:ジョシュの例は、サンプルデータ2でうまく機能します。残念ながら、私が持っているデータは古く、多少異なっているようです。次の 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
4

1 に答える 1

7

編集 (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)))

2 つの詳細を変更する必要があるようです。

まず、*.tif読み込んでいるファイルには 3 つのバンドがあるため、 を使用して読み込む必要がありますstack()。( raster()on を使用すると、単色または「グレースケール」出力を生成する単一のバンド (デフォルトでは最初のバンド) でのみ読み取られます)。

2番目（ここで述べた理由により）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 ファイルを出力します。

ここに画像の説明を入力

于 2015-03-07T18:29:11.997 に答える