为什么 terra 对于同一个栅格会产生不同的值范围？

Question

我一直在使用terra读取然后保存各种全球栅格数据集。今天，我注意到根据我如何看待数据或如何处理数据，值范围有很大不同。在我能解释为什么我会看到我所看到的东西之前，我对进行任何分析感到困惑和紧张。如果你们中有人能让我知道为什么我看到的值范围如此不同，我将不胜感激。

以全球道路密度为例。我使用 terra 导入 GRIP 总密度，所有类型组合 ascii 栅格，如下所示：

roadDir <- c("downloads/GRIP4_density_total/")

tmp <- paste0(roadDir, "grip4_total_dens_m_km2.asc")

roadDensity <- terra::rast(tmp)

如果我查看 roadDensity，我会发现值的范围是 0 到 99445（米每平方公里）：

> roadDensity
class       : SpatRaster 
dimensions  : 2160, 4320, 1  (nrow, ncol, nlyr)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 
source      : grip4_total_dens_m_km2.asc 
name        : grip4_total_dens_m_km2 
min value   :                      0 
max value   :                  99445 

如果我将原始数据导入 QGIS，我会得到与上面代码片段所示相同的像元大小、维度和值范围，因此我将这些值作为数据中的真实值。到目前为止一切顺利。

在以下情况下会出现问题：a) 我查看数据摘要或 b) 将生成的栅格写入 .tif 文件以用于进一步处理。

a) 总结

terra::global(roadDensity, c( "min", "max"), na.rm=TRUE)
                       min    max
grip4_total_dens_m_km2   0 330306

我不明白为什么我现在得到一个非常不同的最大值（330306）。

b) 写入.tif

x <- writeRaster(roadDensity,"test.tif", datatype="INT4U", overwrite=TRUE)

x
class       : SpatRaster 
dimensions  : 2160, 4320, 1  (nrow, ncol, nlyr)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : test.tif 
name        : grip4_total_dens_m_km2 
min value   :                      0 
max value   :                 330306 

所以 terra 的 global 和 writeRaster 产生相同的结果。

我知道，通过将数据类型值更改为其可能性之一（“INT1U”、“INT2U”、“INT4U”），writeRaster 会更改值范围，因为特定数据结构可以容纳的最大可能值。我不明白的是“INT4U”如何能够在单元格大小没有任何变化的情况下“增加”值范围如此显着地超出原始范围。此外，INT*U 值不应更改 global 生成的汇总值。 我确信这是一件非常简单的事情，我应该有一个“啊啊”的时刻，但我会感谢任何指导来帮助我的理解和分析信心。

谢谢你

Answer 1

获取数据：

#setwd(".")
url <- "https://dataportaal.pbl.nl/downloads/GRIP4/GRIP4_density_total.zip"
zip <- basename(url)
if (!file.exists(zip)) {
   download.file(url, zip, mode="wb")
   unzip(zip)
}

ascii 文件不会报告最小值和最大值（因为它们不以文件格式存储）

library(terra)
#terra 1.7.78

r <- rast("grip4_total_dens_m_km2.asc")
#class       : SpatRaster 
#dimensions  : 2160, 4320, 1  (nrow, ncol, nlyr)
#resolution  : 0.08333333, 0.08333333  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#source      : grip4_total_dens_m_km2.asc 
#name        : grip4_total_dens_m_km2 

创建一个新的 SpatRaster 以显示最小值和最大值：

r * 1
#class       : SpatRaster 
#dimensions  : 2160, 4320, 1  (nrow, ncol, nlyr)
#resolution  : 0.08333333, 0.08333333  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#source(s)   : memory
#varname     : grip4_total_dens_m_km2 
#name        : grip4_total_dens_m_km2 
#min value   :                      0 
#max value   :                 330306 

或者像这样计算它们：

minmax(r, compute=TRUE)
#    grip4_total_dens_m_km2
#min                      0
#max                 330306

也许您可以在空文件夹中再试一次。也许之前的工作留下了一个 sidecar 文件（例如，grip4_total_dens_m_km2.xml）来提供这些数字？另请报告您正在使用的 terra 版本。