SMD-Dismo

class: center, middle, inverse, title-slide

.title[
# SMD-Dismo
]
.author[
### Javier Montaño Chiriboga
]
.date[
### 2023-05-18
]

---

background-image: url(https://www.azquotes.com/picture-quotes/quote-all-models-are-wrong-but-some-are-useful-george-e-p-box-53-42-27.jpg)

---
class: center, middle

# Modelos de Distribusion de Especies (SDM)
## con MaXent
---

background-image: url(https://vignette.wikia.nocookie.net/jkenterprises/images/4/48/Golem%2Bof%2BPrague.jpg/revision/latest?cb=20170115084655)

---

background-image: url(https://i.pinimg.com/originals/88/ac/bf/88acbfec6866a3ac0765dfd296e1323c.jpg)

---

## Fuentes de Datos (Dataset)

- Collectas

- Datos de Biodiversidad(ocurrencias)

- Herbarios

- etc..

--
---

##Descarga desde **GBIF**

En esta ocasion vamos a obtener los datos desde los registros almacenados en el GBIF

. Una forma es descargar los datos directamente desde la pagina del GBIF, mendiante un identificardor de objeto digital (DOI)

. Una alternativa es descargar directamente de R desde el paquete *rgbif*

![](https://data-blog.gbif.org/post/2020-10-09-issues-and-flags_files/workflow1.png)

---

###Instalar **rgbif**

Para instalar el paquete se usa la funcion *install.packages* e.g `install.packages("rgbif")`

Para usar el paquete se debe previamente instalar, luego de la instalacion se carga el paquete con la funcion *library*

```r
library(rgbif) #para instalar el paquete install.package("rgbif")
```

---

![](https://docs.ropensci.org/rgbif/logo.png)

---

-Es importante tener cuenta en el gbif para usar sin limites rgbif

![](http://jrsbiodiversity.org/wp-content/uploads/2015/04/GBIF-logo-200x215.jpg)

---

-Para descargar es necesario encontrar el ID de la especie en la base de datos del GBIF, los ID se encuentran en la funcion `name_backbone_checklist()`

```r
spp <- c("Herpothallon rubrocinctum")
gb <- data.frame(Especie = spp)
keys <- 
        gb %>%
        pull("Especie") %>%
        name_backbone_checklist()  %>%
        filter(!matchType == "NONE") %>%
        pull(usageKey) #crea un archivo con los codigos del gbif
```

---

###Descargar datos de ocurrencias

Para descargar sin limite alguno se usa la funcion `occ_download()`, una alternativa a esta funcion es la `occ_search()` donde no es necesario colocar tu cuenta, pero tiene un limite de registros.

```r
keys1 <- occ_download(
        pred_in("taxonKey", keys),
        pred("country", "CO"), #registros de Colombia
        pred("hasCoordinate", TRUE),
        pred("occurrenceStatus","PRESENT"),
        format = "SIMPLE_CSV",
        user = cuenta, pwd = contras, email = correo
) #para descargar los registros

keys2 <- grep("^", keys1, value = T)

keys2
```

```
## [1] "0257166-230224095556074"
```

---

### Datos Espaciales

Mientras se descargan los datos vamos a manejar los datos espaciales.

---

Lo primero que vamos a hacer en conseguir el mapa de Colombia.

```r
colombia1 <- getData('GADM' , country = "COL", level = 1) #shape de Colombia
plot(colombia1)
```

---

datos ambientales que vamos a usar para el modelo

---

La informacion climatica se obtiene en formato raster, basicamente es una matrix, donde cada elemento o pixel tiene un dato y para graficarlo a cada dato se le asigna un color.

![](http://gsp.humboldt.edu/olm_2018/Lessons/GIS/04%20CreatingSpatialData/Images/rasterdatamodel.png)
---

Los raster pueden estar en juntos en diferentes capas o en solo una capa

![](https://raw.githubusercontent.com/NEONScience/NEON-Data-Skills/dev-aten/graphics/raster-general/single_multi_raster.png)

---

![](https://saylordotorg.github.io/text_essentials-of-geographic-information-systems/section_12/2a68ec6d5c214c2518a936e4abf14619.jpg) {width=50% height=50%}
b

---

luego procedemos a cortar los datos de raster para que tengan la misma extencion que nuestor mapa de Colombia.

para eso primero usamos la funcion `crop()` del paquete que corta el raster en un cuadrado de extencion cercana al mapa que tenemos.

Luego con la funcion `mask()` se corta el rater para que tenga una forma cercana a la de nuestro mapa

---

Para tener en cuenta todas las variables, estas se combinan en un raster mutiple.

```r
cb <- list(c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13, c14, c15, c16, c17, c18, c19, c20)

clima <- stack(cb) #se combinan todos los rasters

plot(clima)
```

![](diapo_files/figure-html/combinacion-1.png)

---

Ahora con el mapa debemos cargar los datos de las ocurrencias que descargamos del GBIF a R, para poder modelar la distribucion.

```r
d <- 
        occ_download_get("0252496-230224095556074") %>%
        occ_download_import() #archivo con los datos

d1 <- d[, c(10, 22, 23, 36)] #solo dejar el nombre, las coordenadas y el tipo de observacion
colnames(d1) <- c("especies", "lat", "lon", "basisOfRecord")
knitr::kable(head(d1), format = 'html')
```

<table>
 <thead>
  <tr>
   <th style="text-align:left;"> especies </th>
   <th style="text-align:right;"> lat </th>
   <th style="text-align:right;"> lon </th>
   <th style="text-align:left;"> basisOfRecord </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> Herpothallon rubrocinctum </td>
   <td style="text-align:right;"> 4.615950 </td>
   <td style="text-align:right;"> -74.31527 </td>
   <td style="text-align:left;"> HUMAN_OBSERVATION </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Herpothallon rubrocinctum </td>
   <td style="text-align:right;"> 4.065356 </td>
   <td style="text-align:right;"> -73.59437 </td>
   <td style="text-align:left;"> HUMAN_OBSERVATION </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Herpothallon rubrocinctum </td>
   <td style="text-align:right;"> 6.207769 </td>
   <td style="text-align:right;"> -75.55069 </td>
   <td style="text-align:left;"> HUMAN_OBSERVATION </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Herpothallon rubrocinctum </td>
   <td style="text-align:right;"> 4.780458 </td>
   <td style="text-align:right;"> -71.89709 </td>
   <td style="text-align:left;"> HUMAN_OBSERVATION </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Herpothallon rubrocinctum </td>
   <td style="text-align:right;"> 4.780630 </td>
   <td style="text-align:right;"> -71.89702 </td>
   <td style="text-align:left;"> HUMAN_OBSERVATION </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Herpothallon rubrocinctum </td>
   <td style="text-align:right;"> 4.842261 </td>
   <td style="text-align:right;"> -71.89049 </td>
   <td style="text-align:left;"> HUMAN_OBSERVATION </td>
  </tr>
</tbody>
</table>

---

En este caso vamos a hacer un filtrado leve de los datos, se van a eliminar los duplicados.

```r
d1 <- subset(d1,!is.na(lat) & !is.na(lon)) #filtrar datos sin coordenadas en x o en y
d1 <- data.table(d1) #cambia de data.frame a data.table
d1 <- d1[!especies == ""] #quitar los registros sin nobre cinetifico

pru <- with(d1, data.table(especies, cor = paste(lat, lon))) #la idea es crear pru para filtar duplicaods(registros de las mismas coordenads)
pru <- pru %>% group_by(especies) %>% duplicated(by = key(cor))
d2 <- d1[!pru] #filtro de duplicados
dp <- subset(d2, basisOfRecord == "PRESERVED_SPECIMEN") #filtro de especimes de herbario
with(dp, sort(table(especies))) #numero de registros por especie
```

```
## Herpothallon rubrocinctum 
##                        24
```

```r
spp <- with(d1, unique(especies))
```

---

Para proseguir se debe elegir el sistema de cordenadas en este caso vamos a elegir la WGS 84.

```r
ocu1 <- dp
datos1 <- ocu1[, -c(1, 4)]
ocu1 <- ocu1[,-1]
coordinates(datos1) <- ~lon + lat #espaciales para graficar
src <- CRS("+init=epsg:4326")
crs(datos1) <- src
mapa <- spTransform(colombia1, crs(datos1))
sobrelapan <- over(datos1, mapa)
plot(clima[[1]])
plot(datos1, add = T)
```

![](diapo_files/figure-html/crs-1.png)
---

Ademas de las cordenadas, la proyeccion que se use para visualizar las ocurrencias cobran relevancia al visualizar los resultados

```r
col <- gisco_get_countries(country = "Colombia")
col1 <- col %>% st_transform(crs = "+proj=robin")
col2 <- col %>% st_transform(crs = "+proj=lagrng")
col3 <- col %>% st_transform(crs = "+proj=merc")

pr <- ggplot(col1) + geom_sf(fill = "#078930", col = "white") + theme_minimal()
pl <- ggplot(col2) + geom_sf(fill = "#078930", col = "white") + theme_minimal()
pm <- ggplot(col3) + geom_sf(fill = "#078930", col = "white") + theme_minimal()
```
---
###Robinson

```r
pr
```

![](diapo_files/figure-html/robinson-1.png)

---
###Lagrange

```r
pl
```

![](diapo_files/figure-html/lagrange-1.png)

---
###Mercator

```r
pm
```

![](diapo_files/figure-html/mercator-1.png)
---

Con los datos de ocurrencias y las varibles a analizar podemos emepzar con el SMD.

![](https://www.thenakedscientists.com/sites/default/files/media/media/images/xSDM_framework.jpg.pagespeed.ic.VS6vs-TA0v.jpg)

---

Para esto primero creamos un bufer de 10km en esta ocasion sobre los registros obtenidos

```r
buffer <- buffer(datos1, width = 100000) #bufer de 10km
areas <- crop(clima, extent(buffer)) #un rectangulo del area de los puntos
areas <- mask(areas, buffer) #el area exacta que coge todos los registros
plot(areas[[1]])
```

![](diapo_files/figure-html/buffer-1.png)

---

Teneindo los registros y las variables elegidas, en el bufer, para crear "ruido" en el modelo vamos a crear una suerte de registros falsos es los alrededores de los registros del GBIF.

```r
u <- runif(1)
if (u < .6){bg <- sampleRandom(x = areas, size = 10000, na.rm = T, sp = T)} else {bg <- randomPoints(areas, 10000, p = datos1)} #hacen lo mismo pero diferente

plot(bg)
```

![](diapo_files/figure-html/ruido-1.png)

---

antes de correr el modelo vamos a separar los datos del GBIF en dos grupos, el primer  grupo de datos sera usado en nuestro modelo para el proceso de "*entrenamiento*", el otro grupo de datos que no va a ser tenido en cuenta en el modelo sera para "*testear*" el modelo mas adelante

---

En este caso vamos a separar nuestros datos en tres partes iguales y aleatoriamente se escogera una parte que sera la usada para "testear" y las otras dos en la parte del "entrenamiento"

![](diapo_files/figure-html/partes-1.png)

---

Un ultimo paso antes de correr el modelos, es extraer las condiciones ambientales de cada una de las variables evaluadas en los lugares de los registros. Luego se hace una matrix de prescencia-ausencia de cada una de las variables.

```r
p <- raster::extract(clima, train)
ptest <- raster::extract(clima, test)
a <- raster::extract(clima, bg)
pa <- c(rep(1, nrow(p)), rep(0, nrow(a))) #presente en el test y ausente en el entorno
pder <- as.data.frame(rbind(p, a))
knitr::kable(head(pder), format = 'html')
```

---

ya tenemos todo lo necesario para correr el modelo de distribucion de especies **SMD**, atraves del paquete de dismo con maxima entropia.

```r
mod <- dismo::maxent(x = pder, p = pa, factors = "biome", nbg = 10000, args = c("-J", "-P")) #modelo de distribucion de especies
response(mod) #Ponderacion de la variacion de cada variables
```

![](diapo_files/figure-html/modelo-1.png)

```r
ped1 <- predict(mod, clima) #Extrapolacion del modelo en todo Colombia
```

---

![](diapo_files/figure-html/modelo-grafica-1.png)

---

```r
mod <- dismo::maxent(x = pder, p = pa, factors = "biome", nbg = 10000, args = c("-J", "-P")) #modelo de distribucion de especies
response(mod) #Ponderacion de la variacion de cada variables
```

![](diapo_files/figure-html/testeo-1.png)

```r
ped1 <- predict(mod, clima) #Extrapolacion del modelo en todo Colombia
```

---

Ahora con el modelo hacemos la evaluacion o el "*testeo*" con los datos que separamos previamente

```r
p1 <- predict(mod, data.frame(ptest))
a1 <- predict(mod, data.frame(a))
combinado <- c(p1, a1)
label <- c(rep(1, length(p1)), rep(0, length(a1)))
predic <- prediction(combinado, label)
perfor <- performance(predic, "tpr","fpr") 
auc <- performance(predic, "auc") #auc del modelo
bdf2 <- data.frame("TPR" = perfor@y.values, "FPR" = perfor@x.values, "Cutoff" = perfor@alpha.values)
colnames(bdf2) <- c("TPR", "TFR", "Cutoff")
bdf2 <- data.table(bdf2)
AUC <- function(data = list(ptest, a), i) { 
                        p1 <- predict(mod, data.frame(ptest[i,]))
                        a1 <- predict(mod, data.frame(a))
                        combinado <- c(p1, a1)
                        label <- c(rep(1, length(p1)), rep(0, length(a1)))
                        predic <- prediction(combinado, label)
                        
                        auc <- performance(predic, "auc")
                        result <- list()
                        result[[1]] <- auc@y.values[[1]]
                        result[[2]] <- sd(perfor@y.values[[1]])
                        return(result[[1]])
}
                
bo <- boot(ptest, AUC, 500) #Bootrap del modelo
```

---

---
background-image: url(https://images.fineartamerica.com/images/artworkimages/mediumlarge/3/1-labor-omnia-vincit-improbus-vidddie-publyshd.jpg)

---

background-image: url(http://i2.cdn.turner.com/cnnnext/dam/assets/150126102625-01-auschwitz-liberation-0126-super-169.jpg)