Práctica Maxent
Fabián Santos, Phd.
27/11/2020
Proyecto Maxent para el modelamiento de nichos y distribución de especies
El software MaxEnt es una de las herramientas más populares para la distribución de especies y el modelado de nichos ambientales con más de 1000 solicitudes publicadas desde 2006 (Merow et al. 2013). Gracias a su implementación en R es posible automatizar su uso, así como facilitar la descarga de la información que este require para el modelado. Para esto requerimos las siguientes fuentes de información: * Información bioclimatica. Se la puede obtener de WorldClim https://www.worldclim.org/data/index.html y de la misión SRTM http://srtm.csi.cgiar.org/ pero la librería raster y su función getData hace posible su obtención desde R. * Registros de especies. Existen varias fuentes, e.g. IUCN Red List Toolbox, eBIRD, GBIF, entre otros. Aquí usaremos datos descargados del GBIF https://www.gbif.org/
#Librerías requeridas
library(raster) ## Loading required package: splibrary(dismo)
library(rJava)
library(rgeos)## rgeos version: 0.5-5, (SVN revision 640)
## GEOS runtime version: 3.8.0-CAPI-1.13.1
## Linking to sp version: 1.4-4
## Polygon checking: TRUEArea de estudio
Primero obtenemos la respectiva a los límites administrativos del Ecuador. Notese el uso de “ECU” para refererirse al país y el nivel 1 administrativo (provincias).
#limites administrativos para Ecuador
study.area <- getData("GADM",country="ECU",level=1)
#para extraer solo la región amazonica
head(study.area)## GID_0 NAME_0 GID_1 NAME_1 VARNAME_1 NL_NAME_1 TYPE_1 ENGTYPE_1
## 1 ECU Ecuador ECU.1_1 Azuay <NA> <NA> Provincia Province
## 12 ECU Ecuador ECU.2_1 Bolivar <NA> <NA> Provincia Province
## 18 ECU Ecuador ECU.3_1 Cañar <NA> <NA> Provincia Province
## 19 ECU Ecuador ECU.4_1 Carchi <NA> <NA> Provincia Province
## 20 ECU Ecuador ECU.5_1 Chimborazo <NA> <NA> Provincia Province
## 21 ECU Ecuador ECU.6_1 Cotopaxi <NA> <NA> Provincia Province
## CC_1 HASC_1
## 1 01 EC.AZ
## 12 02 EC.BO
## 18 03 EC.CN
## 19 04 EC.CR
## 20 06 EC.CB
## 21 05 EC.CTstudy.area <- study.area[study.area$NAME_1 %in% c("Orellana","Pastaza","Sucumbios","Morona Santiago","Napo","Zamora Chinchipe"),]
plot(study.area)
Obtención de información sobre altitud
Procedemos ahora a descargar la información de altura. Para esto debemos especificar los puntos referenciales de donde se encuentra nuestra area de estudio. Aquí tomamos los centroides y seleccionamos 2, esto porque la información del SRTM viene recordata en cuadrantes
#derivación de los centroides de las provincias
centroids <- gCentroid(study.area,byid=T)
plot.new()
plot(study.area)
plot(centroids,add=T,col="red")
#observar orden de los puntos y probar con el que esta más al norte
study.area$NAME_1## [1] "Morona Santiago" "Napo" "Orellana" "Pastaza"
## [5] "Sucumbios" "Zamora Chinchipe"pts.ref <- coordinates(centroids[5])[1,]
#obtener SRTM para este y observar
alt.1 <- getData("SRTM",lon=pts.ref[1],lat=pts.ref[2])
plot.new()
plot(alt.1)
plot(study.area,add=T)
plot(centroids[5],add=T,col="red")
#como todavia no cubre asigno arbitrariamente un punto más al norte
pts.ref[2] <- 1
alt.2 <- getData("SRTM",lon=pts.ref[1],lat=pts.ref[2])
plot.new()
plot(alt.2)
plot(study.area,add=T) 
#y junto las capas del SRTM
alt <- list(alt.1,alt.2)
alt <- do.call(merge,alt)
#y lo corto con el area de estudio
alt <- crop(alt,extent(study.area))
plot.new()
plot(alt)
plot(study.area,add=T) 
Obtención de información bioclimática
Usando un procedimiento similar, ahora descargamos la información bioclimática, es decir la del WorldClim. Usaremos en este caso el escenario presente. Considerar que existen la siguientes: BIO1 = Annual Mean Temperature; BIO2 = Mean Diurnal Range (Mean of monthly (max temp – min temp)); BIO3 = Isothermality (BIO2/BIO7) (x 100); BIO4 = Temperature Seasonality (standard deviation x100); BIO5 = Max Temperature of Warmest Month; BIO6 = Min Temperature of Coldest Month; BIO7 = Temperature Annual Range (BIO5-BIO6); BIO8 = Mean Temperature of Wettest Quarter; BIO9 = Mean Temperature of Driest Quarter; BIO10 = Mean Temperature of Warmest Quarter; BIO11 = Mean Temperature of Coldest Quarter; BIO12 = Annual Precipitation; BIO13 = Precipitation of Wettest Month; BIO14 = Precipitation of Driest Month; BIO15 = Precipitation Seasonality (Coefficient of Variation); BIO16 = Precipitation of Wettest Quarter; BIO17 = Precipitation of Driest Quarter; BIO18 = Precipitation of Warmest Quarter; and BIO19 = Precipitation of Coldest Quarter.
#recordar que ya calculamos los centroides y derivamos aqui el que esta mas al norte
pts.ref <- coordinates(centroids[5])[1,]
#descargamos la informacion del WorldClim
bio.1 <- getData("worldclim", var='bio', res=0.5,lon=pts.ref[1],lat=pts.ref[2])
#recordar que son varias capas de informacion pero para comprobar el area, ploteamos solo la primera
plot.new()
plot(bio.1[[1]])
plot(study.area,add=T) 
#como todavia no cubre asigno arbitrariamente un punto más al norte
pts.ref[2] <- 1
bio.2 <- getData("worldclim", var='bio', res=0.5,lon=pts.ref[1],lat=pts.ref[2])
plot.new()
plot(bio.2[[1]])
plot(study.area,add=T) 
#y junto las capas del WorldClim
bio <- list(bio.1,bio.2)
bio <- do.call(merge,bio)
#lo corto con el area de estudio
bio <- crop(bio,extent(study.area))
plot.new()
plot(bio[[1]])
plot(study.area,add=T)
#y ajusto la resolucion de la altitud (es mayor que la de bioclima) para juntarlos en una sola pila
alt <- raster::resample(x=alt, y=bio)
predictors <- stack(bio,alt)
plot.new()
plot(predictors)
Informacion sobre la especie objetivo
El archivo CSV usado aquí, fue descargado desde el GBIF. Aquí usaremos la del mono araña (ateles belzebuth)
#ahora observemos los datos de especies y su resumen
spp <- read.delim("E:/DATA/maxent/spp/0122840-200613084148143.csv",header=T)
summary(spp)## gbifID datasetKey occurrenceID kingdom
## Min. :1.970e+08 Length:41 Length:41 Length:41
## 1st Qu.:1.212e+09 Class :character Class :character Class :character
## Median :2.244e+09 Mode :character Mode :character Mode :character
## Mean :1.896e+09
## 3rd Qu.:2.457e+09
## Max. :2.863e+09
##
## phylum class order family
## Length:41 Length:41 Length:41 Length:41
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## genus species infraspecificEpithet taxonRank
## Length:41 Length:41 Length:41 Length:41
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## scientificName verbatimScientificName verbatimScientificNameAuthorship
## Length:41 Length:41 Length:41
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## countryCode locality stateProvince occurrenceStatus
## Length:41 Length:41 Length:41 Length:41
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## individualCount publishingOrgKey decimalLatitude decimalLongitude
## Min. :1 Length:41 Min. :-4.2830 Min. :-78.91
## 1st Qu.:1 Class :character 1st Qu.:-1.1368 1st Qu.:-77.50
## Median :1 Mode :character Median :-0.8807 Median :-76.51
## Mean :1 Mean :-1.1319 Mean :-76.86
## 3rd Qu.:1 3rd Qu.:-0.5874 3rd Qu.:-76.14
## Max. :1 Max. :-0.2172 Max. :-75.69
## NA's :35 NA's :11 NA's :11
## coordinateUncertaintyInMeters coordinatePrecision elevation
## Min. : 63.2 Mode:logical Min. :900
## 1st Qu.: 31448.0 NA's:41 1st Qu.:900
## Median : 31449.5 Median :900
## Mean : 33158.8 Mean :900
## 3rd Qu.: 31451.0 3rd Qu.:900
## Max. :110737.0 Max. :900
## NA's :13 NA's :40
## elevationAccuracy depth depthAccuracy eventDate
## Mode:logical Mode:logical Mode:logical Length:41
## NA's:41 NA's:41 NA's:41 Class :character
## Mode :character
##
##
##
##
## day month year taxonKey
## Min. : 1.0 Min. : 1.0 Min. :1810 Min. :2436642
## 1st Qu.: 4.5 1st Qu.: 2.5 1st Qu.:2002 1st Qu.:2436642
## Median :20.0 Median : 3.0 Median :2016 Median :2436642
## Mean :17.2 Mean : 5.2 Mean :1992 Mean :2763416
## 3rd Qu.:26.0 3rd Qu.: 7.5 3rd Qu.:2019 3rd Qu.:2436642
## Max. :31.0 Max. :12.0 Max. :2020 Max. :5786076
## NA's :6 NA's :6 NA's :5
## speciesKey basisOfRecord institutionCode collectionCode
## Min. :2436642 Length:41 Length:41 Length:41
## 1st Qu.:2436642 Class :character Class :character Class :character
## Median :2436642 Mode :character Mode :character Mode :character
## Mean :2436642
## 3rd Qu.:2436642
## Max. :2436642
##
## catalogNumber recordNumber identifiedBy dateIdentified
## Length:41 Length:41 Length:41 Length:41
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## license rightsHolder recordedBy typeStatus
## Length:41 Length:41 Length:41 Mode:logical
## Class :character Class :character Class :character NA's:41
## Mode :character Mode :character Mode :character
##
##
##
##
## establishmentMeans lastInterpreted mediaType issue
## Length:41 Length:41 Length:41 Length:41
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## #Observar que en las coordenadas hay 11 valores no existentes (NA). Tomamos solo los validos
spp <- spp[!is.na(spp$decimalLatitude),]
#lo asignamos como una capa espacial de puntos definiendo sus coordenadas como geograficas
spp <- SpatialPointsDataFrame(spp[,c("decimalLongitude","decimalLatitude")],spp,proj4string=study.area@proj4string)
plot.new()
plot(study.area)
plot(spp,add=T,col="green")
Definición de muestras, ejecución del modelo y reporte
Para entrenar el modelo y poder validarlo dividimos en 80% entrenamiento; y 20% validacion. Posterimente se ejecuta el maxent. El resultado es el modelo cartografico de la distribución de esta especie
#para entrenar y validar el modelo, primero divido las muestras (80/20)
fold <- kfold(spp, k=5)
occtest <- spp[fold == 1, ]
occtrain <- spp[fold != 1, ]
#correr maxent y observar reporte
me <- maxent(predictors, occtrain)## This is MaxEnt version 3.4.1me## class : MaxEnt
## variables: layer.1 layer.2 layer.3 layer.4 layer.5 layer.6 layer.7 layer.8 layer.9 layer.10 layer.11 layer.12 layer.13 layer.14 layer.15 layer.16 layer.17 layer.18 layer.19 layer#predecir y mapear
r <- predict(me, predictors)## This is MaxEnt version 3.4.1plot.new()
plot(r)
plot(study.area,add=T)
plot(spp,col="red",add=T)
Validación y almacenamiento del modelo
Finalmente validamos y guardamos el resultado en un archivo TIF y shapefile para observarlos con un software como QGIS
#dado que el modelo contrasta con presencias de fondo (background), creamos un set de puntos al azar sobre el area de estudio
set.seed(666)
bck <- spsample(study.area,500,"random")## Warning in proj4string(obj): CRS object has comment, which is lost in outputplot.new()
plot(study.area)
plot(bck,add=T,col="blue")
plot(spp,add=T,col="green")
#ya que algunos puntos estan cerca de los sitios de observaciones de la especie, eliminamos estos para distancias <10km
spp.buff <- buffer(spp,10000)
bck <- gDifference(bck,spp.buff)
plot.new()
plot(study.area)
plot(bck,add=T,col="blue")
plot(spp,add=T,col="green")
#validar
e1 <- evaluate(me, p=occtest, a=bck, x=predictors)## This is MaxEnt version 3.4.1
## This is MaxEnt version 3.4.1e1## class : ModelEvaluation
## n presences : 6
## n absences : 473
## AUC : 0.6180409
## cor : 0.04352829
## max TPR+TNR at : 0.4309942#grabar modelo para observarlo en QGIS
output.name <- "E:/DATA/maxent/modelo.tif"
writeRaster(r,output.name,overwrite=T)
output.name <- "E:/DATA/maxent/spp.shp"
shapefile(spp,output.name,overwrite=T)## Warning in rgdal::writeOGR(x, filename, layer, driver = "ESRI Shapefile", :
## Field names abbreviated for ESRI Shapefile driver