How to iteratively apply a function using a pair of data rasters, based on their names?

[nikos_geo]

I am having 2 sets of raster data and their names are:

1. ntl_'a number'.tif

2. pop_'a number'.tif

My goal is to create a function that reads the filenames (e.g., ntl_1.tif and pop_1.tif) and then perform the below analysis to raster with the same number in their layer names:

library(raster)
library(DescTools)

#create a data.frame of values from the NTL and pop raster data
ntl = raster("path/ntl_1.tif")
vals_ntl <- as.data.frame(values(ntl))
ntl_coords = as.data.frame(xyFromCell(ntl, 1:ncell(ntl)))
combine <- as.data.frame(cbind(ntl_coords,vals_ntl))

pop<-raster("path/pop_1.tif")
pop = resample(pop, ntl, method = 'bilinear')
vals_pop <- as.data.frame(values(pop))

block.data <- as.data.frame(cbind(combine, vals_pop))

names(block.data)[3] <- "ntl"
names(block.data)[4] <- "pop"

block.data <- na.omit(block.data)

block.data = subset(block.data, select = -c(x, y))

# sort by ntl
block.data <-block.data[order(block.data$ntl),]

ntl_vector <- block.data[ , "ntl"]
pop_vector <- block.data[ , "pop"]

#compute gini index
Gini(ntl_vector, pop_vector, unbiased = FALSE)

The above code is when using a pair of raster. In my case I have hundreds of pairs so I'd like to automate the process and hopefully to get the results (i.e., the gini coefficient) of every pair in my console or, even better, in a data.frame.
I made a lot of tries to create that function so I don't know which of my tries should I post. Here are some dummy data:

new("RasterStack", filename = "", layers = list(new("RasterLayer", 
    file = new(".RasterFile", name = "", datanotation = "FLT4S", 
        byteorder = "little", nodatavalue = -Inf, NAchanged = FALSE, 
        nbands = 1L, bandorder = "BIL", offset = 0L, toptobottom = TRUE, 
        blockrows = 0L, blockcols = 0L, driver = "", open = FALSE), 
    data = new(".SingleLayerData", values = 1:9, offset = 0, 
        gain = 1, inmemory = TRUE, fromdisk = FALSE, isfactor = FALSE, 
        attributes = list(), haveminmax = TRUE, min = 1L, max = 9L, 
        band = 1L, unit = "", names = "layer.1"), legend = new(".RasterLegend", 
        type = character(0), values = logical(0), color = logical(0), 
        names = logical(0), colortable = logical(0)), title = character(0), 
    extent = new("Extent", xmin = -180, xmax = 180, ymin = -90, 
        ymax = 90), rotated = FALSE, rotation = new(".Rotation", 
        geotrans = numeric(0), transfun = function () 
        NULL), ncols = 3L, nrows = 3L, crs = new("CRS", projargs = "+proj=longlat +datum=WGS84 +no_defs"), 
    history = list(), z = list()), new("RasterLayer", file = new(".RasterFile", 
    name = "", datanotation = "FLT4S", byteorder = "little", 
    nodatavalue = -Inf, NAchanged = FALSE, nbands = 1L, bandorder = "BIL", 
    offset = 0L, toptobottom = TRUE, blockrows = 0L, blockcols = 0L, 
    driver = "", open = FALSE), data = new(".SingleLayerData", 
    values = 1:9, offset = 0, gain = 1, inmemory = TRUE, fromdisk = FALSE, 
    isfactor = FALSE, attributes = list(), haveminmax = TRUE, 
    min = 1L, max = 9L, band = 1L, unit = "", names = "layer.2"), 
    legend = new(".RasterLegend", type = character(0), values = logical(0), 
        color = logical(0), names = logical(0), colortable = logical(0)), 
    title = character(0), extent = new("Extent", xmin = -180, 
        xmax = 180, ymin = -90, ymax = 90), rotated = FALSE, 
    rotation = new(".Rotation", geotrans = numeric(0), transfun = function () 
    NULL), ncols = 3L, nrows = 3L, crs = new("CRS", projargs = "+proj=longlat +datum=WGS84 +no_defs"), 
    history = list(), z = list())), title = character(0), extent = new("Extent", 
    xmin = -180, xmax = 180, ymin = -90, ymax = 90), rotated = FALSE, 
    rotation = new(".Rotation", geotrans = numeric(0), transfun = function () 
    NULL), ncols = 3L, nrows = 3L, crs = new("CRS", projargs = "+proj=longlat +datum=WGS84 +no_defs"), 
    history = list(), z = list())

[nikos_geo]

library(purrr)
library(fs)

raster_gini <- function(
    .ntl = "ntl_1.tif", 
    .pop = "pop_1.tif",
    .rdgal = TRUE
) {
  
  if(.rdgal) {
    
    ntl = raster(.ntl)
    vals_ntl <- as.data.frame(values(ntl))
    ntl_coords = as.data.frame(xyFromCell(ntl, 1:ncell(ntl)))
    combine <- as.data.frame(cbind(ntl_coords,vals_ntl))
    
    pop<-raster(.pop)
    pop = resample(pop, ntl, method = 'bilinear')
    vals_pop <- as.data.frame(values(pop))
    
    block.data <- as.data.frame(cbind(combine, vals_pop))
    
    #rename the columns
    names(block.data)[3] <- "ntl"
    names(block.data)[4] <- "pop"
    
    #remove NA values
    block.data <- na.omit(block.data)
    
    #remove the columns x & y
    block.data = subset(block.data, select = -c(x, y))
    
    # sort by ntl
    block.data <-block.data[order(block.data$ntl),]
    
    ntl_vector <- block.data[ , "ntl"]
    pop_vector <- block.data[ , "pop"]
    
    #compute gini index
    gini <- Gini(ntl_vector, pop_vector, unbiased = FALSE)
    
    c(ntl = .ntl, pop = .pop, gini = gini)
    
  } else {
    c(ntl = .ntl, pop = .pop)
  }
  
}


doc_paths_ntl <- fs::dir_ls("path_to_ntl_raster", glob = "*tif*") 

doc_paths_pop <- fs::dir_ls("path_to_pop_raster", glob = "*tif*") 


result_df <- purrr::map2_dfr(.x = doc_paths_ntl, .y = doc_paths_pop, .f = raster_gini)


result_df <- result_df |> 
  dplyr::mutate(ntl = basename(ntl)) |> 
  dplyr::mutate(pop = basename(pop))

result_df