R/aptg-internal.R

In aptg: Automatic Phylogenetic Tree Generator

# Internal helpers for aptg. Not exported.

# ---------------------------------------------------------------------------
# Name resolution: map names to OTT ids and keep only those that are actually
# present in the Open Tree synthetic tree.
#
# Returns a list:
#   matched : the tnrs_match_names data frame
#   ott     : numeric OTT ids, aligned row-for-row with `matched` (NA = no match)
#   keep    : logical, TRUE for rows resolved AND present in the synthetic tree
#   dropped : search strings for rows that were not kept
# ---------------------------------------------------------------------------
.aptg_resolve <- function(species, verbose = TRUE) {
  species <- unique(stats::na.omit(as.character(species)))
  species <- species[nzchar(trimws(species))]
  if (length(species) < 2L) {
    stop("Need at least two distinct, non-empty names to build a tree.",
         call. = FALSE)
  }

  matched <- rotl::tnrs_match_names(species)

  # ott_id() may return a list; coerce to a positional numeric vector aligned
  # row-for-row with `matched` (NA where a name did not resolve).
  ott <- unlist(rotl::ott_id(matched), use.names = FALSE)
  resolved <- !is.na(ott)

  in_tree <- logical(length(ott))
  if (any(resolved)) {
    in_tree[resolved] <- rotl::is_in_tree(ott[resolved])
  }

  keep <- resolved & in_tree
  list(matched = matched,
       ott = ott,
       keep = keep,
       dropped = matched$search_string[!keep])
}

# ---------------------------------------------------------------------------
# Phylum of each OTT id (one batched API call). Returns a character vector
# aligned to `ott_ids`; NA where the Open Tree lineage carries no phylum rank.
# The phylum rank is the safeguard boundary: taxa are grouped by it so that no
# induced subtree ever spans above phylum.
# ---------------------------------------------------------------------------
.aptg_phylum_of <- function(ott_ids) {
  info <- rotl::taxonomy_taxon_info(ott_ids, include_lineage = TRUE)
  lin <- rotl::tax_lineage(info)
  vapply(lin, function(df) {
    hit <- df$name[df$rank == "phylum"]
    if (length(hit)) hit[[1]] else NA_character_
  }, character(1))
}

# ---------------------------------------------------------------------------
# Build an Open Tree topology subtree for a set of OTT ids (all assumed to be
# within one phylum). Open Tree subtrees carry no branch lengths, so unit edges
# are assigned and the returned distance matrix counts edges between tips.
# ---------------------------------------------------------------------------
.aptg_otl_tree <- function(ott_ids, verbose = TRUE) {
  tree <- rotl::tol_induced_subtree(ott_ids = ott_ids, label_format = "name")
  tree$tip.label <- rotl::strip_ott_ids(tree$tip.label,
                                        remove_underscores = TRUE)
  if (is.null(tree$edge.length)) {
    if (verbose) {
      message("Open Tree subtree has no branch lengths; assigning unit ",
              "edges (distances are edge counts, not time).")
    }
    tree <- ape::compute.brlen(tree, 1)
  }
  list(tree = tree, dist = stats::cophenetic(tree))
}

# ---------------------------------------------------------------------------
# Dated tree for ray-finned fishes from the Fish Tree of Life (fishtree).
# fishtree is in Suggests, so guard with requireNamespace(). All of
# Actinopterygii sits within phylum Chordata, so the phylum safeguard is
# satisfied by construction. The chronogram carries real branch lengths, so
# the distance matrix is patristic (time).
# ---------------------------------------------------------------------------
.aptg_fish_tree <- function(species, verbose = TRUE) {
  if (!requireNamespace("fishtree", quietly = TRUE)) {
    stop("source = \"fish\" requires the 'fishtree' package. ",
         "Install it with install.packages(\"fishtree\").", call. = FALSE)
  }
  species <- unique(stats::na.omit(as.character(species)))
  species <- species[nzchar(trimws(species))]
  if (length(species) < 2L) {
    stop("Need at least two distinct species for the fish backend.",
         call. = FALSE)
  }

  tree <- fishtree::fishtree_phylogeny(species = species)

  present <- gsub("_", " ", tree$tip.label)
  dropped <- setdiff(species, present)
  if (verbose && length(dropped)) {
    message("Not in the Fish Tree of Life (dropped): ",
            paste(dropped, collapse = ", "))
  }
  if (length(tree$tip.label) < 2L) {
    stop("Fewer than two of the supplied species are in the Fish Tree of ",
         "Life; cannot build a tree.", call. = FALSE)
  }

  list(tree = tree, dist = stats::cophenetic(tree), dropped = dropped)
}

# ---------------------------------------------------------------------------
# Plot a named list of trees, one panel per phylum, restoring par() on exit.
# ---------------------------------------------------------------------------
.aptg_plot_trees <- function(trees) {
  k <- length(trees)
  if (!k) return(invisible(NULL))
  nc <- ceiling(sqrt(k))
  nr <- ceiling(k / nc)
  op <- graphics::par(mfrow = c(nr, nc))
  on.exit(graphics::par(op), add = TRUE)
  for (nm in names(trees)) {
    ape::plot.phylo(trees[[nm]]$tree)
    graphics::title(main = nm)
  }
  invisible(NULL)
}

# ===========================================================================
# Geographic helpers (region.tree). All GBIF/geosphere access is via
# requireNamespace(), since rgbif and geosphere live in Suggests.
# ===========================================================================

# Normalise a place name to lowercase ASCII letters only, so "Québec",
# "quebec" and "QC " all collapse to the same key.
.aptg_norm <- function(x) {
  x <- tolower(trimws(x))
  x2 <- iconv(x, to = "ASCII//TRANSLIT")
  if (is.na(x2)) x2 <- x
  gsub("[^a-z]", "", x2)
}

# GADM level-1 GIDs for the 13 Canadian provinces and territories, keyed by
# normalised names and common aliases (postal codes etc.). Based on GADM's
# alphabetical level-1 ordering, which GBIF uses for its gadmGid filter.
# `qubec` is included as a safety alias in case iconv transliteration of the
# accent is unavailable on a given platform.
.aptg_canada_gadm <- c(
  alberta = "CAN.1_1",  ab = "CAN.1_1",
  britishcolumbia = "CAN.2_1", bc = "CAN.2_1",
  manitoba = "CAN.3_1", mb = "CAN.3_1",
  newbrunswick = "CAN.4_1", nb = "CAN.4_1",
  newfoundlandandlabrador = "CAN.5_1", newfoundland = "CAN.5_1",
  labrador = "CAN.5_1", nl = "CAN.5_1",
  northwestterritories = "CAN.6_1", nwt = "CAN.6_1", nt = "CAN.6_1",
  novascotia = "CAN.7_1", ns = "CAN.7_1",
  nunavut = "CAN.8_1", nu = "CAN.8_1",
  ontario = "CAN.9_1", on = "CAN.9_1",
  princeedwardisland = "CAN.10_1", pei = "CAN.10_1", pe = "CAN.10_1",
  quebec = "CAN.11_1", qubec = "CAN.11_1", qc = "CAN.11_1",
  saskatchewan = "CAN.12_1", sk = "CAN.12_1",
  yukon = "CAN.13_1", yukonterritory = "CAN.13_1", yt = "CAN.13_1"
)

# Resolve an admin area to a GADM GID. A raw `gadm` GID wins (works anywhere);
# otherwise a Canadian province/territory name/alias is looked up.
.aptg_gadm_gid <- function(province = NULL, gadm = NULL) {
  if (!is.null(gadm) && nzchar(gadm)) return(gadm)
  if (is.null(province) || !nzchar(trimws(province))) {
    stop("Provide either `gadm` (a GADM GID) or `province` (a name).",
         call. = FALSE)
  }
  gid <- unname(.aptg_canada_gadm[.aptg_norm(province)])
  if (is.na(gid)) {
    stop("Unknown Canadian province/territory: '", province, "'. Supported: ",
         "Alberta, British Columbia, Manitoba, New Brunswick, Newfoundland ",
         "and Labrador, Northwest Territories, Nova Scotia, Nunavut, Ontario, ",
         "Prince Edward Island, Quebec, Saskatchewan, Yukon (or their postal ",
         "codes). For areas outside Canada, pass a GADM GID via `gadm`.",
         call. = FALSE)
  }
  gid
}

# Build a Well-Known-Text polygon approximating a geodesic circle of
# `radius_km` around (lat, lon). GBIF requires the exterior ring to be
# counter-clockwise (right-hand rule); geosphere::destPoint walks clockwise as
# the bearing increases, so the ring is reversed. The ring is closed (first
# vertex repeated).
.aptg_wkt_circle <- function(lat, lon, radius_km, n = 64L) {
  if (!requireNamespace("geosphere", quietly = TRUE)) {
    stop("Radius mode requires the 'geosphere' package. ",
         "install.packages(\"geosphere\").", call. = FALSE)
  }
  if (!is.numeric(lat) || !is.numeric(lon) || !is.numeric(radius_km) ||
      abs(lat) > 90 || abs(lon) > 180 || radius_km <= 0) {
    stop("lat in [-90, 90], lon in [-180, 180] and radius_km > 0 required.",
         call. = FALSE)
  }
  bearings <- seq(0, 360, length.out = n + 1L)
  pts <- geosphere::destPoint(c(lon, lat), b = bearings, d = radius_km * 1000)
  ring <- pts[nrow(pts):1L, , drop = FALSE]          # reverse -> CCW
  coords <- paste(sprintf("%.6f %.6f", ring[, 1], ring[, 2]), collapse = ", ")
  sprintf("POLYGON((%s))", coords)
}

# Resolve a taxon name to a GBIF backbone key.
.aptg_gbif_taxonkey <- function(taxon) {
  bb <- rgbif::name_backbone(name = taxon)
  key <- bb$usageKey
  if (is.null(key) || length(key) != 1L || is.na(key)) {
    stop("GBIF could not match the taxon '", taxon, "'.", call. = FALSE)
  }
  key
}

# Distinct species (canonical names) for a taxon within a geographic filter.
# Uses occurrence faceting on speciesKey (no bulk download, no credentials),
# then resolves each key to a name. Only georeferenced records without
# geospatial issues are counted.
.aptg_gbif_species <- function(taxon_key, geometry = NULL, gadm_gid = NULL,
                               max_species = 2000L, verbose = TRUE) {
  args <- list(taxonKey = taxon_key, limit = 0,
               facet = "speciesKey", facetLimit = max_species,
               hasCoordinate = TRUE, hasGeospatialIssue = FALSE)
  if (!is.null(geometry)) args$geometry <- geometry
  if (!is.null(gadm_gid)) args$gadmGid <- gadm_gid

  sr <- do.call(rgbif::occ_search, args)
  fac <- sr$facets$speciesKey
  if (is.null(fac) || !NROW(fac)) {
    stop("GBIF returned no species for this taxon in the requested area.",
         call. = FALSE)
  }
  keys <- as.integer(fac$name)
  if (verbose) {
    message("GBIF returned ", length(keys), " distinct species",
            if (length(keys) >= max_species)
              paste0(" (capped at max_species = ", max_species, ")") else "",
            "; resolving names.")
  }

  nm <- vapply(keys, function(k) {
    d <- tryCatch(rgbif::name_usage(key = k)$data, error = function(e) NULL)
    val <- if (is.null(d)) NULL else d$canonicalName
    if (is.null(val) || !length(val) || is.na(val[[1]])) NA_character_ else val[[1]]
  }, character(1))

  unique(stats::na.omit(nm))
}