R/region.tree.R

In aptg: Automatic Phylogenetic Tree Generator

#' Phylogenetic tree(s) for a taxon within a geographic area
#'
#' Retrieves the species of a given clade recorded in an area from the Global
#' Biodiversity Information Facility (GBIF, via \pkg{rgbif}), then builds the
#' tree(s) with \code{\link{taxa.tree}} (inheriting the per-phylum split and
#' the \code{source} option). The area is specified in one of two ways:
#'
#' \itemize{
#'   \item \strong{Radius mode}: a circle of \code{radius_km} around
#'     (\code{lat}, \code{lon}).
#'   \item \strong{Admin mode}: an administrative area given by \code{gadm}
#'     (a GADM GID, works worldwide) or \code{province} (a Canadian province
#'     or territory name/abbreviation; all 13 are supported).
#' }
#'
#' @param taxon Required. A single clade name to constrain the search
#'   (e.g. \code{"Aves"}, \code{"Mammalia"}, \code{"Salmonidae"}). Without it,
#'   a populated area returns an unusable number of species.
#' @param lat,lon,radius_km Radius mode: latitude, longitude (decimal degrees)
#'   and radius in kilometres. Supply all three.
#' @param province Admin mode: a Canadian province/territory name or postal
#'   code (e.g. \code{"Quebec"}, \code{"QC"}, \code{"Nunavut"}, \code{"BC"}).
#' @param gadm Admin mode: a GADM GID (e.g. \code{"CAN.11_1"} for Quebec, or
#'   any GID worldwide). Takes precedence over \code{province}.
#' @param source Passed to \code{\link{taxa.tree}}: \code{"otl"} (default) or
#'   \code{"fish"}.
#' @param max_species Cap on the number of distinct species retrieved from
#'   GBIF (facet limit). Default 2000.
#' @param plot,verbose Passed to \code{\link{taxa.tree}}.
#'
#' @return Invisibly, the \code{\link{taxa.tree}} result
#'   (\code{list(trees, unmatched)}) with an added \code{species} element: the
#'   distinct species names GBIF returned for the area.
#'
#' @section Requirements:
#'   Needs the suggested \pkg{rgbif} package (both modes) and \pkg{geosphere}
#'   (radius mode only). Both hit the network.
#'
#' @section Caveats:
#'   GBIF returns \emph{occurrence records}, not a curated checklist. The
#'   result is presence-only and sampling-biased (some groups and places are
#'   far better recorded than others), so it reflects what has been observed
#'   and digitised, not a definitive species inventory. Only records with
#'   coordinates and without flagged geospatial issues are counted, but
#'   coordinate precision still varies. Treat the species list as a starting
#'   point, not ground truth.
#'
#' @section Canadian areas:
#'   The \code{province} lookup covers Alberta, British Columbia, Manitoba,
#'   New Brunswick, Newfoundland and Labrador, Northwest Territories, Nova
#'   Scotia, Nunavut, Ontario, Prince Edward Island, Quebec, Saskatchewan and
#'   Yukon (and their postal codes). The GIDs follow GADM's alphabetical
#'   level-1 ordering; the resolved GID is printed when \code{verbose = TRUE}.
#'   If GBIF's GADM version ever disagrees, pass the GID directly via
#'   \code{gadm} (browse them at
#'   \code{https://api.gbif.org/v1/geocode/gadm/browse/CAN}).
#'
#' @seealso \code{\link{taxa.tree}}, \code{\link{downto.tree}}
#'
#' @examples
#' \dontrun{
#' # Birds within 50 km of Montreal:
#' region.tree("Aves", lat = 45.50, lon = -73.57, radius_km = 50)
#'
#' # Mammals of Nunavut (admin mode, Canadian territory):
#' region.tree("Mammalia", province = "Nunavut")
#'
#' # Freshwater fishes of Quebec, as a dated tree:
#' region.tree("Actinopterygii", province = "QC", source = "fish")
#'
#' # Anywhere in the world via a raw GADM GID:
#' region.tree("Reptilia", gadm = "USA.5_1")
#' }
#' @export
region.tree <- function(taxon,
                        lat = NULL, lon = NULL, radius_km = NULL,
                        province = NULL, gadm = NULL,
                        source = c("otl", "fish"),
                        max_species = 2000L, plot = TRUE, verbose = TRUE) {
  source <- match.arg(source)
  if (missing(taxon) || length(taxon) != 1L || !is.character(taxon) ||
      !nzchar(trimws(taxon))) {
    stop("`taxon` (a single clade name to constrain the search) is required.",
         call. = FALSE)
  }
  radius_mode <- !is.null(lat) || !is.null(lon) || !is.null(radius_km)
  admin_mode <- !is.null(province) || !is.null(gadm)
  if (radius_mode && admin_mode) {
    stop("Use either radius mode (lat/lon/radius_km) or admin mode ",
         "(province/gadm), not both.", call. = FALSE)
  }
  if (!radius_mode && !admin_mode) {
    stop("Provide a location: lat/lon/radius_km, or province/gadm.",
         call. = FALSE)
  }
  if (radius_mode && (is.null(lat) || is.null(lon) || is.null(radius_km))) {
    stop("Radius mode needs all of lat, lon and radius_km.", call. = FALSE)
  }

  if (!requireNamespace("rgbif", quietly = TRUE)) {
    stop("region.tree() requires the 'rgbif' package. ",
         "install.packages(\"rgbif\").", call. = FALSE)
  }

  geometry <- NULL
  gadm_gid <- NULL
  if (radius_mode) {
    geometry <- .aptg_wkt_circle(lat, lon, radius_km)
    if (verbose) {
      message("Querying GBIF within ", radius_km, " km of (", lat, ", ",
              lon, ").")
    }
  } else {
    gadm_gid <- .aptg_gadm_gid(province = province, gadm = gadm)
    if (verbose) {
      message("Querying GBIF for GADM area ", gadm_gid,
              if (!is.null(province)) paste0(" (", province, ")") else "", ".")
    }
  }

  taxon_key <- .aptg_gbif_taxonkey(taxon)
  species <- .aptg_gbif_species(taxon_key, geometry = geometry,
                                gadm_gid = gadm_gid,
                                max_species = max_species, verbose = verbose)
  if (length(species) < 2L) {
    stop("Fewer than two species returned for this area; cannot build a tree.",
         call. = FALSE)
  }
  if (verbose) {
    message("Building tree(s) from ", length(species), " species.")
  }

  res <- taxa.tree(species, source = source, plot = plot, verbose = verbose)
  res$species <- species
  invisible(res)
}