writing a raster stack on disk with writeRaster in R changes the range of the values for each layer

Question

I am processing a Landsat 8 scene in R to calculate the NDVI and to run a land cover classification algorithm. I am having problems with the writeRaster function from the raster package, in particular when it comes to write a raster stack on disk.

I start loading the 12 bands of the landsat 8 scene, and I stack them as layers in a raster stack. Since they are delivered as 16-bit images, the range of values for all layers varies from 0 to 65535. After writing the raster stack on the disk, when I reload the newly created file from the disk on the R environment, the value range for all layers is different from the original values. I can't figure out why, and I could not find any solution on the internet.

This is the 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)

When I check the characteristics of the rast.stack object, I can see that the ranges for the values of all bands are 0 -> 65535 Next, I write the raster stack on the disk:

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

When I load this new file in the R environment,

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

the value range for the bands is lower than in the original raster stack. I tried saving the file as .tif and in the original raster .grd format, but that has made no difference. I also tried to specify the data type with the datatype argument as follows:

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

I also tried to write the raster on the disk in chunks using the writeValues function. None of these has solved the problem. Does anyone know what is wrong, and how to fix this? In case you want to run this code, I am using the scene with path 201, row 037 recorded on November 14th 2016, freely available for download on EarthExplorer. Thanks

Answer 1

I'd try to run something like

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

on the original and reloaded stack.

This because I suspect that the differences you are reporting are not "real". In particular, it seems suspect to me that in the original raster all the bands cover the full 0 to 65535 range: that is difficult to be happening (would mean a complete "saturation" of the possible DN range in all bands).

What's in my opinion happening here is that the "ranges" of your "original" stack are only showing the "possible" min-max range given the data type. writeRaster is saving the computed statistics in the XML file associated to the tiff (see R: how to write a raster to disk without auxiliary file?), so that when you read it in again the "ranges" values are "correct".

HTH.