---
title: "Part 2: Ecoregion-Based Range Inference"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Part 2: Ecoregion-Based Range Inference}
  %\VignetteEngine{knitr::rmarkdown_notangle}
  %\VignetteEncoding{UTF-8}
---

## Transparent setup

```{r setup, include = TRUE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  purl = FALSE
)

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)
```

## Scope

This vignette focuses on the core ecological idea behind `gbif.range`: species ranges are inferred from occurrences, but those occurrences are constrained by ecoregions rather than expanded across the landscape as a purely geometric hull.

The two longer examples cover a broad terrestrial workflow using `Panthera tigris` and a finer regional workflow using `Arctostaphylos alpinus` in the European Alps with custom ecoregions derived from environmental rasters.

The relevant functions are:

- `get_range()` for the actual range inference,
- `read_ecoreg()` and `ecoreg_list` for packaged large-scale ecoregions,
- `make_ecoreg()` for custom ecoregions derived from any set of spatial raster layers,
- `cv_range()` and `evaluate_range()` for evaluating the resulting maps.

## How `get_range()` works

`get_range()` is designed for one focal species at a time. The algorithm
originates from Hagen et al. (2019) and can be summarized as four linked
operations:

1. remove isolated spatial outliers from the occurrence table,
2. identify occupied ecoregions and cluster occurrences within them,
3. build convex hulls and buffers for each within-ecoregion cluster,
4. intersect those polygons with the occupied ecoregions to enforce ecological limits.

This means that the final output reflects both the geometry of the occurrences and the eco-geographic structure of the study area.

## A minimal offline example

The example below is deliberately small and fully offline. Two synthetic rasters are used to create a custom ecoregion layer, then a short occurrence table is passed to `get_range()`.

```{r minimal-example}
set.seed(1)

# Create two simple environmental surfaces on a small study area.
r1 <- terra::rast(ncols = 20, nrows = 20, xmin = 0, xmax = 10, ymin = 0, ymax = 10)
terra::values(r1) <- rep(seq(0, 1, length.out = 20), each = 20)
r2 <- terra::rast(r1)
terra::values(r2) <- rep(seq(0, 1, length.out = 20), times = 20)
env <- c(r1, r2)

# Derive four custom ecoregions from the environmental layers.
eco <- make_ecoreg(env = env, nclass = 4)

# Define one species as an occurrence table in the format expected by get_range().
occ <- data.frame(
  input_search = "example_species",
  decimalLongitude = c(0.5, 1.2, 2.4, 3.6, 2.8, 6.0, 7.2, 7.4, 7.6, 8.8),
  decimalLatitude = c(1.0, 0.1, 2.3, 2.5, 2.7, 5.0, 7.1, 7.3, 6.5, 7.7)
)

# Build the range from the occurrences and the ecoregion layer.
range_obj <- get_range(
  occ_coord = occ,
  ecoreg = eco,
  ecoreg_name = "EcoRegion",
  verbose = FALSE
)

# Plot
terra::plot(
  range_obj$rangeOutput,
  col = "#3c8d5a",
  main = "A minimal ecoregion-constrained range"
)
graphics::points(occ$decimalLongitude, occ$decimalLatitude, pch = 4)
```

Three practical points are worth noting here.

First, `get_range()` accepts any set of occurrence coordinates as input via
`occ_coord` — not just a `getGBIF` object returned by `get_gbif()`. Any
`data.frame` with `decimalLongitude`, `decimalLatitude` columns are valid.
This means museum records, citizen science exports, expert-curated datasets, or
coordinates from any other source can be used directly, as long as they follow
that column structure.

Second, custom ecoregions created with `make_ecoreg()` can be passed directly into `get_range()`. This is useful for regional studies where global ecoregions are too coarse.

Third, `get_range()` returns a `getRange` object. The spatial output is stored in `range_obj$rangeOutput`, while the original inputs are kept in `range_obj$init.args`.

## Packaged versus custom ecoregions

Large-scale analyses usually rely on one of the packaged ecoregion layers:
terrestrial (*eco_terra*; Olson et al. 2001; The Nature Conservancy 2009),
marine (*eco_marine*, *eco_hd_marine*; Spalding et al. 2007, 2012;
The Nature Conservancy 2012), and freshwater (*eco_fresh*; Abell et al. 2008).

```{r ecoregions}
vapply(ecoreg_list, `[[`, character(1), "filename")
```

The public pattern for a terrestrial analysis is:

```{r built-in-ecoregions, eval = FALSE}
# Read ecoregion
eco_terra <- read_ecoreg("eco_terra", save_dir = tempdir())

# Construct range (not run)
range <- get_range(
  occ_coord = obs,
  ecoreg = eco_terra,
  ecoreg_name = "ECO_NAME"
)
```

Beyond the packaged layers, any polygon object can serve as the `ecoreg`
argument — a `SpatVector`, `sf`, or other spatial polygon class — provided it
contains a character column whose name is passed to `ecoreg_name`. Habitat
maps, expert-defined biogeographic units, and bioregions derived from species
composition data (Denelle et al. 2025) therefore need no conversion step.

When no suitable polygon layer exists, `make_ecoreg()` builds one by
clustering raster data (Chauvier et al. 2021). Climate layers are the most
common input, but any spatially structured raster works — species richness
surfaces, biodiversity indices, habitat suitability maps, or any combination
of drivers thought to define meaningful spatial units. This is typically the
route for regional analyses, where the packaged global layers are too coarse.
The package ships raster examples under `inst/extdata`, including `rst.tif`:

```{r custom-ecoregions, eval = FALSE}
# Construct ecoregion
rst <- terra::rast(ext_file("rst.tif"))
my_eco <- make_ecoreg(env = rst, nclass = 200)

# Construct range (not run)
range <- get_range(
  occ_coord = obs,
  ecoreg = my_eco,
  ecoreg_name = "EcoRegion",
  res = 0.05
)
```

The key trade-off is ecological detail versus robustness to sampling density. Coarser ecoregions produce broader, smoother ranges. Finer ecoregions can better capture regional structure, but they also make the results more sensitive to local sampling gaps.

## A broad terrestrial workflow: *Panthera tigris*

Tiger occurrences are a good example of the broad terrestrial use case for `gbif.range`. The workflow is straightforward: retrieve occurrences with a coarse precision filter, load the terrestrial ecoregions of the world, and infer the range inside those ecoregional boundaries.

```{r tiger-workflow, eval = FALSE}
# Step 1: download global tiger occurrences with a 100 km precision filter
obs_tiger <- get_gbif(
  sp_name = "Panthera tigris",
  grain = 100
)

# Step 2: load the packaged terrestrial ecoregion layer
eco_terra <- read_ecoreg("eco_terra", save_dir = tempdir())

# Step 3: infer the range with default parameters
range_tiger <- get_range(
  occ_coord = obs_tiger,
  ecoreg = eco_terra,
  ecoreg_name = "ECO_NAME",
  degrees_outlier = 5,
  clust_pts_outlier = 4,
  res = 0.1
)

# Plot
terra::plot(range_tiger$rangeOutput, col = "#238b45")
graphics::points(obs_tiger[, c("decimalLongitude", "decimalLatitude")], pch = 20)
```

```{r fig-tiger-workflow, echo = FALSE, out.width = "50%"}
knitr::include_graphics("../man/figures/Part1_tiger_workflow.png")
```

This is the typical continental to global pattern. The ecoregion layer is externally defined and ecologically interpretable, while the occurrence filter is intentionally conservative. In practice, this is the kind of workflow where `gbif.range` is most useful as a transparent alternative to unconstrained geometric envelopes.

## A finer regional workflow: *Arctostaphylos alpinus* in the Alps

The second workflow is deliberately different. Here the issue is not worldwide coverage, but the opposite: the global terrestrial ecoregions are too coarse for a mountain system with strong local climatic structure. A custom ecoregion layer is therefore built from packaged example CHELSA bioclimatic rasters (Karger et al. 2017).

```{r arcto-workflow, eval = FALSE}
# Step 1: define the Alps study region and retrieve occurrences.
alps_extent <- terra::vect(ext_file("shp_lonlat.shp"))
obs_arcto <- get_gbif(
  sp_name = "Arctostaphylos alpinus",
  geo = alps_extent,
  grain = 1
)

# Step 2: create finer ecoregions from the packaged environmental rasters.
rst <- terra::rast(ext_file("rst.tif"))
eco_alps <- make_ecoreg(env = rst, nclass = 200)

# Step 3: infer the range at finer output resolution (~5x5-km resolution)
range_arcto <- get_range(
  occ_coord = obs_arcto,
  ecoreg = eco_alps,
  ecoreg_name = "EcoRegion",
  res = 0.05,
  format = "SpatRaster"
)

# Plot
countries <- terra::vect(ext_file("world_countries.shp"))
terra::plot(alps_extent, col="#00000030", border = NULL)
plot(countries, col = "#cacee8", add = TRUE)
terra::plot(alps_extent, col="#00000020", border = NULL, add = TRUE)
terra::plot(range_arcto$rangeOutput, col = "#3c8d5a", add = TRUE)
graphics::points(obs_arcto[, c("decimalLongitude", "decimalLatitude")], pch = 20)
```

```{r fig-arcto-workflow, echo = FALSE, out.width = "50%"}
knitr::include_graphics("../man/figures/Part1_arcto_workflow.png")
```

This regional example shows why `make_ecoreg()` is part of the main package workflow rather than an auxiliary convenience. At small spatial extents, the ecological realism of the range map often depends more on the choice of ecoregion layer than on minor changes in buffer parameters.

## Tuning the main range arguments

The most important `get_range()` arguments fall into three groups:

| Group | Arguments | Role |
|---|---|---|
| Outlier removal | `degrees_outlier`, `clust_pts_outlier` | How aggressively isolated occurrences are discarded. Matters most when the input contains obvious anomalies or disjunct stray clusters of observations. |
| Buffering | `buff_width_point`, `buff_incrmt_pts_line`, `buff_width_polygon` | How singletons, linear clusters, and polygon hulls are buffered before the final ecoregion intersection. |
| Output | `format`, `res` | Geometry type and, for rasters, resolution. Vector output suits exploratory work; gridded output is easier to stack across many species. |

Which group needs attention depends on scale.

At broad scale the defaults are usually a reasonable starting point, since the
main goal is simply to remove very isolated records and obvious anomalies.

At regional scale the decisive settings often lie outside `get_range()`
altogether: the coordinate precision filter (`grain`) in `get_gbif()` and the
granularity (`nclass`) of the `make_ecoreg()` layer. Species with narrow,
habitat-constrained distributions are the exception — there, tightening the
outlier and buffer parameters themselves can be necessary, as shown below.

*Paederota bonarota* is a good example: restricted to limestone/dolomite rock
crevices, its Alpine distribution is essentially limited to the Slovenian,
Italian, and Austrian sectors regardless of broader climatic suitability.
Ideally this habitat constraint would be encoded directly in the custom
ecoregions — e.g. by adding a bedrock or geology layer to `make_ecoreg()`.
Since that is not done here (`eco_alps` uses only climate), the outlier and buffer
parameters are tightened instead, as a geographic substitute for the missing
environmental constraint:

```{r paederota-workflow, eval = FALSE}
# Download Paederota bonarota within the same Alps study region
obs_paed <- get_gbif(
  sp_name = "Paederota bonarota",
  geo = alps_extent,
  grain = 1
)

# Build the range with tightened outlier and buffer parameters, reflecting
# the species' narrow, rock-habitat-restricted distribution
range_paed <- get_range(
  occ_coord = obs_paed,
  ecoreg = eco_alps,
  ecoreg_name = "EcoRegion",
  res = 0.05,
  degrees_outlier = 0.5,
  buff_width_point = 0.5,
  buff_incrmt_pts_line = 0.5,
  buff_width_polygon = 0.5,
  format = "sf"
)

# Plot
countries <- terra::vect(ext_file("world_countries.shp"))
terra::plot(alps_extent, col="#00000030", border = NULL)
plot(countries, col = "#cacee8", add = TRUE)
terra::plot(alps_extent, col="#00000020", border = NULL, add = TRUE)
terra::plot(merge_range(range_paed), col = "#3c8d5a", add = TRUE)
graphics::points(obs_paed[, c("decimalLongitude", "decimalLatitude")], pch = 20)
```

```{r fig-paederota-workflow, echo = FALSE, out.width = "50%"}
knitr::include_graphics("../man/figures/Part1_paederota_workflow.png")
```

## Evaluation workflows

The package provides two complementary ways to evaluate range maps.

`cv_range()` performs internal cross-validation of a `getRange` object by repeatedly rebuilding the map from subsets of the original occurrences:

```{r cv-code, eval = FALSE}
cv_res <- cv_range(
  range_object = range_obj,
  cv = "random-cv",
  nfolds = 5
)
```

`evaluate_range()` compares saved range outputs with external validation layers. The package includes a small validation example under `inst/extdata`:

```{r evaluate-code, eval = FALSE}
# Set root
root_dir <- ext_file()

# Evaluate
res_eval <- evaluate_range(
  root_dir = root_dir,
  valData_dir = "SDM",
  ecoRM_dir = "EcoRM",
  print_map = FALSE,
  verbose = FALSE,
  valData_type = "TIFF"
)

# Summary
head(res_eval$df_eval)
```

Together, `cv_range()` and `evaluate_range()` let you tune parameter choices on well-documented taxa before applying the same workflow to data-poor species.

The same logic can be used with the broad and regional examples above. For example, a block cross-validation of the tiger workflow is:

```{r tiger-cv, eval = FALSE}
cv_range(
  range_object = range_tiger,
  cv = "block-cv",
  nfolds = 5,
  nblocks = 2,
  backpoints = 1e4
)
```

That kind of evaluation is particularly useful when deciding how strongly to filter outliers, how fine the ecoregion layer should be, or whether a raster output resolution is appropriate for later biodiversity summaries.

## Take-home message

`get_range()` is best understood as an ecological range-construction tool rather than a pure hull algorithm. The strength of the approach lies in combining occurrence geometry with eco-geographic structure. The rest of the package, including GBIF retrieval and large disk-based workflows, is built to feed that core mapping step in a transparent and reproducible way.

## References

Abell, R., Thieme, M. L., Revenga, C., Bryer, M., Kottelat, M., Bogutskaya, N., … Petry, P. (2008). Freshwater ecoregions of the world: a new map of biogeographic units for freshwater biodiversity conservation. *BioScience*, 58(5), 403–414. https://doi.org/10.1641/B580507

Chauvier, Y., Zimmermann, N. E., Poggiato, G., Bystrova, D., Brun, P., & Thuiller, W. (2021). Novel methods to correct for observer and sampling bias in presence-only species distribution models. *Global Ecology and Biogeography*, 30(11), 2312–2325. https://doi.org/10.1111/geb.13383

Denelle, P., Leroy, B., & Lenormand, M. (2025). Bioregionalization analyses with the bioregion R package. *Methods in Ecology and Evolution*, 16, 496–506. https://doi.org/10.1111/2041-210X.14496

Hagen, O., Vaterlaus, L., Albouy, C., Brown, A., Leugger, F., Onstein, R. E., Novaes de Santana, C., Scotese, C. R., & Pellissier, L. (2019). Mountain building, climate cooling and the richness of cold-adapted plants in the Northern Hemisphere. *Journal of Biogeography*, 46(8), 1792–1807. https://doi.org/10.1111/jbi.13653

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

Spalding, M. D., Fox, H. E., Allen, G. R., Davidson, N., Ferdaña, Z. A., Finlayson, M., … Robertson, J. (2007). Marine ecoregions of the world: a bioregionalization of coastal and shelf areas. *BioScience*, 57(7), 573–583. https://doi.org/10.1641/B570707

Spalding, M. D., Agostini, V. N., Rice, J., & Grant, S. M. (2012). Pelagic provinces of the world: a biogeographic classification of the world's surface pelagic waters. *Ocean & Coastal Management*, 60, 19–30. https://doi.org/10.1016/j.ocecoaman.2011.12.016

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

The Nature Conservancy (2012). Marine Ecoregions and Pelagic Provinces of the World. Cambridge (UK): The Nature Conservancy. https://habitats.oceanplus.org