Range-map as covariate

Calculate expert range layer for each species, to use as a covariate in the model.

• R Libraries

library(patchwork)
library(knitr)
library(units)
library(smoothr)
library(tmap)
library(terra)
library(sf)
sf_use_s2(FALSE)
library(tidyverse)

• Basemaps

equalareaCRS <-  '+proj=laea +lon_0=-73.125 +lat_0=0 +datum=WGS84 +units=m +no_defs'
latam <- st_read('data/latam.gpkg', layer = 'latam', quiet = T)
countries <- st_read('data/latam.gpkg', layer = 'countries', quiet = T)
latam_land <- st_read('data/latam.gpkg', layer = 'latam_land', quiet = T)
latam_raster <- rast('data/latam_raster.tif', lyrs='latam')
latam_countries <- rast('data/latam_raster.tif', lyrs='countries')

• Datasets

data_carnivores_PO <- read_csv('data/data_PO.csv')
data_carnivores_PA <- read_csv('data/data_PA.csv')

Functions

• range_edge_dist_expert_range(): generates a raster with the distance from the edge of expert range map and calculates the values fro the presence-only grid cells and presence-absence blobs (for time1 and time2). The raster will have a value of 0 inside the range, and then a positive value outside as the distance to the range (rescaled from 0 to 1).
• rescale0to1(): rescales the distance to the edge of the expert range map to be a value from 0 to 1, using the formula (value-min)/(max-min). The original code by hannahlowens and their climateStability package.

range_edge_dist_expert_range <- function(expert_range_sf, time1_blobs, time2_blobs, 
                                         Latam.vect, Latam.rast){
  
  area_thresh <- units::set_units(100000, km^2)
  expert_range_sf_filled <- smoothr::fill_holes(expert_range_sf, area_thresh)
  
  dist_expertrange_rast <- rasterize(vect(expert_range_sf_filled), Latam.rast, background=2) %>% 
    mask(vect(Latam.vect)) %>% costDist(1) %>% classify(cbind(NaN, NA)) %>% rescale0to1()
  names(dist_expertrange_rast) <- 'dist_exprt'
  
  dist_expertrange_time1_blobs <- terra::extract(x = dist_expertrange_rast, y = vect(time1_blobs), fun = mean)[,2]
  dist_expertrange_time2_blobs <- terra::extract(x = dist_expertrange_rast, y = vect(time2_blobs), fun = mean)[,2]
  
  return(list(dist_expertrange_rast, dist_expertrange_time1_blobs, dist_expertrange_time2_blobs))
}

rescale0to1 <- function(rasterForCalculation){
  if (!is(rasterForCalculation, "SpatRaster")){
    warning("Supplied argument is not a SpatRaster./n", sep = "")
    return(NULL)
  }
  rasterMinMax <- minmax(rasterForCalculation)
  rescaledRaster <- (rasterForCalculation - rasterMinMax[1])/(rasterMinMax[2] - rasterMinMax[1])
  return(rescaledRaster)
}

Herpailurus yagouaroundi

hyagouaroundi_PA_time1_blobs <- readRDS('data/species_POPA_data/PA_hyagouaroundi_time1_blobs.rds')
hyagouaroundi_PA_time2_blobs <- readRDS('data/species_POPA_data/PA_hyagouaroundi_time2_blobs.rds')
hyagouaroundi_IUCN <- st_read('big_data/hyagouaroundi_IUCN.shp', quiet = T) %>% 
  st_transform(crs=equalareaCRS)


# Calculate raster and blobs
hyagouaroundi_range_edge <- range_edge_dist_expert_range(hyagouaroundi_IUCN, 
                                                         hyagouaroundi_PA_time1_blobs,
                                                         hyagouaroundi_PA_time2_blobs,
                                                         latam, latam_raster)

# grid cell values
data_hyagouaroundi_expert_gridcell <- data.frame(expertCell=hyagouaroundi_range_edge[[1]][])
# blobs in time 1
hyagouaroundi_expert_blob_time1 <- hyagouaroundi_range_edge[[2]]
# blobs in time 2
hyagouaroundi_expert_blob_time2 <- hyagouaroundi_range_edge[[3]]

Plot

tm_shape(hyagouaroundi_range_edge[[1]]) + 
    tm_raster(palette = 'GnBu', n=20, title = 'Distance to the range') +
    tm_shape(hyagouaroundi_IUCN) + 
    tm_borders() +
    tm_shape(latam) + 
    tm_borders(lwd = 0.2) +
    tm_layout(legend.position = c('left', 'bottom'))

Save

terra::writeRaster(hyagouaroundi_range_edge[[1]], 'data/species_POPA_data/hyagouaroundi_expert_gridcell.tif', overwrite=TRUE)
saveRDS(data_hyagouaroundi_expert_gridcell, 'data/species_POPA_data/hyagouaroundi_expert_gridcell.rds')
saveRDS(hyagouaroundi_expert_blob_time1, 'data/species_POPA_data/hyagouaroundi_expert_blob_time1.rds')
saveRDS(hyagouaroundi_expert_blob_time2, 'data/species_POPA_data/hyagouaroundi_expert_blob_time2.rds')

Leopardus pardalis

lpardalis_PA_time1_blobs <- readRDS('data/species_POPA_data/PA_lpardalis_time1_blobs.rds')
lpardalis_PA_time2_blobs <- readRDS('data/species_POPA_data/PA_lpardalis_time2_blobs.rds')
lpardalis_IUCN <- st_read('big_data/lpardalis_IUCN.shp', quiet = T) %>% 
  st_transform(crs=equalareaCRS)


# Calculate raster and blobs
lpardalis_range_edge <- range_edge_dist_expert_range(lpardalis_IUCN, 
                                                         lpardalis_PA_time1_blobs,
                                                         lpardalis_PA_time2_blobs,
                                                         latam, latam_raster)

# grid cell values
data_lpardalis_expert_gridcell <- data.frame(expertCell=lpardalis_range_edge[[1]][])
# blobs in time 1
lpardalis_expert_blob_time1 <- lpardalis_range_edge[[2]]
# blobs in time 2
lpardalis_expert_blob_time2 <- lpardalis_range_edge[[3]]

Plot

tm_shape(lpardalis_range_edge[[1]]) + 
    tm_raster(palette = 'GnBu', n=20, title = 'Distance to the range') +
    tm_shape(lpardalis_IUCN) + 
    tm_borders() +
    tm_shape(latam) + 
    tm_borders(lwd = 0.2) +
    tm_layout(legend.position = c('left', 'bottom'))

Save

terra::writeRaster(lpardalis_range_edge[[1]], 'data/species_POPA_data/lpardalis_expert_gridcell.tif', overwrite=TRUE)
saveRDS(data_lpardalis_expert_gridcell, 'data/species_POPA_data/lpardalis_expert_gridcell.rds')
saveRDS(lpardalis_expert_blob_time1, 'data/species_POPA_data/lpardalis_expert_blob_time1.rds')
saveRDS(lpardalis_expert_blob_time2, 'data/species_POPA_data/lpardalis_expert_blob_time2.rds')

Cerdocyon thous

cthous_PA_time1_blobs <- readRDS('data/species_POPA_data/PA_cthous_time1_blobs.rds')
cthous_PA_time2_blobs <- readRDS('data/species_POPA_data/PA_cthous_time2_blobs.rds')
cthous_IUCN <- st_read('big_data/cthous_IUCN.shp', quiet = T) %>% 
  st_transform(crs=equalareaCRS)


# Calculate raster and blobs
cthous_range_edge <- range_edge_dist_expert_range(cthous_IUCN, 
                                                         cthous_PA_time1_blobs,
                                                         cthous_PA_time2_blobs,
                                                         latam, latam_raster)

# grid cell values
data_cthous_expert_gridcell <- data.frame(expertCell=cthous_range_edge[[1]][])
# blobs in time 1
cthous_expert_blob_time1 <- cthous_range_edge[[2]]
# blobs in time 2
cthous_expert_blob_time2 <- cthous_range_edge[[3]]

Plot

tm_shape(cthous_range_edge[[1]]) + 
    tm_raster(palette = 'GnBu', n=20, title = 'Distance to the range') +
    tm_shape(cthous_IUCN) + 
    tm_borders() +
    tm_shape(latam) + 
    tm_borders(lwd = 0.2) +
    tm_layout(legend.position = c('left', 'bottom'))

Save

terra::writeRaster(cthous_range_edge[[1]], 'data/species_POPA_data/cthous_expert_gridcell.tif', overwrite=TRUE)
saveRDS(data_cthous_expert_gridcell, 'data/species_POPA_data/cthous_expert_gridcell.rds')
saveRDS(cthous_expert_blob_time1, 'data/species_POPA_data/cthous_expert_blob_time1.rds')
saveRDS(cthous_expert_blob_time2, 'data/species_POPA_data/cthous_expert_blob_time2.rds')

Chrysocyon brachyurus

cbrachyurus_PA_time1_blobs <- readRDS('data/species_POPA_data/PA_cbrachyurus_time1_blobs.rds')
cbrachyurus_PA_time2_blobs <- readRDS('data/species_POPA_data/PA_cbrachyurus_time2_blobs.rds')
cbrachyurus_IUCN <- st_read('big_data/cbrachyurus_IUCN.shp', quiet = T) %>% 
  st_transform(crs=equalareaCRS)


# Calculate raster and blobs
cbrachyurus_range_edge <- range_edge_dist_expert_range(cbrachyurus_IUCN, 
                                                         cbrachyurus_PA_time1_blobs,
                                                         cbrachyurus_PA_time2_blobs,
                                                         latam, latam_raster)

# grid cell values
data_cbrachyurus_expert_gridcell <- data.frame(expertCell=cbrachyurus_range_edge[[1]][])
# blobs in time 1
cbrachyurus_expert_blob_time1 <- cbrachyurus_range_edge[[2]]
# blobs in time 2
cbrachyurus_expert_blob_time2 <- cbrachyurus_range_edge[[3]]

Plot

tm_shape(cbrachyurus_range_edge[[1]]) + 
    tm_raster(palette = 'GnBu', n=20, title = 'Distance to the range') +
    tm_shape(cbrachyurus_IUCN) + 
    tm_borders() +
    tm_shape(latam) + 
    tm_borders(lwd = 0.2) +
    tm_layout(legend.position = c('left', 'bottom'))

Save

terra::writeRaster(cbrachyurus_range_edge[[1]], 'data/species_POPA_data/cbrachyurus_expert_gridcell.tif', overwrite=TRUE)
saveRDS(data_cbrachyurus_expert_gridcell, 'data/species_POPA_data/cbrachyurus_expert_gridcell.rds')
saveRDS(cbrachyurus_expert_blob_time1, 'data/species_POPA_data/cbrachyurus_expert_blob_time1.rds')
saveRDS(cbrachyurus_expert_blob_time2, 'data/species_POPA_data/cbrachyurus_expert_blob_time2.rds')

Eira barbara

ebarbara_PA_time1_blobs <- readRDS('data/species_POPA_data/PA_ebarbara_time1_blobs.rds')
ebarbara_PA_time2_blobs <- readRDS('data/species_POPA_data/PA_ebarbara_time2_blobs.rds')
ebarbara_IUCN <- st_read('big_data/ebarbara_IUCN.shp', quiet = T) %>% 
  st_transform(crs=equalareaCRS)


# Calculate raster and blobs
ebarbara_range_edge <- range_edge_dist_expert_range(ebarbara_IUCN, 
                                                         ebarbara_PA_time1_blobs,
                                                         ebarbara_PA_time2_blobs,
                                                         latam, latam_raster)

# grid cell values
data_ebarbara_expert_gridcell <- data.frame(expertCell=ebarbara_range_edge[[1]][])
# blobs in time 1
ebarbara_expert_blob_time1 <- ebarbara_range_edge[[2]]
# blobs in time 2
ebarbara_expert_blob_time2 <- ebarbara_range_edge[[3]]

Plot

tm_shape(ebarbara_range_edge[[1]]) + 
    tm_raster(palette = 'GnBu', n=20, title = 'Distance to the range') +
    tm_shape(ebarbara_IUCN) + 
    tm_borders() +
    tm_shape(latam) + 
    tm_borders(lwd = 0.2) +
    tm_layout(legend.position = c('left', 'bottom'))

Save

terra::writeRaster(ebarbara_range_edge[[1]], 'data/species_POPA_data/ebarbara_expert_gridcell.tif', overwrite=TRUE)
saveRDS(data_ebarbara_expert_gridcell, 'data/species_POPA_data/ebarbara_expert_gridcell.rds')
saveRDS(ebarbara_expert_blob_time1, 'data/species_POPA_data/ebarbara_expert_blob_time1.rds')
saveRDS(ebarbara_expert_blob_time2, 'data/species_POPA_data/ebarbara_expert_blob_time2.rds')

Leopardus wiedii

lwiedii_PA_time1_blobs <- readRDS('data/species_POPA_data/PA_lwiedii_time1_blobs.rds')
lwiedii_PA_time2_blobs <- readRDS('data/species_POPA_data/PA_lwiedii_time2_blobs.rds')
lwiedii_IUCN <- st_read('big_data/lwiedii_IUCN.shp', quiet = T) %>% 
  st_transform(crs=equalareaCRS)


# Calculate raster and blobs
lwiedii_range_edge <- range_edge_dist_expert_range(lwiedii_IUCN, 
                                                         lwiedii_PA_time1_blobs,
                                                         lwiedii_PA_time2_blobs,
                                                         latam, latam_raster)

# grid cell values
data_lwiedii_expert_gridcell <- data.frame(expertCell=lwiedii_range_edge[[1]][])
# blobs in time 1
lwiedii_expert_blob_time1 <- lwiedii_range_edge[[2]]
# blobs in time 2
lwiedii_expert_blob_time2 <- lwiedii_range_edge[[3]]

Plot

tm_shape(lwiedii_range_edge[[1]]) + 
    tm_raster(palette = 'GnBu', n=20, title = 'Distance to the range') +
    tm_shape(lwiedii_IUCN) + 
    tm_borders() +
    tm_shape(latam) + 
    tm_borders(lwd = 0.2) +
    tm_layout(legend.position = c('left', 'bottom'))

Save

terra::writeRaster(lwiedii_range_edge[[1]], 'data/species_POPA_data/lwiedii_expert_gridcell.tif', overwrite=TRUE)
saveRDS(data_lwiedii_expert_gridcell, 'data/species_POPA_data/lwiedii_expert_gridcell.rds')
saveRDS(lwiedii_expert_blob_time1, 'data/species_POPA_data/lwiedii_expert_blob_time1.rds')
saveRDS(lwiedii_expert_blob_time2, 'data/species_POPA_data/lwiedii_expert_blob_time2.rds')

Nasua nasua

nnasua_PA_time1_blobs <- readRDS('data/species_POPA_data/PA_nnasua_time1_blobs.rds')
nnasua_PA_time2_blobs <- readRDS('data/species_POPA_data/PA_nnasua_time2_blobs.rds')
nnasua_IUCN <- st_read('big_data/nnasua_IUCN.shp', quiet = T) %>% 
  st_transform(crs=equalareaCRS)


# Calculate raster and blobs
nnasua_range_edge <- range_edge_dist_expert_range(nnasua_IUCN, 
                                                         nnasua_PA_time1_blobs,
                                                         nnasua_PA_time2_blobs,
                                                         latam, latam_raster)

# grid cell values
data_nnasua_expert_gridcell <- data.frame(expertCell=nnasua_range_edge[[1]][])
# blobs in time 1
nnasua_expert_blob_time1 <- nnasua_range_edge[[2]]
# blobs in time 2
nnasua_expert_blob_time2 <- nnasua_range_edge[[3]]

Plot

tm_shape(nnasua_range_edge[[1]]) + 
    tm_raster(palette = 'GnBu', n=20, title = 'Distance to the range') +
    tm_shape(nnasua_IUCN) + 
    tm_borders() +
    tm_shape(latam) + 
    tm_borders(lwd = 0.2) +
    tm_layout(legend.position = c('left', 'bottom'))

Save

terra::writeRaster(nnasua_range_edge[[1]], 'data/species_POPA_data/nnasua_expert_gridcell.tif', overwrite=TRUE)
saveRDS(data_nnasua_expert_gridcell, 'data/species_POPA_data/nnasua_expert_gridcell.rds')
saveRDS(nnasua_expert_blob_time1, 'data/species_POPA_data/nnasua_expert_blob_time1.rds')
saveRDS(nnasua_expert_blob_time2, 'data/species_POPA_data/nnasua_expert_blob_time2.rds')

Pteronura brasiliensis

pbrasiliensis_PA_time1_blobs <- readRDS('data/species_POPA_data/PA_pbrasiliensis_time1_blobs.rds')
pbrasiliensis_PA_time2_blobs <- readRDS('data/species_POPA_data/PA_pbrasiliensis_time2_blobs.rds')
pbrasiliensis_IUCN <- st_read('big_data/pbrasiliensis_IUCN.shp', quiet = T) %>% 
  st_transform(crs=equalareaCRS)


# Calculate raster and blobs
pbrasiliensis_range_edge <- range_edge_dist_expert_range(pbrasiliensis_IUCN, 
                                                         pbrasiliensis_PA_time1_blobs,
                                                         pbrasiliensis_PA_time2_blobs,
                                                         latam, latam_raster)

# grid cell values
data_pbrasiliensis_expert_gridcell <- data.frame(expertCell=pbrasiliensis_range_edge[[1]][])
# blobs in time 1
pbrasiliensis_expert_blob_time1 <- pbrasiliensis_range_edge[[2]]
# blobs in time 2
pbrasiliensis_expert_blob_time2 <- pbrasiliensis_range_edge[[3]]

Plot

tm_shape(pbrasiliensis_range_edge[[1]]) + 
    tm_raster(palette = 'GnBu', n=20, title = 'Distance to the range') +
    tm_shape(pbrasiliensis_IUCN) + 
    tm_borders() +
    tm_shape(latam) + 
    tm_borders(lwd = 0.2) +
    tm_layout(legend.position = c('left', 'bottom'))

Save

terra::writeRaster(pbrasiliensis_range_edge[[1]], 'data/species_POPA_data/pbrasiliensis_expert_gridcell.tif', overwrite=TRUE)
saveRDS(data_pbrasiliensis_expert_gridcell, 'data/species_POPA_data/pbrasiliensis_expert_gridcell.rds')
saveRDS(pbrasiliensis_expert_blob_time1, 'data/species_POPA_data/pbrasiliensis_expert_blob_time1.rds')
saveRDS(pbrasiliensis_expert_blob_time2, 'data/species_POPA_data/pbrasiliensis_expert_blob_time2.rds')