Getting started with AlphaSDM

Get occurrence records

The GBIF occurrence API is public and needs no account. This helper asks for one year of records inside a box around Tucson, keeping only those with coordinates accurate to 30 m, which suits 10 m embeddings. The API returns at most 300 records per request, so it reads two pages. The models are fitted on the 2022 records; the 2023 records are kept aside to test them.

gbif_records <- function(year) {
  url <- paste0(
    "https://api.gbif.org/v1/occurrence/search?",
    "scientificName=Carnegiea%20gigantea&year=", year,
    "&hasCoordinate=true&hasGeospatialIssue=false",
    "&coordinateUncertaintyInMeters=0,30",
    "&decimalLongitude=-111.4,-110.6&decimalLatitude=31.9,32.6&limit=300")
  do.call(rbind, lapply(c(0, 300), function(offset)
    jsonlite::fromJSON(paste0(url, "&offset=", offset))$results[
      , c("decimalLongitude", "decimalLatitude", "year")]))
}
obs <- gbif_records(2022)
nrow(obs)
#> [1] 438

Format the records

format_data() standardises column names, checks that coordinates are WGS84 longitude and latitude, and drops records outside the years the embeddings cover. With no presence column, every row is a presence.

pres <- format_data(obs, coords = c("decimalLongitude", "decimalLatitude"),
                    year = "year")

Add pseudo-absences

Models need absences too, and where artificial absences go is a modelling decision (Barbet-Massin et al. 2012), so AlphaSDM asks you to choose a strategy. "combined" keeps them away from the presences both geographically and in embedding space, the recipe recommended for the tree models in the default ensemble. aoi = "bbox" draws them inside the bounding box of the presences, and each is read from the same year’s embeddings as the records.

occ <- generate_pseudo_absences(pres, aoi = "bbox", strategy = "combined",
                                n = nrow(pres))
table(occ$present)
#> 
#>   0   1 
#> 282 432

It asked for as many absences as presences but found 282: most of the box looks like saguaro habitat in embedding space, and the function reports a shortfall rather than place absences inside habitat. The exclusion radius (250 m) and envelope threshold it estimated are stored in attr(occ, "pa_settings").

Look at the data

Plot the points on a satellite image before modelling them. This one is a cloud-free Sentinel-2 composite for 2022, made on Earth Engine and downloaded as a small RGB GeoTIFF.

ee <- reticulate::import("ee")
box <- ee$Geometry$Rectangle(c(-111.4, 31.9, -110.6, 32.6))
sentinel2 <- ee$ImageCollection("COPERNICUS/S2_SR_HARMONIZED")$
  filterBounds(box)$
  filterDate("2022-01-01", "2023-01-01")$
  filter(ee$Filter$lt("CLOUDY_PIXEL_PERCENTAGE", 10))$
  median()$
  visualize(bands = list("B4", "B3", "B2"), min = 0, max = 3500)
tif <- tempfile(fileext = ".tif")
utils::download.file(sentinel2$getDownloadURL(list(
  region = box, scale = 100, crs = "EPSG:4326", format = "GEO_TIFF")),
  tif, mode = "wb", quiet = TRUE)

plot(stars::read_stars(tif), rgb = 1:3, reset = FALSE,
     main = "Saguaro records and pseudo-absences")
points(latitude ~ longitude, data = occ[occ$present == 0, ],
       pch = 21, cex = 0.7, bg = "white")
points(latitude ~ longitude, data = occ[occ$present == 1, ],
       pch = 21, cex = 0.8, bg = "gold")
legend("bottomleft", inset = 0.02, bg = "white",
       legend = c("GBIF record", "pseudo-absence"), pch = 21,
       pt.bg = c("gold", "white"))

plot of chunk occurrence-map

The records sit on the desert slopes and foothills around the city, where saguaros grow. The pseudo-absences went where the landscape differs: mostly the forested upper Santa Catalina Mountains to the northeast and the irrigated farmland of the Avra Valley to the west, with a few in town and around the open-pit mine to the south.

Evaluate on the next year’s records

The fairest test uses records the models have never seen. Here that is the 2023 records, scored against background points drawn at random across the study area from the 2023 embeddings. The pseudo-absences are used only for training: scoring against them would reward the rule that placed them, not the models.

pres_2023 <- format_data(gbif_records(2023),
                         coords = c("decimalLongitude", "decimalLatitude"),
                         year = "year")
test <- generate_pseudo_absences(pres_2023, aoi = "bbox", strategy = "random",
                                 n = 2000)

fit <- evaluate_models(data = occ, predict_coords = test)

metrics <- do.call(rbind, lapply(fit$metrics, as.data.frame))
knitr::kable(metrics[, c("auc_roc", "auc_prg", "tss", "cbi")], digits = 3)

The background includes plenty of real saguaro habitat, so AUC cannot reach 1 here even for a perfect model; it measures how well the models separate the 2023 records from the landscape as a whole. The continuous Boyce index (cbi; Hirzel et al. 2006) is the metric designed for presence-only data: it asks whether places the model rates more suitable hold proportionally more of the new records than chance would give, and ranges from -1 to 1. The ensemble scores 0.73.

Adjust a model’s settings

Every model uses Earth Engine’s defaults for its classifier unless you change them (the few exceptions are listed in ?evaluate_models). To change a setting, pass params: a list with one entry per model, using Earth Engine’s argument names.

Earth Engine’s default SVM is a linear classifier. Above, it ranks the 2023 records well but scores poorly on the Boyce index, because its probabilities crowd towards 1. A regression SVM with a radial kernel gives a continuous suitability score instead:

svm_rbf <- list(svmType = "EPSILON_SVR", kernelType = "RBF", cost = 10,
                gamma = 0.05)
fit_svm <- evaluate_models(data = occ, predict_coords = test, methods = "svm",
                           params = list(svm = svm_rbf))

svm <- rbind(`default (linear classifier)` = metrics["svm", ],
             `regression, RBF kernel` = as.data.frame(fit_svm$metrics$svm))
knitr::kable(svm[, c("auc_roc", "auc_prg", "tss", "cbi")], digits = 3)

The Boyce index rises from 0.36 to 0.96. The same params work in generate_map(). In a real study, choose settings using the training records, for example by cross-validation, rather than the test records.

Map suitability

generate_map() refits the models on all the data and writes one GeoTIFF per model plus the ensemble mean. It takes the same params as evaluate_models(). aoi = "bbox" maps the whole area the points cover, about 75 by 78 km. At 30 m that is six map tiles and takes a few minutes; the native 10 m resolution (scale = 10) takes longer.

maps <- generate_map(occ, aoi = "bbox", scale = 30,
                     output_dir = file.path(tempdir(), "saguaro"))
plot(stars::read_stars(maps$ensemble_map), main = "Saguaro suitability, 30 m",
     col = grDevices::hcl.colors(20, "YlGn", rev = TRUE))

plot of chunk map

Suitability follows the slopes and drainages of the Tucson Mountains and the foothills ringing the city, and fades on the valley floors, the farmland and the high Catalinas.

References

Barbet-Massin, M., Jiguet, F., Albert, C. H. & Thuiller, W. (2012). Selecting pseudo-absences for species distribution models: how, where and how many? Methods in Ecology and Evolution 3, 327-338.

Hirzel, A. H., Le Lay, G., Helfer, V., Randin, C. & Guisan, A. (2006). Evaluating the ability of habitat suitability models to predict species presences. Ecological Modelling 199, 142-152.