Vamos a generar modelos usando dos algoritmos -Bioclim y Maxent-

Primero subimos nuestros datos de coordenadas:

library(rgdal)
## Loading required package: sp
## rgdal: version: 1.3-1, (SVN revision 747)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/juvel/Documents/R/win-library/3.5/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/juvel/Documents/R/win-library/3.5/rgdal/proj
##  Linking to sp version: 1.2-7
library(sp)
library(raster)
library(rgeos)
## rgeos version: 0.3-26, (SVN revision 560)
##  GEOS runtime version: 3.6.1-CAPI-1.10.1 r0 
##  Linking to sp version: 1.2-7 
##  Polygon checking: TRUE
library(dismo)


coords <- read.csv("Coordenadas_GBIF_depuradas.csv")
head(coords)
##    X Longitude Latitude
## 1  5 -100.7055 22.49427
## 2 19 -101.0598 22.71030
## 3 20 -101.5221 26.31369
## 4 26 -105.9166 28.65161
## 5 65 -104.8663 26.89161
## 6 66 -100.9707 22.16421

Le quitamos la primera columna, no la necesitamos

coords <- coords[,c(2:3)]
head(coords)
##   Longitude Latitude
## 1 -100.7055 22.49427
## 2 -101.0598 22.71030
## 3 -101.5221 26.31369
## 4 -105.9166 28.65161
## 5 -104.8663 26.89161
## 6 -100.9707 22.16421

Subimos las capas climáticas

clim.act <- list.files("C:/Users/juvel/Dropbox/Cursos_Charlas/Zaragosa_2018/Curso_FES_Zaragosa18/Capas_clima/", pattern="asc", full.names=T)
clim.act <- stack(clim.act) # hacemos un stack

clim.act
## class       : RasterStack 
## dimensions  : 2183, 3799, 8293217, 19  (nrow, ncol, ncell, nlayers)
## resolution  : 0.008333334, 0.008333334  (x, y)
## extent      : -118.3684, -86.71002, 14.5296, 32.72126  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## names       : bio1, bio10, bio11, bio12, bio13, bio14, bio15, bio16, bio17, bio18, bio19, bio2, bio3, bio4, bio5, ...
plot(clim.act[[1]])

Vamos a hacer un buffer a nuestros puntos. Este buffer será nuestra área M del Diagrama BAM y es el área que la especie hipóteticamente ha podido “explorar” durante su historia evolutiva (i.e., dispersión, colonización, extinción):

coordinates(coords) <- ~Longitude+Latitude

library(rgeos)
buffer.sp <- gBuffer(coords, byid=FALSE, id=NULL, width=3.0, quadsegs=5, 
                     capStyle="ROUND", joinStyle="ROUND", mitreLimit=1.0)

plot(buffer.sp)
plot(coords, add=TRUE)

Vamos a cortar las capas climáticas usando este buffer

# Usaremos mask, puede tardar un poco, así que tengan un poco de paciencia.
clim.m <- crop(clim.act, buffer.sp)
plot(clim.m[[1]])

Ahora vamos a calibrar el modelo usando el algoritmo Bioclim

# primero convertimos las coordenadas a dataframe

coordsDF <- as.data.frame(coords)

bioc.modesta <- bioclim(clim.m, coordsDF)

Vamos a ver las gráficas de las variables de respuesta

response(bioc.modesta)

plot(bioc.modesta, a=1, b=2, p=0.95)

Ahora vamos a generar la predicción geográfica:

pred.bioc <- predict(clim.m, bioc.modesta)

plot(pred.bioc) # graficamos

Vamos a generar una predicción binaria (0 ó 1; presencia vs. ausencia). Usaremos el criterio de umbral del “minimum presence training”

coords.bioc <- extract(pred.bioc, coords)
head(coords.bioc)
## [1] 0.12408759 0.00729927 0.30656934 0.18248175 0.15328467 0.03649635
a <- min(coords.bioc)
a
## [1] 0.00729927

Usaremos ese valor para generar la predicción binaria. Por encima de 0.0073 será presencia y por debajo, ausencia.

bin.bioc <- reclassify(pred.bioc,c(-Inf,a,0, a,Inf,1))
bin.bioc
## class       : RasterLayer 
## dimensions  : 1742, 2509, 4370678  (nrow, ncol, ncell)
## resolution  : 0.008333334, 0.008333334  (x, y)
## extent      : -115.61, -94.70169, 18.2046, 32.72126  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : layer 
## values      : 0, 1  (min, max)

Graficamos

plot(bin.bioc)

Vemos que la predicción geográfica ocupa una gran parte del área que se usó para calibrar (entorno, M, background). Es necesario probar con un entorno quizás más grande (valores de 5 ó 6 en el buffer que se hizo para cortar las capas). Pero vamos a hacer el mismo modelo usando la extensión completa de las capas.

bioc.modesta.v2 <- bioclim(clim.act, coordsDF)
pred.bioc.v2 <- predict(clim.act, bioc.modesta.v2)
plot(pred.bioc.v2)

Comparemos las dos predicciones:

par(mfrow=c(1,2))
plot(pred.bioc, main="M específica")
plot(pred.bioc.v2, main="M todo México")

Aquí hemos calibrado un modelo de nicho ecológico usando todos los datos disponibles y por lo tanto no tenemos forma de evaluar que tanto error puede tener nuestra predicción geográfica. Vamos a partir los datos en dos conjuntos. Uno de calibración y uno de validación:

La partición de los datos será aleatoria. Usaremos la función kfold para esto. Usaremos 5 grupos

group <- kfold(coordsDF, 5)
group
##   [1] 5 4 3 1 4 4 3 1 4 2 4 3 2 3 1 2 1 5 2 4 1 4 5 2 5 3 3 3 5 5 4 5 4 4 1
##  [36] 2 3 3 1 5 5 4 5 5 1 2 4 5 4 4 1 3 2 4 3 5 2 3 4 1 1 5 1 2 3 2 1 3 5 1
##  [71] 2 4 1 2 5 2 5 1 1 3 3 4 2 2 3 3 2 3 1 1 1 4 3 2 3 4 4 4 2 4 5 1 5 5 3
## [106] 3 3 2 5 5 3 4 4 2 2 5 5 1 2 4 1 3 1 2 1 1 2 4 2 2 2 1 5 5 5 4 3
length(group)
## [1] 137
pres_train <- coordsDF[group != 1, ]
pres_test <- coordsDF[group == 1, ]

pres_train
##      Longitude Latitude
## 1   -100.70555 22.49427
## 2   -101.05979 22.71030
## 3   -101.52213 26.31369
## 5   -104.86634 26.89161
## 6   -100.97070 22.16421
## 7   -101.00324 27.92715
## 9   -100.23816 24.36905
## 10  -101.10850 23.89352
## 11  -100.84574 25.24488
## 12  -100.66156 26.03307
## 13  -100.88555 23.68940
## 14  -105.00584 30.46371
## 16  -101.89417 22.83517
## 18  -102.62886 28.46958
## 19  -102.73839 28.74661
## 20  -100.70377 25.70094
## 22  -102.80069 28.96800
## 23  -100.95000 28.53361
## 24  -101.67297 25.62062
## 25  -101.28340 24.06499
## 26  -100.94028 24.27906
## 27  -103.37386 25.15272
## 28  -101.27942 25.84200
## 29  -100.06598 24.06585
## 30   -99.87808 23.53440
## 31  -103.15372 26.64569
## 32  -103.35092 26.96283
## 33  -103.48183 27.21056
## 34  -103.10339 26.20875
## 36  -103.67567 28.04439
## 37  -104.43094 28.04197
## 38  -105.33147 29.54167
## 40  -103.69644 26.20253
## 41  -104.16469 27.81803
## 42  -100.24078 22.57936
## 43  -105.26272 29.97431
## 44  -106.08233 30.13389
## 46  -104.22522 28.79339
## 47  -106.37528 31.15117
## 48  -100.41778 23.16293
## 49  -100.50542 22.68775
## 50  -107.80564 31.39231
## 52   -99.23429 24.61317
## 53  -101.07028 26.31472
## 54  -100.15166 22.85944
## 55  -104.51544 26.56981
## 56  -104.00000 27.52000
## 57  -101.36018 24.58038
## 58  -100.45800 25.65000
## 59  -101.54458 24.08967
## 62  -108.06704 31.00012
## 64  -101.67056 22.69444
## 65  -101.62278 22.94639
## 66  -101.24917 22.86833
## 68  -106.17268 28.53829
## 69  -103.68357 25.41701
## 71  -103.20022 23.96040
## 72  -105.50029 28.23356
## 74  -103.76691 25.01703
## 75  -101.39018 24.35039
## 76  -105.20618 29.03603
## 77  -106.36032 29.57018
## 80  -106.50767 29.30628
## 81  -107.61704 31.55010
## 82  -100.85034 27.01700
## 83  -107.30035 30.28348
## 84  -105.78684 30.53603
## 85  -102.37021 27.28027
## 86  -102.27021 26.24031
## 87  -102.55411 29.06306
## 88  -108.15037 30.30014
## 92  -101.42240 25.56363
## 93  -107.31958 31.15539
## 94  -101.22000 22.35000
## 95  -103.19192 25.60550
## 96  -103.24167 23.55833
## 97  -104.94167 27.10833
## 98  -102.15833 23.50833
## 99  -100.05833 25.79167
## 100 -102.02141 29.15362
## 101 -102.87022 25.92032
## 103 -100.70833 25.47500
## 104 -106.40833 29.90833
## 105 -104.87500 27.39167
## 106 -106.23194 28.09008
## 107 -100.02500 21.57500
## 108 -107.47500 28.47500
## 109 -102.15833 25.47500
## 110 -106.17500 29.77500
## 111 -103.62500 25.89167
## 112 -108.29167 31.12500
## 113  -99.70833 22.99167
## 114  -97.70008 21.20054
## 115 -112.60833 31.77500
## 116 -103.15655 27.31756
## 117 -107.64167 30.19167
## 119 -105.08367 29.46339
## 120 -102.89167 25.62500
## 122 -103.22500 25.99167
## 124 -104.45917 24.77750
## 127 -106.52500 30.47500
## 128 -102.50833 23.80833
## 129 -102.72500 26.24167
## 130 -109.67500 29.80833
## 131 -107.77500 29.25833
## 133 -101.52500 27.89167
## 134 -101.54167 25.22500
## 135 -106.24167 30.62500
## 136 -101.04167 25.45833
## 137 -103.25505 27.82208
pres_test
##      Longitude Latitude
## 4   -105.91661 28.65161
## 8    -99.12439 23.64901
## 15  -101.71500 23.85833
## 17  -105.05407 27.65698
## 21  -102.91481 27.19092
## 35  -107.15242 29.98267
## 39  -104.09783 28.11419
## 45  -106.08064 31.30008
## 51  -106.58603 31.37706
## 60  -100.38000 24.09300
## 61  -105.05944 25.07444
## 63  -103.83600 26.73600
## 67  -105.19697 28.04955
## 70  -100.32000 24.84600
## 73  -100.03250 25.24306
## 78  -106.82972 29.86667
## 79  -100.60167 23.43936
## 89  -100.38333 25.28378
## 90  -100.94389 23.34583
## 91  -102.06314 26.98422
## 102 -101.32361 27.06750
## 118 -102.30833 21.90833
## 121 -102.80833 25.34167
## 123 -103.81958 29.15542
## 125 -106.17783 28.84009
## 126 -101.44167 26.79167
## 132  -99.47500 21.50833

Ya tenemos las presencias particionadas en dos conjuntos (calibración o “training” y validación y “testing”). Para validar nuestros modelos necesitamos a fuerza ausencias. No tenemos ausencias, así que vamos a generar unas ausencias aleatorias a partir de la geografía: a este tipo de ausencias se les conoce como “pseudo-ausencias”

ext <- extent(clim.m)
backg <- randomPoints(clim.m, n=1000, ext=ext)
## Warning in couldBeLonLat(mask): CRS is NA. Assuming it is longitude/
## latitude
colnames(backg) = c('Longitude', 'Latitude')
group <- kfold(backg, 5)
backg_train <- backg[group != 1, ]
backg_test <- backg[group == 1, ]

Vamos a graficar los puntos

plot(clim.m[[1]])
points(backg_train, pch='-', cex=1, col='yellow')
points(backg_test, pch='-',  cex=1, col='black')
points(pres_train, pch= '+', col='green')
points(pres_test, pch='+', col='blue')

Ahora vamos a calibrar con el primer conjunto que extrajimos, pres_train:

bc <- bioclim(clim.m, pres_train)
plot(bc, a=1, b=2, p=0.85)

Vamos a validar nuestro modelo de Bioclim

# Usaremos la función evaluate de dismo.
eval.modesta <- evaluate(pres_test, backg_test, bc, clim.m)
eval.modesta
## class          : ModelEvaluation 
## n presences    : 27 
## n absences     : 200 
## AUC            : 0.7464815 
## cor            : 0.3104916 
## max TPR+TNR at : 0.01808182

Vamos a calcular un umbral

tr <- threshold(eval.modesta, 'spec_sens')
tr
## [1] 0.01808182

Vamos a hacer la predicción en la geografía:

pred.bioc.v3 <- predict(clim.m, bc, ext=ext, progress='')
pred.bioc.v3
## class       : RasterLayer 
## dimensions  : 1742, 2509, 4370678  (nrow, ncol, ncell)
## resolution  : 0.008333334, 0.008333334  (x, y)
## extent      : -115.61, -94.70169, 18.2046, 32.72126  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : layer 
## values      : 0, 0.7272727  (min, max)

Graficamos los dos mapas de predicciones, el continuo y el discreto usando el valor de umbral que calculamos antes.

par(mfrow=c(1,2))
plot(pred.bioc.v3, main='Bioclim, valores crudos')
plot(pred.bioc.v3 > tr, main='presencia/ausencia')
points(pres_train, pch='+')

Ahora vamos con Maxent, pero en la plataforma dedicada del programa.

Eso es todo!

Escrito por Julián A. Velasco, 29 de Julio de 2018.