écrire une pile raster sur le disque avec writeRaster dans R modifie la plage des valeurs de chaque couche

Question

Je suis en train de traiter une scène Landsat 8 dans R pour calculer le NDVI et exécuter un algorithme de classification de la couverture terrestre. J'ai des problèmes avec la fonction writeRaster du paquet raster, en particulier quand il s'agit d'écrire une pile raster sur le disque.

Je commence à charger les 12 bandes de la scène Landsat 8, et je les empile en tant que couches dans une pile raster. Comme elles sont livrées en tant qu'images 16 bits, la plage de valeurs de toutes les couches varie de 0 à 65535. Après avoir écrit la pile raster sur le disque, lorsque je recharge le fichier nouvellement créé à partir du disque sur l'environnement R, la plage de valeurs pour toutes les couches est différente des valeurs d'origine. Je n'arrive pas à comprendre pourquoi et je n'ai trouvé aucune solution sur internet.

Voici le code:

library(raster) 
# Load the individual bands of the Landsat scene. 
b01 <- raster(list.files(dirname, pattern = "B1.TIF", full.names = TRUE)) 
b02 <- raster(list.files(dirname, pattern = "B2.TIF", full.names = TRUE)) 
b03 <- raster(list.files(dirname, pattern = "B3.TIF", full.names = TRUE)) 
b04 <- raster(list.files(dirname, pattern = "B4.TIF", full.names = TRUE)) 
b05 <- raster(list.files(dirname, pattern = "B5.TIF", full.names = TRUE)) 
b06 <- raster(list.files(dirname, pattern = "B6.TIF", full.names = TRUE)) 
b07 <- raster(list.files(dirname, pattern = "B7.TIF", full.names = TRUE)) 
b08 <- raster(list.files(dirname, pattern = "B8.TIF", full.names = TRUE)) 
b09 <- raster(list.files(dirname, pattern = "B9.TIF", full.names = TRUE)) 
b10 <- raster(list.files(dirname, pattern = "B10.TIF", full.names = TRUE)) 
b11 <- raster(list.files(dirname, pattern = "B11.TIF", full.names = TRUE)) 
b12 <- raster(list.files(dirname, pattern = "BQA.TIF", full.names = TRUE)) 

# Since the band 8 has a 15m resolution, compared to 30m of all other bands, I 
# need to resample it to match the other bands. 
b08 <- resample(b08, b01) 

allbands <- c(b01, b02, b03, b04, b05, b06, b07, b08, b09, b10, b11, b12) 
rast.stack <- stack(allbands) 

Quand je vérifie les caractéristiques de l'objet rast.stack, je peux voir que les plages pour les valeurs de toutes les bandes sont 0 -> 65535 Ensuite, je vous écris la pile de trame sur le disque:

writeRaster(rast.stack, filename = "LT820103720161114.tif", overwrite = TRUE) 

Quand je charge ce nouveau fichier dans l'environnement R,

rast.stack <- stack("LT820103720161114.tif") 

La plage de valeurs pour les bandes est inférieure à celle de la pile raster d'origine. J'ai essayé d'enregistrer le fichier au format .tif et dans le format raster .grd d'origine, mais cela n'a fait aucune différence. J'ai essayé aussi de spécifier le type de données avec l'argument datatype comme suit:

writeRaster(rast.stack, filename = "LT820103720161114.tif", datatype = "INT2U", overwrite = TRUE) 

J'ai essayé aussi d'écrire la trame sur le disque en morceaux à l'aide de la fonction writeValues. Aucun d'entre eux a résolu le problème. Est-ce que quelqu'un sait ce qui ne va pas, et comment résoudre ce problème? Dans le cas où vous voulez exécuter ce code, j'utilise la scène avec le chemin 201, ligne 037 enregistrée le 14 novembre 2016, disponible gratuitement en téléchargement sur EarthExplorer. Merci

Answer 1

Je vais essayer d'exécuter quelque chose comme

cellStats(rast.stack, mean) cellStats(rast.stack, min) cellStats(rast.stack, max)

sur la pile d'origine et rechargées. Ceci parce que je soupçonne que les différences que vous signalez ne sont pas "réelles". En particulier, il me semble suspect que dans le raster d'origine toutes les bandes couvrent toute la gamme de 0 à 65535: cela est difficile à réaliser (signifierait une "saturation" complète de la gamme de DN possible dans toutes les bandes). À mon avis, ce qui se passe ici, c'est que les «gammes» de votre pile «originale» ne montrent que la plage «possible» min-max étant donné le type de données. writeRaster enregistre les statistiques calculées dans le fichier XML associé au tiff (voir R: how to write a raster to disk without auxiliary file?), de sorte que lorsque vous le relisez, les valeurs "ranges" sont "correctes".

HTH.