library(patchwork)
library(knitr)
library(units)
library(smoothr)
library(tmap)
library(terra)
library(sf)
sf_use_s2(FALSE)
library(tidyverse)
Range-map as covariate
Calculate expert range layer for each species, to use as a covariate in the model.
- R Libraries
- Basemaps
<- '+proj=laea +lon_0=-73.125 +lat_0=0 +datum=WGS84 +units=m +no_defs'
equalareaCRS <- st_read('data/latam.gpkg', layer = 'latam', quiet = T)
latam <- st_read('data/latam.gpkg', layer = 'countries', quiet = T)
countries <- st_read('data/latam.gpkg', layer = 'latam_land', quiet = T)
latam_land <- rast('data/latam_raster.tif', lyrs='latam')
latam_raster <- rast('data/latam_raster.tif', lyrs='countries') latam_countries
- Datasets
<- read_csv('data/data_PO.csv')
data_carnivores_PO <- read_csv('data/data_PA.csv') data_carnivores_PA
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 theirclimateStability
package.
<- function(expert_range_sf, time1_blobs, time2_blobs,
range_edge_dist_expert_range
Latam.vect, Latam.rast){
<- units::set_units(100000, km^2)
area_thresh <- smoothr::fill_holes(expert_range_sf, area_thresh)
expert_range_sf_filled
<- rasterize(vect(expert_range_sf_filled), Latam.rast, background=2) %>%
dist_expertrange_rast mask(vect(Latam.vect)) %>% costDist(1) %>% classify(cbind(NaN, NA)) %>% rescale0to1()
names(dist_expertrange_rast) <- 'dist_exprt'
<- terra::extract(x = dist_expertrange_rast, y = vect(time1_blobs), fun = mean)[,2]
dist_expertrange_time1_blobs <- terra::extract(x = dist_expertrange_rast, y = vect(time2_blobs), fun = mean)[,2]
dist_expertrange_time2_blobs
return(list(dist_expertrange_rast, dist_expertrange_time1_blobs, dist_expertrange_time2_blobs))
}
<- function(rasterForCalculation){
rescale0to1 if (!is(rasterForCalculation, "SpatRaster")){
warning("Supplied argument is not a SpatRaster./n", sep = "")
return(NULL)
}<- minmax(rasterForCalculation)
rasterMinMax <- (rasterForCalculation - rasterMinMax[1])/(rasterMinMax[2] - rasterMinMax[1])
rescaledRaster return(rescaledRaster)
}
Herpailurus yagouaroundi
<- readRDS('data/species_POPA_data/PA_hyagouaroundi_time1_blobs.rds')
hyagouaroundi_PA_time1_blobs <- readRDS('data/species_POPA_data/PA_hyagouaroundi_time2_blobs.rds')
hyagouaroundi_PA_time2_blobs <- st_read('big_data/hyagouaroundi_IUCN.shp', quiet = T) %>%
hyagouaroundi_IUCN st_transform(crs=equalareaCRS)
# Calculate raster and blobs
<- range_edge_dist_expert_range(hyagouaroundi_IUCN,
hyagouaroundi_range_edge
hyagouaroundi_PA_time1_blobs,
hyagouaroundi_PA_time2_blobs,
latam, latam_raster)
# grid cell values
<- data.frame(expertCell=hyagouaroundi_range_edge[[1]][])
data_hyagouaroundi_expert_gridcell # blobs in time 1
<- hyagouaroundi_range_edge[[2]]
hyagouaroundi_expert_blob_time1 # blobs in time 2
<- hyagouaroundi_range_edge[[3]] hyagouaroundi_expert_blob_time2
Plot
Code
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
Code
::writeRaster(hyagouaroundi_range_edge[[1]], 'data/species_POPA_data/hyagouaroundi_expert_gridcell.tif', overwrite=TRUE)
terrasaveRDS(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
<- readRDS('data/species_POPA_data/PA_lpardalis_time1_blobs.rds')
lpardalis_PA_time1_blobs <- readRDS('data/species_POPA_data/PA_lpardalis_time2_blobs.rds')
lpardalis_PA_time2_blobs <- st_read('big_data/lpardalis_IUCN.shp', quiet = T) %>%
lpardalis_IUCN st_transform(crs=equalareaCRS)
# Calculate raster and blobs
<- range_edge_dist_expert_range(lpardalis_IUCN,
lpardalis_range_edge
lpardalis_PA_time1_blobs,
lpardalis_PA_time2_blobs,
latam, latam_raster)
# grid cell values
<- data.frame(expertCell=lpardalis_range_edge[[1]][])
data_lpardalis_expert_gridcell # blobs in time 1
<- lpardalis_range_edge[[2]]
lpardalis_expert_blob_time1 # blobs in time 2
<- lpardalis_range_edge[[3]] lpardalis_expert_blob_time2
Plot
Code
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
Code
::writeRaster(lpardalis_range_edge[[1]], 'data/species_POPA_data/lpardalis_expert_gridcell.tif', overwrite=TRUE)
terrasaveRDS(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
<- readRDS('data/species_POPA_data/PA_cthous_time1_blobs.rds')
cthous_PA_time1_blobs <- readRDS('data/species_POPA_data/PA_cthous_time2_blobs.rds')
cthous_PA_time2_blobs <- st_read('big_data/cthous_IUCN.shp', quiet = T) %>%
cthous_IUCN st_transform(crs=equalareaCRS)
# Calculate raster and blobs
<- range_edge_dist_expert_range(cthous_IUCN,
cthous_range_edge
cthous_PA_time1_blobs,
cthous_PA_time2_blobs,
latam, latam_raster)
# grid cell values
<- data.frame(expertCell=cthous_range_edge[[1]][])
data_cthous_expert_gridcell # blobs in time 1
<- cthous_range_edge[[2]]
cthous_expert_blob_time1 # blobs in time 2
<- cthous_range_edge[[3]] cthous_expert_blob_time2
Plot
Code
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
Code
::writeRaster(cthous_range_edge[[1]], 'data/species_POPA_data/cthous_expert_gridcell.tif', overwrite=TRUE)
terrasaveRDS(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
<- readRDS('data/species_POPA_data/PA_cbrachyurus_time1_blobs.rds')
cbrachyurus_PA_time1_blobs <- readRDS('data/species_POPA_data/PA_cbrachyurus_time2_blobs.rds')
cbrachyurus_PA_time2_blobs <- st_read('big_data/cbrachyurus_IUCN.shp', quiet = T) %>%
cbrachyurus_IUCN st_transform(crs=equalareaCRS)
# Calculate raster and blobs
<- range_edge_dist_expert_range(cbrachyurus_IUCN,
cbrachyurus_range_edge
cbrachyurus_PA_time1_blobs,
cbrachyurus_PA_time2_blobs,
latam, latam_raster)
# grid cell values
<- data.frame(expertCell=cbrachyurus_range_edge[[1]][])
data_cbrachyurus_expert_gridcell # blobs in time 1
<- cbrachyurus_range_edge[[2]]
cbrachyurus_expert_blob_time1 # blobs in time 2
<- cbrachyurus_range_edge[[3]] cbrachyurus_expert_blob_time2
Plot
Code
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
Code
::writeRaster(cbrachyurus_range_edge[[1]], 'data/species_POPA_data/cbrachyurus_expert_gridcell.tif', overwrite=TRUE)
terrasaveRDS(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
<- readRDS('data/species_POPA_data/PA_ebarbara_time1_blobs.rds')
ebarbara_PA_time1_blobs <- readRDS('data/species_POPA_data/PA_ebarbara_time2_blobs.rds')
ebarbara_PA_time2_blobs <- st_read('big_data/ebarbara_IUCN.shp', quiet = T) %>%
ebarbara_IUCN st_transform(crs=equalareaCRS)
# Calculate raster and blobs
<- range_edge_dist_expert_range(ebarbara_IUCN,
ebarbara_range_edge
ebarbara_PA_time1_blobs,
ebarbara_PA_time2_blobs,
latam, latam_raster)
# grid cell values
<- data.frame(expertCell=ebarbara_range_edge[[1]][])
data_ebarbara_expert_gridcell # blobs in time 1
<- ebarbara_range_edge[[2]]
ebarbara_expert_blob_time1 # blobs in time 2
<- ebarbara_range_edge[[3]] ebarbara_expert_blob_time2
Plot
Code
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
Code
::writeRaster(ebarbara_range_edge[[1]], 'data/species_POPA_data/ebarbara_expert_gridcell.tif', overwrite=TRUE)
terrasaveRDS(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
<- readRDS('data/species_POPA_data/PA_lwiedii_time1_blobs.rds')
lwiedii_PA_time1_blobs <- readRDS('data/species_POPA_data/PA_lwiedii_time2_blobs.rds')
lwiedii_PA_time2_blobs <- st_read('big_data/lwiedii_IUCN.shp', quiet = T) %>%
lwiedii_IUCN st_transform(crs=equalareaCRS)
# Calculate raster and blobs
<- range_edge_dist_expert_range(lwiedii_IUCN,
lwiedii_range_edge
lwiedii_PA_time1_blobs,
lwiedii_PA_time2_blobs,
latam, latam_raster)
# grid cell values
<- data.frame(expertCell=lwiedii_range_edge[[1]][])
data_lwiedii_expert_gridcell # blobs in time 1
<- lwiedii_range_edge[[2]]
lwiedii_expert_blob_time1 # blobs in time 2
<- lwiedii_range_edge[[3]] lwiedii_expert_blob_time2
Plot
Code
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
Code
::writeRaster(lwiedii_range_edge[[1]], 'data/species_POPA_data/lwiedii_expert_gridcell.tif', overwrite=TRUE)
terrasaveRDS(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
<- readRDS('data/species_POPA_data/PA_nnasua_time1_blobs.rds')
nnasua_PA_time1_blobs <- readRDS('data/species_POPA_data/PA_nnasua_time2_blobs.rds')
nnasua_PA_time2_blobs <- st_read('big_data/nnasua_IUCN.shp', quiet = T) %>%
nnasua_IUCN st_transform(crs=equalareaCRS)
# Calculate raster and blobs
<- range_edge_dist_expert_range(nnasua_IUCN,
nnasua_range_edge
nnasua_PA_time1_blobs,
nnasua_PA_time2_blobs,
latam, latam_raster)
# grid cell values
<- data.frame(expertCell=nnasua_range_edge[[1]][])
data_nnasua_expert_gridcell # blobs in time 1
<- nnasua_range_edge[[2]]
nnasua_expert_blob_time1 # blobs in time 2
<- nnasua_range_edge[[3]] nnasua_expert_blob_time2
Plot
Code
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
Code
::writeRaster(nnasua_range_edge[[1]], 'data/species_POPA_data/nnasua_expert_gridcell.tif', overwrite=TRUE)
terrasaveRDS(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
<- readRDS('data/species_POPA_data/PA_pbrasiliensis_time1_blobs.rds')
pbrasiliensis_PA_time1_blobs <- readRDS('data/species_POPA_data/PA_pbrasiliensis_time2_blobs.rds')
pbrasiliensis_PA_time2_blobs <- st_read('big_data/pbrasiliensis_IUCN.shp', quiet = T) %>%
pbrasiliensis_IUCN st_transform(crs=equalareaCRS)
# Calculate raster and blobs
<- range_edge_dist_expert_range(pbrasiliensis_IUCN,
pbrasiliensis_range_edge
pbrasiliensis_PA_time1_blobs,
pbrasiliensis_PA_time2_blobs,
latam, latam_raster)
# grid cell values
<- data.frame(expertCell=pbrasiliensis_range_edge[[1]][])
data_pbrasiliensis_expert_gridcell # blobs in time 1
<- pbrasiliensis_range_edge[[2]]
pbrasiliensis_expert_blob_time1 # blobs in time 2
<- pbrasiliensis_range_edge[[3]] pbrasiliensis_expert_blob_time2
Plot
Code
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
Code
::writeRaster(pbrasiliensis_range_edge[[1]], 'data/species_POPA_data/pbrasiliensis_expert_gridcell.tif', overwrite=TRUE)
terrasaveRDS(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')