GBIF_Bioclim
Julián A. Velasco
11/16/2020
Descarga de datos de GBIF
Este ejercicio estará enfocado en descargar datos del repositorio GBIF y usar algunas herramientas sencillas para la limpieza de datos.
# Instalar los paquetes (dismo, raster, sp, spThin, rgdal)
library(sp)
library(raster)
library(dismo)
library(spThin)## Loading required package: spam## Loading required package: dotCall64## Loading required package: grid## Spam version 2.5-1 (2019-12-12) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.##
## Attaching package: 'spam'## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve## Loading required package: fields## See https://github.com/NCAR/Fields for
## an extensive vignette, other supplements and source code## Loading required package: knitrlibrary(rgdal)## rgdal: version: 1.5-18, (SVN revision 1082)
## 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: /usr/share/gdal/2.2
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ shared files: (autodetected)
## Linking to sp version:1.4-4Buscaremos datos de presencia para una especie de lagartija del desierto (Phrynosoma modestum)
sp.gbif <- dismo::gbif(genus="Phrynosoma", species="modestum", geo=TRUE, removeZeros=TRUE, download=TRUE)## Loading required namespace: jsonlite## 3853 records found## 0-300-600-900-1200-1500-1800-2100-2400-2700-3000-3300-3600-3853 records downloadedVamos a indexar solos los datos de las coordenadas geográficas. Debemos saber en que posición están.
names(sp.gbif)## [1] "acceptedNameUsage"
## [2] "acceptedScientificName"
## [3] "acceptedTaxonKey"
## [4] "accessRights"
## [5] "adm1"
## [6] "adm2"
## [7] "associatedReferences"
## [8] "associatedSequences"
## [9] "basisOfRecord"
## [10] "bibliographicCitation"
## [11] "catalogNumber"
## [12] "class"
## [13] "classKey"
## [14] "cloc"
## [15] "collectionCode"
## [16] "collectionID"
## [17] "continent"
## [18] "coordinateUncertaintyInMeters"
## [19] "country"
## [20] "crawlId"
## [21] "created"
## [22] "datasetID"
## [23] "datasetKey"
## [24] "datasetName"
## [25] "dateIdentified"
## [26] "day"
## [27] "depth"
## [28] "depthAccuracy"
## [29] "disposition"
## [30] "dynamicProperties"
## [31] "earliestAgeOrLowestStage"
## [32] "earliestEonOrLowestEonothem"
## [33] "earliestEpochOrLowestSeries"
## [34] "earliestPeriodOrLowestSystem"
## [35] "elevation"
## [36] "elevationAccuracy"
## [37] "endDayOfYear"
## [38] "establishmentMeans"
## [39] "eventDate"
## [40] "eventRemarks"
## [41] "eventTime"
## [42] "family"
## [43] "familyKey"
## [44] "fieldNotes"
## [45] "fieldNumber"
## [46] "formation"
## [47] "fullCountry"
## [48] "gbifID"
## [49] "genericName"
## [50] "genus"
## [51] "genusKey"
## [52] "geodeticDatum"
## [53] "georeferencedBy"
## [54] "georeferencedDate"
## [55] "georeferenceProtocol"
## [56] "georeferenceRemarks"
## [57] "georeferenceSources"
## [58] "georeferenceVerificationStatus"
## [59] "habitat"
## [60] "higherClassification"
## [61] "higherGeography"
## [62] "higherGeographyID"
## [63] "http://rs.tdwg.org/dwc/terms/organismQuantity"
## [64] "http://rs.tdwg.org/dwc/terms/organismQuantityType"
## [65] "http://unknown.org/language"
## [66] "http://unknown.org/nick"
## [67] "http://unknown.org/occurrenceDetails"
## [68] "http://unknown.org/recordId"
## [69] "http://unknown.org/rights"
## [70] "http://unknown.org/rightsHolder"
## [71] "identificationID"
## [72] "identificationQualifier"
## [73] "identificationRemarks"
## [74] "identificationVerificationStatus"
## [75] "identifiedBy"
## [76] "identifier"
## [77] "individualCount"
## [78] "informationWithheld"
## [79] "installationKey"
## [80] "institutionCode"
## [81] "institutionID"
## [82] "ISO2"
## [83] "key"
## [84] "kingdom"
## [85] "kingdomKey"
## [86] "language"
## [87] "lastCrawled"
## [88] "lastInterpreted"
## [89] "lastParsed"
## [90] "lat"
## [91] "latestAgeOrHighestStage"
## [92] "latestEpochOrHighestSeries"
## [93] "latestPeriodOrHighestSystem"
## [94] "license"
## [95] "lifeStage"
## [96] "locality"
## [97] "locationAccordingTo"
## [98] "locationID"
## [99] "locationRemarks"
## [100] "lon"
## [101] "modified"
## [102] "month"
## [103] "municipality"
## [104] "nameAccordingTo"
## [105] "namePublishedIn"
## [106] "namePublishedInYear"
## [107] "nomenclaturalCode"
## [108] "occurrenceID"
## [109] "occurrenceRemarks"
## [110] "occurrenceStatus"
## [111] "order"
## [112] "orderKey"
## [113] "organismID"
## [114] "organismQuantity"
## [115] "organismQuantityType"
## [116] "otherCatalogNumbers"
## [117] "ownerInstitutionCode"
## [118] "parentNameUsage"
## [119] "phylum"
## [120] "phylumKey"
## [121] "preparations"
## [122] "previousIdentifications"
## [123] "projectId"
## [124] "protocol"
## [125] "publishingCountry"
## [126] "publishingOrgKey"
## [127] "recordedBy"
## [128] "recordNumber"
## [129] "references"
## [130] "reproductiveCondition"
## [131] "rights"
## [132] "rightsHolder"
## [133] "samplingProtocol"
## [134] "scientificName"
## [135] "scientificNameID"
## [136] "sex"
## [137] "species"
## [138] "speciesKey"
## [139] "specificEpithet"
## [140] "startDayOfYear"
## [141] "taxonConceptID"
## [142] "taxonID"
## [143] "taxonKey"
## [144] "taxonomicStatus"
## [145] "taxonRank"
## [146] "taxonRemarks"
## [147] "type"
## [148] "typeStatus"
## [149] "typifiedName"
## [150] "verbatimCoordinateSystem"
## [151] "verbatimElevation"
## [152] "verbatimEventDate"
## [153] "verbatimLocality"
## [154] "verbatimSRS"
## [155] "vernacularName"
## [156] "year"
## [157] "downloadDate"sp.coords <- na.omit(sp.gbif[,c(100, 90, 137)]) # lon, lat, species
head(sp.coords)## lon lat species
## 1 -101.0339 25.44163 Phrynosoma modestum
## 2 -100.9753 23.69214 Phrynosoma modestum
## 3 -106.8000 33.55111 Phrynosoma modestum
## 4 -103.8765 26.62552 Phrynosoma modestum
## 5 -109.5921 32.55360 Phrynosoma modestum
## 6 -103.7489 26.67445 Phrynosoma modestumGraficamos
plot(sp.coords$lon, sp.coords$lat)
Aquí vemos que hay muchos puntos muy agrupados. Este tipo de datos puede influir negativamente en la calibración del modelo de nicho y la predicción geográfica. Vamos a tratar de reducir ese efecto adelgazando la muestra de puntos. Lo primero que hacemos es tratar de eliminar puntos a una distancia menor de 25 km. Como tenemos un tamaño de muestra grande (> 3000 registros) podemos hacer esto sin mucho problema. Usaremos el paquete spThin para esto.
sp.thin <-thin(sp.coords, lat.col="lat", long.col="lon", spec.col="species",
thin.par=25, reps=1, locs.thinned.list.return=TRUE, write.files=TRUE,
max.files=1, out.dir=".", out.base = "modestum.train.25km", write.log.file=TRUE, log.file="modestum.log.txt",
verbose=TRUE)## **********************************************
## Beginning Spatial Thinning.
## Script Started at: Thu Nov 19 07:44:44 2020
## lat.long.thin.count
## 397
## 1
## [1] "Maximum number of records after thinning: 397"
## [1] "Number of data.frames with max records: 1"
## [1] "Writing new *.csv files"## Warning in thin(sp.coords, lat.col = "lat", long.col = "lon", spec.col =
## "species", : Writing new *.csv files to output directory: .## [1] "Writing file: ./modestum.train.25km_thin1.csv"Ahora vamos primero a identificar que tipo de objeto se acaba de generar:
class(sp.thin)## [1] "list"str(sp.thin)## List of 1
## $ :'data.frame': 397 obs. of 2 variables:
## ..$ Longitude: num [1:397] -102 -106 -106 -106 -101 ...
## ..$ Latitude : num [1:397] 25.7 31.9 33.8 35.6 26.2 ...head(sp.thin[[1]])## Longitude Latitude
## 10 -102.3471 25.71672
## 20 -106.0402 31.92565
## 23 -106.0220 33.76997
## 57 -106.3340 35.64416
## 58 -100.8209 26.15840
## 63 -107.5259 33.20370Tenemos que tener en cuenta como vienen organizados los datos para luego manipularlos como queramos. Vamos a graficar los puntos adelgazados a una distancia no menor de 25 km.
plot(sp.thin[[1]]$Longitude, sp.thin[[1]]$Latitude)
Si comparamos con el objeto que contenia todos los puntos podemos ver que hemos reducido el sobre-muestreo en algunas regiones y esto nos facilitará la calibración y posterior predicción de nuestro modelo de nicho. Vamos a ver donde caen nuestros puntos dentro de un mapa.
Vamos a convertir nuestros puntos en un objeto espacial.
coords <- sp.thin[[1]]
coordinates(coords)<- ~Longitude+Latitude
plot(coords)
Si solo nos interesan los puntos que se distribuyen en México, podríamos excluir los puntos que caen por fuera del país. Vamos a hacerlo. Primero necesitamos cargar el mapa de México.
library(rgdal)mx <- readOGR(dsn="/media/julian/Behemoth/Behemoth/Dropbox/Taller_R/GBIF_Bioclim/Poligono_Mexico", layer="gadm36_MEX_1")## OGR data source with driver: ESRI Shapefile
## Source: "/media/julian/Behemoth/Behemoth/Dropbox/Taller_R/GBIF_Bioclim/Poligono_Mexico", layer: "gadm36_MEX_1"
## with 32 features
## It has 10 fieldsplot(mx)
Ahora le ponemos los puntos encima al mapa…
plot(coords)
plot(mx, add=TRUE)
Con un clip vamos a quitar los puntos que caen por fuera. Pero primero definimos el sistema de coordenadas para nuestros puntos
mx## class : SpatialPolygonsDataFrame
## features : 32
## extent : -118.3689, -86.71014, 14.53292, 32.71804 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 10
## names : GID_0, NAME_0, GID_1, NAME_1, VARNAME_1, NL_NAME_1, TYPE_1, ENGTYPE_1, CC_1, HASC_1
## min values : MEX, Mexico, MEX.1_1, Aguascalientes, NA, NA, Distrito Federal, Federal District, NA, MX.AG
## max values : MEX, Mexico, MEX.9_1, Zacatecas, NA, NA, Estado, State, NA, MX.ZAcrdref <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
projection(coords)<-crdref
coords## class : SpatialPoints
## features : 397
## extent : -118.1731, -95.86836, 21.20054, 44.71625 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0coords.clip <- coords[mx, ]
coords## class : SpatialPoints
## features : 397
## extent : -118.1731, -95.86836, 21.20054, 44.71625 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0coords.clip## class : SpatialPoints
## features : 169
## extent : -112.6083, -97.70008, 21.20054, 31.775 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0Graficar
plot(coords.clip)
plot(mx, add=TRUE)
Ahora tenemos nuestros puntos únicamente para México y los vamos a convertir en un objeto tipo data.frame con el fin de guardarlos en formato CSVy explorarlos con excel.
coords.clipDF <- as.data.frame(coords.clip)
head(coords.clipDF)## Longitude Latitude
## 10 -102.3471 25.71672
## 58 -100.8209 26.15840
## 148 -102.2906 21.88374
## 211 -101.7094 22.60750
## 214 -101.7068 27.05868
## 235 -101.1094 22.76387write.csv(coords.clipDF, file="Coordenadas_GBIF_depuradas.csv")Modelos de nicho ecológico
BIOCLIM
library(rgdal)
library(sp)
library(raster)
library(rgeos)## rgeos version: 0.5-5, (SVN revision 640)
## GEOS runtime version: 3.6.2-CAPI-1.10.2
## Linking to sp version: 1.4-4
## Polygon checking: TRUElibrary(dismo)
coords <- read.csv("Coordenadas_GBIF_depuradas.csv")
head(coords)## X Longitude Latitude
## 1 10 -102.3471 25.71672
## 2 58 -100.8209 26.15840
## 3 148 -102.2906 21.88374
## 4 211 -101.7094 22.60750
## 5 214 -101.7068 27.05868
## 6 235 -101.1094 22.76387Le quitamos la primera columna, no la necesitamos
coords <- coords[,c(2:3)]
head(coords)## Longitude Latitude
## 1 -102.3471 25.71672
## 2 -100.8209 26.15840
## 3 -102.2906 21.88374
## 4 -101.7094 22.60750
## 5 -101.7068 27.05868
## 6 -101.1094 22.76387Subimos las capas climáticas
clim.act <- list.files("/media/julian/Behemoth/Behemoth/Dropbox/Taller_R/GBIF_Bioclim/Capas_clima", pattern="asc", full.names=T)
clim.act <- stack(clim.act) # hacemos un stack
clim.act## class : RasterStack
## dimensions : 2183, 3799, 8293217, 6 (nrow, ncol, ncell, nlayers)
## resolution : 0.008333334, 0.008333334 (x, y)
## extent : -118.3684, -86.71002, 14.5296, 32.72126 (xmin, xmax, ymin, ymax)
## crs : NA
## names : bio1, bio12, bio15, bio4, bio5, bio6plot(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)
# ajustamos un modelo de bioclim a los datos con las capas climáticas
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.06508876 0.43786982 0.13609467 0.06508876 0.14201183 0.07692308a <- min(coords.bioc)
a## [1] 0.00591716Usaremos 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)
## crs : NA
## source : 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 (e.g., usar valores de 5 ó 6 para generar el buffer). 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] 4 4 1 2 5 1 2 1 2 4 2 3 4 5 5 2 2 2 2 4 2 4 4 3 1 3 2 4 4 5 5 3 4 5 5 5 2
## [38] 2 2 1 3 2 2 3 2 2 4 4 2 4 3 1 2 5 3 3 5 5 3 2 5 1 1 3 1 2 5 5 4 1 3 4 1 4
## [75] 1 5 3 3 3 1 1 1 3 3 4 1 5 4 4 4 3 4 5 5 4 3 4 5 3 3 1 2 2 3 1 4 3 5 1 1 3
## [112] 2 5 5 2 1 3 4 5 3 3 5 4 5 5 2 4 5 1 4 4 3 4 2 2 2 3 1 3 5 2 1 2 5 5 4 1 5
## [149] 1 2 1 3 1 1 2 4 1 5 4 5 3 4 3 1 2 5 1 1 1length(group)## [1] 169pres_train <- coordsDF[group != 1, ]
pres_test <- coordsDF[group == 1, ]
head(pres_train)## Longitude Latitude
## 1 -102.3471 25.71672
## 2 -100.8209 26.15840
## 4 -101.7094 22.60750
## 5 -101.7068 27.05868
## 7 -101.0397 23.55722
## 9 -101.5221 26.31369head(pres_test)## Longitude Latitude
## 3 -102.2906 21.88374
## 6 -101.1094 22.76387
## 8 -106.4816 31.57700
## 25 -100.7743 23.82274
## 40 -104.0409 24.79429
## 52 -104.4999 25.18843Ya tenemos las presencias particionadas en dos conjuntos (calibración o “training” y validación o “testing”). Para validar nuestros modelos necesitamos las “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(x, warnings = warnings): CRS is NA. Assuming it is
## longitude/latitudecolnames(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 : 34
## n absences : 200
## AUC : 0.8193382
## cor : 0.3066443
## max TPR+TNR at : 0.05915926Vamos a calcular un umbral
tr <- threshold(eval.modesta, 'spec_sens')
tr## [1] 0.05915926Vamos 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)
## crs : NA
## source : memory
## names : layer
## values : 0, 0.8222222 (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='+')
Eso es todo!