Created
February 23, 2022 11:04
-
-
Save bniebuhr/9eca21d87691f7d45878f9520054caa3 to your computer and use it in GitHub Desktop.
Example of calculating landscae metrics in multiple buffers, adapted from the landscapemetrics package
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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")) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment