--- title: "Using custom gazetteers" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Using custom gazetteers} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6.5, fig.height = 4, dev = "svglite", fig.ext = "svg", dev.args = list(pointsize = 14)) has_sf <- requireNamespace("sf", quietly = TRUE) has_cc <- requireNamespace("CoordinateCleaner", quietly = TRUE) ``` The tests use three kinds of reference data: - Points marking capitals, centroids, institutions, the GBIF headquarters and artificial hotspots. - Polygons outlining land, urban areas, countries and species ranges. - A table of genus habitats, used together with land polygons by the two habitat tests. OccurrenceCleaner ships fixed versions of these reference layers, except for species ranges, which users supply. Set `reference = "CoordinateCleaner"` to use the corresponding reference tables from [CoordinateCleaner](https://ropensci.github.io/CoordinateCleaner/). To supply your own layer for an individual test, use the `gazetteers` argument of [`flag_occurrences()`](https://gillescolling.com/OccurrenceCleaner/reference/flag_occurrences.html). This serves the same purpose as the `ref` arguments of CoordinateCleaner's `cc_*()` functions. This article shows how to filter and extend the shipped point tables, build polygon layers with [sf](https://r-spatial.github.io/sf/), supply a genus habitat table, and inspect how custom reference layers are recorded with the flags. Installation is covered in [Get started](OccurrenceCleaner.html); the polygon examples also require sf. ## Included reference points `reference_points()` returns the points of a proximity test with their attributes: ```{r points} library(OccurrenceCleaner) caps <- reference_points("capital") head(caps) table(caps$class) inst <- reference_points("institution") head(sort(table(inst$type), decreasing = TRUE)) ``` The institution table joins the GBIF registry of scientific collections (GRSciColl), the botanical gardens, natural history museums, zoos, aquaria and herbaria of Wikidata, the museums of Wikipedia's lists of natural history museums (the English list and six other language editions), the gardens of Botanic Gardens Conservation International (BGCI), three GeoNames places, the herbaria of Index Herbariorum that GRSciColl does not place, the publishing organisations of the GBIF registry, and every further institution of [CoordinateCleaner](https://ropensci.github.io/CoordinateCleaner/)'s list, placed by its name from Wikidata, OpenStreetMap or GeoNames. `source` says which registry each point comes from, `type` what kind of institution it is, and `master_source`, for GRSciColl entries, where GBIF took the entry from (Index Herbariorum, GRSciColl itself, or the GBIF registry); for the Wikipedia museums that the list places only by town and for the rows of CoordinateCleaner's list, it says where the coordinates come from. `extension` is `TRUE` for the points from sources [CoordinateCleaner](https://ropensci.github.io/CoordinateCleaner/)'s own table does not draw on (Wikidata beyond the museums of the English list, and the other language editions); the rows with `extension` `FALSE`, passed as the institution gazetteer, restrict the test to CoordinateCleaner's sources. ## A filtered point gazetteer A point layer is any data frame with longitude and latitude columns (`lon`/`lat`, `longitude`/`latitude` or `decimalLongitude`/`decimalLatitude`, in any case). A plant study, for example, may want to flag records at herbaria only, since a specimen georeferenced to a zoological museum is unlikely in a plant table, and use a wider buffer for them: ```{r herbaria} herbaria <- inst[inst$type %in% "Herbarium", ] nrow(herbaria) near <- rbind(herbaria[1:3, ], inst[inst$type %in% "Museum", ][1:3, ]) x <- data.frame(species = "Quercus robur", decimalLongitude = near$lon + 0.005, decimalLatitude = near$lat) f_all <- flag_occurrences(x, tests = "institution", buffers = list(institution = 1000)) f_herb <- flag_occurrences(x, tests = "institution", buffers = list(institution = 1000), gazetteers = list(institution = herbaria)) data.frame(near = near$name, type = near$type, all = f_all$flag_institution, herbaria_only = f_herb$flag_institution) ``` Each record lies about 0.005 degrees (a few hundred metres) east of an institution: three herbaria, then three museums. The full table flags all six at a 1 km buffer; the herbarium table flags only `r sum(f_herb$flag_institution)`, the records next to the herbaria. ## An extended point gazetteer The two capital tables differ: the shipped table has the Factbook capitals unrounded and adds Natural Earth's alternative and official capitals, while [CoordinateCleaner](https://ropensci.github.io/CoordinateCleaner/)'s has one point per country, rounded to 0.01 degrees. Their union flags records near either: ```{r union, eval = has_cc} both <- rbind(caps[, c("lon", "lat")], reference_points("capital", reference = "CoordinateCleaner")) nrow(both) y <- data.frame(city = c("Cape Town", "Naypyidaw", "Vienna"), decimalLongitude = c(18.42, 96.13, 16.37), decimalLatitude = c(-33.92, 19.75, 48.21)) data.frame(city = y$city, CoordinateCleaner = flag_occurrences(y, tests = "capital", reference = "CoordinateCleaner")$flag_capital, union = flag_occurrences(y, tests = "capital", gazetteers = list(capital = both))$flag_capital) ``` Other point sets fit the same way: regional capitals for a study in which collections were georeferenced to provincial seats, the field stations of a survey, or the coordinates a particular data provider uses as a default. The institution table is one of them. The shipped table holds every institution of CoordinateCleaner's list, at the coordinates of its sources. CoordinateCleaner placed nearly every institution by its own geocoding, so its point can mark a different spot of the same institution. With CoordinateCleaner installed, the union adds its points to the shipped ones at run time: ```{r inst-union, eval = has_cc} inst_both <- rbind(inst[, c("lon", "lat")], reference_points("institution", reference = "CoordinateCleaner")) nrow(inst_both) ``` `gazetteers = list(institution = inst_both)` then flags records near an institution of either table. The Artificial Hotspot Occurrence Inventory (Park et al. 2023) behind the `aohi` test marks each point as a real observation site (`determination = TRUE`, for example a birding reserve) or an artificial aggregate of records (`FALSE`, for example a georeferencing default). A study that keeps records from popular sites can flag only the aggregates: ```{r aohi} ah <- reference_points("aohi") table(ah$taxon, ah$determination) artificial <- ah[!ah$determination, ] v <- data.frame(decimalLongitude = ah$lon[1:4] + 0.02, decimalLatitude = ah$lat[1:4]) data.frame(point = ah$notes[1:4], site = ah$determination[1:4], all = flag_occurrences(v, tests = "aohi")$flag_aohi, aggregates = flag_occurrences(v, tests = "aohi", gazetteers = list(aohi = artificial))$flag_aohi) ``` The four records lie about 2 km from the first four bird hotspots, all of them birding sites, so the full inventory flags them and the aggregates-only table does not. ## Polygon layers Polygon layers are `sf` or `sfc` objects of POLYGON or MULTIPOLYGON geometries in longitude and latitude (EPSG:4326). For `sea`, the polygons are the land: a record outside every polygon is flagged. That makes `sea` with a custom layer a test for "outside the study area": ```{r study-area, eval = has_sf} area <- sf::st_sfc(sf::st_polygon(list(rbind(c(9.5, 46.4), c(17.2, 46.4), c(17.2, 49.0), c(9.5, 49.0), c(9.5, 46.4)))), crs = 4326) z <- data.frame(decimalLongitude = c(16.37, 11.40, 8.54, 19.04), decimalLatitude = c(48.21, 47.26, 47.37, 47.50)) flag_occurrences(z, tests = "sea", gazetteers = list(sea = area))$flag_sea ``` The first two points (Vienna, Innsbruck) lie inside the box, the last two (Zurich, Budapest) outside it. The same applies to `urban`, whose polygons are the areas to flag, for example a city's built-up area at a finer scale than Natural Earth's. For `country_mismatch` the polygons carry codes. Any of the columns `iso_a2`, `iso_a3`, `adm0_a3` or `code` (in any case) is read, and the record's code is matched against all of them. Two territories with their own codes: ```{r countries, eval = has_sf} terr <- sf::st_sf( code = c("XA", "XB"), geometry = sf::st_sfc(sf::st_polygon(list(rbind(c(0, 0), c(10, 0), c(10, 10), c(0, 10), c(0, 0)))), sf::st_polygon(list(rbind(c(10, 0), c(20, 0), c(20, 10), c(10, 10), c(10, 0)))), crs = 4326)) w <- data.frame(decimalLongitude = c(5, 15, 15, 25), decimalLatitude = 5, countryCode = c("XA", "XB", "XA", "XA")) flag_occurrences(w, tests = "country_mismatch", gazetteers = list(country_mismatch = terr))$flag_country_mismatch ``` The third record is in territory XB with code XA, the fourth outside both territories, so both are flagged. ## Species ranges The `range` test flags records outside the range polygons of their species, such as IUCN range maps. The layer is an `sf` object with one or more polygons per species and a species column, named like the record's species column or `species`, `binomial` or `sci_name`. Each species gets an index of its own polygons, so overlapping ranges of different species do not interfere; species without a polygon are not tested: ```{r range, eval = has_sf} ranges <- sf::st_sf( species = c("Bufo bufo", "Bufo bufo", "Lacerta viridis"), geometry = sf::st_sfc(sf::st_polygon(list(rbind(c(-10, 36), c(30, 36), c(30, 60), c(-10, 60), c(-10, 36)))), sf::st_polygon(list(rbind(c(30, 40), c(60, 40), c(60, 60), c(30, 60), c(30, 40)))), sf::st_polygon(list(rbind(c(5, 40), c(30, 40), c(30, 52), c(5, 52), c(5, 40)))), crs = 4326)) r <- data.frame(species = c("Bufo bufo", "Bufo bufo", "Bufo bufo", "Lacerta viridis", "Lacerta viridis", "Rana arvalis"), decimalLongitude = c(10, 45, -20, 15, 15, 100), decimalLatitude = c(50, 50, 50, 45, 55, 10)) data.frame(r, outside = flag_occurrences(r, tests = "range", gazetteers = list(range = ranges))$flag_range, outside_100km = flag_occurrences(r, tests = "range", gazetteers = list(range = ranges), buffers = list(range = 100000))$flag_range) ``` The toad at 45 E lies in the second of its two polygons and is kept; the one at 20 W lies in the Atlantic, outside both. The lizard at 55 N is 3 degrees north of its range and stays flagged with a 100 km buffer. *Rana arvalis* has no polygon and is not tested. ## A genus habitat table The habitat tests `sea_nonmarine` and `land_marine` look up the record's genus in `reference_habitat()`, built from the Interim Register of Marine and Nonmarine Genera. A table with a `genus` column, a logical `marine` column (`TRUE` only marine, `FALSE` not marine) and optionally a `kingdom` column replaces it, for example to class a genus the register leaves open: ```{r habitat} h <- reference_habitat() h[h$genus %in% c("Larus", "Zostera"), ] mine <- rbind(h, data.frame(kingdom = "Plantae", genus = "Zostera", marine = TRUE)) x <- data.frame(kingdom = "Plantae", genus = "Zostera", decimalLongitude = 16.37, decimalLatitude = 48.21) flag_occurrences(x, tests = "land_marine")$flag_land_marine flag_occurrences(x, tests = "land_marine", gazetteers = list(habitat = mine))$flag_land_marine ``` The register marks seagrass (*Zostera*) as marine but leaves its freshwater and terrestrial flags empty, so the shipped table does not class it and its records are not tested; classed as only marine, a record of it in Vienna is flagged. ## Recording the layer A custom layer is indexed once per call, and its size and an MD5 checksum of its coordinates are written into the reference label of the output: ```{r label} f <- flag_occurrences(x, tests = "institution", gazetteers = list(institution = herbaria)) attr(f, "occ_reference") ``` The checksum identifies the exact table, so two analyses that report the same label used the same points. For streamed input the label is stored in the output directory, and a rerun into that directory with a different custom layer stops with an error instead of mixing flags from two gazetteers. ## Comparison with CoordinateCleaner ```{r cc-map, echo = FALSE} knitr::kable(data.frame( CoordinateCleaner = c("cc_cap(ref = ...)", "cc_cen(ref = ...)", "cc_inst(ref = ...)", "cc_sea(ref = ...)", "cc_urb(ref = ...)", "cc_coun(ref = ..., ref_col = ...)", "cc_iucn(range = ...)", "cc_aohi(taxa = ...)", "(none)"), OccurrenceCleaner = c("gazetteers = list(capital = points)", "gazetteers = list(centroid = points)", "gazetteers = list(institution = points)", "gazetteers = list(sea = polygons)", "gazetteers = list(urban = polygons)", "gazetteers = list(country_mismatch = polygons)", "gazetteers = list(range = polygons)", "gazetteers = list(aohi = subset of reference_points(\"aohi\"))", "gazetteers = list(habitat = genera)"), check.names = FALSE)) ``` [CoordinateCleaner](https://ropensci.github.io/CoordinateCleaner/)'s point references need its own column names (`capital.lon`, `centroid.lon`, `decimalLongitude`); here the longitude and latitude columns are found by name in any of the common spellings. Its `cc_coun()` reads one code column chosen with `ref_col`; here every code column present is matched. With the same reference points and polygons the two packages give the same flags, up to the buffer-polygon band described in [Compatibility with CoordinateCleaner](coordinatecleaner.html); the package tests check this for `institution`, `sea` and `range`.