--- title: "Getting started with AlphaSDM" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Getting started with AlphaSDM} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- AlphaSDM fits species distribution models on Google's AlphaEarth satellite embeddings: 64 numbers per 10 m pixel that summarise what the land surface looks like there, for every year since 2017. Sampling, model fitting and mapping all run on Google Earth Engine, so there are no environmental layers to find, download or align. This vignette maps saguaro cactus (*Carnegiea gigantea*) around Tucson, Arizona, from public GBIF records to a habitat-suitability map, and uses each step of the workflow once. **Before you start**, connect to Earth Engine once per machine with `setup_gee(project = "your-cloud-project")`; the README explains the free registration. After that, AlphaSDM connects on its own. ## 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. ``` r 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. ``` r 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. ``` r 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. ``` r 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](AlphaSDM-occurrence-map-1.png) 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. ``` r 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) ``` | | auc_roc| auc_prg| tss| cbi| |:--------|-------:|-------:|-----:|-----:| |svm | 0.807| 0.786| 0.486| 0.355| |rf | 0.766| 0.705| 0.400| 0.911| |gbt | 0.721| 0.632| 0.336| 0.713| |ensemble | 0.762| 0.703| 0.382| 0.730| 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: ``` r 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) ``` | | auc_roc| auc_prg| tss| cbi| |:---------------------------|-------:|-------:|-----:|-----:| |default (linear classifier) | 0.807| 0.786| 0.486| 0.355| |regression, RBF kernel | 0.792| 0.722| 0.463| 0.959| 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. ``` r 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](AlphaSDM-map-1.png) 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.