在 R 中使用 writeRaster 在磁盘上写入栅格堆栈会更改每一层的值范围

Question

我正在 R 中处理 Landsat 8 场景以计算 NDVI 并运行土地覆盖分类算法。 我在使用 raster 包中的 writeRaster 函数时遇到问题，尤其是在磁盘上写入光栅堆栈时。

我开始加载 landsat 8 场景的 12 个波段，并将它们作为图层堆叠在栅格堆栈中。由于它们是作为 16 位图像交付的，因此所有层的值范围从 0 到 65535 不等。在磁盘上写入光栅堆栈后，当我在 R 环境中从磁盘重新加载新创建的文件时，值范围所有层都与原始值不同。我不知道为什么，我在互联网上找不到任何解决方案。

这是代码：

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)

当我检查 rast.stack 对象的特征时，我可以看到所有波段的值的范围是 0 -> 65535 接下来，我将光栅堆栈写入磁盘：

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

当我在 R 环境中加载这个新文件时，

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

波段的值范围低于原始栅格堆栈中的值。 我尝试将文件另存为 .tif 和原始光栅 .grd 格式，但这没有任何区别。我也尝试指定数据类型 datatype 参数如下：

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

我还尝试使用writeValues 函数将光栅分块写入磁盘。这些都没有解决问题。有谁知道出了什么问题，以及如何解决这个问题？如果您想运行此代码，我将使用路径为 201 的场景，第 037 行记录于 2016 年 11 月 14 日，可在 EarthExplorer 上免费下载。 谢谢

Answer 1

我会尝试运行类似的东西

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

在原始和重新加载的堆栈上。

这是因为我怀疑您报告的差异不是“真实的”。特别是，在我看来，在原始栅格中，所有波段都覆盖了 0 到 65535 的整个范围：这很难发生（这意味着所有波段中可能的 DN 范围完全“饱和”）。

我认为这里发生的情况是，您的“原始”堆栈的“范围”仅显示给定数据类型的“可能”最小-最大范围。 writeRaster 将计算的统计信息保存在与 tiff 关联的 XML 文件中（请参阅 R: how to write a raster to disk without auxiliary file?），这样当您再次读取它时，“范围”值是“正确的”。

HTH。