Example of calculating landscae metrics in multiple buffers, adapted from the landscapemetrics package

metrics_multiple_buffers.R

# library(devtools)
# devtools::install_github("bniebuhr/landscapemetrics", ref = "multibuffer")
# devtools::install_github("bniebuhr/landscapetools", ref = "multibuffer")

# packages
library(dplyr)
library(readxl)
library(ggplot2)
library(forcats)

library(raster)
library(sf)
library(landscapetools)
library(landscapemetrics)
library(tmap)
library(tmaptools)

# rasters
forest_1985 <- raster::raster("data/mapbiomas_v06_1985_forest.tif")
forest_2020 <- raster::raster("data/mapbiomas_v06_2020_forest.tif")
# plot(forest_1985)

# CRS
crs_atlantic <- raster::crs(forest_1985)
# resolution of the map in meters
res_orig <- raster::res(forest_1985)[1]
res_m <- 30

# datasets
pol <- readxl::read_xlsx("data/Atlantic_forest_floral_visitor_data_points.xlsx") %>% 
  dplyr::rename(x = longitude, y = latitude) %>% 
  sf::st_as_sf(coords = c("x", "y"), crs = crs_atlantic)
fru <- readxl::read_xlsx("data/Atlantic_forest_frugivory_data_points.xlsx") %>% 
  dplyr::rename(x = longitude, y = latitude) %>% 
  sf::st_as_sf(coords = c("x", "y"), crs = crs_atlantic)

#' # Calculating changes in forest amount and isolation
#' 
#' To calculate the change in forest amount around each sampling point, we
#' followed these steps:
#' 
#' - Create a series of buffers, with extent from 250 m to 10 km .
#' - Calculate the change in absolute area of forest (in hectares) within each buffer, 
#' for each point, between 1985 and 2020.
#' - Calculate the change in proportion of forest within each buffer, for each point, 
#' between 1985 and 2020.
#' - Calculate the change in isolation, both absolute and proportional, 
#' between 1985 and 2020. Isolation was calculated as the mean distance between
#' nearest neighbor forest patches.
#' - Classify points where the changes in forest amount and isolation were 
#' the highest (> 3% of change).
#' 
#' # Pollination networks
#' 
#' ## Calculate metrics and changes

#--- label=pollination, eval=FALSE, echo=TRUE
# pollination

# extract amount of forest - 1985
buffers <- c(250, 750, 1000, 1500, 2000, 3000, 4000, 5000, 7500, 10000)
tictoc::tic()
pol_1985 <- landscapetools::util_extract_multibuffer(forest_1985, pol, buffer_width = buffers,
                                                     rel_freq = TRUE, point_id_text = FALSE) %>%
  dplyr::rename(freq1985 = freq, rel_freq1985 = rel_freq) %>% 
  tibble::as_tibble() %>% 
  dplyr::mutate(id = as.numeric(as.character(id)),
                layer = as.numeric(as.character(layer)))
tictoc::toc() # 17 min, 10 buffer sizes, 200 points

# calculate isolation - 1985
pol_1985_enn <- landscapemetrics::scale_sample(forest_1985, pol[1:2,], shape = "circle",
                                               size = (buffers*res_orig/res_m)[1:2],
                                               what = c("lsm_c_enn_mn"),
                                               verbose = TRUE,
                                               progress = TRUE) %>%
  dplyr::mutate(size = size/res_orig*res_m)
tictoc::tic()
landscape <- terra::rast(forest_1985)
progress <- TRUE

# check points that do not overlap - to be remvoed
y <- fru %>% terra::vect() %>% terra::buffer(width = buffers[1])
not_over <- list()
for(i in 1:length(y)) 
  not_over[[i]] <- try(terra::crop(landscape, y[i]))
to_remove <- grep("Error", not_over)

# these points do not overlap and must be checked
print(to_remove)

# loop over buffer sizes and points
pol_1985_enn <- do.call(rbind, lapply(X = seq_along(buffers), FUN = function(x) {
  
  # create buffers
  y <- pol[-152,] %>% terra::vect() %>% terra::buffer(width = buffers[x]) %>% 
    .[-to_remove]
  
  result_pts <- do.call(rbind, lapply(X = seq_along(y), FUN = function(current_plot) {
    
    if (progress) {
      
      cat("\r> Progress scales: ", sprintf("%02d", x), "/", length(buffers), 
          "; Progress points: ", sprintf("%03d", current_plot), "/", length(y))
    }
    
    # crop sample plot
    landscape_crop <- terra::crop(x = landscape,
                                  y = y[current_plot])
    
    # mask sample plot
    landscape_mask <- terra::mask(x = landscape_crop,
                                  mask = y[current_plot])
    
    # calculate lsm
    result_current_plot <- landscapemetrics::calculate_lsm(landscape = landscape_mask,
                                                           progress = FALSE, 
                                                           what = "lsm_c_enn_mn")
    # add buffer size
    result_current_plot$size <- buffers[x]
    
    # add plot id
    result_current_plot$plot_id <- current_plot
    
    return(result_current_plot)
  }))
  
  return(result_pts)
}))

pol_1985_enn <- pol_1985_enn %>% 
  dplyr::select(id = plot_id, layer = class, metric1985 = metric, 
                value1985 = value, buffer = size)
tictoc::toc()


# join
pol_1985 <- pol_1985 %>% 
  dplyr::left_join(pol_1985_enn,
                   by = c("id", "layer", "buffer"))