# Supplementary: manuscript figures

## Scope

This vignette provides the exact code used to reproduce the figures in Chauvier-Mendes et al. (2025), currently under review, included for transparency and to let readers and reviewers verify how each figure was generated from the functions of this package. All chunks are set to `eval = FALSE`; run them locally to regenerate the figures from GBIF and the bundled example data — the figures shown below are cached, static images.

Note that basemaps here use a low-resolution world map bundled with the package, rather than the `rnaturalearth::ne_countries(scale = 10)` HD dataset used in the manuscript, so coastline detail and ecoregion colors will differ slightly from the published version.

The citation below refers to the current preprint and will be updated once the associated manuscript is published.

> Chauvier-Mendes, Y., Hagen, O., Pinkert, S., Albouy, C., Fopp, F., Brun, P., Descombes, P., Altermatt, F., Pellissier, L., & Csilléry, K. (2025). gbif.range: An R package to generate ecologically-informed species range maps from occurrence data with seamless GBIF integration. Authorea preprint.

    knitr::opts_chunk$set(collapse = TRUE, comment = "#>")

    ext_file <- function(...) {
      path <- system.file("extdata", ..., package = "gbif.range")
      if (nzchar(path)) {
        return(path)
      }
      normalizePath(file.path("..", "inst", "extdata", ...), mustWork = TRUE)
    }

    library(gbif.range)

## Online data retrieval

The figures below share a common set of GBIF downloads and ecoregion layers. To avoid repeating the same network calls in every section, all data used across the figures is retrieved once here. Run this section first if executing chunks individually rather than the whole vignette in order.

    # Study extent, world boundaries, and packaged terrestrial ecoregions
    shp.lonlat <- terra::vect(ext_file("shp_lonlat.shp"))
    countries <- terra::vect(ext_file("world_countries.shp"))
    eco.terra <- read_ecoreg(
      ecoreg_name = "eco_terra", save_dir = tempdir()
    )

    # Custom, high-resolution ecoregions for the Alps study extent
    # ('rst' --> 5x5-km resolution)
    rst <- terra::rast(ext_file("rst.tif"))
    my.eco <- make_ecoreg(env = rst, nclass = 200, format = "SpatVector")

    # Continent extent to keep only terrestrial my.eco
    contExtL.shp <- terra::aggregate(
      terra::crop(countries, terra::ext(rst))
    )
    contExtS.shp <- terra::crop(
      contExtL.shp,
      terra::ext(shp.lonlat)
    )
    my.ecoS <- intersect(my.eco, contExtS.shp)

    # GBIF occurrences for the two focal species used throughout
    obs.arcto <- get_gbif(
      sp_name = "Arctostaphylos alpinus",
      geo = shp.lonlat,
      grain = 1
    )
    obs.pt <- get_gbif(sp_name = "Panthera tigris")

    # Diversity rasters used in Figure 3
    iucn.robin <- terra::rast(ext_file("iucn_div_robin.tif"))
    gf.robin <- terra::rast(ext_file("gf_div_robin.tif"))

## Figure 1a (left) — Custom ecoregion range: Arctostaphylos alpinus

This figure demonstrates a fully custom, high-resolution workflow: ecoregions are derived from CHELSA bioclimatic layers (Karger et al. 2017) via `make_ecoreg()`, and the species’ range is constrained to a small alpine extent.

    # Regional range at ~5x5-km resolution ('my.eco' resolution)
    range.arcto <- get_range(
      occ_coord = obs.arcto,
      ecoreg = my.eco,
      ecoreg_name = "EcoRegion",
      format = "SpatVector"
    )

    # Assign colors to ecoregions
    range.arctoS <- terra::crop(range.arcto$rangeOutput, contExtS.shp)
    col.palette <- grDevices::colorRampPalette(c(
      "#a6cee3", "#1f78b4", "#b2df8a", "#33a02c",
      "#fb9a99", "#e31a1c", "#fdbf6f", "#ff7f00",
      "#cab2d6", "#6a3d9a", "#ffff99", "#b15928"
    ))
    colcol <- col.palette(length(my.ecoS))
    set.seed(7)
    my.ecoS$color <- sample(
      paste0(colcol, ""),
      length(my.ecoS),
      replace = FALSE
    )

    pt.col <- terra::extract(
      x = my.ecoS,
      y = as.data.frame(
        obs.arcto[, c(
          "decimalLongitude", "decimalLatitude"
        )]
      )
    )
    pt.plot <- obs.arcto[
      !is.na(pt.col$color),
      c("decimalLongitude", "decimalLatitude")
    ]
    pt.col2 <- pt.col[!is.na(pt.col$color), "color"]
    pt.col3 <- grDevices::adjustcolor(
      pt.col2, red.f = 0.6, green.f = 0.6,
      blue.f = 0.6
    )

    # Plot
    terra::plot(
      my.ecoS,
      col = paste0(my.ecoS$color, "99"),
      border = NA,
      axes = FALSE
    )
    terra::plot(
      merge_range(range.arctoS),
      border = "black",
      lwd = 1,
      col = "#00000099",
      add = TRUE
    )
    graphics::points(
      pt.plot, col = pt.col2, pch = 16, cex = 1
    )
    graphics::points(
      pt.plot, col = pt.col3, pch = 16, cex = 0.6
    )

## Figure 1a (right) — Packaged ecoregion range: Panthera tigris

This figure demonstrates the complementary broad-scale workflow: a worldwide range built directly from the packaged terrestrial ecoregions (Olson et al. 2001; The Nature Conservancy 2009) via `read_ecoreg()`, requiring no custom ecoregion layer.

    # Global range
    range.tiger <- get_range(
      occ_coord = obs.pt,
      ecoreg = eco.terra,
      ecoreg_name = "ECO_NAME",
      format = "SpatVector"
    )

    # Assign colors to ecoregions
    ext.tiger.eco <- terra::ext(range.tiger$rangeOutput)
    ext.tiger.eco  <- c(
      ext.tiger.eco[1] - 2, ext.tiger.eco[2] + 2,
      ext.tiger.eco[3] - 2, ext.tiger.eco[4] + 2
    )
    eco.local <- terra::crop(eco.terra, ext.tiger.eco)
    col.palette <- grDevices::colorRampPalette(c(
      "#a6cee3", "#1f78b4", "#b2df8a", "#33a02c",
      "#fb9a99", "#e31a1c", "#fdbf6f", "#ff7f00",
      "#cab2d6", "#6a3d9a", "#ffff99", "#b15928"
    ))
    colcol <- col.palette(length(eco.local))
    set.seed(3)
    eco.local$color <- sample(
      paste0(colcol, ""),
      length(eco.local),
      replace = FALSE
    )
    pt.coords <- as.data.frame(
      obs.pt[, c("decimalLongitude", "decimalLatitude")]
    )
    pt.col <- terra::extract(eco.local, pt.coords)
    pt.plot <- obs.pt[
      !is.na(pt.col$color),
      c("decimalLongitude", "decimalLatitude")
    ]
    pt.col2 <- pt.col[!is.na(pt.col$color), "color"]
    pt.col3 <- grDevices::adjustcolor(
      pt.col2, red.f = 0.6, green.f = 0.6,
      blue.f = 0.6
    )
    out.plot <- terra::extract(range.tiger$rangeOutput, pt.coords)
    op.na <- is.na(out.plot[, 2])
    out.plot <- obs.pt[
      op.na, c("decimalLongitude", "decimalLatitude")
    ]

    # Plot
    terra::plot(
      eco.local,
      col = eco.local$color,
      border = NA,
      axes = FALSE
    )
    terra::plot(
      merge_range(range.tiger),
      border = "black",
      lwd = 2,
      col = "#63636399",
      add = TRUE
    )
    graphics::points(
      pt.plot, col = pt.col2, pch = 16, cex = 1.5
    )
    graphics::points(
      pt.plot, col = pt.col3, pch = 16, cex = 0.8
    )
    graphics::points(
      out.plot, col = "black", pch = 4, cex = 1.5,
      lwd = 2
    )

## Figure 2a–b — Cross-validation of ecoregion-constrained ranges

This figure evaluates both ranges from Figure 1 using `cv_range()`’s block cross-validation, reporting mean sensitivity and precision against held-out occurrence folds for the custom-ecoregion (A. alpinus) and packaged-ecoregion (P. tigris) workflows.

    # ------------------------------------------
    # Arctostaphylos alpinus (custom ecoregion)
    # ------------------------------------------

    # Create pseudo-absences
      # First remove observations considered outliers in get_range()
    xy.df <- range.arcto$init.args$occ_coord
    r.ext <- terra::ext(range.arcto$rangeOutput)
    Xrm.cond <- xy.df$decimalLongitude >= r.ext[1] &
      xy.df$decimalLongitude <= r.ext[2]
    Yrm.cond <- xy.df$decimalLatitude >= r.ext[3] &
      xy.df$decimalLatitude <= r.ext[4]
    xy.df <- xy.df[Xrm.cond & Yrm.cond, ]
      # Sample n regular background points over the range extent
    x.interv <- (r.ext[2] - r.ext[1]) / (sqrt(1e4) - 1)
    y.interv <- (r.ext[4] - r.ext[3]) / (sqrt(1e4) - 1)
    lx <- seq(r.ext[1], r.ext[2], x.interv)
    ly <- seq(r.ext[3], r.ext[4], y.interv)
    bp.xy <- expand.grid(
      decimalLongitude = lx,
      decimalLatitude = ly
    )

    # Combine observations with background
    obs.xy <- xy.df[, c("decimalLongitude", "decimalLatitude")]
    all.xy <- rbind(obs.xy, bp.xy)
    all.xy$Pres <- 0
    all.xy[1:nrow(obs.xy), "Pres"] <- 1

    # Run block-cv
    xy.pres <- all.xy$Pres
    cv.strat <- make_blocks(
      nfolds = 5,
      df = all.xy[, c("decimalLongitude", "decimalLatitude")],
      nblocks = 5 * 2,
      pres = xy.pres
    )
    all.xy$bcv <- cv.strat
    all.xy[all.xy$bcv %in% 1, "col"] <- "#e41a1c"
    all.xy[all.xy$bcv %in% 2, "col"] <- "#377eb8"
    all.xy[all.xy$bcv %in% 3, "col"] <- "#4daf4a"
    all.xy[all.xy$bcv %in% 4, "col"] <- "#984ea3"
    all.xy[all.xy$bcv %in% 5, "col"] <- "#ff7f00"
    all.xy[all.xy$Pres %in% 1, "col"] <- substr(
      grDevices::adjustcolor(
        all.xy[all.xy$Pres %in% 1, "col"],
        red.f = 0.5, green.f = 0.5, blue.f = 0.5
      ),
      1, 7
    )

    # Evaluate
    ar.test <- cv_range(
      range_object = range.arcto,
      cv = "block-cv",
      nfolds = 5,
      nblocks = 2
    )

    # Use the broader (L) extent and country borders so
    # the evaluation text fits on the plot
    ext.L <- terra::ext(contExtL.shp)
    world.local.ar <- terra::crop(countries, ext.L)
    world.local.ar <- terra::aggregate(world.local.ar)
    pres <- all.xy[all.xy$Pres %in% 1, ]
    abs.pts <- all.xy[all.xy$Pres %in% 0, ]
    pres_coords <- as.data.frame(
      pres[, c("decimalLongitude", "decimalLatitude")]
    )
    id.in <- terra::extract(
      range.arcto$rangeOutput, pres_coords
    )
    pres <- pres[!is.na(id.in[, 2]), ]

    # Plot
    terra::plot(
      world.local.ar, col = "#bcd1bc",
      axes = FALSE, lwd = 1
    )
    graphics::points(
      abs.pts, col = paste0(abs.pts$col, "50"),
      pch = 16, cex = 0.5
    )
    terra::plot(
      merge_range(range.arctoS),
      border = "black", lwd = 1.7,
      col = "#63636370", add = TRUE
    )
    graphics::points(
      pres, col = paste0(pres$col, "90"),
      pch = 16, cex = 1.3
    )

    # Evaluation text anchored to the L extent's bottom-left corner
    txt.x <- ext.L[1] + 0.02 * (ext.L[2] - ext.L[1])
    txt.y1 <- ext.L[3] + 0.09 * (ext.L[4] - ext.L[3])
    txt.y2 <- ext.L[3] + 0.04 * (ext.L[4] - ext.L[3])

    graphics::text(
      txt.x, txt.y1,
      paste(
        "Mean TSS =",
        round(tail(ar.test[, "TSS"], 1), 2)
      ),
      cex = 1.2, font = 2, adj = 0
    )
    graphics::text(
      txt.x, txt.y2,
      paste(
        "Mean Precision =",
        round(tail(ar.test[, "Precision"], 1), 2)
      ),
      cex = 1.2, font = 2, adj = 0
    )

The tiger panel follows the same block cross-validation, applied instead to the packaged-ecoregion, worldwide range.

    # ------------------------------------------
    # Panthera tigris (packaged ecoregion)
    # ------------------------------------------

    ext.temp.t2 <- terra::ext(range.tiger$rangeOutput)
    ext.temp.t2  <- c(
      ext.temp.t2[1] - 0.2, ext.temp.t2[2] + 0.2,
      ext.temp.t2[3] - 0.2, ext.temp.t2[4] + 0.2
    )

    # Create pseudo-absences
      # First remove observations considered outliers in get_range()
    xy.df <- range.tiger$init.args$occ_coord
    r.ext <- terra::ext(range.tiger$rangeOutput)
    Xrm.cond <- xy.df$decimalLongitude >= r.ext[1] &
      xy.df$decimalLongitude <= r.ext[2]
    Yrm.cond <- xy.df$decimalLatitude >= r.ext[3] &
      xy.df$decimalLatitude <= r.ext[4]
    xy.df <- xy.df[Xrm.cond & Yrm.cond, ]
      # Sample n regular background points over the range extent
    x.interv <- (r.ext[2] - r.ext[1]) / (sqrt(1e4) - 1)
    y.interv <- (r.ext[4] - r.ext[3]) / (sqrt(1e4) - 1)
    lx <- seq(r.ext[1], r.ext[2], x.interv)
    ly <- seq(r.ext[3], r.ext[4], y.interv)
    bp.xy <- expand.grid(
      decimalLongitude = lx,
      decimalLatitude = ly
    )

    # Combine observations with background
    obs.xy <- xy.df[, c("decimalLongitude", "decimalLatitude")]
    all.xy <- rbind(obs.xy, bp.xy)
    all.xy$Pres <- 0
    all.xy[1:nrow(obs.xy), "Pres"] <- 1

    # Run block-cv
    xy.pres <- all.xy$Pres
    cv.strat <- make_blocks(
      nfolds = 5,
      df = all.xy[, c("decimalLongitude", "decimalLatitude")],
      nblocks = 5 * 2,
      pres = xy.pres
    )
    all.xy$bcv <- cv.strat
    all.xy[all.xy$bcv %in% 1, "col"] <- "#e41a1c"
    all.xy[all.xy$bcv %in% 2, "col"] <- "#377eb8"
    all.xy[all.xy$bcv %in% 3, "col"] <- "#4daf4a"
    all.xy[all.xy$bcv %in% 4, "col"] <- "#984ea3"
    all.xy[all.xy$bcv %in% 5, "col"] <- "#ff7f00"
    all.xy[all.xy$Pres %in% 1, "col"] <- substr(
      grDevices::adjustcolor(
        all.xy[all.xy$Pres %in% 1, "col"],
        red.f = 0.5, green.f = 0.5, blue.f = 0.5
      ),
      1, 7
    )

    # Evaluate
    pt.test <- cv_range(
      range_object = range.tiger,
      cv = "block-cv",
      nfolds = 5,
      nblocks = 2
    )

    # Use extent and country borders
    world.local.ti <- terra::crop(countries, ext.temp.t2)
    world.local.ti <- terra::aggregate(world.local.ti)
    pres <- all.xy[all.xy$Pres %in% 1, ]
    abs.pts <- all.xy[all.xy$Pres %in% 0, ]
    pres_coords <- as.data.frame(
      pres[, c("decimalLongitude", "decimalLatitude")]
    )
    id.in <- terra::extract(range.tiger$rangeOutput, pres_coords)
    pres <- pres[!is.na(id.in[, 2]), ]

    # Plot
    terra::plot(
      world.local.ti, col = "#bcd1bc",
      axes = FALSE, lwd = 1
    )
    graphics::points(
      abs.pts, col = paste0(abs.pts$col, "50"),
      pch = 16, cex = 0.6
    )
    terra::plot(
      merge_range(range.tiger),
      border = "black", lwd = 2,
      col = "#63636370", add = TRUE
    )
    graphics::points(
      pres, col = paste0(pres$col, "80"),
      pch = 16, cex = 1.6
    )

    # Evaluation text anchored to the tiger extent's top-left corner
    txt.x.t <- ext.temp.t2[1] + 0.02 * (ext.temp.t2[2] - ext.temp.t2[1])
    txt.y1.t <- ext.temp.t2[4] - 0.05 * (ext.temp.t2[4] - ext.temp.t2[3])
    txt.y2.t <- ext.temp.t2[4] - 0.10 * (ext.temp.t2[4] - ext.temp.t2[3])

    graphics::text(
      txt.x.t, txt.y1.t,
      paste(
        "Mean TSS =",
        round(tail(pt.test[, "TSS"], 1), 2)
      ),
      cex = 1.5, font = 2, adj = 0
    )
    graphics::text(
      txt.x.t, txt.y2.t,
      paste(
        "Mean Precision =",
        round(tail(pt.test[, "Precision"], 1), 2)
      ),
      cex = 1.5, font = 2, adj = 0
    )

## Figure 2c–d — Evaluation against external validation data

This figure evaluates the ranges from Figure 1a (`range.arcto`, `range.tiger`) against independent species distribution model outputs bundled with the package, using `evaluate_range()` to report precision and sensitivity relative to external validation data rather than held-out occurrences.

    # -------------------------------------
    # Arctostaphylos alpinus (custom ecoregion)
    # -------------------------------------

    root.dir <- ext_file()

    res5km <- evaluate_range(
      root_dir = root.dir,
      valData_dir = "SDM",
      ecoRM_dir = "EcoRM",
      verbose = TRUE,
      print_map = FALSE,
      valData_type = "TIFF",
      mask = NULL,
      res_fact = NULL
    )

    # Plot plant
    terra::plot(
      contExtL.shp, col = "#dce8dc",
      axes = FALSE, lwd = 1
    )
    terra::plot(
      merge_range(range.arctoS),
      lwd = 0.1, col = "#63636350", add = TRUE
    )
    terra::plot(
      res5km$overlay_list[[1]],
      col = c(
        "#d1c845", "#e39c59", "#8d60cc", "#6cba6c"
      ),
      breaks = c(-0.5, 0.5, 1.5, 2.5, 3.5),
      axes = FALSE, legend = FALSE,
      las = 1, add = TRUE
    )

    # Evaluation text anchored to the L extent's bottom-left corner
    ext.L <- terra::ext(contExtL.shp)
    txt.x <- ext.L[1] + 0.02 * (ext.L[2] - ext.L[1])
    txt.y1 <- ext.L[3] + 0.09 * (ext.L[4] - ext.L[3])
    txt.y2 <- ext.L[3] + 0.04 * (ext.L[4] - ext.L[3])

    graphics::text(
      txt.x, txt.y1,
      paste(
        "Mean TSS =",
        round(res5km$df_eval[1, "TSS_ecoRM"], 2)
      ),
      cex = 1, font = 2, adj = 0
    )
    graphics::text(
      txt.x, txt.y2,
      paste(
        "Mean Precision =",
        round(res5km$df_eval[1, "Prec_ecoRM"], 2)
      ),
      cex = 1, font = 2, adj = 0
    )
    terra::plot(
      contExtL.shp,
      border = "#383d38", axes = FALSE,
      lwd = 1, add = TRUE
    )

The tiger panel follows the same external-validation evaluation, applied instead to the packaged-ecoregion, worldwide range.

    # -------------------------------------
    # Panthera tigris (packaged ecoregion)
    # -------------------------------------

    # Plot tiger
    terra::plot(
      world.local.ti, col = "#dce8dc",
      axes = FALSE, lwd = 1
    )
    toPlot <- terra::mask(
      res5km$overlay_list[[6]], world.local.ti
    )
    terra::plot(
      toPlot,
      col = c(
        "#d1c845", "#e39c59", "#8d60cc", "#6cba6c"
      ),
      breaks = c(-0.5, 0.5, 1.5, 2.5, 3.5),
      axes = FALSE, legend = FALSE,
      las = 1, add = TRUE
    )

    # Evaluation text anchored to the tiger extent's top-left corner
    txt.x.t <- ext.temp.t2[1] + 0.02 * (ext.temp.t2[2] - ext.temp.t2[1])
    txt.y1.t <- ext.temp.t2[4] - 0.05 * (ext.temp.t2[4] - ext.temp.t2[3])
    txt.y2.t <- ext.temp.t2[4] - 0.10 * (ext.temp.t2[4] - ext.temp.t2[3])

    graphics::text(
      txt.x.t, txt.y1.t,
      paste(
        "Mean TSS =",
        round(res5km$df_eval[6, "TSS_ecoRM"], 2)
      ),
      cex = 1.5, font = 2, adj = 0
    )
    graphics::text(
      txt.x.t, txt.y2.t,
      paste(
        "Mean Precision =",
        round(res5km$df_eval[6, "Prec_ecoRM"], 2)
      ),
      cex = 1.5, font = 2, adj = 0
    )
    terra::plot(
      world.local.ti,
      border = "#383d38", axes = FALSE,
      lwd = 1, add = TRUE
    )
    leg.x <- ext.temp.t2[1] + 0.05 * (ext.temp.t2[2] - ext.temp.t2[1])
    leg.y <- ext.temp.t2[3] + 0.20 * (ext.temp.t2[4] - ext.temp.t2[3])
    graphics::legend(
      leg.x, leg.y,
      legend = c(
        "True Presences", "True Absences",
        "False Presences", "False Absences"
      ),
      fill = c(
        "#6cba6c", "#d1c845", "#e39c59", "#8d60cc"
      ),
      bg = NA, box.col = NA,
      cex = 1.1,
      x.intersp = 0.2
    )

## Figure 3 — Global richness comparison and correlation

This figure compares species richness estimated from `gbif.range`-derived ranges against IUCN expert range maps (IUCN 2024), at the global scale. `iucn.robin` and `gf.robin` were loaded in “Online data retrieval” above.

    # -------------------------------------
    # World maps: IUCN vs gbif.range richness
    # -------------------------------------

    # CRS
    robin <- paste(
      "+proj=robin +lon_0=0 +x_0=0 +y_0=0",
      "+datum=WGS84 +units=m +no_defs +type=crs"
    )

    # Reproject countries
    countries.robin <- terra::project(countries, robin)

    # Boundary box
    bb <- terra::as.polygons(
      terra::ext(-180, 180, -90, 90),
      crs = "EPSG:4326"
    )
    bb <- terra::densify(bb, interval = 1)
    bb.robin <- terra::project(bb, robin)

    # Shared diverging color palette for both maps
    col <- grDevices::colorRampPalette(
      c("#67a9cf", "#f7f7f7", "#ef8a62")
    )

    # Plot IUCN map
    max.iucn <- round(terra::minmax(iucn.robin)[2])
    terra::plot(
      iucn.robin,
      axes = FALSE, legend = FALSE,
      col = col(10), smooth = TRUE,
      mar = c(1, 1, 1, 5)
    )
    terra::plot(countries.robin, add = TRUE, lwd = 2)
    terra::plot(bb.robin, add = TRUE, lwd = 3)
    oldpar <- graphics::par(xpd = NA, lwd = 3)
    cscl(
      colors = col(10),
      crds = c(-9501111, 9734033, -11800000, -13000000),
      zrng = c(0, max.iucn),
      tickle = -0.3, cx = 1.1, lablag = -1.3,
      tria = "b", at = seq(0, max.iucn, 10),
      horiz = TRUE, title = "IUCN richness",
      labs = seq(0, max.iucn, 10), titlag = 2
    )
    graphics::par(oldpar)

    # Plot GBIF.RANGE map
    max.gf <- round(terra::minmax(gf.robin)[2])
    terra::plot(
      gf.robin,
      axes = FALSE, legend = FALSE,
      col = col(10), smooth = TRUE,
      mar = c(1, 1, 1, 5)
    )
    terra::plot(countries.robin, add = TRUE, lwd = 2)
    terra::plot(bb.robin, add = TRUE, lwd = 3)
    oldpar <- graphics::par(xpd = NA, lwd = 3)
    cscl(
      colors = col(10),
      crds = c(-9501111, 9734033, -11800000, -13000000),
      zrng = c(0, max.gf),
      tickle = -0.3, cx = 1.1, lablag = -1.3,
      tria = "b", at = seq(0, max.gf, 10),
      horiz = TRUE, title = "gbif.range richness",
      labs = seq(0, max.gf, 10), titlag = 2
    )
    graphics::par(oldpar)

The scatter plots below summarize richness and range-area agreement between the two products, reporting Lin’s concordance correlation coefficient and Pearson’s r for each.

    # -------------------------------------
    # Scatter plots: richness and area agreement
    # -------------------------------------

    # Load area table
    data(area_data)

    # Extract data to plot
    cor.ras <- terra::rast(list(gf.robin, iucn.robin))
    names(cor.ras) <- c("RANGE", "IUCN")
    set.seed(1)
    samp.div <- terra::spatSample(
      cor.ras, 5000, replace = FALSE, na.rm = TRUE
    )
    dat.plot <- list(samp.div, area_data[, -1])
    strings <- c("richness", "areas")

    # Plot diversities relationship, side by side
    oldpar <- graphics::par(mfrow = c(1, 2), mar = c(5, 5.5, 4, 2))
    lapply(seq_along(dat.plot), function(x) {

      # Extract the data
      xy <- dat.plot[[x]]
      sex <- 2.5
      col <- "#d95f0240"
      add <- ""

      # Log or not depending on the output
      if (strings[x] == "areas") {
        xy <- log(xy + 1)
        sex <- 2.5
        col <- "#d95f0240"
        add <- "(log)"
      }

      # Plot points
      graphics::plot(
        xy[, 2], xy[, 1],
        cex.axis = 1.7, cex = sex, col = col,
        xlim = c(min(xy[, 2]), max(xy[, 2])),
        ylim = c(min(xy[, 1]), max(xy[, 1])),
        pch = 16,
        xlab = paste("IUCN", strings[x], add),
        ylab = paste("gbif.range", strings[x], add),
        cex.lab = 1.7, font.lab = 1
      )

      # Run Linear Models and spearman's correlation
      lm.div <- stats::lm(xy[, 1] ~ xy[, 2])

      # Lin's concordance correlation coefficient
      # (population covariance/variance)
      ccc_lin <- function(x, y) {
        mx <- mean(x)
        my <- mean(y)
        vx <- mean((x - mx)^2)
        vy <- mean((y - my)^2)
        sxy <- mean((x - mx) * (y - my))
        2 * sxy / (vx + vy + (mx - my)^2)
      }
      ccc.div <- ccc_lin(xy[, 2], xy[, 1])
      cor.div <- stats::cor(xy[, 2], xy[, 1])
      adj.r2 <- summary(lm.div)[[9]]

      # Plot text
      text_cor1 <- bquote("ccc" == .(round(ccc.div, 2)))
      text_cor2 <- bquote("r" == .(round(cor.div, 2)))
      fig_label(
        text_cor1,
        region = "plot", pos = "topleft",
        bty = "n", font = 2, col = "#121212",
        cex = 2, margin = 0.02
      )
      fig_label(
        text_cor2,
        region = "plot", pos = "bottomright",
        bty = "n", font = 2, col = "#6f69c2",
        cex = 2, margin = 0.02
      )

      # Plot relationship
      graphics::lines(
        xy[, 2], lm.div$fit, lwd = 7, col = "#7570b3"
      )
      graphics::abline(
        a = 0, b = 1, col = "#252525",
        lwd = 5, lty = 2
      )
    })
    graphics::par(oldpar)

## References

Chamberlain, S., Oldoni, D., & Waller, J. (2022). rgbif: interface to the global biodiversity information facility API. https://doi.org/10.5281/zenodo.6023735

Chauvier-Mendes, Y., Hagen, O., Pinkert, S., Albouy, C., Fopp, F., Brun, P., Descombes, P., Altermatt, F., Pellissier, L., & Csilléry, K. (2025). gbif.range: An R package to generate ecologically-informed species range maps from occurrence data with seamless GBIF integration. Authorea preprint. https://doi.org/10.22541/au.175130858.83083354/v1

IUCN. (2024). The IUCN Red List of Threatened Species. Version 2024-2. https://www.iucnredlist.org. Accessed on 2024-09-25.

Karger, D. N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R. W., Zimmermann, N. E., Linder, H. P., & Kessler, M. (2017). Climatologies at high resolution for the earth’s land surface areas. Scientific Data, 4, 170122. https://doi.org/10.1038/sdata.2017.122

Olson, D. M., Dinerstein, E., Wikramanayake, E. D., Burgess, N. D., Powell, G. V. N., Underwood, E. C., … Kassem, K. R. (2001). Terrestrial ecoregions of the world: a new map of life on Earth. BioScience, 51(11), 933–938. https://doi.org/10.1641/0006-3568(2001)051[0933:TEOTWA]2.0.CO;2

The Nature Conservancy (2009). Global Ecoregions, Major Habitat Types, Biogeographical Realms and The Nature Conservancy Terrestrial Assessment Units. Cambridge (UK): The Nature Conservancy. https://geospatial.tnc.org/datasets/b1636d640ede4d6ca8f5e369f2dc368b/about