Skip to contents

What the package does

mappingAS — Mapping Area of Species. Geographic range metrics (EOO / AOO), land-cover conversion and fire, and protected-area overlap for extinction-risk screening — from raw occurrence points to an interactive app and a written assessment report.

mappingAS is an R package, with an accompanying Shiny application, for screening species against Criterion B of the IUCN Red List. Starting from a set of occurrence points, it estimates each species’ geographic range (EOO / AOO), measures how much of that range has been converted to anthropic land cover and how much has burned, quantifies the overlap with protected areas (global WDPA), and packages everything for mapping, inspection, export and reporting. Land-cover conversion is driven by MapBiomas across South America and by the global Esri / Impact Observatory (Sentinel-2) 10 m land cover elsewhere, so a species anywhere on Earth can be screened with the same workflow — no Google Earth Engine account required.

Screening only. Two categories are produced. provisional_cat rests only on the EOO/AOO size thresholds of Criterion B. category_B combines those thresholds with the sub-criteria the data can support: the estimated number of locations (a), a continuing decline (b) inferred from a significant recent decline of natural habitat or assumed when that cannot be measured, and extreme fluctuation (c) when you supply it. Neither replaces a formal IUCN assessment: the assessor reviews the evidence and decides.

1. Import the points

read_occurrences() accepts .csv/.tsv/.txt, .xlsx/.xls and spatial files (.shp, .gpkg, .geojson, .zip) and tries to detect the species, longitude and latitude columns automatically.

ex  <- system.file("extdata", "example_occurrences.csv", package = "mappingAS")
occ <- read_occurrences(ex)

# if detection fails, name the columns explicitly:
# occ <- read_occurrences("my_data.csv",
#                         species_col = "species",
#                         lon_col = "longitude", lat_col = "latitude")

The result is an sf object of points in WGS84 with a species column.

2. Assess the species

assess_species() runs the whole chain per species: EOO, AOO, subpopulations, locations, conversion, the habitat trend and the Criterion B categories.

res <- assess_species(
  occ,
  year = 2024,          # MapBiomas year for the land-cover snapshot
  collection = 10,      # MapBiomas collection
  backend = "local"     # "local" (no GEE account) or "gee"
)

res            # compact print
res$summary    # one row per species

Key res$summary columns: eoo_km2, aoo_km2, aoo_cells, n_subpop, n_locations, eoo_converted_pct/eoo_natural_pct, aoo_converted_pct/aoo_natural_pct, provisional_cat, category_B, criterion_B_code, decline_basis, habitat_trend_pct_yr and habitat_trend_p. The conversion index uses a terrestrial denominator: converted = anthropic / (anthropic + natural) * 100 (water and “not observed” are excluded by default, change with water_in_denominator = TRUE).

How the main numbers are obtained:

  • EOO and AOO. Following the IUCN Guidelines (sect. 4.10), when the EOO is smaller than the AOO it is set equal to the AOO before the thresholds are applied. The AOO is the minimum count of occupied 2 × 2 km cells over 100 grid offsets laid out on a fixed lattice, so the same data always give the same AOO and number of locations (no random seed involved).
  • Continuing decline (b). When decline is not given, the natural area of the range is measured at six years spanning the last decline_years (default 10). A decline in habitat, b(iii), is inferred only when the linear trend is negative and significant (p < 0.05). The slope (% of the natural area per year) and its p-value are reported. A share of converted land with no ongoing loss is not treated as a decline. When the trend cannot be measured (for example with a single-year src file) the decline is assumed (assume_decline = TRUE, the default) or left unmet (FALSE).
  • Locations and protected areas. With protected = TRUE, each protected area holding records counts as one location. Nested or overlapping designations, such as a national park inside an environmental protection area, are counted once.
  • Code notation. criterion_B_code follows the IUCN style, repeating the sub-criteria for each axis, for example EN B1ab(iii)+2ab(iii).
# Stricter screening: a decline must be measured or documented
res_strict <- assess_species(occ, year = 2024, assume_decline = FALSE)

# Longer window for the habitat trend
res_20y <- assess_species(occ, year = 2024, decline_years = 20)

# Expert input always wins over the inference
res_doc <- assess_species(occ, year = 2024,
                          decline = c("Prepusa connata Gardner" = TRUE,
                                      "Worsleya procera (Lem.) Traub" = FALSE))

3. Map and conversion chart

map_species(res)            # leaflet: points + EOO + AOO
plot_conversion(res)        # bars: natural / altered / water / other

4. Per-class composition

ct <- class_table(res)
head(ct)
# species | range (EOO/AOO) | code | class_pt | class_en | group | area_km2 | pct

5. Land-cover time series

# selected years within the first species' EOO
ts <- timeseries_for_species(res, range = "eoo", by = "class",
                             years = c(2000, 2010, 2020, 2024))
plot_timeseries(ts)         # stacked-area chart, official MapBiomas colours

# or over any geometry:
# ts2 <- cover_timeseries(my_geometry, years = c(1990, 2000, 2010, 2020),
#                         by = "class", backend = "local")

Each year is read separately from MapBiomas, so a full annual series can be slow; the minimum reliable step is 5 years.

6. Red List supporting analyses

These assimilate the sRedList workflow (Cazalis et al. 2024) while keeping the focus on Criterion B, computed from the data the package already produces. Most take the range polygon; take it from calc_eoo().

sp1  <- occ[occ$species == occ$species[1], ]
hull <- calc_eoo(sp1)$hull          # EOO polygon
aoo_lo <- calc_aoo(sp1)$area_km2     # occurrence-based (lower) AOO

# Countries of occurrence + biogeographical realm (needs 'rnaturalearth'):
countries_of_occurrence(range = hull, points = sp1)

# Area of Habitat (suitable land cover within the EOO, optionally masked to an
# elevation band; needs land cover, and 'elevatr' or a DEM for elevation):
# (year defaults to the latest one of the product; large ranges are read
# within `max_pixels` by sampling, which keeps fragmented habitat in the area)
aoh <- calc_aoh(hull, year = 2024, elev_min = 0, elev_max = 1500, points = sp1)
aoh$aoh_km2            # Area of Habitat
aoh$aoo_upper_km2      # habitat-based (upper) AOO bound, at the 2 km scale
aoh$point_prevalence  # share of records inside the habitat (a validation metric)

# Suggested elevation band from the occurrences:
elevation_preferences(sp1)

# Lower/upper AOO bounds and the resulting B2 category range (Brooks et al. 2019):
b <- aoo_bounds(aoo_lower_km2 = aoo_lo, aoh_km2 = aoh$aoo_upper_km2,
                eoo_km2 = calc_eoo(sp1)$area_km2)
b$cat_range
iucn_category_B_range(eoo_km2 = calc_eoo(sp1)$area_km2,
                      aoo_lower_km2 = b$aoo_lower, aoo_upper_km2 = b$aoo_upper)

# Severe-fragmentation screening (Santini et al. 2019): cluster occurrences by a
# dispersal / isolation distance and see how much of the population is in small
# subpopulations. It never sets sub-criterion (a) automatically.
assess_fragmentation(sp1, isolation_km = 20)

7. Narrative report

assessment_report() writes an interpretive Criterion B summary as HTML, plain text or a Word .docx file.

cat(assessment_report(res, output = "text"))

# Word document (needs the 'officer' package):
# assessment_report(res, output = "docx", file = "assessment.docx")

8. Export

export_ranges(res, dir = "out")                 # two shapefiles (EOO, AOO) + CSV
export_ranges(res, dir = "out", format = "gpkg") # one GeoPackage (eoo, aoo layers)
export_ranges(res, dir = "out", zip = TRUE)      # everything in one .zip

The attribute table carries species, eoo_km2/aoo_km2, conv_pct, nat_pct, cat_B1/cat_B2, prov_cat, mb_year and mb_coll; the *_classes.csv file gives each class’s area and percentage.

9. Shiny app

The app has tabs for the Map, Add points, Conversion, Classes, Protected areas, Time series, Fire, Countries, Habitat, Fragmentation, and, at the end of the workflow, Assessment, Report and Methods - with CSV, shapefile/GeoPackage and image (PNG/HTML) downloads.

10. MapBiomas backends

  • "local" (default): reads only the window of the national GeoTIFF over the network with GDAL /vsicurl/. No Earth Engine account and no full-mosaic download. For offline/repeated use, point src= at a local GeoTIFF.
  • "gee": computes the per-class areas server-side on Google Earth Engine (recommended for very large ranges). Requires rgee:
# install.packages("rgee"); rgee::ee_install(); rgee::ee_Initialize()
res <- assess_species(occ, backend = "gee")

Caveats

  • The categories are provisional: do not report them as final categories. The supporting analyses (fragmentation, the habitat trend, the AOO bounds) inform the judgement but the assessor decides.
  • The AOO is reproducible but still depends on sampling effort. Treat it as a lower bound and compare it with the habitat-based upper bound (aoo_bounds()).
  • The habitat trend measures natural land cover, not the specific habitat of the taxon. Check it against what the species actually depends on.
  • MapBiomas Fire covers Brazil only. With initiative = "auto" fire is skipped for species outside Brazil.
  • MapBiomas accuracy varies by class, biome and year.
  • For continental ranges, prefer backend = "gee". Ranges that cross the antimeridian (180°) are handled on the short side of the globe.

Citation

Please also cite MapBiomas (the collection used), the IUCN Red List Guidelines for Criterion B, GeoCAT (Bachman et al. 2011, ZooKeys), and, for the supporting analyses, sRedList (Cazalis et al. 2024, Biological Conservation).