运行 在多核上的 `terra::app()` 中嵌套了用户定义的函数
Running nested user-defined function in `terra::app()` on multiple cores
前言
我的错误是一个简单的“应用于非数字参数的数学运算”错误，但我认为这是由于我如何创建一组用户定义的函数并在 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
经过一夜好眠，我突然意识到答案可能比我意识到的要简单得多。我只需要使用 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
前言
我的错误是一个简单的“应用于非数字参数的数学运算”错误，但我认为这是由于我如何创建一组用户定义的函数并在 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
经过一夜好眠，我突然意识到答案可能比我意识到的要简单得多。我只需要使用 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