R: преобразование однополосного растрового слоя с таблицей цветов в трехканальный RGB rasterStack.

Подобно вопросу, поднятому в R: Обрезка GeoTiff Raster с использованием пакетов rgdal и растр Я пытаюсь обрезать карту Швейцарского федерального топографического управления с помощью пакетов "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 (как указано в комментарии в ранее упомянутом post), который, по-видимому, сохраняет таблицу цветов.

Однако я не знаю, как преобразовать одиночный канал * .tif в 3-полосный RGB rasterStack при сохранении таблицы цветов. Кто-нибудь знает, как можно сделать это преобразование, или у кого-нибудь есть идея получше решить проблему?


r crop raster rgdal geotiff
person Daniel Bernet    schedule 19.03.2015    source источник

Ответы (1)


arrow_upward
2
arrow_downward

Для этого можно использовать 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'))

введите описание изображения здесь

person jbaums    schedule 19.03.2015