-3.39999999999999996e+38 NODATA 值和栅格 R答案

【问题标题】：-3.39999999999999996e+38 NODATA value and raster R-3.39999999999999996e+38 NODATA 值和栅格 R
【发布时间】：2021-04-26 14:42:27
【问题描述】：

NODATA 一些 geotiff 的值 (-3.39999999999999996e+38) 在 R 中被识别为 -Inf。

我使用光栅包读取文件并生成光栅堆栈。

在使用 calc 进行栅格计算后，NODATA 值将转换为 0。

有什么办法可以保留NODATA标志值。

我认为这是一个浮点精度限制。

数据来自https://www.worldclim.org/data/v1.4/cmip5_30s.html

这是我的 R 代码。

library(raster) 
library(rasterVis) 

tiffs <- list.files(".", sprintf("%s.*tif","clip" ), full.names = T)

s <-stack(tiffs)
NAvalue(s)
levelplot(s,margin=FALSE)


composite <- raster::calc(s,sum, na.rm=T)

levelplot(composite,margin=FALSE)

【问题讨论】：

这个问题可能更适合 gis.stackexchange.com .. 但您是否只需将 0 重新分类为 NODATA 还是希望从一开始就保留这些？
不，0 不同于 NODATA 值 (-3.39999999999999996e+38)。我想将其保留为 0，将 NODATA 保留为 -3.39999999999999996e+38
这可能是-3.39e308（不是-3.39e38），它会超出R的浮点最大负值（.Machine$double.xmin）吗？正如@RobertHijmans 所说，我们真的需要minimal reproducible example，然后才能提供很大帮助……这些值是在导入数据后在进行计算之后转换为0，还是仅在计算之后转换为0？你在做什么计算？（例如，1/Inf 或 -1/Inf 在 R 中的结果为零，但转换为 NA（R 的 NODATA 模拟）应该保留缺失...
相关：stackoverflow.com/questions/33446590/…

标签： r r-raster geotiff

【解决方案1】：

混淆源于sum 的工作方式。我将使用 R 附带的示例文件以便于重现

library(raster)
f <- system.file("external/test.grd", package="raster")
r <- raster(f)
s <- stack(r, r*2)

## equivalent to calc(s, sum)
sum(s, na.rm=FALSE)
#class      : RasterLayer 
#dimensions : 115, 80, 9200  (nrow, ncol, ncell)
#resolution : 40, 40  (x, y)
#values     : 416.1212, 5208.174  (min, max)

sum(s, na.rm=TRUE)
#class      : RasterLayer 
#dimensions : 115, 80, 9200  (nrow, ncol, ncell)
#resolution : 40, 40  (x, y)
#values     : 0, 5208.174  (min, max)

注意最小值的差异（0 与 416.1212）。这是因为sum 的工作原理。

sum(c(NA, NA))
#[1] NA
sum(c(NA, NA), na.rm=TRUE)
#[1] 0

无之和为零，而一个或多个NAs 之和为NA。

在您使用 worldclim 数据的情况下，没有理由使用 na.rm=TRUE。否则你可以在calc之后使用mask。

terra 不是这种情况：

library(terra)
x <- rast(s)
sum(x, na.rm=TRUE)
#class       : SpatRaster 
#dimensions  : 115, 80, 1  (nrow, ncol, nlyr)
#resolution  : 40, 40  (x, y)
#extent      : 178400, 181600, 329400, 334000  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +datum=WGS84 +units=m +no_defs 
#source      : memory 
#name        :      sum 
#min value   : 416.1212 
#max value   : 5208.174 

app(x, sum, na.rm=TRUE)
#class       : SpatRaster 
#dimensions  : 115, 80, 1  (nrow, ncol, nlyr)
#resolution  : 40, 40  (x, y)
#extent      : 178400, 181600, 329400, 334000  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +datum=WGS84 +units=m +no_defs 
#source      : memory 
#name        :      sum 
#min value   : 416.1212 
#max value   : 5208.174

甚至认为这与基本 R 有点不一致，但我认为在栅格数据的上下文中这样做会更好。获得与 base-R 和 raster 相同的行为

app(x, function(i) sum(i, na.rm=TRUE))
#class       : SpatRaster 
#dimensions  : 115, 80, 1  (nrow, ncol, nlyr)
#resolution  : 40, 40  (x, y)
#extent      : 178400, 181600, 329400, 334000  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +datum=WGS84 +units=m +no_defs 
#source      : memory 
#name        :    lyr.1 
#min value   :        0 
#max value   : 5208.174

【讨论】：

虽然有用，但它更像是评论而不是答案......
您使用模式"clip.*tif" --- 您是如何生成这些文件的？（worldclim 网站没有类似名称的文件）
我用 QGIS 剪辑了它们。你可以从这里下载剪辑github.com/kokkytos/demo
原来如此简单...谢谢

【解决方案2】：

光栅对象可以保存NA 值。

假设对象被称为x。然后你可以例如将低于阈值的所有值设置为NA。

x[x < -3.39999999999999996e+37] <- NA

【讨论】：

或x[ !is.finite(x)] <- NA ?
@BenBolker 这也会将被视为Inf 的大正值设置为NA。

【解决方案3】：

以 e+38 结尾的数据是 geotiff 文件中的空点

【讨论】：