-3.39999999999999996e+38 NODATA 值和栅格 R

编程问答 2022-05-10

问题描述

一些 geotiff 的

NODATA 值 (-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)

解决方法

混淆源于 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,y)
#values     : 0,max)

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

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

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

就您的 worldclim 数据而言，没有理由使用 na.rm=TRUE。否则，您可以在 mask 之后使用 calc。

terra 不是这种情况：

library(terra)
x <- rast(s)
sum(x,na.rm=TRUE)
#class       : SpatRaster 
#dimensions  : 115,1  (nrow,nlyr)
#resolution  : 40,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,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,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
,

光栅对象可以保存 NA 值。

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

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

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

geotiffrrr-raster