05_Overview.Rmd

# Descripción general
Ejemplo de flujo de análisis que demuestra brevemente algunas de las tareas que abordaremos durante el taller.  
  
## Preparación
Asegurémonos de que tenemos los paquetes que necesitamos, sino los instalamos desde CRAN:  
  
```{r pkgs1, message = FALSE}
pkgs <- c("dplyr", "ncdf4", "terra", "remotes", "robis", "worrms")
pkgs <- setdiff(pkgs, installed.packages()[, 1])
if (length(pkgs) > 0) install.packages(pkgs)
```

También podemos instalarlos desde GitHub:

```{r pkgs2}
## paquetes necesarios
pkgs <- c("SCAR/antanym" = NA, "AustralianAntarcticDivision/blueant" = NA,
          "AustralianAntarcticDivision/SOmap" =  "0.5.1")
for (pkg in names(pkgs)) {
    if (!basename(pkg) %in% installed.packages()[, 1] ||
        (!is.na(pkgs[[pkg]]) && packageVersion(basename(pkg)) < pkgs[[pkg]])) {
        remotes::install_github(pkg)
    }
}

remotes::install_github("ropensci/antanym")
```

## Llamando paquetes relevantes
  
```{r}
#Manipulacion de datos
library(dplyr)
#Taxonomia
library(worrms)
#Datos de ocurrencia desde OBIS
library(robis)
#Graficando mapas
library(SOmap)
#Datos ambientales
library(blueant)
#Modelos aditivos generalizados
library(mgcv)
#Nombres de lugares
library(antanym)
```


## Taxonomía

```{r tax1, message = FALSE}
#Especie de interes
especie_interes <- "Euphausia crystallorophias"

#Revisar records en WoRMS
tax <- wm_records_names(name = especie_interes)

# Aphia ID (identificación taxonómica) de nuestras especies
my_aphia_id <- tax |> 
  as.data.frame() |> 
  pull(valid_AphiaID)
```

## Occurrencias

Obtengamos todos los datos disponibles de la familia *Euphausiacea* de Antártica desde el conjunto de datos obtenidos en las expediciones "recursos marinos antárticos alemanes": [Antarctic Euphausiacea occurence data from "German Antarctic Marine Living Resources" (GAMLR) Expeditions](https://ipt.biodiversity.aq/resource?r=gamlr)  

```{r occ1, message = FALSE}
x <- occurrence(datasetid = "cb16377b-56a8-4d95-802d-4eec02466773")
```

Creemos el gráfico de la distribución completa de las muestras de la familia *Euphasia* en negro y de *Euphausia crystallorophias* en verde:

```{r map1, message = FALSE}
SOmap_auto(x$decimalLongitude, x$decimalLatitude, input_lines = FALSE, 
           pcol = "black", pcex = 0.2)
with(x |>
       dplyr::filter(aphiaID == my_aphia_id), 
     SOplot(decimalLongitude, decimalLatitude, pch = 19, cex = 0.2, 
            col = "green"))
```
  
También podemos crear un mapa polar estereográfico:
  
```{r map2, message = FALSE}
basemap <- SOmap(trim = ceiling(max(x$decimalLatitude))+1, 
                 bathy_legend = FALSE)
plot(basemap)
SOplot(x$decimalLongitude, x$decimalLatitude, pch = 19, cex = 0.2, 
       col = "black")
with(x |>
       dplyr::filter(aphiaID == my_aphia_id), 
     SOplot(decimalLongitude, decimalLatitude, pch = 19, cex = 0.2, 
            col = "green"))
```
  
Vamos a ajustar un modelo de distribución de especies (presencia/ausencia), así que primero reorganizaremos nuestros datos de presencia/ausencia por sitio de muestreo:  
  
```{r dat1}
xfit <- x |>
  dplyr::rename(lon = "decimalLongitude", lat = "decimalLatitude") |>
  group_by(lon, lat) |>
  dplyr::summarize(present = any(my_aphia_id %in% aphiaID))
```

## Datos ambientales

```{r envdat1, message = FALSE}
## añadimos los datos en un directorio temporal
my_data_directory <- tempdir()

## agregamos el origen de los datos que queremos
data_source <- sources_sdm("Southern Ocean marine environmental data")

## buscamos los datos
status <- bb_get(data_source, local_file_root = my_data_directory, 
                 verbose = TRUE)
```

```{r envdat2, message = FALSE}
nc_files <- Filter(function(z) grepl("\\.nc$", z), status$files[[1]]$file)

## creamos un raster apilado con todas las capas
env_stack <- rast(nc_files)

## veamos las primeras capas
head(names(env_stack))
```

Seleccionemos solo las capas de `depth` (Profundidad) y `ice_cover_mean` (cobertura de hielo promedio); y extraigamos sus valores en nuestros sitios de muestreo:

```{r envdat3}
env_stack <- subset(env_stack, c("depth", "ice_cover_mean"))

temp <- as.data.frame(terra::extract(env_stack, xfit[, c("lon", "lat")]))

xfit <- bind_cols(xfit, temp)

head(xfit)
```

## Ajustando el modelo

Tenemos datos de presencia/ausencia, así que ajustaremos un modelo binomial simple. La probabilidad de presencia de *Euphausia crystallorophias* se ajusta como una función suave de la profundidad (smooth function; s()) y la cobertura media de hielo marino:

```{r mod1, message = FALSE}
fit <- gam(present ~ s(depth) + s(ice_cover_mean), 
           family = binomial,
           data = xfit)

summary(fit)
```

Los ajustes para la profundidad y la cobertura de hielo:

```{r mod2, message = FALSE}
plot(fit, pages = 1, shade = TRUE)
```

Esto sugiere que es más probable encontrar *Euphausia crystallorophias* en áreas poco profundas con alta cobertura anual de hielo, lo cual parece plausible dado que esta especie se encuentra típicamente en aguas costeras y de la plataforma continental.

## Predicciones basadas en el modelo

```{r pred1}
xpred <- expand.grid(lon = seq(from = floor(min(xfit$lon)), 
                               to = ceiling(max(xfit$lon)), by = 0.25),
                     lat = seq(from = floor(min(xfit$lat)), 
                               to = ceiling(max(xfit$lat)), by = 0.25))

xpred <- bind_cols(as.data.frame(xpred), 
                   as.data.frame(terra::extract(env_stack, 
                                                xpred[, c("lon", "lat")])))

xpred$predicted <- predict(fit, newdata = xpred, type = "response")

## create raster
pr <- rast(xpred[, c("lon", "lat", "predicted")], type = "xyz")
projection(pr) <- "+proj=longlat +datum=WGS84"

my_cmap <- if(getRversion() >= "3.6"){
  hcl.colors(21, "Geyser")
  }else{
    c("#008585", "#359087", "#539B8A", "#6DA590", "#85AF97", "#9BBAA0",
      "#AEC4AA", "#BED0B0", "#D0DCB5", "#E5E7BC", "#FBF2C4", "#F3E3B2",
      "#EDD59F", "#E7C68C", "#E3B77A", "#DEA868", "#DA9857", "#D58847",
      "#D1773A", "#CC6530", "#C7522B")
  }
```

Gráfica:

```{r predmap2}
p <- SOmap_auto(x = pr, bathy = pr)
p$bathy_palette <- my_cmap
p
```

```{r predmap}
plot(basemap)
SOplot(pr, col = my_cmap)
```

## Detalles extra

### Nombres de lugares

```{r placenames1}
## obtengamos datos completos del diccionario geográfico de SCAR
xn <- an_read(cache = "session", sp = TRUE)

## reduzcamos a un nombre por lugar
xn <- an_preferred(xn, origin = "United Kingdom")

## solicitemos sugerencias en nuestra región de estudio para añadir a nuestro
# mapa
xns <- an_suggest(xn, map_scale = 20e6, map_extent = extent(pr))

## transformemos a nuestra proyección, convirtimos datos a dataframe (marco de 
# datos), mostramos las primeros 10 filas
xns <- as_tibble(SOproj(xns, target = basemap$projection)) |>
  head(10)
```

Añadimos nombres al mapa:

```{r predmap3}
plot(basemap)
SOplot(pr, col = my_cmap)

## indicando la ubicación de los nombres
with(xns, points(x = coords.x1, y = coords.x2, pch = 16, cex = 0.4))
## añadir etiquetas
with(xns, text(x = coords.x1, y = coords.x2, labels = place_name, 
               cex = 0.75, pos = 2 + 2*(seq_len(nrow(xns)) %% 2)))
```