R: Convertir un rasterLayer à une seule bande avec une table de couleurs en un rasterStack RVB 3 bandes

Question

Similaire à la question soulevée dans R: Crop GeoTiff Raster using packages "rgdal" and "raster" J'essaie de recadrer une carte de l'Office fédéral de topographie avec les paquets "rgdal" et " raster "tout en préservant la table des couleurs d'origine. Pour un fichier * .tif à bande unique, l'image recadrée perd les informations de la table de couleurs et n'est donc pas affichée correctement (l'image résultante est presque noire).

Le fichier d'entrée peut être téléchargé here et doit être extrait dans le dossier "C:/files". Voici le code:

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

La solution présentée dans le post dont nous avons parlé seulement travaillé pour un * .tif 3 bandes, mais pas pour une 1 bande * .tif (tels que le fichier exemple).

Une solution qui devrait fonctionner est de convertir le rasterLayer simple bande qui inclut une table de couleurs dans un rasterStack RVB à 3 bandes (comme indiqué dans un commentaire dans le post précédemment mentionné) qui conserve apparemment la table de couleurs.

Cependant, je ne sais pas comment convertir une bande * .tif simple en une bande raster RGB à 3 bandes tout en préservant la table de couleurs. Est-ce que quelqu'un sait comment cette conversion peut être faite ou quelqu'un a-t-il une meilleure idée pour résoudre le problème?

Answer 1

Vous pouvez utiliser gdalUtils::gdalwarp pour cela:

library(raster) 
library(gdalUtils) 

Téléchargement des données:

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

Recadrage avec gdalwarp:

cropped <- gdalwarp(
    file.path(tempdir(), 'PK25_KOMB_20L_2004_1_1.tif'), 
    'cropped.tif', te=c(600500, 196500, 601500, 197500), output_Raster = TRUE) 

Notez que l'étendue doit être c(xmin, ymin, xmax, ymax), ce qui est différent de l'ordre utilisé pour raster::extent.

Confirmant que cela a fonctionné:

plot(raster('cropped.tif'))