05_Overview.Rmd

# Overview

An example analysis that very briefly demonstrates some of the tasks that we'll tackle during the workshop.

## Preparation

Make sure we have the packages we need from CRAN:

```{r pkgs1, message = FALSE}
pkgs <- c("dplyr", "ncdf4", "raster", "remotes", "robis", "worrms")
pkgs <- setdiff(pkgs, installed.packages()[, 1])
if (length(pkgs) > 0) install.packages(pkgs)

library(dplyr)
```

And from GitHub:

```{r pkgs2}
## packages with required minimum version
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)
    }
}
```

## Taxonomy

```{r tax1, message = FALSE}
library(worrms)

my_species <- "Euphausia crystallorophias"
tax <- wm_records_names(name = my_species)
tax

## the Aphia ID (taxonomic ID) of our species
my_aphia_id <- tax[[1]]$valid_AphiaID
```

## Occurrences

Get all data from the [Antarctic Euphausiacea occurence data from "German Antarctic Marine Living Resources" (GAMLR) Expeditions](https://ipt.biodiversity.aq/resource?r=gamlr) data set:

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

Plot the complete distribution of samples in black, and *Euphausia crystallorophias* in green:

```{r map1, message = FALSE}
library(SOmap)
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"))
```

Or as a polar stereo map:

```{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"))
```

We're going to fit a species (presence/absence) distribution model, so first let's reorganise our data into presence/absence by sampling site:

```{r dat1}
xfit <- x %>% dplyr::rename(lon = "decimalLongitude", lat = "decimalLatitude") %>%
    group_by(lon, lat) %>%
    dplyr::summarize(present = any(my_aphia_id %in% aphiaID))
```

## Environmental data

```{r envdat1, message = FALSE}
library(blueant)
## put the data into a temporary directory
my_data_directory <- tempdir()

## the data source we want
data_source <- sources_sdm("Southern Ocean marine environmental data")

## fetch the data
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)

## create a raster stack of all layers
env_stack <- stack(nc_files)

## the first few layers
head(names(env_stack))
```

Select just the `depth` and `ice_cover_mean` layers and extract their values at our sampling locations:

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

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

xfit <- bind_cols(xfit, temp)
head(xfit)

```

## Fit model

We have presence/absence data, so we'll fit a simple binomial model. The probability of presence of *Euphausia crystallorophias* is fitted as a smooth function of depth and mean sea ice cover:

```{r mod1, message = FALSE}
library(mgcv)

fit <- gam(present ~ s(depth) + s(ice_cover_mean), family = binomial, data = xfit)

summary(fit)
```

The fits to depth and ice cover:

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

This suggests that we are more likely to find *Euphausia crystallorophias* in shallow areas with high annual ice cover, which seems plausible given that this species is typically found in coastal/shelf waters.

## Predict from model

```{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(raster::extract(env_stack, xpred[, c("lon", "lat")])))

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

## create raster
pr <- rasterFromXYZ(xpred[, c("lon", "lat", "predicted")])
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")
```

Plot it:

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

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

## Other bits and pieces

### Place names

```{r placenames1}
library(antanym)

## get full SCAR gazetteer data
xn <- an_read(cache = "session", sp = TRUE)

## reduce to one name per feature
xn <- an_preferred(xn, origin = "United Kingdom")

## ask for suggestions in our region to show on our map
xns <- an_suggest(xn, map_scale = 20e6, map_extent = extent(pr))

## transform to our map projection, convert to data frame, take the top 10
xns <- as_tibble(SOproj(xns, target = basemap$projection)) %>% head(10)
```

Add them to the map:

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

## placename points
with(xns, points(x = longitude, y = latitude, pch = 16, cex = 0.4))
## and labels
with(xns, text(x = longitude, y = latitude, labels = place_name, 
               cex = 0.75, pos = 2 + 2*(seq_len(nrow(xns)) %% 2)))
```