Pull request

.gitignore

@@ -1 +1,6 @@
 .DS_Store
+*~
+.Rproj.user
+.Rhistory
+data-library
+_bookdown_files/

01_Schedule.Rmd

@@ -0,0 +1,58 @@
+# Course schedule
+
+## Monday
+
+Time  Item
+----  ----
+9:00  [Course introduction][Course introduction], housekeeping, presentation of participants
+10:00 Break
+10:15 [Introduction to R and the Tidyverse]
+11:15 Break
+11:30 [Open data and data cleaning]
+12:30 Lunch
+13:00 Finding and downloading data through the rOpenSci packages
+14:00 Break
+14:15 [Mapping]
+15:15 Break
+15:30 [Data visualisation]
+17:30 End day 1
+
+
+
+## Tuesday
+Time  Item
+----  ----
+9:00  [What are species distribution models?][SDM introduction and review]
+10:00 Break
+10:15 [Recent and past application of SDMs]
+11:15 Break
+11:30 Preparation for exercises
+12:30 Lunch
+13:30 [Run a simple SDM, understand and generate SDM outputs][SDM outputs]
+15:00 Break
+15:15 [SDM calibration: what you should know before running a model][Calibration]
+17:30 End day 2
+
+## Wednesday
+Time  Item
+----  ----
+9:30  [Introduction to different SDM algorithms][SDM algorithms]
+10:30 Break
+10:45 Questions, complete previous courses if necessary
+12:30 Lunch
+13:30 Personal work
+17:30 End day 3
+
+
+## Thursday
+Time  Item
+----  ----
+9:30  Personal work (with breaks)
+17:30 End Day 4
+
+
+## Friday
+Time  Item
+----  ----
+9:30  Personal work, final presentations, and discussion (with breaks)
+17:00 End Day 5

02_Course_Introduction.Rmd

@@ -0,0 +1,3 @@
+# Course Introduction
+
+[Presentation (PDF)](./course_material/AVDP_01_INTRODUCTION.pdf)

05_Overview.Rmd

@@ -0,0 +1,196 @@
+# 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)))
+```