Luis Verde Arregoitia bio photo

Luis Verde Arregoitia

Personal research page - ecology . evolution . conservation . biogeography

Twitter Publons ResearchGate Google Scholar

This post is available in English here.

Riqueza de especies y otros procesos espaciales en R

Para calcular y visualizar de manera resumida la riqueza de especies en una determinada región, tenemos que calcular el número de registros o de polígonos que resumen la distribución de una especie con respecto a los elementos de otra capa espacial (superpuesta sobre la misma). Esta es una tarea bastante común que se puede hacer con diferentes programas de SIG. He estado trabajando en cuestiones espaciales usando R y el paquete sf, y quería compartir esta forma de cuantificar y graficar mapas de riqueza de especies.

El código en este ejemplo es 100% reproducible, siempre y cuando tengamos instalados los paquetes necesarios. Para este ejemplo escogí trabajar con Costa Rica. En vez de descargar datos de distribución reales de especies que sí existen, aquí vamos a generar puntos al azar, y lo haremos iterativamente gracias al paquete purrr. De todas maneras es fácil importar archivos .shp y usarlos en R con el paquete sf y la función st_read. Aquí vamos a calcular riqueza de especies usando una cuadrícula y dos diferentes formas de representar la distribución de una especie: con datos puntuales y con polígonos derivados de una envoltura que agrupa a los puntos.

Estos son los pasos principales:

Preparación

  • filtrar un subconjunto de datos a partir de un mapamundi para trabajar dentro de un solo país
  • generar un número aleatorio de puntos aleatorios dentro del país de estudio para n diferentes especies
  • generar envolturas que agrupan a los puntos, y ‘suavizar’ sus esquinas

Análisis espaciales

  • generar una cuadrícula que cubra nuestra área de estudio
  • hacer una intersección y unión espacial entre las distribuciones de las especies con cada celda de la cuadrícula

Mapas

  • visualizar la riqueza de especies con colores nítidos, fregones, y adecuados gracias a ggplot, scico y sf

Así van quedando los mapas con los distintos elementos:

División política:

Puntos generados al azar para n especies

Envolturas

Cuadrícula

Riqueza de especies calculada a partir de datos puntuales

Riqueza de especies calculada a partir de polígonos

Ojo:

La función rerun del paquete purrr está muy buena y es un reemplazo útil para no tener que escribir ‘for loops’. No la conocía y la recomiendo ampliamente.

Al dibujar los mapas, usamos un ‘bounding box’ y la función st_touches para que el mapa quede acotado a nuestra área de estudio y para dibujar los países aledaños, sin necesidad de especificar todo ésto a mano.

Acepto cualquier duda o sugerencia.

El script:

# cargar paquetes
library(sf)
library(dplyr)
library(ggplot2)
library(scico)
library(rnaturalearth)
library(purrr)
library(smoothr)

# mapamundi
worldMap <- ne_countries(scale = "medium", type = "countries", returnclass = 'sf')

# subconjunto
CRpoly <- worldMap %>% filter(sovereignt=="Costa Rica")

# puntos al azar, con nombres para los elementos de la lista
sp_occ <- rerun(12,st_sample(CRpoly,sample(3:20,1)))
names(sp_occ) <- paste0("sp_",letters[1:length(sp_occ)])

# pasar a objeto sf
sflisss <- 
  map(sp_occ,st_sf) %>% 
      map2(.,names(.),~mutate(.x,id=.y)) #%>% 
sp_occ_sf <- sflisss %>% reduce(rbind)

# convertir a multipunto
sp_occ_sf <- sp_occ_sf %>% group_by(id) %>% summarise()

# acotar a nuestra área de estudio
limsCR <- st_buffer(CRpoly,dist = 0.7) %>% st_bbox()
# paises aledaños
adjacentPolys <- st_touches(CRpoly,worldMap)
neighbours <- worldMap %>% slice(pluck(adjacentPolys,1))

# división política
divpolPlot <- 
  ggplot()+
  geom_sf(data=neighbours,color="white")+
  geom_sf(data=CRpoly)+
  coord_sf(xlim = c(limsCR["xmin"], limsCR["xmax"]), 
           ylim = c(limsCR["ymin"], limsCR["ymax"]))+
  scale_x_continuous(breaks = c(-84))+
  theme(plot.background = element_rect(fill="#f1f2f3"),
        panel.background = element_rect(fill="#2F4051"),
        panel.grid = element_blank(),
        line = element_blank(), 
        rect = element_blank())

# plot points
spPointsPlot <- 
  ggplot()+
  geom_sf(data=neighbours,color="white")+
  geom_sf(data=CRpoly)+
  geom_sf(data=sp_occ_sf,aes(fill=id),pch=21)+
  scale_fill_scico_d(palette = "davos",direction=-1,end=0.9,guide=FALSE)+
  coord_sf(xlim = c(limsCR["xmin"], limsCR["xmax"]), 
           ylim = c(limsCR["ymin"], limsCR["ymax"]))+
  scale_x_continuous(breaks = c(-84))+
  theme(plot.background = element_rect(fill="#f1f2f3"),
        panel.background = element_rect(fill="#2F4051"),
        panel.grid = element_blank(),
        line = element_blank(), 
        rect = element_blank())
spPointsPlot
ggsave("curso_ammac/002_pointss.png",spPointsPlot,width = 4,height = 3.5, units = "in",dpi = 200)

# smoothed convex hulls
spEOOs <- st_convex_hull(sp_occ_sf) %>% smooth()

# plot hulls
hullsPlot <- 
  ggplot()+
  geom_sf(data=neighbours,color="white")+
  geom_sf(data=CRpoly)+
  geom_sf(data=spEOOs,aes(fill=id),alpha=0.7)+
  scale_fill_scico_d(palette = "davos",direction=-1,end=0.9,guide=FALSE)+
  coord_sf(xlim = c(limsCR["xmin"], limsCR["xmax"]), 
           ylim = c(limsCR["ymin"], limsCR["ymax"]))+
  scale_x_continuous(breaks = c(-84))+
  theme(plot.background = element_rect(fill="#f1f2f3"),
        panel.background = element_rect(fill="#2F4051"),
        panel.grid = element_blank(),
        line = element_blank(), 
        rect = element_blank())
hullsPlot
ggsave("curso_ammac/003_hulls.png",hullsPlot,width = 4,height = 3.5, units = "in",dpi = 200)


# grid
CRGrid <- CRpoly %>%
  st_make_grid(cellsize = 0.2) %>%
  st_intersection(CRpoly) %>%
  st_cast("MULTIPOLYGON") %>%
  st_sf() %>%
  mutate(cellid = row_number())

# calculate n per grid square for points
richness_grid <- CRGrid %>%
  st_join(sp_occ_sf) %>%
  group_by(cellid) %>%
  summarize(num_species = n())
# calculate n per grid square for hulls
richness_gridEOO <- CRGrid %>%
  st_join(spEOOs) %>%
  group_by(cellid) %>%
  summarize(num_species = n())

# blank grid
gridPlot <- 
  ggplot()+
  geom_sf(data=neighbours,color="white")+
  geom_sf(data=CRpoly)+
  geom_sf(data=CRGrid)+
  coord_sf(xlim = c(limsCR["xmin"], limsCR["xmax"]), 
           ylim = c(limsCR["ymin"], limsCR["ymax"]))+
  scale_x_continuous(breaks = c(-84))+
  theme(plot.background = element_rect(fill="#f1f2f3"),
        panel.background = element_rect(fill="#2F4051"),
        panel.grid = element_blank(),
        line = element_blank(), 
        rect = element_blank())
gridPlot
ggsave("curso_ammac/004_grid.png",gridPlot,width = 4,height = 3.5, units = "in",dpi = 200)


# plot gridded richness
gridRichCR <- 
  ggplot(richness_grid)+
  geom_sf(data=neighbours,color="white")+
  geom_sf(data=CRpoly,fill="grey",size=0.1)+
  geom_sf(aes(fill=num_species),color=NA)+
  scale_fill_scico(palette = "davos",direction=-1,end=0.9)+
  coord_sf(xlim = c(limsCR["xmin"], limsCR["xmax"]), 
           ylim = c(limsCR["ymin"], limsCR["ymax"]))+
  scale_x_continuous(breaks = c(-84))+
  theme(plot.background = element_rect(fill="#f1f2f3"),
        panel.background = element_rect(fill="#2F4051"),
        panel.grid = element_blank(),
        line = element_blank(), 
        rect = element_blank())+labs(fill="richness")
gridRichCR
ggsave("curso_ammac/005_gridRichpts.png",gridRichCR,width = 4,height = 3.5, units = "in",dpi = 200)

# richness based on hulls
gridRichCR_eoo <- 
  ggplot(richness_gridEOO)+
  geom_sf(data=neighbours,color="white")+
  geom_sf(data=CRpoly,fill="grey",size=0.1)+
  geom_sf(aes(fill=num_species),color=NA)+
  scale_fill_scico(palette = "davos",direction=-1,end=0.9)+
  coord_sf(xlim = c(limsCR["xmin"], limsCR["xmax"]), 
           ylim = c(limsCR["ymin"], limsCR["ymax"]))+
  scale_x_continuous(breaks = c(-84))+
  theme(plot.background = element_rect(fill="#f1f2f3"),
        panel.background = element_rect(fill="#2F4051"),
        panel.grid = element_blank(),
        line = element_blank(), 
        rect = element_blank())+labs(fill="richness")