R y GRASS: Modelado de nicho

En esta última sesión, vamos a demostrar y ejemplificar el uso combinado de GRASS y R para modelar la distribución de Aedes aegypti en la provincia de Córdoba en función de variables ambientales. Algunas de estas variables provienen de los ejercicios realizados en la unidad de series de tiempo y otras serán generadas durante este ejercicio.

Antes de empezar y para ganar tiempo, conectemos nuestro drive e instalemos GRASS y los paquetes de R que vamos a usar en esta sesión.

# import drive from google colab
from google.colab import drive
# define mounting point for drive
dmp = "/content/drive"
# mount drive
drive.mount(dmp)

%%bash
DEBIAN_FRONTEND=noninteractive 
sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable 
apt update 
apt install grass subversion grass-dev
apt remove libproj22

!grass --config path

Instalamos y cargamos al entorno el paquete de python que nos permite hacer interfaz con R dentro de una Jupyter notebook.

!pip install rpy2==3.5.1
%reload_ext rpy2.ipython

Chequeamos nuestra sesión de R.

%R sessionInfo()

Instalamos todos los paquetes necesarios para este ejercicio. This might take a while…

%%R
install.packages("rgrass")
install.packages("terra")
install.packages("raster")
install.packages("sf")
install.packages("biomod2")
install.packages("dismo")
install.packages("usdm")
install.packages("SDMtune")
install.packages("glmnet")
install.packages("zeallot")
install.packages("ggpubr")
install.packages("tmap")
install.packages("tmaptools")

rgrass

El paquete de R que hace las veces de interfaz entre R y GRASS GIS se denomina rgrass. Este paquete ha sido desarrollado y sigue siendo mantenido por Bivand (2022) y se puede encontrar en: https://github.com/rsbivand/rgrass/.

Nota

Pueden ver la viñeta del paquete para más explicaciones, algo de contexto y ejemplos: https://rsbivand.github.io/rgrass/.

Las principales funciones que podemos encontrar en el paquete rgrass son las siguientes:

• initGRASS(): inicia una sesión de GRASS GIS desde R, similar a gs.setup.init() o gj.init().
• execGRASS(): ejecuta comandos de GRASS, similar a gs.run_command()
• gmeta(): muestra los metadatos de localización de GRASS.
• read_VECT() y read_RAST(): leen mapas vectoriales y raster desde la base de datos de GRASS a objetos SpatVect y SpatRast del paquete terra de R.
• write_VECT() y write_RAST(): escriben objetos del paquete terra en la base de datos de GRASS GIS.

Nota

El paquete terra es desarrollado y mantenido por Hijmans (2022) y eventualmente reemplazará a raster. Detalles sobre el paquete terra pueden encontrarse en: https://rspatial.github.io/terra/reference/terra-package.html y https://rspatial.org/spatial/index.html.

Cómo usamos rgrass?

GRASS GIS y R se pueden utilizar juntos de dos maneras:

A. Usar R dentro de una sesión de GRASS GIS, es decir, iniciar R (o RStudio) desde una sesión de GRASS

• escribimos R o rstudio & en la terminal GRASS GIS o en la pestaña Consola de la interfaz gráfica
• una vez en R (o RStudio) cargamos el paquete rgrass (previo haberlo instalado)
• usamos read_VECT(), read_RAST() para leer datos de GRASS en R
• accedemos a los módulos y la base de datos de GRASS GIS a través de execGRASS()
• escribimos datos resultantes en la base de datos de GRASS con write_VECT() y write_RAST()

B. Iniciar y usar GRASS GIS dentro de una sesión de R, es decir, nos conectamos a la base de datos de GRASS GIS desde R (o RStudio).

• Primero cargamos el paquete rgrass
• Necesitamos iniciar GRASS GIS con initGRASS() desde R y para ello, necesitamos especificar el ejecutable de GRASS y las ubicaciones de la base de datos, el proyecto (location) y mapset
• Accedemos a los módulos GRASS GIS a través de execGRASS()
• usamos read_VECT(), read_RAST(), write_VECT() y write_RAST() para leer datos desde y hacia la base de datos GRASS.

Nota

Originalmente, rgrass estaba destinado a aplicar funciones de GRASS en datos fuera de la base de datos de GRASS; de ahí que algunos prefieran crear proyectos (i.e., locations) desechables o temporarios. Por ejemplo:

# No ejecutar, solo a modo de ejemplo
library(terra)

f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)

library(rgrass)
initGRASS(home=tempdir(), SG=r, override=TRUE)

SDM - modelado de la distribución de especies

El flujo de trabajo comúnmente utilizado para el modelado de la distribución de especies (SDM por su sigla en Inglés) consiste en la recopilación de sitios de ocurrencia de la especie o evento de interés (ej., casos de una enfermedad, focos de fuego, ietc.), así como de variables ambientales que puedan ser relevantes para explicar su distribución. Los datos de ocurrencia y las variables predictoras son las entradas de algún
algoritmo determinado o modelo (ej., GLM, Random Forest, MaxEnt, etc.). Una vez estamos “contentos” con el ajuste del modelo, el próximo paso es la predicción. En este caso, hablamos de prediccón ambiental cuando nos interesa explicar e interpretar la relación entre las variables y la probabilidad de ocurrencia (curvas de respuesta) y predicción geográfica cuando nos interesa ver plasmadas en el espacio geografico el resultado de esas relaciones (mapas).

Existen varios paquetes para realizar SDM. En este caso usaremos SDMtune desarrollado y mantenido por Vignali et al. (2020). Este paquete proporciona funciones que cubren todo el flujo de trabajo de SDM, desde la preparación de datos hasta la selección de variables, la optimización y evaluación de los modelos. Echen un vistazo a los artículos en el sitio web del paquete para obtener más detalles y tutoriales: https://consbiol-unibern.github.io/SDMtune/index.html.

Manos a la obra

Cargamos los paquetes de R

%%R
library(rgrass)
library(sf)
library(terra)
library(biomod2)
library(dismo)
library(usdm)
library(SDMtune)
library(zeallot)
library(tmap)
library(tmaptools)

Iniciamos GRASS

Usaremos la opción B, es decir, iniciamos GRASS GIS desde R en un proyecto y mapset existentes. Notar las similitudes con gj.init() de grass.jupyter y gs.setup.init() de grass.script.

%%R
# path to GRASS binaries (run `grass --config path`)
grassbin <- system("grass --config path", intern = TRUE)
# path to GRASS database in GDrive
grassdata <- "/content/drive/MyDrive/curso_grass_2023/grassdata"
# path to project
project <- "posgar2007_4_cba"
# path to mapset
mapset <- "PERMANENT"

# start GRASS GIS from R
initGRASS(gisBase = grassbin, 
          gisDbase = grassdata, 
          location = project, 
          mapset = mapset, 
          override = TRUE,
          remove_GISRC= TRUE)

Exploramos los mapas raster y vectoriales disponibles en el mapset PERMANENT:

%%R
r <- execGRASS("g.list", 
               parameters = list(type = "raster",
                                 mapset = "."))
v <- execGRASS("g.list", 
               parameters = list(type = "vector",
                                 mapset = "."))

Tenemos alguna máscara?

%%R
execGRASS("r.mask", flags="r")

Lo siguiente antes de comenzar es establecer la región computacional.

%%R
execGRASS("g.region", 
          parameters = list(vector="provincia_cba",
                            align="elevation"),
          flags="p")

Datos de presencia/ocurrencia

Para este ejercicio vamos a generar los puntos de ocurrencia del mosquito Aedes aegypti como un sub-conjunto al azar de las localidades (áreas edificadas) de la provincia de Córdoba que estén a una altura de menos de 850 m. Notar que estamos trabajando en la base de datos de GRASS GIS, pero desde R.

%%R
# Extract centroids from built-up area polygons
execGRASS("v.extract", 
          parameters = list(input="area_edificada_cba", 
                            type="centroid", 
                            output="area_edificada_cba_centroides", 
                            random=250))

# convert centroids to points
execGRASS("v.type", 
          parameters = list(input="area_edificada_cba_centroides", 
                            output="area_edificada_cba_puntos", 
                            from_type="centroid", 
                            to_type="point"))

# extract elevation data for points
execGRASS("v.what.rast", 
          parameters = list(map="area_edificada_cba_puntos", 
                            raster="elevation", 
                            column="elevation"))

# filter points by elevation, keep those <= 850m
execGRASS("v.extract",
          parameters = list(input="area_edificada_cba_puntos",
                            where="elevation <= 850",
                            output="aedes_aegypti"))

Ahora leemos desde GRASS los datos de ocurrencia que generamos, los convertimos en un objeto sf, y los graficamos.

%%R
# Read vector layers
presence <- st_as_sf(read_VECT("aedes_aegypti"))

%%R
# Display presence vector
plot(presence)

Podemos usar la función plot() para visualizar la geometría o la geometría y los atributos (o una selección de los mismos) de los objetos sf.

%%R
# Plot only the geometry 
plot(st_geometry(presence))

# Plot geometry + attr
plot(presence["elevation"])

Datos de background

El algoritmo MaxEnt que vamos a usar en este ejercicio, requiere contrastar las variables ambientales en los sitios de ocurrencia con el resto del ambiente disponible para la especie, el background. Entonces, necesitamos generar puntos de background para caracterizar este ambiente disponible.

Una opción es generar puntos al azar sobre nuestra área de estudio. No obstante, no es cierto que toda la extensión del área de estudio está disponible para Aedes aegypti. Esta especie de mosquito no cría en aguas abiertas y tampoco sobre las salinas. Entonces, vamos a enmascarar primero esas áreas del ambiente disponible.

%%R
# Check region
gmeta()

%%R
# Generar máscara a partir del raster de LULC
expression <- 
  "no_water = if(landcover_2018 == 7 || landcover_2018 == 8, null(), landcover_2018)"

execGRASS("r.mapcalc",
          parameters = list(expression=expression))

Luego, leemos a R la máscara creada para visualizarla

%%R
# Import mask
no_water <- raster(read_RAST("no_water"))
plot(no_water)

Como la región tiene una resolución de 30 m, vamos a llevarla a 1 km, para asegurarnos una mejor separación de los puntos de background, y volvemos a leer la máscara dentro de R.

%%R
# Change resolution
execGRASS("g.region",
          parameters = list(res="1000"),
          flags = c("a","p"))

# Upscale
execGRASS("r.resamp.stats",
          parameters = list(input="no_water",
                            output="MASK",
                            method="mode"))

# Leer en R la mascara_spp que vive en GRASS
mask <- raster(read_RAST("MASK"))

Ahora sí, generamos los puntos de background utilizando una función del paquete dismo(Hijmans et al. 2023) y los convertimos a sf para luego visualizarlos.

%%R
# Generate random points within mask
set.seed(123)
background <- randomPoints(mask = mask, 
                           n = 500)

# Convert to sf to plot
background_sf <- st_as_sf(as.data.frame(background), 
                          coords = c(1,2), 
                          crs = st_crs(mask))

Pregunta

Con qué funcion de GRASS podríamos haber hecho algo similar?

Visualizamos el mapa que usamos como máscara junto con los vectores de presencia y background usando la librería tmap (Tennekes 2018).

%%R
# Aux data
bbox <- st_bbox(mask)

fig_puntos <- 
  tm_shape(mask, 
           bbox = bbox) +
  tm_raster(title = "Classes") +
  tm_shape(presence) +
  tm_dots(size = 0.02) +
  tm_layout(main.title = "Aedes aegypti",
            main.title.fontface = "italic",
            main.title.size = 0.7,
            main.title.position = "left",
            legend.show = TRUE,
            legend.outside = TRUE)

tmap_save(fig_puntos, 
          filename = "fig_puntos_y_mascara.png", 
          width = 1000, height = 1300)

Mapa generado con tmap

Variables ambientales

Antes de leer las variables ambientales que obtuvimos a partir de las series de tiempo de LST y NDVI, vamos a generar dos mapas ráster que representan la distancia a fuentes de agua y rutas y caminos, respectivamente. Para eso, vamos a usar mapas ya disponibles en el mapset PERMANENT.

%%R
# Patch water lines + water bodies
execGRASS("v.patch",
          parameters = list(input="lineas_aguas_continentales,areas_aguas_continentales,embalses",
                            output="lineas_y_cuerpos_de_agua_cba"))

# Convert to raster
execGRASS("v.to.rast",
          parameters = list(input="lineas_y_cuerpos_de_agua_cba",
                            output="lineas_y_cuerpos_de_agua_cba",
                            use="val"))

# Distance to water and roads
execGRASS("r.grow.distance",
          parameters = list(input="lineas_y_cuerpos_de_agua_cba",
                            distance="distancia_agua"))

%%R
# Patch primary + secondary roads
execGRASS("v.patch",
          parameters = list(input="vial_primaria,vial_secundaria",
                            output="red_vial_cba"))
# Convert to raster
execGRASS("v.to.rast",
          parameters = list(input="red_vial_cba",
                            output="red_vial_cba",
                            use="val"))

# Distance to roads
execGRASS("r.grow.distance",
          parameters = list(input="red_vial_cba",
                            distance="distancia_caminos"))

Leemos los mapas generados y los visualizamos con plot(), que en este caso reconoce los objetos SpatRast.

%%R

distancia_agua <- read_RAST("distancia_agua")
distancia_caminos <- read_RAST("distancia_caminos")

plot(c(distancia_agua,distancia_caminos), 
     main=c("Distancia agua", "Distancia rutas"))

También podemos usar tmap y su función tm_facets()

%%R
distancia <- read_RAST(c("distancia_agua","distancia_caminos"))

fig_raster_facet <- 
  tm_shape(distancia, 
           bbox = bbox) +
  tm_raster(style = "cont",
            palette = "magma",
            legend.show = FALSE) +
  tm_facets()

fig_raster_facet

Lectura de datos ráster de otros mapsets

Para leer mapas de otros mapsets, necesitamos agregar esos mapsets a la lista de mapsets accesibles en el mapset donde estamos ahora.

%%R
execGRASS("g.mapsets",
          parameters = list(mapset="modis_lst",
                            operation="add"))
execGRASS("g.mapsets",
          parameters = list(mapset="modis_ndvi",
                            operation="add"))
execGRASS("g.mapsets",
          flags = "p")

Leemos ahora algunas de las variables que derivamos de las series temporales de LST y NDVI. Primero, necesitamos aplicar la máscara de los límites de la provincia.

%%R 
execGRASS("r.mask", 
          parameters = list(vector="provincia_cba"))

%%R
# List rasters to import (modify as you want)
to_import <- c("LST_Day_minimum",
               "LST_Day_maximum",
               "LST_Day_average",
               "ndvi_maximum",
               "ndvi_minimum")

# Read raster layers
predictors <- list()
for (i in to_import){ 
  predictors[i] <- read_RAST(i) }

Agrupamos todos los rasters de variables ambientales.

%%R
# Stack rasters
predictors_r <- rast(c(predictors,distancia))
capas <- c("LST_Day_minimum", 
           "LST_Day_maximum", 
           "LST_Day_average",
           "ndvi_maximum",
           "ndvi_minimum", 
           "distancia_agua",
           "distancia_caminos")
names(predictors_r) <- capas

Preparación de los datos

Ahora que hemos creado y/o importado los registros de presencia, los puntos de background y las variables predictivas, necesitamos preparar los datos en un formato llamado samples with data (SWD). Éste es básicamente una tabla con coordenadas de presencia y background más los valores correspondientes a las variables predictoras para cada punto.

%%R
# Variables for models
sp <- "Aedes aegypti"
presence_coords <- st_coordinates(presence)
background_coords <- background
env <- predictors_r

# Prepare data: SWD
data_sp <- prepareSWD(species = sp, 
                      p = presence_coords, 
                      a = background_coords, 
                      env = env)

data_sp

Definición de parámetros

Aquí definimos algunos de los valores de entrada necesarios para el flujo de trabajo:

%%R
seed=123
perc_test = 0.2
k = 4
method="Maxnet" # open source implementation of MaxEnt
cor_th=0.7
perm=10
imp_th=10

Datos de entrenamiento y evaluación

Entrenaremos el modelo con un 80% de muestras de presencia, y dejaremos el 20% restante para la evaluación al final.

%%R
# Create training and test sets
c(train_sp, test_sp) %<-% 
  trainValTest(data_sp, 
               test = perc_test,
               only_presence = TRUE, 
               seed = seed)

%%R train_sp

%%R test_sp

Creación de folds para validación cruzada

Como usaremos validación cruzada durante el entrenamiento del modelo, creamos los folds con anticipación. En este caso utilizamos folds aleatorios, pero existen otros métodos de determinarlos. Como estamos limitados por la cantidad de registros de presencia, crearemos solo 4 folds o subconjuntos. El algoritmo utilizará iterativamente 3 subconjuntos para entrenar y 1 para validar, pero siempre dentro del entrenamiento.

%%R
# Create folds 
ran_folds <- randomFolds(train_sp, 
                         k = k,
                         only_presence = TRUE, 
                         seed = seed)

Entrenamiento con validación cruzada

Primero entrenaremos un llamado modelo completo, es decir, un modelo con todos los predictores, y de allí eliminaremos aquellos que estén altamente correlacionados y cuya contribución a la predicción no sea importante.

%%R
# Train a full model
full_model_sp <- train(method = method,
                       data = train_sp, 
                       folds = ran_folds)

full_model_sp

Veamos las predicciones geográficas del modelo completo o full model

%%R
pred_full_model <- predict(full_model_sp,
                           data = env,
                           type = "cloglog")

plot(pred_full_model)

Selección de variables

Remover variables altamente correlacionadas

Luego procedemos a eliminar los predictores correlacionados ya que proporcionan información altamente redundante y pueden afectar el rendimiento de los modelos, es decir, como con todos los modelos, queremos que sea simple y del mayor rendimiento posible. Usaremos el área bajo la curva ROC (AUC) como métrica de rendimiento y eliminaremos las variables correlacionadas solo si el AUC disminuye si las mantenemos.

%%R
# Prepare background locations to test correlation
bg_sp <- prepareSWD(species = sp, 
                    a = background_coords,
                    env = env)

# Remove variables with correlation higher than 0.7 
# while accounting for the AUC
vs_sp <- varSel(full_model_sp,
                metric = "auc", 
                bg4cor = bg_sp, 
                cor_th = cor_th,
                permut = perm,
                interactive = FALSE)

Exploremos el objeto de salida

%%R vs_sp@data

Remover variables de menor importancia

Después de descartar las variables correlacionadas, también eliminaremos las variables que tengan una contribución porcentual o una importancia inferior al 10%, considerando como su mantenimiento o remoción afecta al AUC.

%%R
# remove less important variables only if auc does not decrease
reduc_var_sp <- reduceVar(vs_sp,
                          th = imp_th, 
                          metric = "auc", 
                          test = TRUE, 
                          permut = perm, 
                          use_jk = TRUE,
                          interactive = FALSE)

Exploremos el objeto resultante

%R reduc_var_sp

Ahora necesitamos recrear el objeto SWD y los conjuntos de datos de entrenamiento y evaluación, pero solo con las variables seleccionadas, para poder ejecutar el modelo final y hacer predicciones.

%%R
# Get only relevant variables from the reduced model
retained_varnames <- names(reduc_var_sp@models[[1]]@data@data)

# Subset stack
env <- terra::subset(env, retained_varnames)

# SWD with the selected vars
subset_train_sp <- prepareSWD(species = sp, 
                              p = presence_coords,
                              a = background_coords,
                              env = env)

c(train_sp, test_sp) %<-% 
  trainValTest(subset_train_sp, 
               test = perc_test, 
               only_presence = TRUE, 
               seed = seed)

Predicciones con el modelo seleccionado

Ahora entrenamos el modelo final con el conjunto de entrenamiento completo, ya no necesitamos los folds en este punto. Tengan en cuenta que también utilizamos las feature classes (fc) y la regularización (reg) del mejor modelo obtenido anteriormente. En este caso, son solo valores predeterminados, pero si también realizamos una optimización de hiperparámetros, pueden diferir.

%%R
final_model_sp <- train(method = method, 
                        data = train_sp,
                        fc = reduc_var_sp@models[[1]]@model@fc,
                        reg = reduc_var_sp@models[[1]]@model@reg)

Nota

Si les interesa conocer opciones de optimización de hiperparámetros en el contexto de los SDM, puedes chequear el siguiente artículo: https://consbiol-unibern.github.io/SDMtune/articles/tune-hyperparameters.html

Hagamos las predicciones en el espacio geográfico y exploremos el resultado

%%R
map_sp <- predict(final_model_sp,
                  data = env, 
                  type = "cloglog")

plot(map_sp)

Guardamos la predicción en GRASS

Ahora podemos escribir el ráster con las predicciones del modelo final en la base de datos de GRASS.

%%R
write_RAST(map_sp, 
           paste0("Aedes_aegypti_",method), 
           flags = c("o","overwrite"))

Corroboramos que el mapa creado esté allí

%%R
execGRASS("g.list", 
          parameters = list(type="raster",
                            pattern="Aedes*"))

Evaluación del modelo

Queremos saber qué tan bueno es nuestro modelo, por eso en este paso usamos el conjunto de datos de evaluación que separamos al principio. Un AUC de 0,5 significaría que el modelo funciona equivalentemente a lanzar una moneda al aire. AUC es lo que llamamos una métrica de evaluación independiente de umbral.

%%R
# AUC
auc <- auc(final_model_sp, test = test_sp)
auc

Normalmente el resultado del SDM se convierte en mapas de presencia/ausencia. Para determinar qué umbral utilizar, realizamos evaluaciones dependientes del umbral.

%%R
# Threshold dependent evaluation
thres <- thresholds(final_model_sp, 
                    type = "cloglog", 
                    test = test_sp)

knitr::kable(thres, format = 'html', digits = 2)

Elegimos un umbral y creamos un mapa binario, i.e., de presencia y ausencia

%%R
p = map_sp >= 0.5
a = map_sp < 0.5
map_sp[p] <- 1
map_sp[a] <- 0

plot(map_sp)

Importancia de las variables

La importancia de las variables es un indicador de la contribución variable a la predicción.

%%R
vi_model_sp <- varImp(final_model_sp)
vi_model_sp

%%R 
plotVarImp(vi_model_sp)

Curvas de respuesta

Las curvas de respuesta nos dan una idea de la relación entre las variables predictoras y la probabilidad de ocurrencia del evento de interés.

%%R
my_rp <- function(i){
  plotResponse(reduc_var_sp, i)
}

plotlist <- lapply(retained_varnames, my_rp)
labels <- LETTERS[1:length(retained_varnames)]
ggpubr::ggarrange(plotlist = plotlist, labels = labels)

Cerramos el mapset y terminamos :)

%%R
# close the mapset
unlink_.gislock()

Disclaimer

Recordar que éste es sólo un ejemplo sencillo para hacer SDM y sólo el comienzo… Hay:

• otros modelos para probar
• ajuste/optimización de hiperparámetros
• ensemble de modelos
• evaluación de la incertidumbre, i.e., dónde podemos predecir con confianza
• muchos otros paquetes relevantes:
  • dismo, sdm, kuenm, caret, CAST, etc.

Referencias

Bivand, R. (2022), rgrass: Interface Between ’GRASS’ Geographical Information System and ’R’.

Hijmans, R. J. (2022), terra: Spatial Data Analysis.

Hijmans, R. J., Phillips, S., Leathwick, J., y Elith, J. (2023), dismo: Species Distribution Modeling.

Tennekes, M. (2018), «tmap: Thematic Maps in R», Journal of Statistical Software, 84, 1-39. https://doi.org/10.18637/jss.v084.i06.

Vignali, S., Barras, A. G., Arlettaz, R., y Braunisch, V. (2020), «SDMtune: An R package to tune and evaluate species distribution models», Ecology and Evolution, 00, 1-18. https://doi.org/10.1002/ece3.6786.