Spatial overlays and raster context in biofetchR

library(biofetchR)

Purpose

Occurrence records become more useful once their coordinates are linked to the spatial units needed for an analysis. The right spatial framework depends on the question being asked: administrative boundaries support reporting, ecoregions support ecological interpretation, hydrological layers support freshwater context, marine regions support offshore attribution, and raster layers can add local environmental or anthropogenic-pressure information.

biofetchR supports these workflows by adding spatial context to imported occurrence records. This vignette focuses on how to choose vector overlays, add raster context, inspect spatial-attribution outputs, and diagnose common problems such as missing joins, boundary effects or unsuitable overlay choices.

The aim is not to add every available layer to every dataset. Start with the spatial framework that directly matches the analysis question, then add additional overlays or raster context only when they provide a clear interpretive benefit.

Using the examples

The examples use a mixture of local and live code. Local chunks create small example tables, inspect registry objects, or show workflow templates. These can be read and adapted without contacting external services.

Live chunks may submit GBIF requests, retrieve spatial layers, create provider caches, or write example outputs. They are controlled by the run_live parameter and are skipped by default.

To render a live version locally, use:

rmarkdown::render(
  "vignettes/biofetchR-spatial-overlays.Rmd",
  params = list(run_live = TRUE)
)

Use the default version to inspect the overlay workflow structure. Use the live version when you want to test the examples with GBIF credentials, spatial-provider access and local cache folders.

Vector overlays and raster context

biofetchR can add spatial context to occurrence records in two main ways: vector overlays and raster extraction.

Vector overlays assign occurrence points to spatial features with boundaries or geometries. These may include administrative units, terrestrial ecoregions, freshwater ecoregions, protected areas, Ramsar wetlands, EEZs, marine regions, rivers, lakes, basins, reservoirs, dams, wastewater-treatment sites or mining layers. The result is usually a feature name, identifier or category attached to each record.

Raster context extracts values from gridded datasets at the occurrence location. These values may describe environmental conditions, land cover, soil properties, anthropogenic pressure or other cell-based information.

The distinction is important because the outputs answer different spatial questions. A vector overlay asks which mapped feature a record falls within, intersects or is nearest to. Raster context asks what grid-cell value is associated with the record coordinates. Both can be useful, but they should be chosen according to the spatial meaning needed for the analysis.

Choosing a spatial framework

Choose the spatial framework according to the unit that should be attached to each occurrence record. The table below gives common matches between analysis goals, suitable spatial context and the fields users should expect to inspect in the output.

#>                                         analysis_goal
#> 1           Summarise records by administrative units
#> 2             Describe terrestrial ecological setting
#> 3              Describe freshwater ecological setting
#> 4    Add freshwater or aquatic-infrastructure context
#> 5               Add protected-area or wetland context
#> 6           Assign offshore records to marine regions
#> 7 Attach local environmental or human-pressure values
#>                                                     suitable_context
#> 1                                                               GADM
#> 2                                               TEOW or RESOLVE 2017
#> 3                                     FEOW, basins, rivers, or lakes
#> 4 Rivers, lakes, basins, GDW barriers, HydroWASTE, or GDW reservoirs
#> 5                                                     WDPA or Ramsar
#> 6                                             EEZ, LME, MEOW, or IHO
#> 7                                                     Raster context
#>                                     typical_output
#> 1             Administrative identifiers and names
#> 2    Terrestrial ecoregion, biome, or realm labels
#> 3 Freshwater ecoregion or hydrological identifiers
#> 4    Hydrological or infrastructure-context fields
#> 5         Protected-area or wetland-context fields
#> 6              Marine region identifiers and names
#> 7                 Raster-derived covariate columns

A good overlay choice should match the scale and meaning of the analysis. For example, GADM is appropriate when records need to be summarised by political or administrative units, while FEOW or basin layers are more suitable when the freshwater setting is the main unit of interpretation.

For terrestrial and freshwater workflows, start with the framework that will be used for the main summary or grouping variable. Add secondary overlays only when they answer a separate question, such as whether records also fall in protected areas, wetlands, river basins or infrastructure-linked features.

Example occurrence table

Spatial joins require occurrence records with coordinate columns that can be converted into point geometry. The example below is a minimal table with a taxon name, longitude, latitude, country code and year.

example_occurrences <- data.frame(
  species = c(
    "Xenopus laevis",
    "Xenopus laevis",
    "Carcinus maenas"
  ),
  decimalLongitude = c(2.35, 1.44, -5.10),
  decimalLatitude = c(48.86, 43.60, 50.15),
  countryCode = c("FR", "FR", "GB"),
  year = c(2019, 2021, 2020),
  stringsAsFactors = FALSE
)

example_occurrences
#>           species decimalLongitude decimalLatitude countryCode year
#> 1  Xenopus laevis             2.35           48.86          FR 2019
#> 2  Xenopus laevis             1.44           43.60          FR 2021
#> 3 Carcinus maenas            -5.10           50.15          GB 2020

In a full biofetchR workflow, point tables are usually created after GBIF records have been imported and processed. This local example is deliberately small so that the overlay examples below focus on the spatial-join step rather than the download workflow.

Cleaning and thinning before spatial overlays

Spatial overlays should be applied only after the coordinate columns have been checked. Records with missing coordinates, impossible longitude or latitude values, exact 0,0 coordinates, or very high coordinate uncertainty can be assigned to the wrong spatial unit. These records can then affect summaries by administrative unit, ecoregion, basin, protected area or marine region.

biofetchR provides direct access to these checks through thin_spatial_points(). The same cleaning and thinning controls are also exposed inside the high-level terrestrial/freshwater and marine pipelines. Using the helper directly is useful when users want to inspect coordinate processing before running a full pipeline, or when they already have an occurrence table that only needs spatial preparation before overlay assignment.

The main options are:

#>                             option                            argument
#> 1             Coordinate screening                           automatic
#> 2 Coordinate-uncertainty filtering filter_uncertain, max_uncertainty_m
#> 3        High-uncertainty flagging                  warn_uncertainty_m
#> 4               Cleaning-only mode                         dist_km = 0
#> 5          Distance-based thinning                         dist_km > 0
#> 6           Density-based thinning                         max_density
#> 7         Priority-based retention                        priority_col
#> 8                Diagnostic output                   return_all = TRUE
#>                                                                                  purpose
#> 1                Remove missing, non-finite, impossible, or exact zero-zero coordinates.
#> 2                   Remove records with coordinate uncertainty above a chosen threshold.
#> 3                                Flag retained records with high coordinate uncertainty.
#> 4                                            Apply coordinate cleaning without thinning.
#> 5         Reduce spatial clustering by retaining points separated by a minimum distance.
#> 6                               Limit the approximate number of retained points per km².
#> 7                         Preferentially retain higher-priority records during thinning.
#> 8 Return all post-cleaning rows with a `.thinned` flag showing which rows were retained.

The example below deliberately includes several common coordinate problems: valid records, nearby clustered records, an impossible longitude, an exact 0,0 point, a missing coordinate and one record with high coordinate uncertainty. It also includes individualCount as an example priority field that can be used during thinning.

# Example occurrence table for demonstrating cleaning and thinning.
# This is deliberately small and includes several common coordinate problems.
cleaning_example <- data.frame(
  species = c(
    "Xenopus laevis",
    "Xenopus laevis",
    "Xenopus laevis",
    "Xenopus laevis",
    "Xenopus laevis",
    "Xenopus laevis",
    "Xenopus laevis"
  ),

  # Longitude values include valid coordinates, an impossible longitude,
  # an exact zero-zero point, and a missing value.
  decimalLongitude = c(
    2.3500,
    2.3505,
    2.3520,
    181.0000,
    0.0000,
    NA,
    1.4400
  ),

  # Latitude values include valid coordinates and an exact zero-zero point.
  decimalLatitude = c(
    48.8600,
    48.8605,
    48.8620,
    48.8600,
    0.0000,
    43.6000,
    43.6000
  ),

  # Coordinate uncertainty is optional in GBIF-style data.
  # Here, one otherwise valid record has high uncertainty.
  coordinateUncertaintyInMeters = c(
    30,
    50,
    12000,
    20,
    10,
    100,
    300
  ),

  # individualCount is used here as an example priority column.
  # Higher values are retained first during thinning when priority_col is used.
  individualCount = c(
    10,
    2,
    1,
    1,
    1,
    1,
    5
  ),

  year = c(2019, 2019, 2020, 2020, 2021, 2021, 2022),

  stringsAsFactors = FALSE
)

cleaning_example
#>          species decimalLongitude decimalLatitude coordinateUncertaintyInMeters
#> 1 Xenopus laevis           2.3500         48.8600                            30
#> 2 Xenopus laevis           2.3505         48.8605                            50
#> 3 Xenopus laevis           2.3520         48.8620                         12000
#> 4 Xenopus laevis         181.0000         48.8600                            20
#> 5 Xenopus laevis           0.0000          0.0000                            10
#> 6 Xenopus laevis               NA         43.6000                           100
#> 7 Xenopus laevis           1.4400         43.6000                           300
#>   individualCount year
#> 1              10 2019
#> 2               2 2019
#> 3               1 2020
#> 4               1 2020
#> 5               1 2021
#> 6               1 2021
#> 7               5 2022

Cleaning only

Use cleaning-only mode when records should be checked for coordinate quality but their spatial clustering should be left unchanged. This is useful before assigning records to an overlay when the analysis needs all retained locations, rather than a thinned subset.

Set dist_km = 0 to apply the coordinate checks without distance-based thinning. In the example below, records with missing, impossible or exact 0,0 coordinates are removed. Records with coordinate uncertainty greater than max_uncertainty_m are also removed because filter_uncertain = TRUE.

# Apply coordinate cleaning only.
# Records with missing, impossible, or exact zero-zero coordinates are removed.
# Records above max_uncertainty_m are also removed because filter_uncertain = TRUE.
cleaned_only <- thin_spatial_points(
  cleaning_example,
  dist_km = 0,
  filter_uncertain = TRUE,
  max_uncertainty_m = 10000,
  quiet = TRUE
)

# Drop geometry for compact display in the vignette.
sf::st_drop_geometry(cleaned_only)
#>          species decimalLongitude decimalLatitude coordinateUncertaintyInMeters
#> 1 Xenopus laevis           2.3500         48.8600                            30
#> 2 Xenopus laevis           2.3505         48.8605                            50
#> 7 Xenopus laevis           1.4400         43.6000                           300
#>   individualCount year .high_uncertainty
#> 1              10 2019             FALSE
#> 2               2 2019             FALSE
#> 7               5 2022             FALSE

The returned table keeps the records that passed the coordinate-quality checks. This is the safest option when the user wants cleaner coordinates for overlay assignment but does not want to remove nearby records only because they are spatially clustered.

Keeping high-uncertainty records but flagging them

High coordinate uncertainty does not always mean a record should be removed. In some workflows, users may want to retain uncertain records so they can be reviewed later, compared with other evidence, or excluded only in a sensitivity analysis.

Set filter_uncertain = FALSE to keep these records. Records above warn_uncertainty_m are retained, but flagged in the .high_uncertainty column so they can be inspected before overlay results are interpreted.

# Retain high-uncertainty records but flag them for review.
cleaned_with_uncertainty_flags <- thin_spatial_points(
  cleaning_example,
  dist_km = 0,
  filter_uncertain = FALSE,
  warn_uncertainty_m = 5000,
  quiet = TRUE
)

# Inspect the retained records and the uncertainty flag.
sf::st_drop_geometry(cleaned_with_uncertainty_flags)[
  ,
  intersect(
    c(
      "species",
      "decimalLongitude",
      "decimalLatitude",
      "coordinateUncertaintyInMeters",
      ".high_uncertainty"
    ),
    names(cleaned_with_uncertainty_flags)
  ),
  drop = FALSE
]
#>          species decimalLongitude decimalLatitude coordinateUncertaintyInMeters
#> 1 Xenopus laevis           2.3500         48.8600                            30
#> 2 Xenopus laevis           2.3505         48.8605                            50
#> 3 Xenopus laevis           2.3520         48.8620                         12000
#> 7 Xenopus laevis           1.4400         43.6000                           300
#>   .high_uncertainty
#> 1             FALSE
#> 2             FALSE
#> 3              TRUE
#> 7             FALSE

This mode is useful when uncertain records should remain visible rather than being removed automatically. The flag allows users to decide later whether those records are suitable for the selected overlay, should be checked manually, or should be excluded from a stricter analysis.

Distance-based thinning

Spatial thinning reduces clustering by retaining points that are at least a chosen distance apart. This is useful when heavily sampled locations would otherwise dominate maps, summaries, or models. It can also reduce sampling-driven spatial autocorrelation in downstream analyses, although it should not be treated as a substitute for checking spatial structure during model diagnostics.

# Apply coordinate cleaning and then distance-based thinning.
# return_all = TRUE keeps all cleaned records and adds a `.thinned` diagnostic
# column showing which rows were retained by thinning.
distance_thinned <- thin_spatial_points(
  cleaning_example,
  dist_km = 5,
  priority_col = "individualCount",
  return_all = TRUE,
  filter_uncertain = FALSE,
  warn_uncertainty_m = 5000,
  quiet = TRUE
)

# Inspect thinning diagnostics.
sf::st_drop_geometry(distance_thinned)[
  ,
  intersect(
    c(
      "species",
      "decimalLongitude",
      "decimalLatitude",
      "individualCount",
      "coordinateUncertaintyInMeters",
      ".high_uncertainty",
      ".thinned"
    ),
    names(distance_thinned)
  ),
  drop = FALSE
]
#>          species decimalLongitude decimalLatitude individualCount
#> 1 Xenopus laevis           2.3500         48.8600              10
#> 7 Xenopus laevis           1.4400         43.6000               5
#> 2 Xenopus laevis           2.3505         48.8605               2
#> 3 Xenopus laevis           2.3520         48.8620               1
#>   coordinateUncertaintyInMeters .high_uncertainty .thinned
#> 1                            30             FALSE     TRUE
#> 7                           300             FALSE     TRUE
#> 2                            50             FALSE    FALSE
#> 3                         12000              TRUE    FALSE

Distance-based thinning

Distance-based thinning retains a subset of records so that kept points are separated by at least the chosen distance. This is useful when many records come from the same locality or from very densely sampled areas, and those clusters would otherwise dominate maps, overlay summaries or spatial analyses. In downstream modelling, strong spatial clustering can also contribute to spatial autocorrelation, so thinning can help reduce one source of sampling-driven dependence before model fitting.

In this example, dist_km = 5 applies a 5 km thinning distance after coordinate cleaning. The priority_col argument is used to retain higher-priority records first when nearby records compete for the same retained location. Here, individualCount is used only as an example priority field.

# Apply coordinate cleaning and then distance-based thinning.
# return_all = TRUE keeps all cleaned records and adds a `.thinned` diagnostic
# column showing which rows were retained by thinning.
distance_thinned <- thin_spatial_points(
  cleaning_example,
  dist_km = 5,
  priority_col = "individualCount",
  return_all = TRUE,
  filter_uncertain = FALSE,
  warn_uncertainty_m = 5000,
  quiet = TRUE
)

# Inspect thinning diagnostics.
sf::st_drop_geometry(distance_thinned)[
  ,
  intersect(
    c(
      "species",
      "decimalLongitude",
      "decimalLatitude",
      "individualCount",
      "coordinateUncertaintyInMeters",
      ".high_uncertainty",
      ".thinned"
    ),
    names(distance_thinned)
  ),
  drop = FALSE
]
#>          species decimalLongitude decimalLatitude individualCount
#> 1 Xenopus laevis           2.3500         48.8600              10
#> 7 Xenopus laevis           1.4400         43.6000               5
#> 2 Xenopus laevis           2.3505         48.8605               2
#> 3 Xenopus laevis           2.3520         48.8620               1
#>   coordinateUncertaintyInMeters .high_uncertainty .thinned
#> 1                            30             FALSE     TRUE
#> 7                           300             FALSE     TRUE
#> 2                            50             FALSE    FALSE
#> 3                         12000              TRUE    FALSE

Because return_all = TRUE, the output keeps all cleaned records and adds the .thinned column. Rows marked TRUE are the records retained by thinning; rows marked FALSE passed coordinate cleaning but were not selected for the thinned subset. This makes it possible to inspect what thinning changed before using the retained records in an overlay workflow.

Thinning is a sampling-bias reduction step, not a complete solution to spatial autocorrelation. If the retained records are used in models, users should still check spatial structure during model diagnostics and consider whether additional spatial terms, blocking, or sensitivity analyses are needed.

Cleaning and thinning inside the pipelines

The same cleaning and thinning controls can be used directly inside the high-level pipelines. This is usually the simplest route when records are being downloaded, processed and assigned to overlays in one workflow.

Use apply_cleaning = TRUE to run coordinate-quality checks before spatial attribution. Use apply_thinning = TRUE with dist_km to apply distance-based thinning after cleaning.

process_gbif_terrestrial_freshwater_pipeline(
  df = data.frame(
    species = "Xenopus laevis",
    iso2c = "FR",
    stringsAsFactors = FALSE
  ),
  output_dir = file.path(output_root, "example_outputs"),
  user = Sys.getenv("GBIF_USER"),
  pwd = Sys.getenv("GBIF_PWD"),
  email = Sys.getenv("GBIF_EMAIL"),

  region_source = "gadm",
  gadm_unit = 1,
  cache_dir_gadm = file.path(output_root, "cache", "gadm"),

  # Apply coordinate quality control.
  apply_cleaning = TRUE,

  # Apply spatial thinning after cleaning.
  apply_thinning = TRUE,

  # Retain records separated by at least 5 km where possible.
  dist_km = 5,

  export_summary = TRUE,
  store_in_memory = FALSE
)

After the run, inspect the record-count fields in gbif_summary.csv. These show how many records were returned from GBIF, how many remained after coordinate cleaning, and how many remained after thinning. Check these values before interpreting the overlay fields, especially when few records remain or when thinning removes a large proportion of the cleaned dataset.

Primary and secondary spatial context

A primary spatial framework is the main unit used to organise the output. This is usually the unit that records will be grouped by, summarised across, mapped or used in later analyses. For example, a terrestrial workflow may use GADM level 1 as the primary framework when the main output needs to be interpreted by subnational administrative units.

A secondary overlay adds another layer of context to the same records. It should not be treated as a replacement for the primary framework. For example, a freshwater workflow may retain GADM level 1 as the reporting unit while also adding FEOW identifiers to describe the freshwater ecoregion associated with each record.

The key point is to keep these fields interpretable. Administrative units, ecoregions, hydrological features, protected areas, marine regions and raster values describe different spatial concepts. Keeping them as separate fields makes it easier to check the output, compare alternative summaries and decide which spatial unit is appropriate for each part of the analysis.

#>                            workflow                          primary_context
#> 1 Subnational terrestrial reporting                             GADM level 1
#> 2      Freshwater invasion analysis               GADM level 1 or basin unit
#> 3          Protected-area screening                GADM or ecological region
#> 4    Marine jurisdictional analysis                                      EEZ
#> 5  Environmental-covariate analysis GADM, ecoregion, basin, or marine region
#>                       secondary_context
#> 1                  TEOW or RESOLVE 2017
#> 2 FEOW, rivers, lakes, dams, reservoirs
#> 3                        WDPA or Ramsar
#> 4          LME, MEOW, IHO, or high seas
#> 5                        Raster context
#>                                                                                   reason
#> 1                Administrative outputs can be interpreted alongside ecological regions.
#> 2            Reporting units and freshwater ecological units answer different questions.
#> 3 Protected-area context is usually an additional attribute, not the main grouping unit.
#> 4                  Marine records need marine spatial units rather than land boundaries.
#> 5          Raster values describe point-level covariates rather than polygon membership.

Use the primary framework for the main grouping structure, then add secondary fields only when they support a specific interpretation. This avoids forcing different spatial concepts into one combined unit too early.

Common overlay types

The table below provides a quick reference to overlay names commonly used in biofetchR workflows. It is intended as a guide to the main spatial meaning of each overlay, not as a complete description of every provider layer.

Availability can depend on optional packages, internet access, provider availability and local cache settings. Before scaling up, test the selected overlay on a small example and check that the output contains the spatial names or identifiers you need.

overlay_reference <- data.frame(
  overlay = c(
    "gadm",
    "teow",
    "feow",
    "lakes",
    "rivers",
    "basins",
    "ne_admin1",
    "resolve2017",
    "wdpa",
    "ramsar",
    "gdw_barriers",
    "hydrowaste",
    "gdw_reservoirs",
    "global_mining",
    "eez",
    "lme",
    "meow",
    "iho"
  ),
  typical_use = c(
    "Administrative attribution for terrestrial and freshwater records",
    "Terrestrial ecoregion context",
    "Freshwater ecoregion context",
    "Lake-associated occurrence context",
    "River-associated occurrence context",
    "Hydrological basin context",
    "Natural Earth administrative-unit context",
    "RESOLVE ecoregion and biome context",
    "Protected-area context",
    "Ramsar wetland context",
    "Barrier and dam context",
    "Wastewater-treatment context",
    "Reservoir context",
    "Mining or extraction-pressure context",
    "Marine jurisdictional attribution",
    "Large Marine Ecosystem attribution",
    "Marine ecoregion attribution",
    "Sea-area attribution"
  ),
  broad_domain = c(
    "terrestrial/freshwater",
    "terrestrial",
    "freshwater",
    "freshwater",
    "freshwater",
    "freshwater",
    "terrestrial/freshwater",
    "terrestrial",
    "terrestrial/freshwater/marine depending on source",
    "wetland/freshwater/coastal",
    "freshwater",
    "freshwater/infrastructure",
    "freshwater",
    "terrestrial/infrastructure",
    "marine",
    "marine",
    "marine",
    "marine"
  ),
  stringsAsFactors = FALSE
)

overlay_reference
#>           overlay
#> 1            gadm
#> 2            teow
#> 3            feow
#> 4           lakes
#> 5          rivers
#> 6          basins
#> 7       ne_admin1
#> 8     resolve2017
#> 9            wdpa
#> 10         ramsar
#> 11   gdw_barriers
#> 12     hydrowaste
#> 13 gdw_reservoirs
#> 14  global_mining
#> 15            eez
#> 16            lme
#> 17           meow
#> 18            iho
#>                                                          typical_use
#> 1  Administrative attribution for terrestrial and freshwater records
#> 2                                      Terrestrial ecoregion context
#> 3                                       Freshwater ecoregion context
#> 4                                 Lake-associated occurrence context
#> 5                                River-associated occurrence context
#> 6                                         Hydrological basin context
#> 7                          Natural Earth administrative-unit context
#> 8                                RESOLVE ecoregion and biome context
#> 9                                             Protected-area context
#> 10                                            Ramsar wetland context
#> 11                                           Barrier and dam context
#> 12                                      Wastewater-treatment context
#> 13                                                 Reservoir context
#> 14                             Mining or extraction-pressure context
#> 15                                 Marine jurisdictional attribution
#> 16                                Large Marine Ecosystem attribution
#> 17                                      Marine ecoregion attribution
#> 18                                              Sea-area attribution
#>                                         broad_domain
#> 1                             terrestrial/freshwater
#> 2                                        terrestrial
#> 3                                         freshwater
#> 4                                         freshwater
#> 5                                         freshwater
#> 6                                         freshwater
#> 7                             terrestrial/freshwater
#> 8                                        terrestrial
#> 9  terrestrial/freshwater/marine depending on source
#> 10                        wetland/freshwater/coastal
#> 11                                        freshwater
#> 12                         freshwater/infrastructure
#> 13                                        freshwater
#> 14                        terrestrial/infrastructure
#> 15                                            marine
#> 16                                            marine
#> 17                                            marine
#> 18                                            marine

Use this table to narrow the overlay choice before running a pipeline. For example, choose gadm when the output needs administrative units, feow or basin layers when freshwater context is central, wdpa or ramsar when conservation-site context is needed, and marine overlays when records should be interpreted in offshore or coastal marine space.

GADM administrative attribution

GADM is used when occurrence records need to be assigned to administrative boundaries. In biofetchR, this is most useful when the analysis needs outputs by country, first-level administrative unit, or finer subnational subdivision.

Administrative attribution is different from ecological or hydrological context. It is useful when records need to be summarised for reporting units, management regions, national or subnational comparisons, or first-detection summaries within political boundaries.

In the high-level pipeline, GADM attribution is requested with region_source = "gadm". The gadm_unit argument sets the administrative level: level 0 refers to the country, level 1 to the first subnational unit, and level 2 to the next subdivision below that.

# Load GADM level 1 polygons for Portugal.
# This is useful for checking that GADM can be retrieved and cached before
# running a larger GBIF workflow.
gadm_prt_l1 <- run_live_safely(
  load_gadm(
    iso2c = "PT",
    level = 1,
    cache_dir = file.path(output_root, "cache", "gadm"),
    quiet = FALSE
  )
)

# Inspect the returned GADM object if loading succeeded.
if (!is.null(gadm_prt_l1)) {
  gadm_prt_l1
}

Choose the GADM level according to the scale of the question and the density of the occurrence data. Level 0 is suitable for country-level summaries, but it can hide important subnational variation. Level 2 provides finer spatial detail, but may produce sparse or uneven outputs when records are limited.

GADM level 1 is often a practical default for subnational work. It provides more detail than country-level summaries while usually retaining enough records per unit for maps, tables and comparative analyses.

GADM inside the terrestrial and freshwater pipeline

GADM can also be used directly inside the terrestrial and freshwater pipeline. The example below runs a small species-country workflow, processes the GBIF records, and assigns the retained occurrences to GADM level 1 units.

# Run the live GADM pipeline example only when GBIF credentials are available.
if (!gbif_ready) {

  # Skip cleanly rather than failing the vignette.
  skip_live_message("set GBIF_USER, GBIF_PWD, and GBIF_EMAIL to run this example")

} else {

  # Run a minimal terrestrial/freshwater pipeline with GADM attribution.
  gadm_pipeline_result <- run_live_safely(
    process_gbif_terrestrial_freshwater_pipeline(

      # Species-country input table.
      df = data.frame(
        species = "Xenopus laevis",
        iso2c = "FR",
        stringsAsFactors = FALSE
      ),

      # Output folder for this example.
      output_dir = file.path(output_root, "gadm_pipeline"),

      # GBIF credentials read from environment variables.
      user = gbif_user,
      pwd = gbif_pwd,
      email = gbif_email,

      # Use GADM level 1 as the primary spatial framework.
      region_source = "gadm",
      gadm_unit = 1,
      cache_dir_gadm = file.path(output_root, "cache", "gadm"),

      # Clean coordinates and thin retained records.
      apply_cleaning = TRUE,
      apply_thinning = TRUE,
      dist_km = 5,

      # Keep the example small and write outputs to disk.
      batch_size = 1,
      export_summary = TRUE,
      store_in_memory = FALSE,
      return_all_results = FALSE,

      # Show progress messages.
      quiet = FALSE
    )
  )

  gadm_pipeline_result
}

After a live run, inspect gbif_summary.csv before opening the exported occurrence file. The summary shows whether the row completed, how many records were retained after processing, and where the output file was written.

# Path to the summary file written by the GADM example.
summary_path <- file.path(
  output_root,
  "gadm_pipeline",
  "gbif_summary.csv"
)

# Read the summary if it exists.
if (file.exists(summary_path)) {
  gadm_summary <- read.csv(summary_path, stringsAsFactors = FALSE)
  gadm_summary
} else {
  skip_live_message("summary file was not available from the GADM live example")
}

For GADM workflows, check that the exported records contain the expected administrative identifiers or names for the selected level. Rows without a GADM assignment should be reviewed before summaries are produced, especially when records fall near borders, coastlines or small administrative units.

Combining administrative and ecological context

Some workflows need two spatial descriptions for the same occurrence records. For example, GADM level 1 may be used for subnational reporting, while FEOW adds the freshwater ecoregion associated with each retained record.

These fields should usually remain separate because they represent different spatial concepts. GADM identifies administrative location; FEOW identifies freshwater ecological context. Keeping both fields visible makes it easier to summarise records by the main reporting unit while still retaining ecological information for interpretation.

The example below uses GADM as the primary administrative framework and FEOW as a secondary ecological overlay. Setting overlay_mode = "dual_separate" keeps the two joins separate in the exported output.

# Run the live GADM + FEOW example only when GBIF credentials are available.
if (!gbif_ready) {

  skip_live_message("set GBIF_USER, GBIF_PWD, and GBIF_EMAIL to run this example")

} else {

  # Run a pipeline with both administrative and freshwater ecological context.
  feow_result <- run_live_safely(
    process_gbif_terrestrial_freshwater_pipeline(
      df = data.frame(
        species = "Xenopus laevis",
        iso2c = "FR",
        stringsAsFactors = FALSE
      ),

      output_dir = file.path(output_root, "feow_pipeline"),

      user = gbif_user,
      pwd = gbif_pwd,
      email = gbif_email,

      # Use GADM and FEOW together.
      region_source = c("gadm", "feow"),

      # Keep the two spatial contexts as separate output fields.
      overlay_mode = "dual_separate",

      # GADM and FEOW cache directories.
      gadm_unit = 1,
      cache_dir_gadm = file.path(output_root, "cache", "gadm"),
      feow_cache_dir = file.path(output_root, "cache", "feow"),

      apply_cleaning = TRUE,
      apply_thinning = TRUE,
      dist_km = 5,

      batch_size = 1,
      export_summary = TRUE,
      store_in_memory = FALSE,
      return_all_results = FALSE,
      quiet = FALSE
    )
  )

  feow_result
}

After the run, inspect the exported occurrence table and gbif_summary.csv together. The occurrence table shows the spatial fields added to each retained record, while the summary file shows whether the row completed successfully and where the output was written.

If expected overlay fields are missing or mostly empty, check the selected overlay, cache folders and records without spatial assignments before using the output in summaries.

Direct provider loaders

Some spatial layers can be loaded directly before they are used in a pipeline. This is useful when the layer itself needs to be inspected, cached or reused across more than one workflow.

Direct loading separates provider access from occurrence processing. If a layer does not download, cache or open correctly, users can fix that issue before running a GBIF workflow. It also allows users to inspect the attributes and coverage of a layer before deciding whether it is suitable for their records.

Direct provider loaders are most useful when:

# Load Natural Earth administrative level 1 context for Great Britain and Ireland.
ne_admin1 <- run_live_safely(
  bf_load_ne_admin1(
    iso2c = c("GB", "IE"),
    cache_dir = file.path(output_root, "cache", "ne_admin1"),
    quiet = FALSE
  )
)

# Load a RESOLVE 2017 terrestrial ecoregion subset.
resolve2017 <- run_live_safely(
  bf_load_resolve2017(
    iso2c = c("GB", "IE"),
    cache_dir = file.path(output_root, "cache", "resolve2017"),
    quiet = FALSE
  )
)

if (!is.null(ne_admin1)) {
  ne_admin1
}

if (!is.null(resolve2017)) {
  resolve2017
}

After loading a layer directly, check that it contains the expected records, attribute columns and geometry. If the layer looks correct, it can then be used with greater confidence inside a pipeline or in a separate spatial-join workflow.

Using overlays with an existing occurrence table

Some projects already have cleaned occurrence records and do not need to submit a new GBIF download. In those cases, users can load an overlay directly and join it to their own point table.

This approach gives more control, but also moves more responsibility to the analysis script. Users need to make sure that coordinates are converted to point geometry, coordinate reference systems are compatible, joins are checked, and unmatched records are reviewed.

The template below is not evaluated by default. It shows the basic structure of a custom point-in-polygon join using an existing occurrence table and a polygon overlay.

# Example only: start from an existing occurrence table with coordinates.
existing_occurrences <- example_occurrences

# Convert the occurrence table to an sf point object.
occ_sf <- sf::st_as_sf(
  existing_occurrences,
  coords = c("decimalLongitude", "decimalLatitude"),
  crs = 4326,
  remove = FALSE
)

# Load a provider layer.
resolve2017 <- bf_load_resolve2017(
  iso2c = c("GB", "FR"),
  cache_dir = file.path(output_root, "cache", "resolve2017"),
  quiet = FALSE
)

# Join points to polygons.
# In a real workflow, inspect unmatched points and check CRS/geometry validity.
occ_with_resolve <- sf::st_join(
  occ_sf,
  resolve2017,
  left = TRUE
)

occ_with_resolve

This pattern is appropriate for polygon overlays. Other spatial resources may need different join logic. For example, line or point-like layers may require nearest-feature snapping or a distance threshold rather than a simple point-in-polygon join.

Custom joins are most useful when the occurrence table has already been reviewed or when the project needs a spatial operation that sits outside the high-level pipelines. For standard GBIF-based workflows, the pipelines are usually simpler because download, cleaning, thinning, attribution and summary outputs are handled together.

Raster context

Raster context is used when retained occurrence records need values from gridded datasets. Instead of assigning a point to a named polygon or feature, the workflow extracts the raster value at the record location and appends it to the occurrence table.

This is useful for variables such as land-cover class, soil properties, elevation, human footprint or other environmental and human-pressure layers. The added columns can then be used as model covariates, habitat or land-use filters, or descriptive summaries of the conditions associated with retained records.

Raster extraction is best added after the main occurrence and vector-overlay workflow is working. Raster layers can be large, and extraction may involve provider access, local caching, cropping, reading gridded files and checking coordinate reference systems. A small one-species test makes it easier to confirm that the selected layers are available and that non-missing values are being returned.

# Run the raster-context example only when GBIF credentials are available.
if (!gbif_ready) {

  skip_live_message("set GBIF_USER, GBIF_PWD, and GBIF_EMAIL to run this example")

} else {

  # Run a small pipeline with GADM attribution and raster context.
  raster_context_result <- run_live_safely(
    process_gbif_terrestrial_freshwater_pipeline(
      df = data.frame(
        species = "Xenopus laevis",
        iso2c = "FR",
        stringsAsFactors = FALSE
      ),

      output_dir = file.path(output_root, "raster_context"),

      user = gbif_user,
      pwd = gbif_pwd,
      email = gbif_email,

      # Keep GADM as the primary vector framework.
      region_source = "gadm",
      gadm_unit = 1,
      cache_dir_gadm = file.path(output_root, "cache", "gadm"),

      # Add gridded environmental or human-pressure context.
      raster_context = c("worldcover", "soilgrids", "human_footprint"),
      raster_cache_dir = file.path(output_root, "cache", "raster"),

      apply_cleaning = TRUE,
      apply_thinning = TRUE,
      dist_km = 5,

      batch_size = 1,
      export_summary = TRUE,
      store_in_memory = FALSE,
      return_all_results = FALSE,
      quiet = FALSE
    )
  )

  raster_context_result
}

After the run, inspect the exported occurrence file for the raster-derived columns. If the columns are missing or contain mostly missing values, check the selected raster source, cache directory, coordinate columns and whether the retained records fall within the raster extent.

# Run the raster-context example only when GBIF credentials are available.
if (!gbif_ready) {

  skip_live_message("set GBIF_USER, GBIF_PWD, and GBIF_EMAIL to run this example")

} else {

  # Run a small pipeline with GADM attribution and raster context.
  raster_context_result <- run_live_safely(
    process_gbif_terrestrial_freshwater_pipeline(
      df = data.frame(
        species = "Xenopus laevis",
        iso2c = "FR",
        stringsAsFactors = FALSE
      ),

      output_dir = file.path(output_root, "raster_context"),

      user = gbif_user,
      pwd = gbif_pwd,
      email = gbif_email,

      # Keep GADM as the primary vector framework.
      region_source = "gadm",
      gadm_unit = 1,
      cache_dir_gadm = file.path(output_root, "cache", "gadm"),

      # Add gridded environmental or human-pressure context.
      raster_context = c("worldcover", "soilgrids", "human_footprint"),
      raster_cache_dir = file.path(output_root, "cache", "raster"),

      apply_cleaning = TRUE,
      apply_thinning = TRUE,
      dist_km = 5,

      batch_size = 1,
      export_summary = TRUE,
      store_in_memory = FALSE,
      return_all_results = FALSE,
      quiet = FALSE
    )
  )

  raster_context_result
}

After the run, inspect the exported occurrence file for the raster-derived columns. If the columns are missing or contain mostly missing values, check the selected raster source, cache directory, coordinate columns and whether the retained records fall within the raster extent.

Protected-area and wetland overlays

WDPA and Ramsar overlays are used when occurrence records need conservation-site or designated-wetland context.These overlays are useful for questions such as whether retained records overlap protected area, whether they fall within a Ramsar-listed wetland, or which designated feature is associated with a record.

WDPA and Ramsar are external resources with licensing and attribution requirements, so the examples are opt-in. Users should enable these chunks only after checking that their intended use is consistent with the relevant provider terms.

# WDPA and Ramsar loading is opt-in because these datasets have
# external licensing and attribution requirements.
#
# To run this chunk locally, set the option below in the console before
# rendering the vignette:
#
# options(biofetchR.wdpa_opt_in = TRUE)

protected_overlay_opt_in <- isTRUE(
  getOption("biofetchR.wdpa_opt_in", FALSE)
)

if (!protected_overlay_opt_in) {

  skip_live_message(
    "set options(biofetchR.wdpa_opt_in = TRUE) after checking provider terms"
  )

} else {

  # Load protected-area polygons for Great Britain.
  wdpa_gb <- run_live_safely(
    bf_load_wdpa(
      iso2c = "GB",
      cache_dir = file.path(output_root, "cache", "wdpa"),
      opt_in = TRUE,
      quiet = FALSE
    )
  )

  # Load Ramsar wetland polygons for Great Britain.
  ramsar_gb <- run_live_safely(
    bf_load_ramsar(
      iso2c = "GB",
      cache_dir = file.path(output_root, "cache", "ramsar"),
      opt_in = TRUE,
      quiet = FALSE
    )
  )

  list(
    wdpa = wdpa_gb,
    ramsar = ramsar_gb
  )
}

Inside the terrestrial/freshwater pipeline, WDPA and Ramsar can be added as secondary overlays. The example below keeps GADM as the administrative framework and adds WDPA as protected-area context.

# Run the live protected-area overlay example only when GBIF credentials
# are available and WDPA opt-in has been confirmed.
if (!gbif_ready) {

  skip_live_message("set GBIF_USER, GBIF_PWD, and GBIF_EMAIL to run this example")

} else if (!isTRUE(getOption("biofetchR.wdpa_opt_in", FALSE))) {

  skip_live_message(
    "set options(biofetchR.wdpa_opt_in = TRUE) after checking provider terms"
  )

} else {

  protected_overlay_result <- run_live_safely(
    process_gbif_terrestrial_freshwater_pipeline(
      df = data.frame(
        species = "Xenopus laevis",
        iso2c = "GB",
        stringsAsFactors = FALSE
      ),

      output_dir = file.path(output_root, "protected_overlay_pipeline"),

      user = gbif_user,
      pwd = gbif_pwd,
      email = gbif_email,

      # Keep GADM as the main administrative framework and add WDPA.
      region_source = c("gadm", "wdpa"),
      overlay_mode = "dual_separate",
      gadm_unit = 1,
      cache_dir_gadm = file.path(output_root, "cache", "gadm"),

      # WDPA is opt-in because it is an externally licensed dataset.
      wdpa_cache_dir = file.path(output_root, "cache", "wdpa"),
      wdpa_opt_in = TRUE,
      wdpa_exclude_marine = TRUE,

      apply_cleaning = TRUE,
      apply_thinning = TRUE,
      dist_km = 5,

      batch_size = 1,
      export_summary = TRUE,
      store_in_memory = FALSE,
      quiet = FALSE
    )
  )

  protected_overlay_result
}

When interpreting these outputs, treat WDPA and Ramsar fields as site-overlap information. A missing match does not mean that a record is unimportant or that the surrounding habitat has no conservation value. It only means that the record was not assigned to the selected protected-area or wetland layer under the settings used in that run.

Inspecting spatial outputs

After a spatial workflow, start with the summary file. This gives a compact view of row status, record retention and exported file paths before the occurrence tables are opened.

spatial_summary_example <- data.frame(
  species = c(
    "Xenopus laevis",
    "Xenopus laevis",
    "Carcinus maenas"
  ),
  workflow = c(
    "gadm",
    "gadm_plus_feow",
    "marine_eez"
  ),
  primary_spatial_unit = c(
    "GADM level 1",
    "GADM level 1",
    "EEZ"
  ),
  secondary_spatial_context = c(
    NA,
    "FEOW",
    NA
  ),
  n_total = c(1340, 1340, 980),
  n_cleaned = c(1309, 1309, 944),
  n_thinned = c(861, 861, 702),
  status = c("success", "success", "success"),
  output_file = c(
    "Xenopus_laevis_FR_gadm1.csv",
    "Xenopus_laevis_FR_gadm1_feow.csv",
    "Carcinus_maenas_eez.csv"
  ),
  stringsAsFactors = FALSE
)

spatial_summary_example
#>           species       workflow primary_spatial_unit secondary_spatial_context
#> 1  Xenopus laevis           gadm         GADM level 1                      <NA>
#> 2  Xenopus laevis gadm_plus_feow         GADM level 1                      FEOW
#> 3 Carcinus maenas     marine_eez                  EEZ                      <NA>
#>   n_total n_cleaned n_thinned  status                      output_file
#> 1    1340      1309       861 success      Xenopus_laevis_FR_gadm1.csv
#> 2    1340      1309       861 success Xenopus_laevis_FR_gadm1_feow.csv
#> 3     980       944       702 success          Carcinus_maenas_eez.csv

The values above are illustrative. In a real run, read gbif_summary.csv from the relevant output directory and inspect the columns most useful for checking spatial processing.

summary_path <- file.path("path", "to", "gbif_summary.csv")

summary_tbl <- read.csv(summary_path, stringsAsFactors = FALSE)

summary_tbl[
  ,
  intersect(
    c(
      "species",
      "region_id",
      "region_type",
      "n_total",
      "n_cleaned",
      "n_thinned",
      "status",
      "fail_stage",
      "fail_reason",
      "output_file"
    ),
    names(summary_tbl)
  ),
  drop = FALSE
]

Focus on three checks: whether each row completed, whether record counts changed as expected during cleaning and thinning, and whether each successful row has an exported file path. Rows with failed status, missing paths or unexpectedly low record counts should be reviewed before the spatial outputs are summarised.

Troubleshooting spatial joins

When a spatial workflow gives unexpected output, use the summary file to identify where the problem occurred. Most issues fall into one of four areas: records were removed during coordinate processing, the spatial layer did not load, the join did not assign many records, or the output was not written where expected.

#>                               symptom
#> 1           No records after cleaning
#> 2   Many records have no spatial unit
#> 3        Provider layer fails to load
#> 4                Raster context fails
#> 5          Marine records are missing
#> 6 Output files are not where expected
#> 7               Very large memory use
#>                                                                                              possible_cause
#> 1                     Coordinate cleaning removed records with invalid, missing, or suspicious coordinates.
#> 2          Coordinates may fall outside the overlay extent, offshore, near boundaries, or in the wrong CRS.
#> 3               Provider service, URL, dependency, cache, or licence/opt-in requirement may be unavailable.
#> 4        Raster file may be large, unavailable, in an unexpected format, or incompatible with local memory.
#> 5                   Marine species may have been submitted through a terrestrial/country-filtered workflow.
#> 6                      The pipeline may have written split outputs by spatial unit or failed before export.
#> 7 Large GBIF downloads, high-resolution overlays, raster extraction, or in-memory storage may be too heavy.
#>                                                                               first_check
#> 1 Inspect n_total, n_cleaned, n_thinned, fail_stage, and fail_reason in gbif_summary.csv.
#> 2    Check coordinates, overlay choice, CRS, coastal/offshore points, and unmatched rows.
#> 3                Try the provider loader directly with a small country or overlay subset.
#> 4                              Test one species-region pair and inspect raster_cache_dir.
#> 5 Use process_gbif_marine_pipeline() and a marine overlay such as EEZ, LME, MEOW, or IHO.
#> 6           Inspect output_dir, output_file, status, and fail_reason in gbif_summary.csv.
#> 7               Set store_in_memory = FALSE and add overlays or raster context in stages.

For a failed or confusing run, reduce the workflow to the smallest useful test: one taxon, one spatial framework and one output directory. Once that works, add the next overlay, raster layer or evidence setting and check the summary again.

Provider caches and citations

biofetchR can download or cache third-party spatial resources for use in an analysis. These resources remain provider datasets; they are not converted into package data. For any workflow that will be reported, shared or rerun, keep enough information to identify the exact spatial products used.

Useful details to save include:

Provider citation needs depend on the selected layers. Common examples include:

Use the citation supplied by the provider page, download record or metadata file for the exact product used. This is especially important for datasets that are versioned, updated regularly, licensed separately, or distributed through more than one provider route.

Final checklist

Before using spatially enriched occurrence outputs, check that the spatial processing decisions are clear and that the exported files can be traced back to their inputs and provider layers.

Useful checks include:

A good spatial workflow should be easy to inspect later. The key files are the input table, gbif_summary.csv, exported occurrence files, provider-cache notes and any citation or licence information needed for the spatial layers used.

mirror server hosted at Truenetwork, Russian Federation.