Running nested user-defined function in `terra::app()` on multiple cores

Question

PREFACE My error is a simple "Mathematical operation applied to non-numeric argument" error, but I think this arise from how I create a suite of user defined functions and use them within the terra::app() function. I am going to describe the full workflow to elucidate what I am after, so please bare with me.

ISSUE I am trying to apply a statistical-empirical topographic correction to Sentinel 2A data in R. To apply the topographic correction, I have appended rasters of solar azimuth, solar zenith, and slope and aspect to the multiband scene. I am first taking a random sample of each band within a scene to get its intensity value as well as the corresponding solar zenith, aspect, and slope values. From there, I calculate the cosine of the solar incident angle for each sampled cell using the zenith, azimuth, slope and aspect using a user-define function. Then, I run a linear regression between the cosine of solar incident angle and the respective intensity value for each band. I then apply a user-defined function that calls the solar incident function above using the terra::app() function to finalize the topographic information from these linear regressions. This works on one core just fine with fake data, but is painfully slow with real Sentinel data, so I want it to work on multiple cores. When I try to run on multiple cores I get the error:

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

Reading the documentation in terra::app() I found that to export a function to multiple cores, I have to have a user-defined function in the fun= argument of terra::app(). I did this for the final function, but I suspect I am getting this error because I previously defined the cosine of solar incident angle function. I am not quite sure how to resolve this, and would greatly appreciate any advice

Below is my reproducible example with fake data:

##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

It dawned on me after a good night's sleep that the answer might have been much simpler than I realized. I just needed to run the cos_i() function using terra::app(), but everything else would work fine and quickly using standard raster algebra provided by the terra:: package. Therefore, I could eliminate the RAST_CORR() function and make the extra steps basic raster algebra.

##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