R: Обрезать GeoTiff Raster, используя пакеты "rgdal" и "raster"
Я бы хотел обрезать растровые файлы GeoTiff, используя два упомянутых пакета: "rgdal" и "raster". Все работает отлично, за исключением того, что качество полученного результата очень низкое и в оттенках серого, а не в цвете. Исходные данные представляют собой высококачественные растровые карты из швейцарского федерального офиса Topography, примеры файлов можно скачать здесь.
Это мой код:
## 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/". Как ни странно, с использованием данных примера качество обрезанного изображения в порядке, но все же в оттенках серого.
Я играл с опциями "тип данных", "формат", но с этим ничего не получалось. Кто-нибудь может указать на решение? Должен ли я предоставить больше информации входных данных?
РЕДАКТИРОВАТЬ: пример Джош отлично работает с образцами данных 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
1 ответ
Изменить (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() он будет считываться только в одной полосе (первая по умолчанию), получая монохроматический или "полутоновый" вывод).
Второе ( по причинам, указанным здесь) 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:
