+
+
+```{r echo = FALSE}
+showtbl(tribble(~"decimal places", ~"decimal degrees", ~DMS, ~"Object that can be unambiguously recognized at this scale", ~"N/S or E/W at equator", ~"E/W at 67N/S",
+0, "1.0", "1° 00' 0\"", "country or large region", "111.32 km", "43.496 km",
+1, "0.1", "0° 06' 0\"", "large city or district", "11.132 km", "4.3496 km",
+2, "0.01", "0° 00' 36\"", "town or village", "1.1132 km", "434.96 m",
+3, "0.001", "0° 00' 3.6\"", "neighborhood, street", "111.32 m", "43.496 m",
+4, "0.0001", "0° 00' 0.36\"", "individual street, land parcel", "11.132 m", "4.3496 m",
+5, "0.00001", "0° 00' 0.036\"", "individual trees, door entrance", "1.1132 m", "434.96 mm",
+6, "0.000001", "0° 00' 0.0036\"", "individual humans", "111.32 mm", "43.496 mm",
+7, "0.0000001", "0° 00' 0.00036\"", "practical limit of commercial surveying", "11.132 mm", "4.3496 mm",
+8, "0.00000001", "0° 00' 0.000036\"", "specialized surveying (e.g. tectonic plate mapping)", "1.1132 mm", "434.96 µm"))
+```
+
+### When
+
+#### [eventDate](https://dwc.tdwg.org/terms/#dwc:eventDate)
+
+Additional Terms
+
+*[year](https://dwc.tdwg.org/terms/#dwc:year)*,
+*[month](https://dwc.tdwg.org/terms/#dwc:month)*,
+*[Day](https://dwc.tdwg.org/terms/#dwc:Day)*,
+*[Time](https://dwc.tdwg.org/terms/#dwc:Time)*,
+*[Timezone](https://dwc.tdwg.org/terms/#dwc:Timezone)*
+
+
+
+**[eventDate](https://dwc.tdwg.org/terms/#dwc:eventDate)** The date and time when the occurrence/event was recorded. Can be the date-time or interval during which an event occurred. This term uses the ISO 8601 standard. OBIS recommends using the extended ISO 8601 format with hyphens.
+
+
+
+For time intervals ISO 8601 uses / as a separator.
+Date and time are separated by T.
+Times can have a time zone indicator at the end, if this is not the case then the time is assumed to be local time.
+When a time is UTC, a Z is added. Some examples of ISO 8601 dates are:
+
+1973-02-28T15:25:00
+2005-08-31T12:11+12
+1993-01-26T04:39+12/1993-01-26T05:48+12
+2008-04-25T09:53
+1948-09-13
+1993-01/02
+1993-01
+1993
+
+It is an option to split out the *eventDate* into its seperate components. This can be a good safeguard against software like Excel reformatting your dates and times.
+
+*year*,
+*month*,
+*Day*,
+*Time*,
+*Timezone*
+
+*verbatimEventDate* You can put the original date format here, not really needed but might be useful, in case where you are digitising old literature records.
+
+
+### Who
+
+*[identifiedBy](https://dwc.tdwg.org/terms/#dwc:identifiedBy)*, *[dateIdentified](https://dwc.tdwg.org/terms/#dwc:dateIdentified)*
+
+
+*[identifiedBy](https://dwc.tdwg.org/terms/#dwc:identifiedBy)*
+
+A list (concatenated and separated) of names of people, groups, or organizations who assigned the Taxon to the subject. Recommended best practice is to separate the values in a list with space vertical bar space ( | )
+
+*[dateIdentified](https://dwc.tdwg.org/terms/#dwc:dateIdentified)*
+
+This one is very seldom filled out but it is actually an important one and it should be specific. It should be the person that did the actual identification not the suprvisor.
+
+In case of a reindentification it is often useful to be able to contact the person who did the initial identification.
+
+### How
+*[samplingProtocol](https://dwc.tdwg.org/terms/#dwc:samplingProtocol)*
+
+
+This is a tricky one it is quite relevent but this field often doesn't provide enough space to describe what you did. For describing a specific gear you can use [BODC vocab : L22](https://www.bodc.ac.uk/resources/vocabularies/vocabulary_search/L22/)
+
diff --git a/30_Biological_Data.Rmd b/30_Biological_Data.Rmd
new file mode 100644
index 0000000..df9d933
--- /dev/null
+++ b/30_Biological_Data.Rmd
@@ -0,0 +1,22 @@
+# Biological data
+
+## Taxonomy
+
+The Register of Antarctic Marine Species (RAMS) is the authoritative taxonomic database for Antarctic marine organisms. RAMS is part of the World Register of Marine Species (WoRMS).
+
+- [worrms](https://cran.r-project.org/package=worrms) client for the [WoRMS](http://www.marinespecies.org/) API. Contains mostly taxonomic data, but also trait data.
+- [taxize](https://cran.r-project.org/package=taxize) provides access to 20ish sources of taxonomic data sources, including WoRMS. It has consistent data outputs and function interfaces across the different data sources so that you don't need to tailor your code to different taxonomic data providers.
+
+RAMS is currently being extended to cover non-marine taxa, which will become the Register of Antarctic Species (RAS). Hopefully this will remain covered by `worrms` and the server-side infrastructure hosted by VLIZ.
+
+For more detail on R packages dealing with taxonomy in general, see the [rOpenSci taxonomy task view](https://github.com/ropensci/taxonomy).
+
+## Other packages
+
+Tracking of animals using satellite, GPS, or light-level geolocation tags is common, and there are many R packages that can help with this. See the [spatiotemporal task view](https://cloud.r-project.org/web/views/SpatioTemporal.html) for a more complete list. Also of interest may be:
+
+- [TwilightFree](https://github.com/ABindoff/TwilightFree) provides a method for processing light-level geolocation data that is robust to noise (sensor shading and obscuration) and may be particularly suitable for Southern Ocean applications.
+
+- [traipse](https://github.com/Trackage/traipse) provides tools for tracking data, for common metrics of distance, direction, and speed.
+
+- [foieGras](https://github.com/ianjonsen/foieGras) fits a continuous-time model (RW or CRW) in state-space form to filter Argos satellite location data.
diff --git a/40_Environmental_Data.Rmd b/40_Environmental_Data.Rmd
new file mode 100644
index 0000000..4e766f5
--- /dev/null
+++ b/40_Environmental_Data.Rmd
@@ -0,0 +1,189 @@
+# Environmental data
+
+## Bowerbird/blueant
+
+Very commonly, we want to know about the environmental conditions at our points of interest. For the remote and vast Southern Ocean these data typically come from satellite or model sources. Some data centres provide extraction tools that will pull out a subset of data to suit your requirements, but often it makes more sense to cache entire data collections locally first and then work with them from there.
+
+[bowerbird](https://github.com/ropensci/bowerbird) provides a framework for downloading data files to a local collection, and keeping it up to date. The companion [blueant](https://github.com/AustralianAntarcticDivision/blueant) package provides a suite of definitions for Southern Ocean and Antarctic data sources that can be used with `bowerbird`. It encompasses data such as sea ice, bathymetry and land topography, oceanography, and atmospheric reanalysis and weather predictions, from providers such as NASA, NOAA, Copernicus, NSIDC, and Ifremer.
+
+Why might you want to maintain local copies of entire data sets, instead of just fetching subsets of data from providers as needed?
+
+- many analyses make use of data from a variety of providers (in which case there may not be dynamic extraction tools for all of them),
+- analyses might need to crunch through a whole collection of data in order to calculate appropriate statistics (temperature anomalies with respect to a long-term mean, for example),
+- different parts of the same data set are used in different analyses, in which case making one copy of the whole thing may be easier to manage than having different subsets for different projects,
+- a common suite of data are routinely used by a local research community, in which case it makes more sense to keep a local copy for everyone to use, rather than multiple copies being downloaded by different individuals.
+
+In these cases, maintaining local copies of a range of data from third-party providers can be extremely beneficial, especially if that collection is hosted with a fast connection to local compute resources (virtual machines or high-performance computational facilities).
+
+
+Install from GitHub:
+
+```{r ed_install_pkg, eval = FALSE}
+remotes::install_github("AustralianAntarcticDivision/blueant")
+```
+
+And load the package before use.
+
+```{r ed1, cache = FALSE}
+library(blueant)
+```
+### Available data sets
+
+First, we can see the available data sets via the `sources` function.
+
+```{r ed2}
+srcs <- blueant::sources()
+## the names of the first few
+head(srcs$name)
+## the full details of the first one
+srcs[1, ]
+```
+
+### Usage
+
+Choose a directory into which to download the data. Usually this would be a persistent directory on your machine so that data sets downloaded in one session would remain available for use in later sessions, and not need re-downloading. A persistent directory could be something like `c:\data\` (on Windows), or you could use the `rappdirs` package (the `user_cache_dir` function) to suggest a suitable directory (cross-platform).
+
+```{r ed0hidden, echo = FALSE, results = "asis", cache = FALSE}
+if (grepl("ben_ray", tempdir())) {
+ cat("Here we'll use the `c:/data/cache` directory:\n")
+ my_data_dir <- "/data/cache"
+} else {
+ cat("Here we'll use a temporary directory:\n")
+ my_data_dir <- tempdir()
+}
+```
+
+```{r ed0a, include = grepl("ben_ray", tempdir()), eval = FALSE, cache = FALSE}
+my_data_dir <- "/data/cache"
+```
+```{r ed0b, include = !grepl("ben_ray", tempdir()), eval = FALSE, cache = FALSE}
+my_data_dir <- tempdir()
+```
+
+Select the data source that we want:
+
+```{r ed4}
+data_source <- sources("Southern Ocean marine environmental data")
+
```
+
+Note that it's a good idea to check the dataset size before downloading it, as some are quite large! (Though if you are running the download interactively, it will ask you before downloading a large data set).
+
+```{r ed5}
+data_source$collection_size ## size in GB
+```
+
+And fetch the data:
+```{r ed6}
+result <- bb_get(data_source, local_file_root = my_data_dir, verbose = TRUE)
+```
+
+Now we have a local copy of our data. The sync can be run daily so that the local collection is always up to date - it will only download new files, or files that have changed since the last download. For more information on `bowerbird`, see the [package vignette](https://ropensci.github.io/bowerbird/articles/bowerbird.html).
+
+The `result` object holds information about the data that we downloaded:
+
+```{r ed7}
+result
+```
+
+The `result$files` element tells us about the files:
+```{r ed8}
+head(result$files[[1]])
+```
+
+These particular files are netCDF, and so could be read using e.g. the `raster` or `ncdf4` packages. However, different data from different providers will be different in terms of grids, resolutions, projections, variable-naming conventions, and other facets, which tends to complicate these operations. In the next section we'll look at the `raadtools` package, which provides a set of tools for doing common operations on these types of data.
+
+
+## RAADtools
+
+The [`raadtools`](https://github.com/AustralianAntarcticDivision/raadtools) package provides a consistent interface to a range of environmental and similar data, and tools for working with them. It is designed to work data with collections maintained by the `bowerbird`/`blueant` packages, and builds on R's existing ecosystem of packages for working with spatial, raster, and multidimensional data.
+
+Here we'll use two different environmental data sets: sea ice and water depth. Water depth does not change with time but sea ice is provided at daily time resolution.
+
+First download daily sea ice data (from 2013 only), and the ETOPO2 bathymetric data set. ETOPO2 is somewhat dated and low resolution compared to more recent data, but will do as a small dataset for demo purposes. This may take a few minutes, depending on your connection speed:
+
+```{r raadt1, eval = FALSE}
+src <- bind_rows(
+ sources("NSIDC SMMR-SSM/I Nasateam sea ice concentration", hemisphere = "south", time_resolutions = "day",
+ years = 2013),
+ sources("ETOPO2 bathymetry"))
+result <- bb_get(src, local_file_root = my_data_dir, clobber = 0, verbose = TRUE, confirm = NULL)
+```
+
+```{r raadt2, echo = FALSE, message = FALSE}
+## code not shown in output: capture the output and trim it down a bit
+src <- bind_rows(
+ sources("NSIDC SMMR-SSM/I Nasateam sea ice concentration", hemisphere = "south", time_resolutions = "day", years = 2013),
+ sources("ETOPO2 bathymetry"))
+op <- capture.output(result <- bb_get(src, local_file_root = my_data_dir, clobber = 0, verbose = TRUE, confirm = NULL))
+op[4] <- ""
+op[5] <- " [... output truncated]"
+for (oo in op[1:5]) cat(oo, "\n")
+```
+
+
+Now load the `raadtools` package and tell it where our data collection has been stored:
+
+```{r raadt3, results = "hold", cache = FALSE}
+library(raadtools)
+set_data_roots(my_data_dir)
+```
+
+Let's say that we have some points of interest in the Southern Ocean — perhaps a ship track, or some stations where we took marine samples, or as we'll use here, the [track of an elephant seal](http://www.meop.net/) as it moves from the Kerguelen Islands to Antarctica and back again (Data from IMOS 2018[^1], provided as part of the `SOmap` package).
+
+```{r get_track_data, message = FALSE, warning = FALSE, cache = FALSE}
+data("SOmap_data", package = "SOmap")
+ele <- SOmap_data$mirounga_leonina %>% dplyr::filter(id == "ct96-05-13")
+```
+
+Define our spatial region of interest and extract the bathymetry data from this region, using the ETOPO2 files we just downloaded:
+
+```{r raadt4}
+roi <- round(c(range(ele$lon), range(ele$lat)) + c(-2, 2, -2, 2))
+bx <- readtopo("etopo2", xylim = roi)
+```
+
+And now we can make a simple plot of our our track superimposed on the bathymetry:
+
+```{r raadt5}
+plot(bx)
+lines(ele$lon, ele$lat)
+```
+
+The real power of `raadtools` comes from its extraction functions. We can extract the depth values along our track using the `raadtools::extract()` function. We pass it the data-reader function to use (`readtopo`), the data to apply it to (`ele[, c("lon", "lat")]`), and any other options to pass to the reader function (in this case, specifying the topographic data source `topo = "etopo2"`):
+
+```{r raadt6}
+ele$depth <- raadtools::extract(readtopo, ele[, c("lon", "lat")], topo = "etopo2")
+```
+
+Plot the histogram of depth values, showing that most of the track points are located in relatively shallow waters:
+
+```{r raadt7}
+with(ele, hist(depth, breaks = 20))
+```
+
+This type of extraction will also work with time-varying data — for example, we can extract the sea-ice conditions along our track, based on each track point's location and time:
+
+```{r raadt8, results = "none", message = FALSE, warning = FALSE}
+ele$ice <- raadtools::extract(readice, ele[, c("lon", "lat", "date")])
+```
+```{r raadt9}
+## points outside the ice grid will have missing ice values, so fill them with zeros
+ele$ice[is.na(ele$ice)] <- 0
+with(ele, plot(date, ice, type = "l"))
+```
+
+## Other useful packages
+
+- the [PolarWatch](https://polarwatch.noaa.gov/) project aims to enable data discovery and broader use of high-latitude ocean remote sensing data sets. The dedicated ERDDAP server (https://polarwatch.noaa.gov/erddap) is accessible to R users with [rerddap](https://cran.r-project.org/package=rerddap).
+
+- [rsoi](https://cran.r-project.org/package=rsoi) downloads the most up to date Southern Oscillation Index, Oceanic Nino Index, and North Pacific Gyre Oscillation data.
+
+- satellite reflectance data are a common basis for estimating chlorophyll-*a* and other phytoplankton parameters at ocean-basin scales. Global products are widely available; however, Southern-Ocean specific algorithms are likely to provide better estimates in these regions. [croc](https://github.com/sosoc/croc) implements the Johnson et al. (2013) Southern Ocean algorithm.
+
+- more broadly, [oce](https://cran.r-project.org/package=oce) provides a wide range of tools for reading, processing, and displaying oceanographic data, including measurements from Argo floats and CTD casts, sectional data, sea-level time series, and coastline and topographic data.
+
+- [fda.oce](https://github.com/EPauthenet/fda.oce) provides functional data analysis of oceanographic profiles for front detection, water mass identification, unsupervised or supervised classification, model comparison, data calibration, and more.
+
+- [distancetocoast](https://github.com/mdsumner/distancetocoast) provides "distance to coastline" data for longitude and latitude coordinates.
+
+- [geodist](https://github.com/hypertidy/geodist) for very fast calculation of geodesic distances.
diff --git a/50_Mapping.Rmd b/50_Mapping.Rmd
new file mode 100644
index 0000000..eefe464
--- /dev/null
+++ b/50_Mapping.Rmd
@@ -0,0 +1,605 @@
+# Mapping
+
+## Maps in R
+
+The oldest and most *core* general mapping package in R is the `maps` package. It has a simple whole-world coastline data set for immediate use.
+
+```{r maps-package}
+maps::map()
+```
+
+The data underlying this live map is available by capturing the output as an actual object. Notice that the coastline for Antarctica does not extend
+to the south pole, and that parts of Russia that are east of 180 longitude are not in the western part of the map.
+
+```{r maps-assign}
+m <- maps::map(plot = FALSE)
+
+lonlat <- cbind(m$x, m$y)
+
+plot(lonlat, pch = "+", cex = 0.4, axes = FALSE)
+lines(lonlat, col = "dodgerblue")
+
+abline(h = c(-90, 90), v = c(-180, 180))
+
+```
+
+A very similar and more modern data set is available in the `maptools` package.
+
+```{r wrld-simpl}
+data("wrld_simpl", package = "maptools")
+library(sp)
+plot(wrld_simpl)
+```
+
+This data set aligns exactly to the conventional -180/180 -90/90 extent of the longitude/latitude projection.
+
+```{r lon180-lat90}
+plot(0, type = "n", axes = FALSE, xlab = "", ylab = "", xlim = c(-180, 180), ylim = c(-90, 90))
+rect(xleft = -180, ybottom = -90, xright = 180, ytop = 90, border = "darkred", lwd = 4, lty = 2)
+plot(wrld_simpl, add = TRUE)
+```
+
+### Exercises
+
+1. How can we find the longitude and latitude ranges of the maps data `m` and the maptools data `wrld_simpl`?
+2. Can we draw polygons with a fill colour with the maps package?
+
+
+Answer 1: `range(m$x, na.rm = TRUE)` `range(m$y, na.rm = TRUE`) also `m$range`
+
+Answer 2: `polygon(lonlat, col = "grey")` does not work, and `map(mp, fill = TRUE, col = "grey")` does not work, but `maps::map(fill = TRUE, col = "grey")` does seem to work.
+
+
+What's going on? Look at the very south-eastern corner of the map. The "coastline" has been extended to the very south boundary of the available area.
+
+
+
+```{r south-east}
+plot(0, type = "n", axes = FALSE, xlab = "", ylab = "", xlim = c(-150, 180), ylim = c(-90, -60))
+plot(wrld_simpl, add = TRUE, col = "grey")
+rect(xleft = -180, ybottom = -90, xright = 180, ytop = 90, border = "darkred", lwd = 4, lty = 2)
+
+maps::map(add = TRUE, col = "dodgerblue", lwd = 3)
+```
+
+When we add the old maps coastline see that it does not extend to 90S and it does not traverse the southern boundary.
+
+One reason for this is that if we choose a projection where the east and west edges of the Antarctic coastline meet then we get what looks a fairly clean join.
+
+```{r laea}
+## scale factor
+f <- 3e6
+plot(rgdal::project(lonlat, "+proj=laea +lat_0=-90 +datum=WGS84"), asp = 1, type = "l",
+ xlim = c(-1, 1) * f, ylim = c(-1, 1) * f, xlab = "", ylab = "")
+
+```
+
+If we try the same with `wrld_simpl` it's not as neat. We have a strange "seam" that points exactly to the south pole (our projection is centred on longitude = 0, and latitude = -90.
+
+
+
+```{r wrld}
+plot(sp::spTransform(wrld_simpl, "+proj=laea +lat_0=-90 +datum=WGS84"), asp = 1,
+ xlim = c(-1, 1) * f, ylim = c(-1, 1) * f, xlab = "", ylab = "", lwd = 3)
+abline(v = 0, h = 0, lty = 2, col = "grey")
+```
+
+
+### Let's use the maps data!
+
+In `m` we have the maps data structure, and this looks promising.
+
+```{r maps-data}
+str(m)
+
+mp <- m
+pxy <- rgdal::project(lonlat, "+proj=laea +lat_0=-90 +datum=WGS84")
+mp$x <- pxy[,1]
+mp$y <- pxy[,2]
+mp$range <- c(range(mp$x,na.rm = TRUE), range(mp$y, na.rm = TRUE))
+mp$range
+plot(c(-1, 1) * f, c(-1, 1) * f, type = "n", asp = 1)
+maps::map(mp, add = TRUE)
+
+
+## but it doesn't take much to go awry
+plot(c(-1, 1) * f, c(-1, 1) * f, type = "n", asp = 1)
+maps::map(mp, add = TRUE, fill = TRUE, col = "grey")
+
+```
+
+The problem is that the maps database has enough internal structure to join lines correctly, with `NA` gaps between different connected linestrings, but not enough to draw these things as polygons. A similar problem occurs in the default projection. While `wrld_simpl` has been extend by placing two dummy coordinates at the east and west versions of the south pole, this data set does not have those.
+
+We have to look quite carefully to understand what is happening, but this is wrapping around overlapping itself and so close to the southern bound we barely notice.
+
+
+```{r no-pole}
+plot(0, type = "n", axes = FALSE, xlab = "", ylab = "", xlim = c(-180, -110), ylim = c(-90, -60))
+
+rect(xleft = -180, ybottom = -90, xright = 180, ytop = 90, border = "darkred", lwd = 4, lty = 2)
+
+maps::map(add = TRUE,col = "grey", fill = TRUE)
+
+maps::map(col = "grey", fill = TRUE)
+
+
+mpmerc <- m
+pxy <- rgdal::project(lonlat, "+proj=merc +datum=WGS84")
+mpmerc$x <- pxy[,1]
+mpmerc$y <- pxy[,2]
+mpmerc$range <- c(range(mpmerc$x,na.rm = TRUE), range(mpmerc$y, na.rm = TRUE))
+mpmerc$range
+
+## the catastrophe made a little clearer
+plot(0, xlim = range(mpmerc$range[1:2]), ylim = c(mpmerc$range[1], 0))
+maps::map(mpmerc, fill = TRUE, col = "grey", add = TRUE)
+
+
+```
+
+## SOmap
+
+The `SOmap` package is intended to solve some of these problems, and provide an easier way to produce nice-looking maps of Antarctica and the Southern Ocean. It is primarily focused on maps in polar stereographic projection (although the `SOmap_auto` function extends this to other projections). `SOmap` won't necessarily get you exactly the map you want in every circumstance, but the idea is that in most cases it should get you close enough, and if need be you can make modifications to suit your exact purposes.
+
+Please bear in mind that `SOmap` is still in development, and so its functionality (function parameters and/or behaviour) may change.
+
+By default, `SOmap` works with base graphics (and associated functionality from packages such as `raster` and `sp`). It is also possible to work with `ggplot2`-based graphics, as described below.
+
+Start by installing the `SOmap` package if you haven't done so already:
+
+```{somap_install eval = FALSE}
+remotes::install_github("AustralianAntarcticDivision/SOmap")
+
```
+
+Then load the package:
+
+```{r somap_init, cache = FALSE}
+library(SOmap)
+## also define a colour map to use for some examples
+my_cmap <- colorRampPalette(c("#4D4140", "#596F7E", "#168B98",
+ "#ED5B67", "#E27766", "#DAAD50", "#EAC3A6"))(51)
+```
+
+```{r somap_pkgs, echo = FALSE, message = FALSE, warning = FALSE, cache = FALSE}
+## ensure these are loaded, library calls below are cached and might not get evaluated every run
+library(sp)
+library(raster)
+library(ggplot2)
+```
+
+### Circumpolar maps
+
+A basic circumpolar map in polar stereographic projection:
+
+```{r somap1
}
+SOmap()
+```
+
+`SOmanagement()` provides a number of contextual layers such as MPA boundaries and management zones.
+
+```{r somap2}
+SOmap(trim = -40) ## plot to 40S
+## add the exclusive economic zones management layer
+SOmanagement(eez = TRUE)
+```
+
+#### Adding points
+
+```{r somap_pts1}
+## some longitude/latitude data
+library(sp)
+my_points_ll <- data.frame(lon = seq(0, 350, by = 10), lat = -55, z = runif(36))
+coordinates(my_points_ll) <- c("lon", "lat")
+projection(my_points_ll) <- "+proj=longlat +datum=WGS84"
+```
+
+Our data need to be reprojected to match our map before plotting. The `SOproj` function does this:
+
+```{r somap_pts2}
+## reproject to our SOmap projection
+my_points <- SOproj(my_points_ll)
+## and plot
+SOmap()
+plot(my_points, col = "blue", add = TRUE)
+```
+
+Or use `SOplot` to reproject and plot in one step:
+
+```{r somap_pts3}
+SOmap()
+SOplot(my_points_ll, col = "blue")
+```
+
+#### Adding raster layers
+
+First let's construct some artificial raster data (in longitude-latitude space) for demonstration purposes:
+
+```{r somap_raster1, message = FALSE, warning = FALSE}
+library(raster)
+temp <- as.data.frame(expand.grid(lon = seq(100, 140, by = 0.25),
+ lat = seq(-65, -45, by = 0.1)))
+temp$val <- sqrt((temp$lon - 120)^2/3 + (temp$lat - -40)^2/5)
+## create raster object
+xr <- rasterFromXYZ(temp)
+projection(xr) <- "+proj=longlat +datum=WGS84"
+```
+
+`SOplot` will reproject and plot this for us:
+
+```{r somap_raster2}
+SOmap()
+SOplot(xr)
+```
+
+The legend is out of character with the rest of the map. We can use `SOleg` to fix that:
+
+```{r somap_raster3}
+## draw the base map
+SOmap()
+## add our raster
+SOplot(xr, legend = FALSE, col = my_cmap)
+## add the legend
+SOleg(xr, position = "topright", col = my_cmap, ticks = 6,
+ type = "continuous", label = "My variable")
+```
+
+OK, well that worked but clearly the labels need tidying up. The easiest way is probably to set the number of decimal places in the label values via the `rnd` parameter:
+
+
+```{r somap_raster4rnd}
+SOmap()
+SOplot(xr, legend = FALSE, col = my_cmap)
+SOleg(xr, position = "topright", col = my_cmap, ticks = 6, rnd = 2,
+ type = "continuous", label = "My variable")
+```
+
+Alternatively, we could explicitly set the colour range and labels.
+
+```{r somap_raster4}
+## draw the base map
+SOmap()
+## add our raster, controlling the colour range to span the values 0 to 30
+colour_breaks <- seq(0, 30, length.out = length(my_cmap) + 1)
+SOplot(xr, legend = FALSE, col = my_cmap, breaks = colour_breaks)
+## add the legend, again controlling the colour range
+label_breaks <- seq(0, 30, length.out = 7)
+SOleg(position = "topright", col = my_cmap, breaks = label_breaks,
+ type = "continuous", label = "My variable")
+```
+
+Note that if we *don't* want to show the bathymetric legend, we may run into problems:
+
+```{r somap_legend1}
+SOmap(bathy_legend = FALSE) ## suppress the bathy legend
+SOleg(position = "topright", col = my_cmap, breaks = label_breaks,
+ type = "continuous", label = "My variable")
+```
+
+The legend has been chopped off because the layout has not left enough space around the map for the curved legend. Currently, the best solution is probably to generate the `SOmap` object *with* the bathymetric legend, but then remove it before plotting (see the [Modifying map objects][Modifying map objects (advanced usage)] section for more details on this):
+
+```{r somap_legend2}
+temp <- SOmap()
+temp$bathy_legend <- NULL ## remove the bathy legend
+plot(temp)
+SOleg(position = "topright", col = my_cmap, breaks = label_breaks,
+ type = "continuous", label = "My variable")
+```
+
+
+Multiple rasters:
+
+```{r somap_raster5}
+xr2 <- raster::shift(xr, -70) ## offset in longitude
+SOmap()
+SOplot(xr, legend = FALSE, col = my_cmap)
+SOplot(xr2, legend = FALSE, col = my_cmap)
+```
+
+
+### Non-circumpolar maps
+
+The `SOmap_auto` function will take your input data and make a guess at an appropriate projection and extent to use. Note that this is not always going to guess the *best* projection and extent, so you should view it as a starting point from which you can generate a map to your exact requirements.
+
+Use the elephant seal track data bundled with the package:
+
+```{r soauto1}
+ellie <- SOmap_data$mirounga_leonina
+## construct and plot the map
+SOmap_auto(ellie$lon, ellie$lat)
+```
+
+Just a blank map to which you could add other things:
+```{r soauto2}
+SOmap_auto(ellie$lon, ellie$lat, input_points = FALSE, input_lines = FALSE)
+```
+
+You can pass a raster as input data, but note that it won't plot the raster (it uses its extent to infer an appropriate extent for the map):
+```{r soauto3}
+SOmap_auto(xr)
+```
+
+But we can add the raster if we wish:
+```{r soauto4}
+SOmap_auto(xr)
+SOplot(xr, col = my_cmap)
+```
+
+We can force a particular projection:
+
+```{r soauto5}
+SOmap_auto(xr, target = "laea", centre_lon = 147, centre_lat = -42)
+SOplot(xr, col = my_cmap)
+```
+
+Same but by supplying a full proj4 string to `target`:
+
+```{r soauto6}
+SOmap_auto(xr, target = "+proj=laea +lat_0=-42 +lon_0=147")
+SOplot(xr, col = my_cmap)
+```
+
+See [the SOmap_auto vignette](https://australianantarcticdivision.github.io/SOmap/articles/many-automap-examples.html) for more examples.
+
+
+### Plotting via ggplot2
+
+The `SOmap` and `SOmap_auto` functions do their plotting using base graphics. If you are more comfortable working with `ggplot2`, this is also possible. The `SOgg` function takes an object created by one of those functions (using base graphics) and converts it to use `ggplot2` graphics instead. As with other `SOmap` functions, this returns an object (of class `SOmap_gg` or `SOmap_auto_gg`) that contains all of the information needed to generate the map. Printing or plotting this object will cause it to construct a `ggplot` object. Printing or plotting *that* object will cause it to be drawn to the graphics device, just like any other `ggplot` object.
+
+```{r somap3}
+myplot <- SOmap()
+myplotgg <- SOgg(myplot) ## creates a SOmap_gg object
+class(myplotgg)
+my_ggplot <- plot(myplotgg) ## creates a ggplot object
+class(my_ggplot)
+plot(my_ggplot) ## plot it
+```
+
+Or in one step (this will cause myplot to be converted to SOmap's internal gg format, then a ggplot object constructed from that, then that object will be plotted):
+
+```{r somap4}
+SOgg(myplot)
+```
+
+### Modifying map objects (advanced usage)
+
+The goal of `SOmap` is to make it fairly easy to produce a fairly-good-looking map that will be adequate for most mapping requirements. It will never be possible to automatically produce a *perfect* map in *every circumstance*, but the aim is to have a low-effort way of getting fairly close most of the time.
+
+This section describes some approaches to modifying a map to get it closer to your particular needs. Be warned: depending on the exact modifications needed, this might get you pretty close to the crumbling edge of `SOmap` development. In particular, anything that requires modifying the internal structure of an `SOmap` object may change in the future (with luck, we'll make this sort of thing easier - but we're not there yet.)
+
+#### Modifying base graphics maps
+
+Calls to `SOmap()`, `SOmanagement()`, `SOmap_auto()` return an object of class `SOmap`, `SOmap_management`, or `SOmap_auto`. These objects contain all of the data and plotting instructions required to draw the map. Calling `print()` or `plot()` on one of these objects will cause that code to be executed, and the object to be drawn in the current graphics device. Hence, calling `SOmap()` directly *without* assigning the result to a variable will make it appear in the graphics device, because the returned object is being printed to the console (and thus triggering the `print` method). But you can also assign the result to a variable, e.g. `myplot <- SOmap()` and then explicitly plot the object with `plot(myplot)`. The advantage of this is that you can potentially manipulate the `myplot` object to make changes to the map before plotting it.
+
+Note, this is likely to be fragile. Proceed at your own risk!
+
+```{r somap2a}
+mymap <- SOmap()
+names(mymap)
+```
+
+The object contains a `plot_sequence` component, which defines the order in which each part of the plot is drawn. The other components of the object contain the code required to draw each part. Take e.g. the ice component (this is the ice shelves, glacier tongues, etc). This is a list (in this case with only one element). Each element of the list specifies a function to run along with arguments to pass to it:
+
+```{r somap2b}
+str(mymap$ice)
+```
+
+We can modify the function and/or its arguments:
+
+```{r somap2c}
+mymap$ice[[1]]$plotargs$col <- "green"
+```
+
+```{r somap2d}
+plot(mymap)
+```
+
+We can remove entire components:
+```{r somap2e}
+temp <- mymap
+temp$coastline <- NULL
+temp$ice <- NULL
+plot(temp)
+```
+
+But note that some elements are required. In particular, the bathymetry layer can't currently be removed because the code that draws this is also the code that creates the plot page via `plot.new()`. The code below would fail outright if there was no existing existing plot. If there was an existing plot in the graphics device, this code would run but give unpredictable results because it would draw on top of the previously-setup plot:
+
+```{r somap2f, eval = FALSE}
+## code not run here
+temp <- mymap
+temp$bathy <- NULL
+plot(temp)
+```
+
+One way around this would be to simply replace all of the bathymetric data values with `NA`s. The plotting code would still have the extent of the bathymetric layer that it needs in order to set up the plot, but no data would be shown:
+
+```{r somap2h2}
+temp <- mymap
+## the bathy data is held in temp$bathy[[1]]$plotargs$x
+## and it's a raster, so we can set its values to NA with
+raster::values(temp$bathy[[1]]$plotargs$x) <- NA_real_
+temp$bathy_legend <- NULL
+plot(temp)
+
+```
+
+We could also replace the bathymetry data with another raster object. Note that we do need to be careful about the extent and projection of this raster. For example, replacing the bathymetry raster with the ice raster (which has the same polar stereographic projection but smaller extent) gives:
+
+```{r somap2g}
+temp <- mymap
+temp$bathy[[1]]$plotargs$x <- ice
+temp$bathy_legend <- NULL
+plot(temp)
+```
+
+It's chopped off because the extent of the ice raster is being used to set the plot extent. But if we extend the ice raster to match the map extent:
+
+```{r somap2h}
+temp <- mymap
+temp$bathy[[1]]$plotargs$x <- raster::extend(ice, mymap$target)
+temp$bathy_legend <- NULL
+plot(temp)
+```
+
+#### Modifying ggplot maps
+
+We can modify `ggplot2`-based maps at two levels.
+
+##### Modifying the `ggplot` object.
+
+Remember that printing or plotting a `SOmap_gg` object produces a `ggplot` object. This can be modified by adding e.g. layers or themes just like a normal `ggplot`. Remember to load the `ggplot2` library now that we are using `ggplot2` functions directly.
+
+```{r somap6, message = FALSE, warning = FALSE}
+library(ggplot2)
+my_ggplot + geom_point(data = as.data.frame(my_points), aes(lon, lat, colour = z), size = 3) +
+ scale_colour_distiller(palette = "Spectral")
+```
+
+
Multiple rasters or multiple sets of points gets tricky if they are on different scales, because `ggplot2` is only designed to work with a single colour scale per geometry type. You can try your luck with the `ggnewscale` or `relayer` packages, although both are in a fairly experimental stage of development.
+
+```{r somap5f, warning = FALSE, message = FALSE}
+## remotes::install_github("clauswilke/relayer")
+library(relayer)
+plot(SOgg(SOmap(straight = TRUE))) +
+ rename_geom_aes(geom_raster(data = as.data.frame(SOproj(xr), xy = TRUE),
+ aes(x = x, y = y, fill2 = val)), new_aes = c(fill = "fill2")) +
+ scale_fill_gradientn(aesthetics = "fill2", colors = my_cmap, na.value = NA,
+ name = "My variable", guide = "legend")
+```
+
+```{r somap5g, warning = FALSE, message = FALSE}
+## remotes::install_github("eliocamp/ggnewscale")
+library(ggnewscale)
+plot(SOgg(SOmap(straight = TRUE))) +
+ new_scale_fill() +
+ geom_raster(data = as.data.frame(SOproj(xr), xy = TRUE),
+ aes(x = x, y = y, fill = val)) +
+ scale_fill_gradientn(colors = my_cmap, na.value = NA, name = "My variable")
+```
+
+##### Modifying the `SOmap_gg` object
+
+`SOmap_gg` objects are similar in structure to `SOmap` objects, in that they contain all of the data and plotting instructions required to draw the map:
+
+```{r somap5}
+names(myplotgg)
+```
+
+However, instead of base plotting functions, `SOmap_gg` objects use `ggplot2` function calls, e.g.:
+
+```{r somap5a}
+myplotgg$ice[[1]]$plotfun
+```
+
+We can modify these functions and/or arguments in a similar manner to `SOmap` objects.
+
+```{r somap5b}
+myplotgg$ice[[1]]$plotargs$fill <- "green"
+```
+
+```{r somap5c}
+plot(myplotgg)
+```
+
+Or remove the bathymetric raster layer:
+
+```{r somap5d}
+temp <- myplotgg
+temp$bathy <- NULL
+temp$bathy_legend <- NULL
+plot(temp)
+```
+
+Or replace it with a different raster (use the `ice` raster as an example):
+
+```{r somap5e, warning = FALSE, message = FALSE}
+temp <- myplotgg
+## convert ice raster to suitable data.frame
+ice_raster_as_df <- raster::as.data.frame(SOproj(ice), xy = TRUE)
+names(ice_raster_as_df)[3] <- "ice"
+## add this to our object in place of bathy
+temp$bathy <- SO_plotter(plotfun = "ggplot2::geom_raster",
+ plotargs = list(data = ice_raster_as_df,
+ mapping = aes_string(fill = "ice"))
)
+## change the colour scale
+temp$scale_fill[[1]]$plotargs <- list(colours = my_cmap, na.value = "#FFFFFF00", guide = FALSE)
+## remove the bathy legend
+temp$bathy_legend <- NULL
+plot(temp)
+```
+
+### Other SOmap gotchas
+
+Some other things worth noting.
+
+#### Automatic printing and for-loops
+
+If you type a variable/object name directly into the console then it triggers that object's `print` method automatically. Typing `SOmap()` at the console returns an object of class `SOmap`, and because it's happening at the console that object's `print` method is called, which causes the map to be plotted in the current graphics device.
+
+However, R turns off automatic printing inside `for` loops and functions. So this code:
+
+```{r eval = FALSE}
+for (i in 1:5) {
+ SOmap_auto()
+}
+```
+
+won't produce anything, because the `print` method never gets called. If you are generating maps using loops, you will need to explicitly call the `print` method:
+
+```{r eval = FALSE}
+for (i in 1:5) {
+ print(SOmap_auto())
+}
+```
+
+## Supporting data for maps
+
+When constructing maps, we commonly want to show features like oceanographic fronts, ice extent, coastline, place names, and MPA boundaries. There are a few sources of such data:
+
+- some layers are bundled into `SOmap`, see the `SOmap::SOmap_data` object
+- `antanym` provides access to the SCAR Composite Gazetteer of place names
+- the `quantarcticR` package provides access to [Quantarctica](http://quantarctica.npolar.no/) data layers.
+
+### quantarcticR
+
+Note, this package is still in development, so the usage as shown here might change in later versions. Install if needed:
+
+```{r eval = FALSE}
+remotes::install_github("SCAR-sandpit/quantarcticR")
+```
+
+Example usage:
+
+```{r qR1}
+library(quantarcticR)
+```
+
+```{r qR1hidden, echo = FALSE, results = "asis", cache = FALSE}
+if (grepl("ben_ray", tempdir())) {
+ quantarcticR::qa_cache_dir("c:/data/Quantarctica3")
+}
+```
+
+```{r qR1b}
+ds <- qa_datasets() ## all available layers
+head(ds)
+
+## more info about a particular layer
+my_layer <- qa_dataset("Median sea ice extent 1981-2010")
+my_layer
+
+## fetch the actual data for that layer
+layer_data <- qa_get(my_layer)
+
+## plot it
+plot(layer_data[layer_data$MONTH == "October", ])
+
+## or add to a SOmap
+SOmap(trim = -50, border_width = 0.5)
+SOplot(layer_data[layer_data$MONTH == "October", ], col = "red")
+
+
+```
+
+### antanym
+
+See the [overview](overview.html) article.
diff --git a/55_Data_Visualisation.Rmd b/55_Data_Visualisation.Rmd
new file mode 100644
index 0000000..f3aafab
--- /dev/null
+++ b/55_Data_Visualisation.Rmd
@@ -0,0 +1,3 @@
+# Data visualisation
+
+[Presentation (PDF)](./course_material/AVDP_03_DATAVis_GRANT.pdf)
diff --git a/57_SDM_History.Rmd b/57_SDM_History.Rmd
new file mode 100644
index 0000000..98f39f0
--- /dev/null
+++ b/57_SDM_History.Rmd
@@ -0,0 +1,3 @@
+# Recent and past application of SDMs
+
+[Presentation (PDF)](./course_material/C2T_Recent and past applications SDM.pdf)
diff --git a/60_Species_Distribution_Modelling.Rmd b/60_Species_Distribution_Modelling.Rmd
new file mode 100644
index 0000000..3ce51a5
--- /dev/null
+++ b/60_Species_Distribution_Modelling.Rmd
@@ -0,0 +1,1255 @@
+# Species distribution modelling
+
+Exercise scripts are at the bottom of this page.
+
+## SDM introduction and review
+
+What are species distribution models?
+
+[Presentation (PDF)](./course_material/Course 1 Tuesday What_are_SDM.pdf)
+
+### References
+
+- [Araujo and Guisan 2006](./course_material/literature/Araujo _ Guisan2006.pdf)
+- [Elith 2006](./course_material/literature/Elith 2006.pdf)
+- [Elith and Leathwick 2009](./course_material/literature/Elith and Leathwick 2009.pdf)
+- [Guisan and Zimmermann 2000](./course_material/literature/guisan and zimmermann 2000.pdf)
+- [Peterson 2011](./course_material/literature/Peterson 2011.pdf)
+- [Soberon and Peterson 2005](./course_material/literature/Soberon Peterson 2005.pdf)
+
+
+## SDM outputs
+
+Understand and generate SDM outputs.
+
+[Presentation (PDF)](./course_material/Course 3 Tuesday_ SDM outputs.pdf)
+
+### References
+
+ - [Elith et al. 2010](./course_material/literature/Elith_et_al-2010.pdf)
+ - [Guillaumot et al. 2018 Ecol & Evol](./course_material/literature/Guillaumot et al 2018-Ecology_and_Evolution.pdf)
+ - [Liu et al. 2013](./course_material/literature/liu 2013.pdf)
+ - [Torres 2015](./course_material/literature/Torres 2015.PDF)
+
+
+## Calibration
+
+What you should know before running a model.
+
+[Presentation (PDF)](./course_material/Course 5 Calibration- all you should think about before running a model.pdf)
+
+### References
+
+ - [Appendix S1 Elith 2008](./course_material/literature/Appendix S1 Elith 2008.pdf)
+ - [Barbet-Massin et al. 2012](./course_material/literature/Barbet-Massin_et_al-2012.pdf)
+ - [Barve 2011](./course_material/literature/Barve 2011.pdf)
+ - [Boria 2014](./course_material/literature/Boria 2014.pdf)
+ - [Brotons 2004](./course_material/literature/Brotons 2004.pdf)
+ - [Elith et al. 2008](./course_material/literature/Elith_et_al-2008.pdf)
+ - [Elith et al. 2010](./course_material/literature/Elith_et_al-2010.pdf)
+ - [Guillaumot et al. 2018 Ecol & Evol](./course_material/literature/Guillaumot et al 2018-Ecology_and_Evolution.pdf)
+ - [Guillaumot et al. 2019](./course_material/literature/Guillaumot et al. 2019_ Progress in Oceano.pdf)
+ - [Guillaumot et al. 2018 MEPS](./course_material/literature/Guillaumot et al 2018 MEPS final version.pdf)
+ - [Muscarella et al. 2014](./course_material/literature/Muscarella_et_al-2014-Methods_in_Ecology_and_Evolution.pdf)
+ - [Phillips et al. 2009](./course_material/literature/Phillips 2009.pdf)
+ - [Wisz and Guisan 2009](./course_material/literature/Wisz Guisan 2009.pdf)
+
+## Exercise scripts
+
+**Exercise - Run your first SDM**
+
+[SDM_EXERCISE.pdf](./course_material/exercises/SDM_EXERCISE.pdf)
+
+**Final practical work**
+
+Run models, compare different model calibrations and compare your model outputs!
+
+[SDM_FINAL.pdf](./course_material/exercises/SDM_FINAL.pdf)
+
+**Other scripts**
+
+[calibration_BRT.R](./course_material/exercises/scripts/calibration_BRT.R)
+
+```{r eval = FALSE}
+
+#### Open environmental data and occurrences (run_yOur_SDM.R file)
+#---------------------------------------------------------
+# compute there the MATRIX_OCC_ENVI variable (1 set)
+# in the run_yOur_SDM.R file, run until the MyFold definition
+# for set.seed(1) (j=1)
+# run the model and compare the predictive deviance according to the set of parameters chosen
+#=============================================================
+source("scripts/Function_gbm.R")
+tc=4 # tree complexity
+lr=0.011 # learning rate
+bf=0.8 # bag fraction
+
+model.res<- gbm.step_v2 (data=MATRIX_OCC_ENVI,
+ gbm.x = 2:ncol(MATRIX_OCC_ENVI),
+ gbm.y = 1,
+ family = "bernoulli",
+ n.folds=4,
+ fold.vector = MyFold,
+ tree.complexity = tc,
+ learning.rate = lr,
+ bag.fraction =bf)
+
+model1<-model.res # tc3 lr0.001 bf=0.8 nbt 10000
+model2<-model.res # tc4 lr0.001 bf=0.8 nbt 10000
+model3<-model.res # tc4 lr0.005 bf=0.8 nbt 10000
+model4<-model.res # tc3 lr0.005 bf=0.8 nbt 10000
+model5<-model.res # tc4 lr0.01 bf=0.8 nbt 10000
+model6<-model.res # tc4 lr0.008 bf=0.8 nbt 10000
+model7<-model.res # tc4 lr0.005 bf=0.9 nbt 10000
+model8<-model.res # tc4 lr0.011 bf=0.8 nbt 10000
+
+
+# GENERATE THE PLOTS
+tree.list1 <- seq(100, model1$gbm.call$best.trees, by = 100)
+tree.list2 <- seq(100, model2$gbm.call$best.trees, by = 100)
+tree.list3 <- seq(100, model3$gbm.call$best.trees, by = 100)
+tree.list4 <- seq(100, model4$gbm.call$best.trees, by = 100)
+tree.list5 <- seq(100, model5$gbm.call$best.trees, by = 100)
+tree.list6 <- seq(100, model6$gbm.call$best.trees, by = 100)
+tree.list7 <- seq(100, model7$gbm.call$best.trees, by = 100)
+tree.list8 <- seq(100, model8$gbm.call$best.trees, by = 100)
+
+pred1 <- predict.gbm(model1,MATRIX_OCC_ENVI, n.trees = tree.list1, "response")
+pred2 <- predict.gbm(model2,MATRIX_OCC_ENVI, n.trees = tree.list2, "response")
+pred3 <- predict.gbm(model3,MATRIX_OCC_ENVI, n.trees = tree.list3, "response")
+pred4 <- predict.gbm(model4,MATRIX_OCC_ENVI, n.trees = tree.list4, "response")
+pred5 <- predict.gbm(model5,MATRIX_OCC_ENVI, n.trees = tree.list5, "response")
+pred6 <- predict.gbm(model6,MATRIX_OCC_ENVI, n.trees = tree.list6, "response")
+pred7 <- predict.gbm(model7,MATRIX_OCC_ENVI, n.trees = tree.list7, "response")
+pred8 <- predict.gbm(model8,MATRIX_OCC_ENVI, n.trees = tree.list8, "response")
+
+graphe.deviance1 <- rep(0,max(tree.list1)/100)
+for (i in 1:length(graphe.deviance1)) {
+ graphe.deviance1 [i] <- calc.deviance(MATRIX_OCC_ENVI$id, pred1[,i],calc.mean=T)
+}
+graphe.deviance2 <- rep(0,max(tree.list2)/100)
+for (i in 1:length(graphe.deviance2)) {
+ graphe.deviance2 [i] <- calc.deviance(MATRIX_OCC_ENVI$id, pred2[,i],calc.mean=T)
+}
+graphe.deviance3 <- rep(0,max(tree.list3)/100)
+for (i in 1:length(graphe.deviance3)) {
+ graphe.deviance3 [i] <- calc.deviance(MATRIX_OCC_ENVI$id, pred3[,i],calc.mean=T)
+}
+graphe.deviance4 <- rep(0,max(tree.list4)/100)
+for (i in 1:length(graphe.deviance4)) {
+ graphe.deviance4 [i] <- calc.deviance(MATRIX_OCC_ENVI$id, pred4[,i],calc.mean=T)
+}
+graphe.deviance5 <- rep(0,max(tree.list5)/100)
+for (i in 1:length(graphe.deviance5)) {
+ graphe.deviance5 [i] <- calc.deviance(MATRIX_OCC_ENVI$id, pred5[,i],calc.mean=T)
+}
+graphe.deviance6 <- rep(0,max(tree.list6)/100)
+for (i in 1:length(graphe.deviance6)) {
+ graphe.deviance6 [i] <- calc.deviance(MATRIX_OCC_ENVI$id, pred6[,i],calc.mean=T)
+}
+graphe.deviance7 <- rep(0,max(tree.list7)/100)
+for (i in 1:length(graphe.deviance7)) {
+ graphe.deviance7 [i] <- calc.deviance(MATRIX_OCC_ENVI$id, pred7[,i],calc.mean=T)
+}
+graphe.deviance8 <- rep(0,max(tree.list8)/100)
+for (i in 1:length(graphe.deviance8)) {
+ graphe.deviance8 [i] <- calc.deviance(MATRIX_OCC_ENVI$id, pred8[,i],calc.mean=T)
+}
+
+par(mai=c(0.5,0.5,0.5,0.5))
+plot(tree.list1,graphe.deviance1,xlim = c(-100,4000), ylim=c(0,1),type='l', xlab = "number of trees",ylab = "predictive deviance", cex.lab = 1.5)
+lines(tree.list2,graphe.deviance2,col="blue")
+lines(tree.list3,graphe.deviance3,col="red")
+lines(tree.list4,graphe.deviance4,col="green")
+lines(tree.list5,graphe.deviance5,col="orange")
+lines(tree.list6,graphe.deviance6,col="pink")
+lines(tree.list7,graphe.deviance7,col="deepskyblue")
+lines(tree.list8,graphe.deviance8,col="purple")
+legend("topright",legend=c("tc3 lr0.001 bf=0.8","tc4 lr0.001 bf=0.8","tc4 lr0.005 bf=0.8 ","tc3 lr0.005 bf=0.8","tc4 lr0.01 bf=0.8","tc4 lr0.008 bf=0.8","tc4 lr0.005 bf=0.9","tc4 lr0.011 bf=0.8"),col=c("black","blue","red","green","orange","pink","deepskyblue","purple"),pch=16,cex=0.7)
+
+```
+
+[clock4_crossValidation.R](./course_material/exercises/scripts/clock4_crossValidation.R)
+
+```{r eval = FALSE}
+clock4 <- function(occ, bg.coords){
+ # the spatial sampling aims at defining areas containing either test or training data
+ # for the CLOCK-4 method, 4 triangle areas are defined, cutting the Southern Ocean circle into 4 equal areas. Three of these 4 areas contain the data that will be used to train the model, the other remaining one, the data that will be used to test the model after predictions
+ # the areas used to train and test the model are randomly defined at each loop step
+
+ #initialise vectors that are used to generate the spatial sampling structure
+ random_longA <- seq(0,179,1)
+ random_longB <- seq(-180,-1,1)
+ random_longC <- seq(0,180,1)
+ random_longD <- seq(-179,0,1)
+ random_long <- c(random_longA,random_longB,random_longC,random_longD)
+
+ # sample a number between -180 and 180 to define the random sampling transect
+ tirage <- sample(seq(181,541,1),1)
+ random_long_tirage <- random_long[tirage]
+ random_long_tirage_A1 <- random_long[tirage+90]
+ random_long_tirage_A2 <- random_long[tirage+180]
+ random_long_tirage_A3 <- random_long[tirage-90]
+ #random_long_tirage_A4 <- random_long[tirage-180]
+
+
+ ## define training and test groups (composed of both presence and background data)
+ occ.grp <- rep(NA, nrow(occ))
+ bg.coords.grp <- rep(NA, nrow(bg.coords))
+
+ ## define training and test groups (composed of both presence and background data)
+ presence_tot <- occ
+ background_data <- bg.coords
+ training_presences.occ1 <- NA;training_presences.occ3 <-NA
+ training_presences.occ1_part1 <- NA;training_presences.occ3_part1 <-NA
+ training_presences.occ1_part2 <- NA;training_presences.occ3_part2 <-NA
+ training_backgr.occ1 <- NA; training_backgr.occ3 <- NA
+ training_backgr.occ1_part1 <- NA; training_backgr.occ3_part1 <- NA
+ training_backgr.occ1_part2 <- NA; training_backgr.occ3_part2 <- NA
+ training_presences.occ2 <- NA;training_presences.occ4 <- NA
+ training_presences.occ2_part1 <- NA;training_presences.occ4_part1 <- NA
+ training_presences.occ2_part2 <- NA;training_presences.occ4_part2 <- NA
+ training_backgr.occ2 <- NA;training_backgr.occ4 <- NA
+ training_backgr.occ2_part1 <- NA;training_backgr.occ4_part1 <- NA
+ training_backgr.occ2_part2 <- NA;training_backgr.occ4_part2 <- NA
+
+
+ #---------------------
+ # TRAINING PRESENCES
+ #---------------------
+
+ ### ZONE 1 ####
+ if (random_long_tirage <= 0 & random_long_tirage_A1 <=0){ # imply that A2 >= 0 and A3>=0
+ training_presences.occ1 <- which(presence_tot[,1] > random_long_tirage & presence_tot[,1] < random_long_tirage_A1)
+ training_backgr.occ1 <- which(background_data[,1] > random_long_tirage & background_data[,1] < random_long_tirage_A1)
+ }
+ if (random_long_tirage >= 0 & random_long_tirage_A1 >=0){
+ training_presences.occ1 <- which(presence_tot[,1] > random_long_tirage & presence_tot[,1] < random_long_tirage_A1)
+ training_backgr.occ1 <- which(background_data[,1] > random_long_tirage & background_data[,1] < random_long_tirage_A1)
+ }
+
+ if (random_long_tirage >= 0 & random_long_tirage_A1 <=0){
+ training_presences.occ1_part1 <- which(presence_tot[,1] > random_long_tirage & presence_tot[,1] < 180 )
+ training_presences.occ1_part2 <- which(presence_tot[,1] > -180 & presence_tot[,1] < random_long_tirage_A1 )
+ training_backgr.occ1_part1 <- which(background_data[,1] > random_long_tirage & background_data[,1] < 180 )
+ training_backgr.occ1_part2 <- which(background_data[,1] > -180 & background_data[,1] < random_long_tirage_A1 )
+ }
+
+ if (random_long_tirage <= 0 & random_long_tirage_A1 >=0){
+ training_presences.occ1_part1 <- which(presence_tot[,1] > random_long_tirage & presence_tot[,1] < 0 )
+ training_presences.occ1_part2 <- which(presence_tot[,1] > 0 & presence_tot[,1] < random_long_tirage_A1 )
+ training_backgr.occ1_part1 <- which(background_data[,1] > random_long_tirage & background_data[,1] < 0 )
+ training_backgr.occ1_part2 <- which(background_data[,1] > 0 & background_data[,1] < random_long_tirage_A1 )
+ }
+
+ ### ZONE 2 ####
+
+ if (random_long_tirage_A1 <= 0 & random_long_tirage_A2 <=0){
+ training_presences.occ2 <- which(presence_tot[,1] > random_long_tirage_A1 & presence_tot[,1] < random_long_tirage_A2 )
+ training_backgr.occ2 <- which(background_data[,1] > random_long_tirage_A1 & background_data[,1] < random_long_tirage_A2)
+ }
+ if (random_long_tirage_A1 >= 0 & random_long_tirage_A2 >=0){
+ training_presences.occ2 <- which(presence_tot[,1] > random_long_tirage_A1 & presence_tot[,1] < random_long_tirage_A2 )
+ training_backgr.occ2 <- which(background_data[,1] > random_long_tirage_A1 & background_data[,1] < random_long_tirage_A2)
+ }
+
+ if (random_long_tirage_A1 >= 0 & random_long_tirage_A2 <=0){
+ training_presences.occ2_part1 <- which(presence_tot[,1] > random_long_tirage_A1 & presence_tot[,1] < 180)
+ training_presences.occ2_part2 <- which(presence_tot[,1] > -180 & presence_tot[,1] < random_long_tirage_A2)
+ training_backgr.occ2_part1 <- which(background_data[,1] > random_long_tirage_A1 & background_data[,1] < 180)
+ training_backgr.occ2_part2 <- which(background_data[,1] > -180 & background_data[,1] < random_long_tirage_A2)
+ }
+
+ if (random_long_tirage_A1 <= 0 & random_long_tirage_A2 >=0){
+ training_presences.occ2_part1 <- which(presence_tot[,1] > random_long_tirage_A1 & presence_tot[,1] < 0 )
+ training_presences.occ2_part2 <- which(presence_tot[,1] > 0 & presence_tot[,1] < random_long_tirage_A2 )
+ training_backgr.occ2_part1 <- which(background_data[,1] > random_long_tirage_A1 & background_data[,1] < 0 )
+ training_backgr.occ2_part2 <- which(background_data[,1] > 0 & background_data[,1] < random_long_tirage_A2 )
+ }
+
+ ### ZONE 3 ####
+ if (random_long_tirage_A2 <= 0 & random_long_tirage_A3 <=0){
+ training_presences.occ3 <- which(presence_tot[,1] > random_long_tirage_A2 & presence_tot[,1] < random_long_tirage_A3 )
+ training_backgr.occ3 <- which(background_data[,1] > random_long_tirage_A2 & background_data[,1] < random_long_tirage_A3)
+ }
+ if (random_long_tirage_A2 >= 0 & random_long_tirage_A3 >=0){
+ training_presences.occ3 <- which(presence_tot[,1] > random_long_tirage_A2 & presence_tot[,1] < random_long_tirage_A3 )
+ training_backgr.occ3 <- which(background_data[,1] > random_long_tirage_A2 & background_data[,1] < random_long_tirage_A3)
+ }
+
+ if (random_long_tirage_A2 >= 0 & random_long_tirage_A3 <=0){
+ training_presences.occ3_part1 <- which(presence_tot[,1] > random_long_tirage_A2 & presence_tot[,1] < 180 )
+ training_presences.occ3_part2 <- which(presence_tot[,1] > -180 & presence_tot[,1] < random_long_tirage_A3 )
+ training_backgr.occ3_part1 <- which(background_data[,1] > random_long_tirage_A2 & background_data[,1] < 180 )
+ training_backgr.occ3_part2 <- which(background_data[,1] > -180 & background_data[,1] < random_long_tirage_A3 )
+ }
+
+ if (random_long_tirage_A2 <= 0 & random_long_tirage_A3 >=0){
+ training_presences.occ3_part1 <- which(presence_tot[,1] > random_long_tirage_A2 & presence_tot[,1] < 0 )
+ training_presences.occ3_part2 <- which(presence_tot[,1] > 0 & presence_tot[,1] < random_long_tirage_A3 )
+ training_backgr.occ3_part1 <- which(background_data[,1] > random_long_tirage_A2 & background_data[,1] < 0 )
+ training_backgr.occ3_part2 <- which(background_data[,1] > 0 & background_data[,1] < random_long_tirage_A3 )
+ }
+
+ ### set the groups
+ training_presence_grp1_all <- as.vector(na.omit(c(training_presences.occ1,training_presences.occ1_part1,training_presences.occ1_part2)))
+ training_presence_grp2_all <- as.vector(na.omit(c(training_presences.occ2,training_presences.occ2_part1,training_presences.occ2_part2)))
+ training_presence_grp3_all <- as.vector(na.omit(c(training_presences.occ3,training_presences.occ3_part1,training_presences.occ3_part2)))
+
+
+ for (i in training_presence_grp1_all){
+ occ.grp[i] <- 1
+ }
+ for (i in training_presence_grp2_all){
+ occ.grp[i] <- 2
+ }
+ for (i in training_presence_grp3_all){
+ occ.grp[i] <- 3
+ }
+ occ.grp[which(is.na(occ.grp))]=4
+
+
+ training_backgr_grp1_all <- as.vector(na.omit(c(training_backgr.occ1,training_backgr.occ1_part1,training_backgr.occ1_part2)))
+ training_backgr_grp2_all <- as.vector(na.omit(c(training_backgr.occ2,training_backgr.occ2_part1,training_backgr.occ2_part2)))
+ training_backgr_grp3_all <- as.vector(na.omit(c(training_backgr.occ3,training_backgr.occ3_part1,training_backgr.occ3_part2)))
+
+ for (i in training_backgr_grp1_all){
+ bg.coords.grp[i] <- 1
+ }
+ for (i in training_backgr_grp2_all){
+ bg.coords.grp[i] <- 2
+ }
+ for (i in training_backgr_grp3_all){
+ bg.coords.grp[i] <- 3
+ }
+ bg.coords.grp[which(is.na(bg.coords.grp))]=4
+
+
+ out <- list(occ.grp=occ.grp,bg.coords.grp= bg.coords.grp, tirage=tirage)
+ return(out)
+}
+
+
+```
+
+[Function_gbm.R](./course_material/exercises/scripts/Function_gbm.R)
+
+```{r eval = FALSE}
+### gbm.stepv2
+
+.roc <-function (obsdat, preddat) {
+ # code adapted from Ferrier, Pearce and Watson's code, by J.Elith
+ #
+ # see:
+ # Hanley, J.A. & McNeil, B.J. (1982) The meaning and use of the area
+ # under a Receiver Operating Characteristic (ROC) curve.
+ # Radiology, 143, 29-36
+ #
+ # Pearce, J. & Ferrier, S. (2000) Evaluating the predictive performance
+ # of habitat models developed using logistic regression.
+ # Ecological Modelling, 133, 225-245.
+ # this is the non-parametric calculation for area under the ROC curve,
+ # using the fact that a MannWhitney U statistic is closely related to
+ # the area
+ #
+
+ # in dismo, this is used in the gbm routines, but not elsewhere (see evaluate).
+
+ if (length(obsdat) != length(preddat)) {
+ stop("obs and preds must be equal lengths")
+ }
+ n.x <- length(obsdat[obsdat == 0])
+ n.y <- length(obsdat[obsdat == 1])
+ xy <- c(preddat[obsdat == 0], preddat[obsdat == 1])
+ rnk <- rank(xy)
+ wilc <- ((n.x * n.y) + ((n.x * (n.x + 1))/2) - sum(rnk[1:n.x]))/(n.x * n.y)
+ return(round(wilc, 4))
+}
+
+
+.calibration <- function(obs, preds, family = "binomial") {
+ #
+ # j elith/j leathwick 17th March 2005
+ # calculates calibration statistics for either binomial or count data
+ # but the family argument must be specified for the latter
+ # a conditional test for the latter will catch most failures to specify
+ # the family
+ #
+
+ if (family == "bernoulli") {
+ family <- "binomial"
+ }
+
+ pred.range <- max(preds) - min(preds)
+
+ if(pred.range > 1.2 & family == "binomial") {
+ print(paste("range of response variable is ", round(pred.range, 2)), sep = "", quote = F)
+ print("check family specification", quote = F)
+ return()
+ }
+
+ if(family == "binomial") {
+ pred <- preds + 1e-005
+ pred[pred >= 1] <- 0.99999
+ mod <- glm(obs ~ log((pred)/(1 - (pred))), family = binomial)
+ lp <- log((pred)/(1 - (pred)))
+ a0b1 <- glm(obs ~ offset(lp) - 1, family = binomial)
+ miller1 <- 1 - pchisq(a0b1$deviance - mod$deviance, 2)
+ ab1 <- glm(obs ~ offset(lp), family = binomial)
+ miller2 <- 1 - pchisq(a0b1$deviance - ab1$deviance, 1)
+ miller3 <- 1 - pchisq(ab1$deviance - mod$deviance, 1)
+ } else if(family == "poisson") {
+ mod <- glm(obs ~ log(preds), family = poisson)
+ lp <- log(preds)
+ a0b1 <- glm(obs ~ offset(lp) - 1, family = poisson)
+ miller1 <- 1 - pchisq(a0b1$deviance - mod$deviance, 2)
+ ab1 <- glm(obs ~ offset(lp), family = poisson)
+ miller2 <- 1 - pchisq(a0b1$deviance - ab1$deviance, 1)
+ miller3 <- 1 - pchisq(ab1$deviance - mod$deviance, 1)
+ }
+
+ calibration.result <- c(mod$coef, miller1, miller2, miller3)
+ names(calibration.result) <- c("intercept", "slope", "testa0b1", "testa0|b1", "testb1|a")
+ return(calibration.result)
+}
+
+## VI.
+gbm.step_v2 <- function (data, gbm.x, gbm.y, offset = NULL, fold.vector = NULL,
+ tree.complexity = 1, learning.rate = 0.01, bag.fraction = 0.75,
+ site.weights = rep(1, nrow(data)), var.monotone = rep(0,
+ length(gbm.x)),
+ n.folds = 10, prev.stratify = TRUE, family = "bernoulli",
+ n.trees = 50, step.size = n.trees, max.trees = 10000, tolerance.method = "auto",
+ tolerance = 0.001, plot.main = TRUE, plot.folds = FALSE,
+ verbose = TRUE, silent = FALSE, keep.fold.models = FALSE,
+ keep.fold.vector = FALSE, keep.fold.fit = FALSE, ...)
+{
+ if (!requireNamespace("gbm")) {
+ stop("you need to install the gbm package to run this function")
+ }
+ requireNamespace("splines")
+ if (silent)
+ verbose <- FALSE
+ z1 <- Sys.time()
+ x.data <- data[, gbm.x, drop = FALSE]
+ y.data <- data[, gbm.y]
+ sp.name <- names(data)[gbm.y]
+ if (family == "bernoulli") {
+ prevalence <- mean(y.data)
+ }
+ n.cases <- nrow(data)
+ n.preds <- length(gbm.x)
+ if (!silent) {
+ cat("\n", "\n", "GBM STEP - version 2.9", "\n", "\n")
+ cat("Performing cross-validation optimisation of a boosted regression tree model \n")
+ cat("for", sp.name, "and using a family of", family,
+ "\n")
+ cat("Using", n.cases, "observations and", n.preds, "predictors \n")
+ }
+ if (is.null(fold.vector)) {
+ if (prev.stratify & family == "bernoulli") {
+ presence.mask <- data[, gbm.y] == 1
+ absence.mask <- data[, gbm.y] == 0
+ n.pres <- sum(presence.mask)
+ n.abs <- sum(absence.mask)
+ selector <- rep(0, n.cases)
+ temp <- rep(seq(1, n.folds, by = 1), length = n.pres)
+ temp <- temp[order(runif(n.pres, 1, 100))]
+ selector[presence.mask] <- temp
+ temp <- rep(seq(1, n.folds, by = 1), length = n.abs)
+ temp <- temp[order(runif(n.abs, 1, 100))]
+ selector[absence.mask] <- temp
+ }
+ else {
+ selector <- rep(seq(1, n.folds, by = 1), length = n.cases)
+ selector <- selector[order(runif(n.cases, 1, 100))]
+ }
+ } else {
+ if (length(fold.vector) != n.cases) {
+ stop("supplied fold vector is of wrong length")
+ }
+ cat("loading user-supplied fold vector \n")
+ selector <- fold.vector
+ }
+ pred.values <- rep(0, n.cases)
+ cv.loss.matrix <- matrix(0, nrow = n.folds, ncol = 1)
+ training.loss.matrix <- matrix(0, nrow = n.folds, ncol = 1)
+ trees.fitted <- n.trees
+ model.list <- list(paste("model", 1:n.folds, sep = ""))
+ if (is.null(offset)) {
+ gbm.call <- paste("gbm::gbm(y.subset ~ .,data=x.subset, n.trees = n.trees, interaction.depth = tree.complexity, shrinkage = learning.rate, bag.fraction = bag.fraction, weights = weight.subset, distribution = as.character(family), var.monotone = var.monotone, verbose = FALSE)",
+ sep = "")
+ } else {
+ gbm.call <- paste("gbm::gbm(y.subset ~ . + offset(offset.subset), data=x.subset, n.trees = n.trees, interaction.depth = tree.complexity, shrinkage = learning.rate, bag.fraction = bag.fraction, weights = weight.subset, distribution = as.character(family), var.monotone = var.monotone, verbose = FALSE)",
+ sep = "")
+ }
+ n.fitted <- n.trees
+ y_i <- y.data
+ u_i <- sum(y.data * site.weights)/sum(site.weights)
+ u_i <- rep(u_i, length(y_i))
+ total.deviance <- calc.deviance(y_i, u_i, weights = site.weights,
+ family = family, calc.mean = FALSE)
+ mean.total.deviance <- total.deviance/n.cases
+ tolerance.test <- tolerance
+ if (tolerance.method == "auto") {
+ tolerance.test <- mean.total.deviance * tolerance
+ }
+ if (!silent) {
+ cat("creating", n.folds, "initial models of", n.trees,
+ "trees", "\n")
+ if (prev.stratify & family == "bernoulli") {
+ cat("\n", "folds are stratified by prevalence", "\n")
+ } else {
+ cat("\n", "folds are unstratified", "\n")
+ }
+ cat("total mean deviance = ", round(mean.total.deviance,
+ 4), "\n")
+ cat("tolerance is fixed at ", round(tolerance.test, 4),
+ "\n")
+ if (tolerance.method != "fixed" & tolerance.method !=
+ "auto") {
+ stop("invalid argument for tolerance method - should be auto or fixed")
+ }
+ }
+ if (verbose) {
+ cat("ntrees resid. dev.", "\n")
+ }
+
+ for (i in 1:n.folds) {
+ model.mask <- selector != i
+ pred.mask <- selector == i
+ y.subset <- y.data[model.mask]
+ x.subset <- x.data[model.mask, , drop = FALSE]
+ weight.subset <- site.weights[model.mask]
+ if (!is.null(offset)) {
+ offset.subset <- offset[model.mask]
+ }
+ else {
+ offset.subset <- NULL
+ }
+ model.list[[i]] <- eval(parse(text = gbm.call))
+ fitted.values <- model.list[[i]]$fit
+ if (!is.null(offset)) {
+ fitted.values <- fitted.values + offset[model.mask]
+ }
+ if (family == "bernoulli") {
+ fitted.values <- exp(fitted.values)/(1 + exp(fitted.values))
+ }
+ else if (family == "poisson") {
+ fitted.values <- exp(fitted.values)
+ }
+ pred.values[pred.mask] <- gbm::predict.gbm(model.list[[i]],
+ x.data[pred.mask, , drop = FALSE], n.trees = n.trees)
+ if (!is.null(offset)) {
+ pred.values[pred.mask] <- pred.values[pred.mask] +
+ offset[pred.mask]
+ }
+ if (family == "bernoulli") {
+ pred.values[pred.mask] <- exp(pred.values[pred.mask])/(1 +
+ exp(pred.values[pred.mask]))
+ }
+ else if (family == "poisson") {
+ pred.values[pred.mask] <- exp(pred.values[pred.mask])
+ }
+ y_i <- y.subset
+ u_i <- fitted.values
+ weight.fitted <- site.weights[model.mask]
+ training.loss.matrix[i, 1] <- calc.deviance(y_i, u_i,
+ weight.fitted, family = family)
+ y_i <- y.data[pred.mask]
+ u_i <- pred.values[pred.mask]
+ weight.preds <- site.weights[pred.mask]
+ cv.loss.matrix[i, 1] <- calc.deviance(y_i, u_i, weight.preds,
+ family = family)
+ }
+ delta.deviance <- 1
+ cv.loss.values <- apply(cv.loss.matrix, 2, mean)
+ if (verbose) {
+ cat(n.fitted, " ", round(cv.loss.values, 4), "\n")
+ }
+ if (!silent) {
+ cat("now adding trees...", "\n")
+ }
+ j <- 1
+ while (delta.deviance > tolerance.test & n.fitted < max.trees) {
+ training.loss.matrix <- cbind(training.loss.matrix, rep(0,
+ n.folds))
+ cv.loss.matrix <- cbind(cv.loss.matrix, rep(0, n.folds))
+ n.fitted <- n.fitted + step.size
+ trees.fitted <- c(trees.fitted, n.fitted)
+ j <- j + 1
+ for (i in 1:n.folds) {
+ model.mask <- selector != i
+ pred.mask <- selector == i
+ y.subset <- y.data[model.mask]
+ x.subset <- x.data[model.mask, , drop = FALSE]
+ weight.subset <- site.weights[model.mask]
+ if (!is.null(offset)) {
+ offset.subset <- offset[model.mask]
+ }
+ model.list[[i]] <- gbm::gbm.more(model.list[[i]],
+ weights = weight.subset, step.size)
+ fitted.values <- model.list[[i]]$fit
+ if (!is.null(offset)) {
+ fitted.values <- fitted.values + offset[model.mask]
+ }
+ if (family == "bernoulli") {
+ fitted.values <- exp(fitted.values)/(1 + exp(fitted.values))
+ }
+ else if (family == "poisson") {
+ fitted.values <- exp(fitted.values)
+ }
+ pred.values[pred.mask] <- gbm::predict.gbm(model.list[[i]],
+ x.data[pred.mask, , drop = FALSE], n.trees = n.fitted)
+ if (!is.null(offset)) {
+ pred.values[pred.mask] <- pred.values[pred.mask] +
+ offset[pred.mask]
+ }
+ if (family == "bernoulli") {
+ pred.values[pred.mask] <- exp(pred.values[pred.mask])/(1 +
+ exp(pred.values[pred.mask]))
+ }
+ else if (family == "poisson") {
+ pred.values[pred.mask] <- exp(pred.values[pred.mask])
+ }
+ y_i <- y.subset
+ u_i <- fitted.values
+ weight.fitted <- site.weights[model.mask]
+ training.loss.matrix[i, j] <- calc.deviance(y_i,
+ u_i, weight.fitted, family = family)
+ u_i <- pred.values[pred.mask]
+ y_i <- y.data[pred.mask]
+ weight.preds <- site.weights[pred.mask]
+ cv.loss.matrix[i, j] <- calc.deviance(y_i, u_i, weight.preds,
+ family = family)
+ }
+ cv.loss.values <- apply(cv.loss.matrix, 2, mean)
+ if (j < 5) {
+ if (cv.loss.values[j] > cv.loss.values[j - 1]) {
+ if (!silent) {
+ cat("restart model with a smaller learning rate or smaller step size...")
+ }
+ return()
+ }
+ }
+ if (j >= 20) {
+ test1 <- mean(cv.loss.values[(j - 9):j])
+ test2 <- mean(cv.loss.values[(j - 19):(j - 9)])
+ delta.deviance <- test2 - test1
+ }
+ if (verbose) {
+ cat(n.fitted, " ", round(cv.loss.values[j], 4), "\n")
+ flush.console()
+ }
+ }
+ training.loss.values <- apply(training.loss.matrix, 2, mean)
+ cv.loss.ses <- rep(0, length(cv.loss.values))
+ cv.loss.ses <- sqrt(apply(cv.loss.matrix, 2, var))/sqrt(n.folds)
+ y.bar <- min(cv.loss.values)
+ target.trees <- trees.fitted[match(TRUE, cv.loss.values ==
+ y.bar)]
+ if (plot.main) {
+ y.min <- min(cv.loss.values - cv.loss.ses)
+ y.max <- max(cv.loss.values + cv.loss.ses)
+ if (plot.folds) {
+ y.min <- min(cv.loss.matrix)
+ y.max <- max(cv.loss.matrix)
+ }
+ plot(trees.fitted, cv.loss.values, type = "l", ylab = "holdout deviance",
+ xlab = "no. of trees", ylim = c(y.min, y.max), ...)
+ abline(h = y.bar, col = 2)
+ lines(trees.fitted, cv.loss.values + cv.loss.ses, lty = 2)
+ lines(trees.fitted, cv.loss.values - cv.loss.ses, lty = 2)
+ if (plot.folds) {
+ for (i in 1:n.folds) {
+ lines(trees.fitted, cv.loss.matrix[i, ], lty = 3)
+ }
+ }
+ abline(v = target.trees, col = 3)
+ title(paste(sp.name, ", d - ", tree.complexity, ", lr - ",
+ learning.rate, sep = ""))
+ }
+ cv.deviance.stats <- rep(0, n.folds)
+ cv.roc.stats <- rep(0, n.folds)
+ cv.cor.stats <- rep(0, n.folds)
+ cv.calibration.stats <- matrix(0, ncol = 5, nrow = n.folds)
+
+
+ if (family == "bernoulli") {
+ threshold.stats <- rep(0, n.folds)
+
+
+
+ ############################
+ ## JA editions / 08/02/19 ##
+ ############################
+
+ cv.corr.class <- rep(0, n.folds)
+ cv.length <- rep(0, n.folds)
+ tss.stat <- rep(0, n.folds)
+
+ ####################
+ ## End / 08/02/19 ##
+ ####################
+
+
+ }
+ fitted.matrix <- matrix(NA, nrow = n.cases, ncol = n.folds)
+ fold.fit <- rep(0, n.cases)
+ for (i in 1:n.folds) {
+ pred.mask <- selector == i
+ model.mask <- selector != i
+ fits <- gbm::predict.gbm(model.list[[i]], x.data[model.mask,
+ , drop = FALSE], n.trees = target.trees)
+ if (!is.null(offset)) {
+ fits <- fits + offset[model.mask]
+ }
+ if (family == "bernoulli") {
+ fits <- exp(fits)/(1 + exp(fits))
+ } else if (family == "poisson") {
+ fits <- exp(fits)
+ }
+ fitted.matrix[model.mask, i] <- fits
+ fits <- gbm::predict.gbm(model.list[[i]], x.data[pred.mask,
+ , drop = FALSE], n.trees = target.trees)
+ if (!is.null(offset))
+ fits <- fits + offset[pred.mask]
+ fold.fit[pred.mask] <- fits
+ if (family == "bernoulli") {
+ fits <- exp(fits)/(1 + exp(fits))
+ }else if (family == "poisson") {
+ fits <- exp(fits)
+ }
+ fitted.matrix[pred.mask, i] <- fits
+ y_i <- y.data[pred.mask]
+ u_i <- fitted.matrix[pred.mask, i]
+ weight.preds <- site.weights[pred.mask]
+ cv.deviance.stats[i] <- calc.deviance(y_i, u_i, weight.preds,
+ family = family)
+ cv.cor.stats[i] <- cor(y_i, u_i)
+ if (family == "bernoulli") {
+ cv.roc.stats[i] <- .roc(y_i, u_i)
+ cv.calibration.stats[i, ] <- .calibration(y_i, u_i,
+ "binomial")
+ threshold.stats[i] <- approx(ppoints(u_i), sort(u_i,
+ decreasing = T), prevalence)$y
+
+ ############################
+ ## JA editions / 08/02/19 ##
+ ############################
+
+ c_i <- c(1, 0)[as.factor(u_i <= threshold.stats[i])]
+ tab1 <- table(c_i, y_i)
+
+ # Correctly classified test data
+ cv.corr.class[i] <- sum(diag(tab1)) / sum(tab1)
+
+ # % cv / total data
+ cv.length[i] <- table(y_i)["1"] / table(y.data)["1"]
+
+ # TSS
+ # Sensitivity = TP/(TP + FN)
+ # Specificity = TN/(TN + FP)
+ tp <- tab1["1", "1"]
+ fn <- tab1["0", "1"]
+ tn <- tab1["0", "0"]
+ fp <- tab1["1", "0"]
+
+ sen <- tp/(tp + fn)
+ spe <- tn/(tn + fp)
+
+ # TSS
+ tss.stat[i] <- spe + sen -1
+
+
+ ####################
+ ## End / 08/02/19 ##
+ ####################
+
+
+ }
+ if (family == "poisson") {
+ cv.calibration.stats[i, ] <- .calibration(y_i, u_i,
+ "poisson")
+ }
+ }
+ fitted.vars <- apply(fitted.matrix, 1, var, na.rm = TRUE)
+ cv.dev <- mean(cv.deviance.stats, na.rm = TRUE)
+ cv.dev.se <- sqrt(var(cv.deviance.stats))/sqrt(n.folds)
+ cv.cor <- mean(cv.cor.stats, na.rm = TRUE)
+ cv.cor.se <- sqrt(var(cv.cor.stats, use = "complete.obs"))/sqrt(n.folds)
+ cv.roc <- 0
+ cv.roc.se <- 0
+ if (family == "bernoulli") {
+ cv.roc <- mean(cv.roc.stats, na.rm = TRUE)
+ cv.roc.se <- sqrt(var(cv.roc.stats, use = "complete.obs"))/sqrt(n.folds)
+ cv.threshold <- mean(threshold.stats, na.rm = T)
+ cv.threshold.se <- sqrt(var(threshold.stats, use = "complete.obs"))/sqrt(n.folds)
+ }
+ cv.calibration <- 0
+ cv.calibration.se <- 0
+ if (family == "poisson" | family == "bernoulli") {
+ cv.calibration <- apply(cv.calibration.stats, 2, mean)
+ cv.calibration.se <- apply(cv.calibration.stats, 2, var)
+ cv.calibration.se <- sqrt(cv.calibration.se)/sqrt(n.folds)
+ }
+ if (is.null(offset)) {
+ gbm.call <- paste("gbm::gbm(y.data ~ .,data=x.data, n.trees = target.trees, interaction.depth = tree.complexity, shrinkage = learning.rate, bag.fraction = bag.fraction, weights = site.weights, distribution = as.character(family), var.monotone = var.monotone, verbose = FALSE)",
+ sep = "")
+ }
+ else {
+ gbm.call <- paste("gbm::gbm(y.data ~ . + offset(offset),data=x.data, n.trees = target.trees, interaction.depth = tree.complexity, shrinkage = learning.rate, bag.fraction = bag.fraction, weights = site.weights, distribution = as.character(family), var.monotone = var.monotone, verbose = FALSE)",
+ sep = "")
+ }
+ if (!silent) {
+ message("fitting final gbm model with a fixed number of ",
+ target.trees, " trees for ", sp.name)
+ }
+ gbm.object <- eval(parse(text = gbm.call))
+ best.trees <- target.trees
+ gbm.summary <- summary(gbm.object, n.trees = target.trees,
+ plotit = FALSE)
+ fits <- gbm::predict.gbm(gbm.object, x.data, n.trees = target.trees)
+ if (!is.null(offset))
+ fits <- fits + offset
+ if (family == "bernoulli") {
+ fits <- exp(fits)/(1 + exp(fits))
+ }
+ else if (family == "poisson") {
+ fits <- exp(fits)
+ }
+ fitted.values <- fits
+ y_i <- y.data
+ u_i <- fitted.values
+ resid.deviance <- calc.deviance(y_i, u_i, weights = site.weights,
+ family = family, calc.mean = FALSE)
+ self.cor <- cor(y_i, u_i)
+ self.calibration <- 0
+ self.roc <- 0
+ if (family == "bernoulli") {
+ deviance.contribs <- (y_i * log(u_i)) + ((1 - y_i) *
+ log(1 - u_i))
+ residuals <- sqrt(abs(deviance.contribs * 2))
+ residuals <- ifelse((y_i - u_i) < 0, 0 - residuals, residuals)
+ self.roc <- .roc(y_i, u_i)
+ self.calibration <- .calibration(y_i, u_i, "binomial")
+
+
+
+ ############################
+ ## JA editions / 08/02/19 ##
+ ############################
+
+ self.th <- approx(ppoints(u_i), sort(u_i,
+ decreasing = T), prevalence)$y
+
+ ####################
+ ## End / 08/02/19 ##
+ ####################
+
+ }
+ if (family == "poisson") {
+ deviance.contribs <- ifelse(y_i == 0, 0, (y_i * log(y_i/u_i))) -
+ (y_i - u_i)
+ residuals <- sqrt(abs(deviance.contribs * 2))
+ residuals <- ifelse((y_i - u_i) < 0, 0 - residuals, residuals)
+ self.calibration <- .calibration(y_i, u_i, "poisson")
+ }
+ if (family == "gaussian" | family == "laplace") {
+ residuals <- y_i - u_i
+ }
+ mean.resid.deviance <- resid.deviance/n.cases
+ z2 <- Sys.time()
+ elapsed.time.minutes <- round(as.numeric(z2 - z1)/60, 2)
+ if (verbose) {
+ cat("\n")
+ cat("mean total deviance =", round(mean.total.deviance,
+ 3), "\n")
+ cat("mean residual deviance =", round(mean.resid.deviance,
+ 3), "\n", "\n")
+ cat("estimated cv deviance =", round(cv.dev, 3), "; se =",
+ round(cv.dev.se, 3), "\n", "\n")
+ cat("training data correlation =", round(self.cor, 3),
+ "\n")
+ cat("cv correlation = ", round(cv.cor, 3), "; se =",
+ round(cv.cor.se, 3), "\n", "\n")
+ if (family == "bernoulli") {
+ cat("training data AUC score =", round(self.roc,
+ 3), "\n")
+ cat("cv AUC score =", round(cv.roc, 3), "; se =",
+ round(cv.roc.se, 3), "\n", "\n")
+ }
+ cat("elapsed time - ", round(elapsed.time.minutes, 2),
+ "minutes", "\n")
+ }
+ if (n.fitted == max.trees & !silent) {
+ cat("\n", "########### warning ##########", "\n", "\n")
+ cat("maximum tree limit reached - results may not be optimal",
+ "\n")
+ cat(" - refit with faster learning rate or increase maximum number of trees",
+ "\n")
+ }
+ gbm.detail <- list(dataframe = data, gbm.x = gbm.x, predictor.names = names(x.data),
+ gbm.y = gbm.y, response.name = sp.name, offset = offset,
+ family = family, tree.complexity = tree.complexity, learning.rate = learning.rate,
+ bag.fraction = bag.fraction, cv.folds = n.folds, prev.stratification = prev.stratify,
+ max.fitted = n.fitted, n.trees = target.trees, best.trees = target.trees,
+ train.fraction = 1, tolerance.method = tolerance.method,
+ tolerance = tolerance, var.monotone = var.monotone, date = date(),
+ elapsed.time.minutes = elapsed.time.minutes)
+ training.stats <- list(null = total.deviance, mean.null = mean.total.deviance,
+ resid = resid.deviance, mean.resid = mean.resid.deviance,
+ correlation = self.cor, discrimination = self.roc, calibration = self.calibration)
+ cv.stats <- list(deviance.mean = cv.dev, deviance.se = cv.dev.se,
+ correlation.mean = cv.cor, correlation.se = cv.cor.se,
+ discrimination.mean = cv.roc, discrimination.se = cv.roc.se,
+ calibration.mean = cv.calibration, calibration.se = cv.calibration.se)
+ if (family == "bernoulli") {
+ cv.stats$cv.threshold <- cv.threshold
+ cv.stats$cv.threshold.se <- cv.threshold.se
+ }
+ gbm.object$gbm.call <- gbm.detail
+ gbm.object$fitted <- fitted.values
+ gbm.object$fitted.vars <- fitted.vars
+ gbm.object$residuals <- residuals
+ gbm.object$contributions <- gbm.summary
+ gbm.object$self.statistics <- training.stats
+ gbm.object$cv.statistics <- cv.stats
+ gbm.object$weights <- site.weights
+ gbm.object$trees.fitted <- trees.fitted
+ gbm.object$training.loss.values <- training.loss.values
+ gbm.object$cv.values <- cv.loss.values
+ gbm.object$cv.loss.ses <- cv.loss.ses
+ gbm.object$cv.loss.matrix <- cv.loss.matrix
+ gbm.object$cv.roc.matrix <- cv.roc.stats
+
+ ############################
+ ## JA editions / 08/02/19 ##
+ ############################
+
+ gbm.object$cv.cor.matrix <- cv.cor.stats
+ gbm.object$cv.th.matrix <- threshold.stats
+ gbm.object$cv.corr.class <- cv.corr.class
+ gbm.object$self.th <- self.th
+ gbm.object$cv.length <- cv.length
+ gbm.object$tss.cv <- tss.stat
+
+ ####################
+ ## End / 08/02/19 ##
+ ####################
+
+ if (keep.fold.models) {
+ gbm.object$fold.models <- model.list
+ }
+ else {
+ gbm.object$fold.models <- NULL
+ }
+ if (keep.fold.vector) {
+ gbm.object$fold.vector <- selector
+ }
+ else {
+ gbm.object$fold.vector <- NULL
+ }
+ if (keep.fold.fit) {
+ gbm.object$fold.fit <- fold.fit
+ }
+ else {
+ gbm.object$fold.fit <- NULL
+ }
+ return(gbm.object)
+}
+
+
+## II.
+## https://github.com/cran/ENMeval/tree/master/R
+
+#########################################################
+######### MAKE BLOCK EVALUATION GROUPS #############
+#########################################################
+
+get.block <- function(occ, bg.coords){
+ # SPLIT OCC POINTS INTO FOUR SPATIAL GROUPS
+ noccs <- nrow(occ)
+ n1 <- ceiling(nrow(occ)/2)
+ n2 <- floor(nrow(occ)/2)
+ n3 <- ceiling(n1/2)
+ n4 <- ceiling(n2/2)
+ grpA <- occ[order(occ[, 2]),][1:n1,]
+ grpB <- occ[rev(order(occ[, 2])),][1:n2,]
+ grp1 <- grpA[order(grpA[, 1]),][1:(n3),]
+ grp2 <- grpA[!rownames(grpA)%in%rownames(grp1),]
+ grp3 <- grpB[order(grpB[, 1]),][1:(n4),]
+ grp4 <- grpB[!rownames(grpB)%in%rownames(grp3),]
+
+ # SPLIT BACKGROUND POINTS BASED ON SPATIAL GROUPS
+ bvert <- mean(max(grp1[, 1]), min(grp2[, 1])) + 0.5
+ tvert <- mean(max(grp3[, 1]), min(grp4[, 1])) + 0.5
+ horz <- mean(max(grpA[, 2]), min(grpB[, 2])) + 0.5
+ bggrp1 <- bg.coords[bg.coords[, 2] <= horz & bg.coords[, 1]
+
Anton Van de Putte (Royal Belgian Institute for Natural Sciences, KULeuven) manages the SCAR Antarctic Biodiversity Portal and is a strong advocate of FAIR data and reproducible science. You can ask him questions about biological open data standards and data cleaning.
+
+
Charlène Guillaumot (ULB, Bruxelles) is an ecological modeller.She studies the potential and limits of species distribution models and dynamic energy budget models on Southern Ocean case studies, and develops/adapts methodologies according to the encountered issues.
+
+
Grant Humphries (HiDef Aerial Surveying Ltd) is a seabird biologist and data scientist who has worked on terrestrial and marine ecosystems with a focus on applying machine learning to complex ecosystems. You can chat with Grant about creating decision support tools in R and Python, and what his favourite programming language is.
+
+
Ryan Reisinger (CEBC-CNRS) is an ecologist who models the movement behaviour, habitat use and distribution of marine predators, using telemetry and biologging data.
+
+
Ben Raymond (Australian Antarctic Division and the University of Tasmania) is a data-focused research scientist working across a range of Antarctic and Southern Ocean science. He is a co-proponent of the SCAR/rOpenSci R community initiative (https://scar.github.io/ropensci/)
+
+
Huw Griffiths (British Antarctic Survey) is a marine biologist who studies large-scale biogeographic and ecological patterns in space and time. He uses selected groups of animals to investigate trends in distribution in the Polar Regions, utilising large biological and physical databases to understand key processes.
+
+
+## Course Length
+
+The first part of the course (3 days) will focus on introducing various R packages and overall workflows. The second part of the course (2 days) will be a hands-on opportunity for participants to work on their own data with the assistance of the course instructors. You can choose to only register for the first part of the course if you prefer.
+The course is targeting Polar Early Career scientists, but is open to any other research students or scientists that require these skills for their work. There is space for a maximum of 20 participants.
+
+The course will generally run from 9am to 5:30pm each day, with breaks throughout the day for refreshments and discussion. See the [course schedule](https://scar.github.io/EGABIcourse19/course-schedule.html).
+
+## Venue
+
+The course will take place at the Meeting Room of the Collegium Veteranorum (Cove 02.10)
+
+Sint-Michielstraat 2-4
+
+3000 Leuven
## Requirements
-Some experience in the computer language R is desirable, and experience with GIS is helpful. All participants must bring their own personal laptop.
-If you are not sure if your experience in R is sufficient, you should follow the one of the free online course below:
+
+Some experience in the computer language R is desirable, and experience with GIS is helpful.
+
+**For the course particpipants should bring**
+- their own personal laptop
+- Headphones
+
+**Please make sure you have installed**
+- R https://cran.r-project.org
+- Rstudio https://www.rstudio.com
+
+*Feel free to also install*
+QGIS https://www.qgis.org/en/site/ and Quantarctica http://quantarctica.npolar.no
+
+GGIS is a free, cross-platform, open-source software that allows to create, edit, visualise, analyse and publish geospatial information/ Quantarctica is a collection of Antarctic geographical datasets which works with theQGIS software.
+
+
+If you are not sure if your experience in R is sufficient, you should follow the one of the free online courses below:
https://swcarpentry.github.io/r-novice-gapminder/01-rstudio-intro/
https://obis.org/manual/intror/
-## Language
+## Language
+
English
-## Fees and Funding
-There is a 75 euro registration fee for the course.
-The fee include full participation in the course, and digital copies of the course material, as well as Coffee breaks and a small lunch.
-Participants will have to arrange their own travel and accommodation, but limited travel support will be available in collaboration with Ent-Eco.
+## Code of Conduct
+
+To ensure a safe, enjoyable, and friendly experience for everyone who participates, we have a strict [code of conduct](coc.md) that all participants are expected to follow.
-## Registration
-Application for the course are closed.
-• Deadline for payment is Friday 7th July 2019 where after your place will be offered to someone else.
-If you have questions or sugestion you can create an issue at www.github.com/biodiversity-aq/EGABIcourse19
+## Registration
-Alternatively you can contact datascience.biodiversity.aq( at ) gmail.com
+Applications for the course are closed.
+If you have questions or suggestions you can [create an issue here](https://github.com/SCAR/EGABIcourse19/issues). Alternatively you can contact datascience.biodiversity.aq (at) gmail.com
-## About rOpenSci
+----
-
-
-rOpenSci is a non-profit organisation that maintains a number of widely used R packages and is very active in promoting a community spirit around the R-world. Their core values are to have open and reproducible research, shared data and easy-to-use tools and to make all this accessible to a large number of people.
+## About rOpenSci
+
+8 Biological data
+8.1 Taxonomy
+The Register of Antarctic Marine Species (RAMS) is the authoritative taxonomic database for Antarctic marine organisms. RAMS is part of the World Register of Marine Species (WoRMS).
+-
+
- worrms client for the WoRMS API. Contains mostly taxonomic data, but also trait data. +
- taxize provides access to 20ish sources of taxonomic data sources, including WoRMS. It has consistent data outputs and function interfaces across the different data sources so that you don’t need to tailor your code to different taxonomic data providers. +
RAMS is currently being extended to cover non-marine taxa, which will become the Register of Antarctic Species (RAS). Hopefully this will remain covered by worrms and the server-side infrastructure hosted by VLIZ.
For more detail on R packages dealing with taxonomy in general, see the rOpenSci taxonomy task view.
+8.2 Other packages
+Tracking of animals using satellite, GPS, or light-level geolocation tags is common, and there are many R packages that can help with this. See the spatiotemporal task view for a more complete list. Also of interest may be:
+-
+
TwilightFree provides a method for processing light-level geolocation data that is robust to noise (sensor shading and obscuration) and may be particularly suitable for Southern Ocean applications.
+traipse provides tools for tracking data, for common metrics of distance, direction, and speed.
+foieGras fits a continuous-time model (RW or CRW) in state-space form to filter Argos satellite location data.
+