PREFACE我的错误是一个简单的"应用于非数字参数的数学运算"；错误，但我认为这源于我如何创建一套用户定义的函数并在terra::app()函数中使用它们。我将描述完整的工作流程，以阐明我想要什么，所以请与我一起展示。
ISSUE我正在尝试对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
睡了一夜之后，我突然意识到答案可能比我意识到的要简单得多。我只需要使用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