MigClim_Hippos
Julián Velasco
10/11/2021
La mayoría de estudios que han usado ENM/SDM para generar proyecciones a futuro bajo escenarios de cambio climático o de uso de suelo ignoran aspectos de historia de vida de las especies, particularmente limitaciones por dispersión y aspectos demográficos. Dependiendo de la fragmentación del hábitat y las capacidades de dispersión de las especies, estos supuestos suelen ser inadecuados y generan proyecciones con mucha incertidumbre.
En este ejercicio usaremos un modelo que fue conceptualizado en el 2009 por Robin Engler y Antoine Guisan llamado MigClim. El modelo fue utilizado para evaluar la distribución potencial futura de los hipopótamos invasores en Colombia que se encuentran en la cuenca del río Magdalena. El modelo permite integrar la tasa a la que el clima cambia, fragmentación del paisaje, capacidades de dispersión de las especies, y potencial reproductivo. MigClim está basado en autómatas celulares (von Neumann circa 1950), que es un modelo matemático y computacional de evolución en una serie de pasos discretos. En el caso particular de la modelación de nicho, los valores de favorabilidad (“suitability”) de cada celda evolucionan en el tiempo de acuerdo a una serie de reglas (e.g., cambio climático) y que están influenciadas por el vecindario de celdas. El modelo simula dispersión, colonización, crecimiento y extinción de una especie en el paisaje.
Para simular la distribución potencial de los hipopótamos invasores en Colombia necesitamos varios parámetros:
- Frecuencia de eventos de dispersión (años). Este indica el tiempo entre dos eventos de dispersión sucesivos.
- Distancia de dispersión (DispDist) (m).
- Máxima distancia de dispersión de larga distancia (LDD) (m).
- Frecuencia de LDD.
- Potencial reproductivo del pixel. Indica cuál es la probabilidad de que un pixel colonizado sea un pixel fuente e incremente con el tiempo.
- Barreras geográficas.
- Filtros (e.g., áreas urbanas o innaccesibles).
Antes de comenzar hay que instalar el paquete y es un poco tedioso. La versión que se usó para el artículo de Castelblanco-Martínez et al. (2020) es la 1.6.
# Instalar paquetes
# primero instalar SDMTools
# install.packages("remotes")
# remotes::install_version("SDMTools", "1.1-221")
# luego instalar MigClim versión 1.6
# véase https://cran.r-project.org/src/contrib/Archive/MigClim/
# install_version("MigClim", version = "1.6", repos = "http://cran.us.r-project.org")
library(MigClim)## Loading required package: SDMTools## Loading required package: raster## Loading required package: splibrary(raster)
# Funciones para ver cuáles funciones y la ayuda de cada paquete
lsf.str("package:MigClim")## MigClim.genClust : function (hsMap = "hsMap", barrier = "barrier", nrClusters = 4, nrIterations = 1,
## threshold = 445, outFile = "out", initFile = "")
## MigClim.migrate : function (iniDist = "InitialDist", hsMap = "HSmap", rcThreshold = 0, envChgSteps = 1,
## dispSteps = 1, dispKernel = c(1, 1), barrier = "", barrierType = "strong",
## iniMatAge = 1, propaguleProd = c(1), lddFreq = 0, lddMinDist = NULL,
## lddMaxDist = NULL, simulName = "MigClimTest", replicateNb = 1, overWrite = FALSE,
## testMode = FALSE, fullOutput = FALSE, keepTempFiles = FALSE)
## MigClim.plot : function (asciiFile, outDir = "", fileFormat = "jpeg", fullOutput = FALSE)
## MigClim.userGuide : function ()
## MigClim.validate : function (validateFile = "Validation.txt", nrPoints = 0, simFile = "out1.asc",
## nrClusters = 4)help(package = MigClim)Vamos a cargar los datos pre-fabricados para ahorrar tiempo. Esta tabla de datos muestra para cada pixel en la región geográfica la presencia/ausencia (0-1) de los hipopótamos y los valores de favorabilidad (suitability) obtenidos de un ensamble de modelos. Los valores de favorabilidad fueron estandarizados de 0 a 1000. Hay una columna al final que indica si en la región un pixel determinado es accessible (1) o no (0).
Para ajustar el modelo hay que primero obtener los parámetros de los hipopótamos de la literatura.
Descargue los datos de aquí: https://www.dropbox.com/sh/btys7ucruafoemd/AAARkElKrvUqUeQbiaXlnP6na?dl=0
# Cargar los datos pre-fabricados
MigClim.Data.v1 <- read.csv("./Datos/Dataset_Hippos_MigClim.csv", row.names = 1)
rounds <- round(MigClim.Data.v1[c(-1,-2)], digits = 0)
MigClim.Data.v1 <- cbind(MigClim.Data.v1[,c(1,2)], rounds)
# Explorar los datos
head(MigClim.Data.v1)## X Y InitialDist hsmap1 hsmap2 hsmap3 hsmap4 hsmap5 hsmap6
## 1 -74.85417 11.10417 0 303 284 302 257 127 367
## 2 -74.81250 11.06250 0 378 390 305 257 128 362
## 3 -74.77083 11.06250 0 392 385 272 257 128 287
## 4 -74.89583 11.02083 0 418 183 236 325 99 224
## 5 -74.85417 11.02083 0 450 181 236 257 99 219
## 6 -74.81250 11.02083 0 600 189 236 257 99 219
## hsmap7 hsmap8 hsmap9 hsmap10 hsmap11 hsmap12 hsmap13 hsmap14 hsmap15 hsmap16
## 1 253 199 222 222 100 337 185 268 179 284
## 2 258 155 286 225 130 341 165 269 183 223
## 3 263 155 286 225 130 341 165 269 183 223
## 4 260 158 290 175 132 226 137 268 126 226
## 5 191 155 286 172 130 272 134 268 126 222
## 6 191 155 286 145 130 223 147 268 126 222
## hsmap17 hsmap18 hsmap19 hsmap20 hsmap21 hsmap22 hsmap23 hsmap24 hsmap25
## 1 102 376 101 363 307 257 122 102 171
## 2 132 396 102 285 388 298 81 85 190
## 3 132 396 102 285 388 260 96 81 209
## 4 102 327 79 363 387 173 60 81 171
## 5 132 355 79 359 388 178 85 80 171
## 6 117 508 79 359 388 171 85 80 171
## hsmap26 hsmap27 hsmap28 hsmap29 hsmap30 hsmap31 hsmap32 hsmap33 hsmap34
## 1 103 333 206 226 279 140 288 303 162
## 2 91 337 208 263 285 134 295 238 186
## 3 92 337 239 263 286 111 253 238 186
## 4 51 260 161 232 374 109 227 238 140
## 5 69 290 208 229 455 116 224 234 139
## 6 71 263 208 150 475 111 224 299 141
## hsmap35 hsmap36 hsmap37 hsmap38 hsmap39 hsmap40 hsmap41 hsmap42 hsmap43
## 1 87 391 324 117 268 204 193 309 190
## 2 88 396 440 119 272 265 195 360 192
## 3 88 396 440 119 272 265 195 359 192
## 4 68 311 280 119 242 265 195 269 153
## 5 68 305 336 119 249 265 195 268 192
## 6 68 309 344 131 253 266 195 225 192
## hsmap44 hsmap45 hsmap46 hsmap47 hsmap48 hsmap49 hsmap50 hsmap51 hsmap52
## 1 277 247 204 193 284 324 323 187 108
## 2 360 288 209 224 285 377 281 146 125
## 3 361 250 206 195 285 328 324 146 126
## 4 429 253 162 227 285 289 258 146 101
## 5 460 251 159 224 322 290 260 146 101
## 6 490 276 204 196 281 292 263 187 102
## hsmap53 hsmap54 hsmap55 hsmap56 hsmap57 hsmap58 hsmap59 hsmap60 hsmap61
## 1 158 145 236 273 150 275 273 146 326
## 2 183 145 272 276 117 277 351 174 308
## 3 183 133 272 276 117 282 351 174 308
## 4 143 147 245 217 119 279 351 183 268
## 5 142 145 275 240 117 275 352 164 269
## 6 182 145 254 330 117 275 352 229 271
## hsmap62 hsmap63 hsmap64 hsmap65 hsmap66 hsmap67 hsmap68 hsmap69 hsmap70
## 1 136 187 133 201 174 355 199 394 348
## 2 137 210 126 199 143 358 201 483 300
## 3 137 210 126 229 143 360 201 398 326
## 4 137 252 110 242 144 362 243 399 448
## 5 167 243 126 243 176 371 246 394 584
## 6 167 244 110 244 179 356 246 394 736
## hsmap71 hsmap72 hsmap73 hsmap74 hsmap75 hsmap76 hsmap77 hsmap78 hsmap79
## 1 196 186 152 229 153 82 161 220 252
## 2 226 178 125 220 127 78 133 209 207
## 3 197 178 143 220 127 83 136 214 207
## 4 196 155 128 191 130 82 166 208 208
## 5 196 155 126 190 158 82 169 263 255
## 6 196 155 154 194 165 82 180 320 255
## hsmap80 hsmap81 hsmap82 hsmap83 hsmap84 hsmap85 hsmap86 hsmap87 hsmap88
## 1 170 176 212 264 183 246 104 138 112
## 2 161 144 178 223 225 234 86 113 113
## 3 172 165 178 225 224 234 86 113 114
## 4 170 175 203 160 210 235 108 141 115
## 5 170 176 208 237 225 235 110 141 118
## 6 171 178 227 301 227 211 114 143 127
## hsmap89 hsmap90 hsmap91 hsmap92 hsmap93 hsmap94 hsmap95 Barrier
## 1 134 186 294 197 87 109 192 0
## 2 127 153 232 199 107 110 193 0
## 3 127 153 232 200 107 110 193 0
## 4 113 186 295 197 87 113 192 0
## 5 113 186 295 199 104 114 195 0
## 6 138 187 297 225 106 118 207 0PRIMER PASO:
- Crear un kernel de dispersión (dispKernel)
Distribución de probabilidad de la dispersión para un periodo de 95 años Uniforme - Aplanado. Asume que la dispersión es constante a través del tiempo.
dispKernel.flat <- runif(95, min=1, max=1)
plot(dispKernel.flat) # n1
Decaimiento exponencial en la probabilidad de dispersión a través de los pasos de dispersión (dispSteps) en un periodo de 95 años. Primero creamos una función que simule una distribución exponencial que vaya cayendo a través del tiempo (Y = Yo * Exp(-beta*t)).
# decaimiento exponencial en la probabilidad de dispersión
alpha <- 1 # dispersión inicial
beta <- -0.05 # tasa de decaimiento
n <- 95 # años
x <- seq(n) # secuencia de años
dispKernel.expdecay <- alpha * exp(beta * x)
plot(dispKernel.expdecay) # n2
Producción de propágulos (plantas) o crías (animales)
propaguleProdProb <- c(0.02, 0.1, 0.3, 0.5, 0.9)
# Información reproductiva tomada de Tabla 1 (Castelblanco-Martínez et al 2020)
initMatAge=5 # Edad inicial de madurez sexual (hembras)
maxMatAge=43 # Edad máxima de madurez sexual (hembras)Crearemos una distribución uniforme de la probabilidad de reproducción (entre 0.9 y 1) de los hipopótamos en función de la duración de la etapa reproductiva (43 años) una vez han colonizado una celda en particular.
# distribución aplanada o uniforme
pp1 <- runif(43, min=0.9, max=1)
plot(pp1, ylim=c(0,1))
Crearemos una distribución exponencial de la probabilidad de reproducción de los hipopótamos en función de la duración de la etapa reproductiva (43 años) una vez han colonizado una celda en particular.
# distribución exponencial
alpha <- 0.1 # probabilidad de reproducción apenas coloniza una celda
beta <- 0.05 # tasa de crecimiento
n <- 43 # años
x <- seq(n) # secuencia de años
pp2 <- alpha * exp(beta * x) # n2
plot(pp2)
Ahora vamos a simular las probabilidades de dispersión, colonización, y extinción de los hipopótamos en la región de la cuenca del Madgalena. Se crearon cuatro escenarios. Vamos a generar uno solo y les queda de tarea los otros tres.
PRIMER ESCENARIO:
Probabilidad de dispersión uniforme a través del tiempo (dispKernel.flat) y probabilidad de reproducción exponencial de los hipopótamos una vez colonicen una celda (pp2), y una barrera geográfica fuerte.
Explicación breve de los parámetros:
iniDist: Distribución inicial de la especies (presencia-ausencia)
hsMap: Favorabilidad (suitability) para cada paso del tiempo. Como son 95 años, estos valores de favorabilidad se obtuvieron de proyecciones de un ensamble de modelos (recuerden el ejercicio pasado con el ratón del género Neotoma).
rcThreshold: Si el hsMap es continuo, entonces este parámetro indica a partir de cuál valor de favorabilidad los hipopótamos pueden colonizar una celda en particular. Recuerden el ejercicio donde se selecciona el umbral y se binarizan los mapas.
envChgSteps: Número de pasos ambientales. 95 pasos (años) para este caso.
dispSteps: Indica cuántos eventos de dispersión habrá por año (5 eventos) en un periodo de 95 años, lo que indica que serán 950 pasos de dispersión en la simulación completa.
dispKernel: El kernel de dispersión es una función que indica la probabilidad que un propágulo/cría se disperse a cierta distancia de la fuente. Debe ser expresado en unidades de celda y por lo tanto da la probabilidad que se disperse a 1, 2, 3, 4, 5, 6 celdas.
barrier: Indica cuáles celdas puede o no colonizar una especie (binario: 0-1)
lddFreq: Frecuencia de long distance dispersal (LDD): rango entre 0 y 1. Un valor de 0.01 indica que, en promedio, 1 celda entre 100 generarán un evento de LDD definido por su distancia mínima y máxima.
lddMinDist: 96 (celdas de 5 km2; 2.23 x 2.23 km; 96 x 2.23 = 480 km) # basado en datos de mvto.
lddMaxDist: 1000 (celdas de 5 km2: 2.23 x 2.23 km; 96 x 2.23 = 2230 km) # hipótetico.
propaguleProd: Probabilidad que una celda ocupa produzca propágulos/crías en función del tiempo desde que fue colonizada.
iniMatAge: Edad inicial de madurez sexual (hembras)
Correr la simulación con 1 repetición:
n1 <- MigClim.migrate (iniDist=MigClim.Data.v1[,1:3],
hsMap=MigClim.Data.v1[,4:98], rcThreshold=64,
envChgSteps=95, dispSteps=95, dispKernel=dispKernel.flat,
barrier=MigClim.Data.v1[,99], barrierType="strong",
iniMatAge=5, propaguleProd=pp2,
lddFreq=0.01, lddMinDist=96, lddMaxDist=1000,
simulName="SimA", replicateNb=1, overWrite=TRUE,
testMode=FALSE, fullOutput=TRUE, keepTempFiles=FALSE)## Converting data to ascii grid format...
## Starting simulation for SimA ...
## Running MigClim simulation SimA.
## 1...
## 2...
## 3...
## 4...
## 5...
## 6...
## 7...
## 8...
## 9...
## 10...
## 11...
## 12...
## 13...
## 14...
## 15...
## 16...
## 17...
## 18...
## 19...
## 20...
## 21...
## 22...
## 23...
## 24...
## 25...
## 26...
## 27...
## 28...
## 29...
## 30...
## 31...
## 32...
## 33...
## 34...
## 35...
## 36...
## 37...
## 38...
## 39...
## 40...
## 41...
## 42...
## 43...
## 44...
## 45...
## 46...
## 47...
## 48...
## 49...
## 50...
## 51...
## 52...
## 53...
## 54...
## 55...
## 56...
## 57...
## 58...
## 59...
## 60...
## 61...
## 62...
## 63...
## 64...
## 65...
## 66...
## 67...
## 68...
## 69...
## 70...
## 71...
## 72...
## 73...
## 74...
## 75...
## 76...
## 77...
## 78...
## 79...
## 80...
## 81...
## 82...
## 83...
## 84...
## 85...
## 86...
## 87...
## 88...
## 89...
## 90...
## 91...
## 92...
## 93...
## 94...
## 95...
## All dispersal steps completed. Final output in progress...
## Simulation SimA completed successfully. Outputs stored in /home2/Dropbox/Curso/Ejercicios_R/MigClim/SimAVamos a generar la figura del raster final
```r
# Descargue un mapa base
colombia <- raster::getData('GADM', country="COL", level=1)
# cargar simulación final
run1 <- raster("./SimA/SimA_raster.asc")
# crear matriz de reclasificación
reclass <- c(-Inf, 0, 0, 0, 1, 1, 1, Inf, 2)
reclass_m <- matrix(reclass, ncol = 3, byrow = TRUE)
reclass2 <- c(-Inf, 0, 0, 0, 1, 1)
reclass_m2 <- matrix(reclass2, ncol = 2, byrow = TRUE)
# graficar salida cruda
plot(run1)
plot(colombia, add=T, lwd=0.2)
# Definir el esquema de colores
cols <- c("gray89", "blue", "green")
# reclasificar
run1.r <- reclassify(run1,reclass_m)
run1.r2 <- reclassify(run1,reclass_m2)Primer mapa
# graficar
plot(run1.r, col=cols, legend=TRUE, main="Corrida 1")
plot(colombia, add=T, lwd=0.2)
Segundo mapa
# graficar
plot(run1.r2, legend=TRUE, main="Corrida 1")
plot(colombia, add=T, lwd=0.2)
La tarea consiste en probar varias combinaciones de parámetros:
Kernels de dispersión: 1) dispKernel.flat 2) dispKernel.expdecay
Probabilidad de producción de propágulos/crías 1) pp1: uniforme 2) pp2: exponencial
Frecuencia de eventos de dispersión de larga distancia 1) 0.01 2) 0.1 3) 0.5
Tipos de barreras de dispersión 1) fuerte (strong) 2) débil (weak)
A jugar!