包含栅格计算的函数在 R 中返回错误

问题描述

我正在尝试使用其中的栅格计算创建自己的函数。但是当我调用该函数时，它会返回一些错误。代码如下

BT <- function(B_1,B_2){
  tir_10 <- raster::raster(B_1)
  tir_11 <- raster::raster(B_2)
  
  l_lambda_10 <- 3.3420*10^-4*tir_10 + 0.1
  l_lambda_11 <- 3.3420*10^-4*tir_11 + 0.1
  
  BT_10 <- (1321.0789/(log(774.8853/l_lambda_10 + 1)))
  BT_11 <- (1201.1442/(log(480.8883/l_lambda_11 + 1)))
  
  return(list(BT_10 = BT_10,BT_11 = BT_11))
}

现在我正在尝试使用以下代码调用该函数

#Read the bands
library(raster)

a <- raster(ncol=100,nrow=100)
set.seed(2)     
values(a) = runif(10000,min=27791,max=30878)

b <- raster(ncol=100,nrow=100)
set.seed(2)     
values(b) = runif(10000,min=25686,max=28069)

BT(B_1 = a,B_2 = b)

但它返回以下错误

3.342 * 10^-4 * tir_10 中的错误：hasValues(e2) 不是 TRUE 如何解决？

解决方法

您必须在函数内部的输入中使用来自 RasterLayers B_1 和 B_2 的 min 和 max 来运行 runif，如下所示：

BT <- function(B_1,B_2){
  tir_10 <- raster::raster(B_1)
  set.seed(2)     
  values(tir_10) = runif(10000,min=round(min(getValues(B_1)),0),max=round(max(getValues(B_1)),0))
  
  tir_11 <- raster::raster(B_2)
  set.seed(2)     
  values(tir_11) = runif(10000,min=round(min(getValues(B_2)),max=round(max(getValues(B_2)),0))
  
  l_lambda_10 <- 3.3420*10^-4*tir_10 + 0.1
  l_lambda_11 <- 3.3420*10^-4*tir_11 + 0.1
  
  BT_10 <- (1321.0789/(log(774.8853/l_lambda_10 + 1)))
  BT_11 <- (1201.1442/(log(480.8883/l_lambda_11 + 1)))
  
  return(list(BT_10 = BT_10,BT_11 = BT_11))
}

这两行是错误的，因为 tir_10 和 tir_11 没有值。

tir_10 <- raster::raster(B_1)
tir_11 <- raster::raster(B_2)

你可以改为：

BT <- function(tir_10,tir_11){ 
  l_lambda_10 <- 3.3420*10^-4*tir_10 + 0.1
  l_lambda_11 <- 3.3420*10^-4*tir_11 + 0.1
  
  BT_10 <- (1321.0789/(log(774.8853/l_lambda_10 + 1)))
  BT_11 <- (1201.1442/(log(480.8883/l_lambda_11 + 1)))
  stack(list(BT_10=BT_10,BT_11=BT_11))
}

library(raster)
a <- raster(ncols=100,nrows=100)
set.seed(2)     
values(a) = runif(10000,min=27791,max=30878)
b <- setValues(a,runif(10000,min=25686,max=28069))

x <- BT(a,b)
x    
#class      : RasterStack 
#dimensions : 100,100,10000,2  (nrow,ncol,ncell,nlayers)
#resolution : 3.6,1.8  (x,y)
#extent     : -180,180,-90,90  (xmin,xmax,ymin,ymax)
#crs        : +proj=longlat +datum=WGS84 +no_defs 
#names      :    BT_10,BT_11 
#min values : 298.5254,297.9016 
#max values : 305.6353,304.4019

r r r-raster