在多核上运行 `terra::app()` 中的嵌套用户定义函数

Question

前言 我的错误是一个简单的“数学运算应用于非数字参数”错误，但我认为这是由于我创建一套用户定义函数并在 terra::app() 中使用它们而引起的功能。我将描述完整的工作流程来阐明我的目的，所以请跟我说清楚。

问题 我正在尝试对 R 中的 Sentinel 2A 数据应用统计经验地形校正。为了应用地形校正，我已将太阳方位角、太阳天顶以及坡度和方位角的栅格附加到多波段场景。我首先对场景中的每个波段进行随机采样，以获得其强度值以及相应的太阳天顶、方位和坡度值。从那里，我使用用户定义函数使用天顶角、方位角、斜率和方位角计算每个采样单元的太阳入射角的余弦。然后，我在太阳入射角的余弦与每个波段各自的强度值之间进行线性回归。然后，我使用 terra::app() 函数应用一个用户定义的函数调用上面的太阳入射函数，以最终确定这些线性回归的地形信息。这在一个核心上可以很好地处理假数据，但在处理真实的 Sentinel 数据时速度非常慢，所以我希望它可以在多个核心上工作。当我尝试在多核上运行时，出现错误：

Error: [app] cannot use this function
Error in cos(zen): non-numeric argument to mathematical function

阅读 terra::app() 中的文档我发现要将函数导出到多核，我必须有一个用户定义的函数在 terra::app() 的 fun= 参数中。我为最终函数执行了此操作，但我怀疑我收到此错误，因为我之前定义了太阳入射角函数的余弦。我不太确定如何解决这个问题，并且非常感谢任何建议，

下面是我使用虚假数据的可重现示例：

##Loading Necessary Packages##
library(terra)
library(tidyverse)

##For reproducibility##
set.seed(84)

##Creating a Fake multi-band raster##
RAS<-rast(nrows= 200, ncols=200, nlyrs = 7, ymin=45.1, ymax=45.2, xmin=-120.9, xmax=-120.8, crs="EPSG:2992")
b1<-runif(40000,0, 1000)
b2<-runif(40000,50, 2500)
b3<-runif(40000,1500, 3000)
slope<-runif(40000,0, 0.5*pi)
aspect<-runif(40000,0, 1.99*pi)
azimuth<-runif(40000,0, 1.99*pi)
zenith<-runif(40000,0, 1.99*pi)
values(RAS)<-c(b1,b2,b3,slope,aspect,azimuth,zenith)
names(RAS)<-c("band_1", "band_2", "band_3", "slope", "aspect", "azimuth", "zenith")

##Random Sample from raster bands##
SMP<-spatSample(RAS, size=999, xy=TRUE, as.df=TRUE, na.rm=TRUE)

##Function to calculate the cosine of the solar incident angle##
cos_i<-function(azm, zen, slope, aspect){
  slope_angle<-slope*(pi*0.25)
  out<-cos(slope_angle)*cos(zen)+sin(slope_angle)*sin(zen)*cos(azm - aspect) 
  return(out)
}

##Function to run linear regression on each band and output a dataframe of slopes and intercepts##
TOPO_lm<-function(df){
  df[,"X"]<-cos_i(azm=df$azimuth, zen = df$zenith, slope=df$slope, aspect=df$aspect)  
  models <- df %>% 
    pivot_longer(
      cols = c(3:5),
      names_to = "y_name",
      values_to = "y_value"
    ) %>% 
    split(.$y_name) %>% 
    map(~lm(y_value ~ X, data = .)) %>% 
    tibble(
      dvsub = names(.),
      untidied = .
    ) %>%
    mutate(tidy = map(untidied, broom::tidy)) %>%
    unnest(tidy) %>% 
    pivot_wider(id_cols="dvsub",
                names_from="term",
                values_from="estimate")
  out<-as.data.frame(models)
  colnames(out)<-c("band", "Beta_0", "Beta_1")
  return(out)
}

LM_DF<-TOPO_lm(SMP)

##Function to calculate mean intensity value from sampled data## 
L_bar_fxn<-function(df){
  df2<-df %>% summarize(across(.cols = c(3:5), mean)) %>% 
    pivot_longer(cols=everything(),
                 names_to="band",
                 values_to="intensity")
  out<-as.data.frame(df2)
  return(out)
}

MEAN_DF<-L_bar_fxn(SMP)

##Creating dataframe for topographic correction 
CORR_MTX<-merge(MEAN_DF, LM_DF, by = "band")

## Function to do the topographic correction ##
RAST_CORR<-function(df, SOLAR){
  Step1<- cos_i(azm=SOLAR[["azimuth"]], 
                zen=SOLAR[['zenith']], 
                slope=SOLAR[["slope"]], 
                aspect=SOLAR[["aspect"]])*df$Beta_1 - df$Beta_0+ df$intensity
  return(Step1)
  #out<- X - Step1
  #return(out)
}

##Applying the topographic correction to the intensity bands##
TEST<-app(RAS, function(i, ff, df) ff(i, df), ff=RAST_CORR, df=CORR_MTX, cores=5)#Throws error
TEST<-app(RAS, RAST_CORR, df=CORR_MTX)#Works
FINAL<-RAS[[1:3]]-TEST

Answer 1

经过一夜良好的睡眠后，我突然意识到答案可能比我想象的要简单得多。我只需要使用 terra::app() 运行 cos_i() 函数，但使用 提供的标准栅格代数，其他一切都可以正常快速地工作terra:: 包。因此，我可以消除 RAST_CORR() 函数并使额外的步骤成为基本的栅格代数。

##Loading Necessary Packages##
library(terra)
library(tidyverse)

##For reproducibility##
set.seed(84)

##Creating a Fake multi-band raster##
RAS<-rast(nrows= 200, ncols=200, nlyrs = 7, ymin=45.1, ymax=45.2, xmin=-120.9, xmax=-120.8, crs="EPSG:2992")
b1<-runif(40000,0, 1000)
b2<-runif(40000,50, 2500)
b3<-runif(40000,1500, 3000)
slope<-runif(40000,0, 0.5*pi)
aspect<-runif(40000,0, 1.99*pi)
azimuth<-runif(40000,0, 1.99*pi)
zenith<-runif(40000,0, 1.99*pi)
values(RAS)<-c(b1,b2,b3,slope,aspect,azimuth,zenith)
names(RAS)<-c("band_1", "band_2", "band_3", "slope", "aspect", "azimuth", "zenith")

##Random Sample from raster bands##
SMP<-spatSample(RAS, size=999, xy=TRUE, as.df=TRUE, na.rm=TRUE)

##Function to calculate the cosine of the solar incident angle##
cos_i<-function(azm, zen, slope, aspect){
  slope_angle<-slope*(pi*0.25)
  out<-cos(slope_angle)*cos(zen)+sin(slope_angle)*sin(zen)*cos(azm - aspect) 
  return(out)
}

##Function to run linear regression on each band and output a dataframe of slopes and intercepts##
TOPO_lm<-function(df){
  df[,"X"]<-cos_i(azm=df$azimuth, zen = df$zenith, slope=df$slope, aspect=df$aspect)  
  models <- df %>% 
    pivot_longer(
      cols = c(3:5),
      names_to = "y_name",
      values_to = "y_value"
    ) %>% 
    split(.$y_name) %>% 
    map(~lm(y_value ~ X, data = .)) %>% 
    tibble(
      dvsub = names(.),
      untidied = .
    ) %>%
    mutate(tidy = map(untidied, broom::tidy)) %>%
    unnest(tidy) %>% 
    pivot_wider(id_cols="dvsub",
                names_from="term",
                values_from="estimate")
  out<-as.data.frame(models)
  colnames(out)<-c("band", "Beta_0", "Beta_1")
  return(out)
}

LM_DF<-TOPO_lm(SMP)

##Function to calculate mean intensity value from sampled data## 
L_bar_fxn<-function(df){
  df2<-df %>% summarize(across(.cols = c(3:5), mean)) %>% 
    pivot_longer(cols=everything(),
                 names_to="band",
                 values_to="intensity")
  out<-as.data.frame(df2)
  return(out)
}

MEAN_DF<-L_bar_fxn(SMP)

##Creating dataframe for topographic correction 
CORR_MTX<-merge(MEAN_DF, LM_DF, by = "band")

##Adapting the cos_i() function for a raster##
cosI<-function(SOLAR){
  slope_angle<-SOLAR[["slope"]]*(pi*0.25)
  zen<-SOLAR[["zenith"]]
  azm<-SOLAR[["azimuth"]]
  aspect<-SOLAR[["aspect"]]
  out<-cos(slope_angle)*cos(zen)+sin(slope_angle)*sin(zen)*cos(azm - aspect) 
  return(out)
}

##Applying the topographic correction to the intensity bands##
TEST<-app(RAS, function(i, ff) ff(i), ff=cosI, cores=5)#Works now
FINAL<-RAS[[1:3]] - TEST*CORR_MTX$Beta_1-CORR_MTX$Beta_0 + CORR_MTX$intensity