Distribución Potencial de Especies Utilizando Random Forest

Cargar librearías

require(pacman)
pacman::p_load(rnaturalearthdata, rnaturalearth, cptcity, SSDM, ggspatial,
               raster, terra, sf, fs, glue, tidyverse, rgbif, geodata)

Fijar directorio de trabajo

setwd('C:/Users/Victor/OneDrive - Universidad Arturo Prat/Escritorio/LU_24_RANDOM_FOREST')

Limpiar memoria

g <- gc(reset = T)
rm(list = ls())
options(scipen = 999, warn = -1)

Cargar datos

spce <- 'Araucaria araucana' # Pehuen. Gymnosperma cuyo fruto es comestible (piñón)
occr <- occ_data(scientificName = spce, limit = 2e5, hasCoordinate = T, hasGeospatialIssue = F)

Seleccionar el dataframe desde el objeto

occr <- occr[[2]]
colnames(occr)
unique(occr$continent) %>% sort()
occr <- filter(occr, continent %in% c('SOUTH_AMERICA'))
occr

Shapefiles

wrld <- ne_countries(returnclass = 'sf', scale = 50)
class(wrld)
plot(st_geometry(wrld))

Obtener shapefile desde el objeto wrld

sudam <- wrld %>% filter(continent == "South America")
ggplot() + geom_sf(data = sudam, fill="cyan2") 

Dibujar mapa de sudamerica con ocurrencias de A. araucana

plot(sudam$geometry)
points(occr$decimalLongitude, occr$decimalLatitude, pch = 16, col = 'red')

Obtener raster bioclimáticos desde “Worldclim”

bioc <- geodata::worldclim_global(var='bioc', res=2.5,
                                  path='tmpr',version="2.1")

Cortar el raster multicapa bioc usando como máscara “sudam”

bioc <- terra::crop(bioc, sudam) %>% terra::mask(., sudam)
names(bioc) <- glue('bioc{1:19}')

bioc

Se crea un data frame con puntos de ocurrencias y datos climaticos

occr <- dplyr::select(occr, x = decimalLongitude, y = decimalLatitude)
vles <- terra::extract(bioc, occr[,c('x', 'y')])
occr <- cbind(occr[,c('x', 'y')], vles[,-1])
occr <- as_tibble(occr)
occr <- mutate(occr, pb = 1)
occr

Generar puntos de background (pseudoausencias)

cell <- terra::extract(bioc[[1]], occr[,1:2], cells = T)$cell
duplicated(cell)
mask <- bioc[[1]] * 0
mask[cell] <- NA
back <- terra::as.data.frame(mask, xy = T) %>% as_tibble()
back <- sample_n(back, size = nrow(occr) * 2, replace = FALSE)
colnames(back)[3] <- 'pb'
back <- mutate(back, pb = 0)
back <- cbind(back, terra::extract(bioc, back[,c(1, 2)])[,-1])
back <- as_tibble(back)

Unir df “occr” con los puntos de background

tble <- rbind(occr, back)

Aplicar modelo Random forest

bioc <- stack(bioc)
tble <- as.data.frame(tble)
sdrf <- modelling(algorithm = 'RF', Env = bioc, Occurrences = tble,
                  Pcol = 'pb', Xcol = 'x', cv.parm = c(0.75, 0.25), 
                  Ycol = 'y', metric = 'TSS', select.metric = 'AUC')

Gráficos de proyección continua y binaria

plot(sdrf@projection)
plot(sdrf@binary)
sdrf@parameters

Data frame

sdrf@name
sdrf@variable.importance
as.numeric(sdrf@variable.importance) %>% sum()

Proyección (rslt)

rstr <- sdrf@projection
rstr
rstr <- terra::rast(rstr)
plot(rstr)
writeRaster(rstr, filename='proyeccion.tif', overwrite=TRUE)
rslt <- terra::as.data.frame(rstr, xy = T) %>% as_tibble()
head(rslt)

Hacer el mapa

windowsFonts(georg = windowsFont('Georgia'))
gmap <- ggplot() + 
  geom_tile(data = rslt, aes(x = x, y = y, fill = layer)) +
  scale_fill_gradientn(colors = cpt(pal = 'mpl_inferno', n = 8, rev = TRUE)) +
  geom_sf(data = wrld, fill = NA, col = 'grey40', lwd = 0.2) + 
  geom_sf(data = st_as_sf(sudam), fill = NA, col = 'grey40', lwd = 0.3) + 
  coord_sf(xlim = ext(sudam)[1:2], ylim = ext(sudam)[3:4]) + 
  labs(x = 'Lon', y = 'Lat', fill = 'Puntaje de idoneidad') + 
  ggtitle(label = 'Idoneidad para Araucaria araucana en sudamérica', 
          subtitle = 'Modelo Random Forest') +
  theme_bw() + 
  theme(text = element_text(family = 'georg', color = 'grey50'), 
        legend.position = 'bottom', 
        plot.title = element_text(hjust = 0.5, face = 'bold', color = 'grey30'),
        plot.subtitle = element_text(hjust = 0.5, face = 'bold', color = 'grey30'),
        # legend.key.width = unit(3, 'line'),
        panel.border = element_rect(color = 'grey80')) +
  guides(fill = guide_legend( 
    direction = 'horizontal',
    keyheight = unit(1.15, units = "mm"),
    keywidth = unit(15, units = "mm"),
    title.position = 'top',
    title.hjust = 0.5,
    label.hjust = .5,
    nrow = 1,
    byrow = T,
    reverse = F,
    label.position = "bottom"
  )) + 
  annotation_scale(location =  "bl", width_hint = 0.5, text_family = 'georg', text_col = 'grey60', bar_cols = c('grey60', 'grey99'), line_width = 0.2) +
  annotation_north_arrow(location = "tr", which_north = "true", 
                         pad_x = unit(0.1, "in"), pad_y = unit(0.2, "in"), 
                         style = north_arrow_fancy_orienteering(text_family = 'georg', text_col = 'grey40', line_col = 'grey60', fill = c('grey60', 'grey99'))) 

gmap

Grabar este mapa en el disco duro

ggsave(plot = gmap, filename = 'mapa_rf.tif', device='tiff',units = 'in', width = 9, height = 7, dpi = 300)