R/get_tiles.R

In terrainr: Landscape Visualizations in R and 'Unity'

#' A user-friendly way to get USGS National Map data tiles for an area
#'
#' This function splits the area contained within a bounding box into a set of
#' tiles, and retrieves data from the USGS National map for each tile. As of
#' version 0.5.0, the method for lists has been deprecated.
#'
#' @param data An `sf` or `SpatRaster` object; tiles will be downloaded for the
#' full extent of the provided object.
#' @param output_prefix The file prefix to use when saving tiles.
#' @param side_length The length, in meters, of each side of tiles to download.
#' If \code{NULL}, defaults to the maximum side length permitted by the least
#' permissive service requested.
#' @param resolution How many meters are represented by each pixel? The default
#' value of 1 means that 1 pixel = 1 meter, while a value of 2 means that
#' 1 pixel = 2 meters, and so on.
#' @param services A character vector of services to download data from. Current
#' options include "3DEPElevation", "USGSNAIPPlus", and "nhd". Users can also
#' use short codes to download a specific type of data without specifying the
#' source; current options for short codes include "elevation" (equivalent to
#' "3DEPElevation"), "ortho" (equivalent to "USGSNAIPPlus), and "hydro" ("nhd").
#' Short codes are
#' not guaranteed to refer to the same source across releases. Short codes are
#' converted to their service name and then duplicates are removed, so any given
#' source will only be queried once per tile.
#' @param verbose Logical: should tile retrieval functions run in verbose mode?
#' @param georeference Logical: should tiles be downloaded as PNGs without
#' georeferencing, or should they be downloaded as georeferenced TIFF files?
#' This option does nothing when only elevation data is being downloaded.
#' @param projected Logical: is `data` in a projected coordinate reference
#' system? If `NULL`, the default, inferred from [sf::st_is_longlat].
#' @param ... Additional arguments passed to [hit_national_map_api].
#' These can be used to change default query parameters or as additional options
#' for the National Map services. See below for more details.
#'
#' @section Available Datasources:
#' The following services are currently available
#' (with short codes in parentheses where applicable). See links for API
#' documentation.
#'
# nolint start
#'
#' * [3DEPElevation](https://elevation.nationalmap.gov/arcgis/rest/services/3DEPElevation/ImageServer)
#'   (short code: elevation)
#' * [USGSNAIPPlus](https://imagery.nationalmap.gov/arcgis/rest/services/USGSNAIPPlus/ImageServer/exportImage)
#'   (short code: ortho)
# * [HRO](https://imagery.nationalmap.gov/arcgis/rest/services/HRO/ImageServer)
#' * [USGSNAIPImagery](https://imagery.nationalmap.gov/arcgis/rest/services/USGSNAIPImagery/ImageServer)
#' * [nhd](https://hydro.nationalmap.gov/arcgis/rest/services/nhd/MapServer)
#'   (short code: hydro)
#' * [govunits](https://carto.nationalmap.gov/arcgis/rest/services/govunits/MapServer)
#' * [contours](https://carto.nationalmap.gov/arcgis/rest/services/contours/MapServer)
#' * [geonames](https://carto.nationalmap.gov/arcgis/rest/services/geonames/MapServer)
#' * [NHDPlus_HR](https://hydro.nationalmap.gov/arcgis/rest/services/NHDPlus_HR/MapServer)
#' * [structures](https://carto.nationalmap.gov/arcgis/rest/services/structures/MapServer)
#' * [transportation](https://carto.nationalmap.gov/arcgis/rest/services/transportation/MapServer)
#' * [wbd](https://hydro.nationalmap.gov/arcgis/rest/services/wbd/MapServer)
#'   ("short code": watersheds)
#' * [ecosystems](https://www.usgs.gov/centers/geosciences-and-environmental-change-science-center/science/global-ecosystems)
#' * [USGSTopo](https://basemap.nationalmap.gov/arcgis/rest/services/USGSTopo/MapServer)
#' * [USGSShadedReliefOnly](https://basemap.nationalmap.gov/arcgis/rest/services/USGSShadedReliefOnly/MapServer)
#' * [USGSImageryOnly](https://basemap.nationalmap.gov/arcgis/rest/services/USGSImageryOnly/MapServer)
#' * [USGSHydroCached](https://basemap.nationalmap.gov/arcgis/rest/services/USGSHydroCached/MapServer)
#' * [USGSTNMBlank](https://basemap.nationalmap.gov/arcgis/rest/services/USGSTNMBlank/MapServer)
#'
#'
# nolint end
#'
#' @section Additional Arguments:
#' The `...` argument can be used to pass additional arguments to the
#' National Map API or to edit the hard-coded defaults used by this function.
#' More information on common arguments to change can be found in
#' [hit_national_map_api]. Note that `...` can also be used to change
#' the formats returned by the server, but that doing so while using this
#' function will likely cause the function to error (or corrupt the output
#' data). To download files in different formats, use [hit_national_map_api].
#'
#' @family data retrieval functions
#'
#' @return A list of the same length as the number of unique services requested,
#' containing named vectors of where data files were saved to. Returned
#' invisibly.
#'
#' @examples
#' \dontrun{
#' simulated_data <- data.frame(
#'   id = seq(1, 100, 1),
#'   lat = runif(100, 44.04905, 44.17609),
#'   lng = runif(100, -74.01188, -73.83493)
#' )
#'
#' simulated_data <- sf::st_as_sf(simulated_data, coords = c("lng", "lat"))
#'
#' get_tiles(simulated_data, tempfile())
#' }
#'
#' @rdname get_tiles
#' @md
#' @export
get_tiles <- function(data,
                      output_prefix = tempfile(),
                      side_length = NULL,
                      resolution = 1,
                      services = "elevation",
                      verbose = FALSE,
                      georeference = TRUE,
                      projected = NULL,
                      ...) {
  UseMethod("get_tiles")
}

#' @rdname get_tiles
#' @export
get_tiles.sf <- function(data,
                         output_prefix = tempfile(),
                         side_length = NULL,
                         resolution = 1,
                         services = "elevation",
                         verbose = FALSE,
                         georeference = TRUE,
                         projected = NULL,
                         ...) {
  dots <- list(...)

  if (is.null(projected)) {
    projected <- !sf::st_is_longlat(data)
    if (is.na(projected)) {
      rlang::warn(c(
        "Assuming geographic CRS.",
        i = "Set 'projected' to TRUE if projected."
      ))
      projected <- FALSE
    }
  }

  if (!any(names(dots) == "bboxSR")) {
    bboxSR <- handle_bboxSR(data, projected) # nolint
  } else {
    bboxSR <- dots[["bboxSR"]]
  } # nolint

  if (!any(names(dots) == "imageSR")) {
    imageSR <- bboxSR # nolint
  } else {
    imageSR <- dots[["imageSR"]]
  } # nolint

  data <- sf::st_bbox(data)

  get_tiles_internal(
    data,
    output_prefix = output_prefix,
    side_length = side_length,
    resolution = resolution,
    services = services,
    verbose = verbose,
    georeference = georeference,
    bboxSR = bboxSR,
    imageSR = imageSR,
    projected = projected,
    ...
  )
}

#' @rdname get_tiles
#' @export
get_tiles.sfc <- function(data,
                          output_prefix = tempfile(),
                          side_length = NULL,
                          resolution = 1,
                          services = "elevation",
                          verbose = FALSE,
                          georeference = TRUE,
                          projected = NULL,
                          ...) {
  data <- sf::st_as_sf(data)

  get_tiles(data,
    output_prefix = output_prefix,
    side_length = side_length,
    resolution = resolution,
    services = services,
    verbose = verbose,
    georeference = georeference,
    projected = projected,
    ...
  )
}

#' @rdname get_tiles
#' @export
get_tiles.Raster <- function(data,
                             output_prefix = tempfile(),
                             side_length = NULL,
                             resolution = 1,
                             services = "elevation",
                             verbose = FALSE,
                             georeference = TRUE,
                             projected = NULL,
                             ...) {
  data <- terra::rast(data)
  get_tiles.SpatRaster(data,
                       output_prefix = output_prefix,
                       side_length = side_length,
                       resolution = resolution,
                       services = services,
                       verbose = verbose,
                       georeference = georeference,
                       projected = projected,
                       ...)
}


#' @rdname get_tiles
#' @export
get_tiles.SpatRaster <- function(data,
                               output_prefix = tempfile(),
                               side_length = NULL,
                               resolution = 1,
                               services = "elevation",
                               verbose = FALSE,
                               georeference = TRUE,
                               projected = NULL,
                               ...) {
  dots <- list(...)
  if (is.null(projected)) {
    projected <- !sf::st_is_longlat(data)
    if (is.na(projected)) {
      rlang::warn(c(
        "Assuming geographic CRS.",
        i = "Set 'projected' to TRUE if projected."
      ))
      projected <- FALSE
    }
  }
  if (!any(names(dots) == "bboxSR")) {
    bboxSR <- handle_bboxSR(data, projected) # nolint
  } else {
    bboxSR <- dots[["bboxSR"]]
  } # nolint

  if (!any(names(dots) == "imageSR")) {
    imageSR <- bboxSR # nolint
  } else {
    imageSR <- dots[["imageSR"]]
  } # nolint

  data <- as.vector(terra::ext(data))
  data <- data.frame(
    lng = c(data[[1]], data[[2]]),
    lat = c(data[[3]], data[[4]])
  )
  data <- sf::st_as_sf(data, coords = c("lng", "lat"))
  data <- sf::st_bbox(data)

  get_tiles_internal(
    data,
    output_prefix = output_prefix,
    side_length = side_length,
    resolution = resolution,
    services = services,
    verbose = verbose,
    georeference = georeference,
    projected = projected,
    bboxSR = bboxSR,
    imageSR = imageSR,
    ...
  )
}


#' @rdname get_tiles
#' @export
get_tiles.list <- function(data,
                           output_prefix = tempfile(),
                           side_length = NULL,
                           resolution = 1,
                           services = "elevation",
                           verbose = FALSE,
                           georeference = TRUE,
                           projected = NULL,
                           ...) {
  .Deprecated(
    "get_tiles.list",
    "terrainr",
    msg = paste("'get_tiles.list' is deprecated as of terrainr 0.5.0.",
      "Convert your list to an sf object instead.",
      sep = "\n"
    )
  )
  projected <- FALSE

  bbox <- terrainr_bounding_box(data[[1]], data[[2]])
  bbox <- sf::st_as_sfc(terrainr_st_bbox(bbox))
  get_tiles(
    bbox,
    output_prefix = output_prefix,
    side_length = side_length,
    resolution = resolution,
    services = services,
    verbose = verbose,
    georeference = georeference,
    projected = projected,
    ...
  )
}

get_tiles_internal <- function(data,
                               output_prefix = tempfile(),
                               side_length = NULL,
                               resolution = 1,
                               services = "elevation",
                               verbose = FALSE,
                               georeference = TRUE,
                               projected = NULL,
                               ...) {

  # short codes are assigned as names; we'll cast them into the full name later
  # full names are from the API URL, hence capitalization woes
  list_of_services <- c(
    "elevation" = "3DEPElevation",
    "ortho" = "USGSNAIPPlus",
    "hydro" = "nhd",
    "govunits",
    "contours",
    "geonames",
    "NHDPlus_HR",
    "structures",
    "transportation",
    "watersheds" = "wbd",
    "ecosystems",
    "USGSTopo",
    "USGSShadedReliefOnly",
    "USGSImageryOnly",
    "USGSHydroCached",
    "USGSTNMBlank"
  )

  stopifnot(all(services %in% list_of_services |
    services %in% names(list_of_services)))

  tif_files <- c("3DEPElevation")
  png_files <- list_of_services[!(list_of_services %in% tif_files)]

  services <- stats::setNames(services, services)

  if (any(services %in% names(list_of_services))) { # cast short codes now
    replacements <- which(services %in% names(list_of_services))
    services[replacements] <- as.vector(
      list_of_services[services[replacements]]
    )
  }

  # duplicated instead of unique to preserve names
  services <- services[!duplicated(services)]

  if (is.null(side_length)) {
    if (any(services %in% png_files)) {
      side_length <- 4096
    } else {
      side_length <- 8000
    }
  }

  if (any(services %in% png_files) && side_length > 4096) {
    rlang::abort(c(
      "These tiles have a maximum side length of 4096.",
      i = "Set `side_length` to 4096 or less"
    ))
  }
  if (("3DEPElevation" %in% services) && side_length > 8000) {
    rlang::abort(c(
      "3DEPElevation tiles have a maximum side length of 8000.",
      i = "Set `side_length` to 8000 or less"
    ))
  }

  bbox_splits <- split_bbox(data, side_length, resolution, projected)
  tile_boxes <- bbox_splits[[1]]
  x_tiles <- bbox_splits[[2]]
  y_tiles <- bbox_splits[[3]]

  if (requireNamespace("progressr", quietly = TRUE)) { # nocov start
    p <- progressr::progressor(steps = x_tiles * y_tiles * length(services))
  } # nocov end

  for (i in seq_len(x_tiles)) {
    for (j in seq_len(y_tiles)) {
      current_box <- tile_boxes[tile_boxes$x_tiles == i &
        tile_boxes$y_tiles == j, ]
      current_bbox <- data.frame(
        lat = c(current_box$min_y, current_box$max_y),
        lng = c(current_box$min_x, current_box$max_x)
      )
      current_bbox <- sf::st_as_sf(current_bbox, coords = c("lng", "lat"))
      current_bbox <- sf::st_bbox(current_bbox)
      for (k in seq_along(services)) {
        if (requireNamespace("progressr", quietly = TRUE)) { # nocov start
          p(message = sprintf(
            "Retrieving %s tile (%d, %d)",
            services[[k]],
            i,
            j
          ))
        } # nocov end

        if (services[[k]] %in% tif_files | georeference) {
          fileext <- ".tif"
        } else if (services[[k]] %in% png_files) {
          fileext <- ".png"
        }

        final_path <- paste0(
          output_prefix,
          "_",
          services[[k]],
          "_",
          i,
          "_",
          j,
          fileext
        )

        if (georeference && services[[k]] != "3DEPElevation") {
          rm_path <- TRUE
          cur_path <- tempfile(fileext = ".png")
        } else {
          rm_path <- FALSE
          cur_path <- final_path
        }

        counter <- 0
        while ((!file.exists(cur_path) ||
          file.size(cur_path) == 0) &&
          counter < 5) {
          img_bin <- hit_national_map_api(
            current_bbox,
            current_box[["img_width"]],
            current_box[["img_height"]],
            services[[k]],
            verbose = verbose,
            ...
          )
          if (is.raw(img_bin$imageData)) {
            writeBin(img_bin$imageData, cur_path)
          } else {
            outconn <- file(cur_path, "wb")
            base64enc::base64decode(
              what = img_bin$imageData,
              output = outconn
            )
            close(outconn)
          }
        }

        if (georeference && services[[k]] != "3DEPElevation") {
          cur_raster <- terra::rast(png::readPNG(cur_path))
          terra::crs(cur_raster) <- paste0(
            "+init=EPSG:",
            img_bin$extent$spatialReference$wkid
          )
          terra::ext(cur_raster) <- c(
            img_bin$extent$xmin,
            img_bin$extent$xmax,
            img_bin$extent$ymin,
            img_bin$extent$ymax
          )

          terra::writeRaster(cur_raster,
            final_path,
            overwrite = TRUE
          )
          if (rm_path) unlink(cur_path)
        }
      }
    }
  }

  res <- vector("list")
  for (i in seq_along(services)) {
    if (services[[i]] %in% tif_files || georeference) {
      fileext <- ".tif"
    } else if (services[[i]] %in% png_files) {
      fileext <- ".png"
    }
    res[[i]] <- paste0(
      output_prefix,
      "_",
      services[[i]],
      "_",
      outer(1:x_tiles, 1:y_tiles, paste, sep = "_"),
      fileext
    )
    names(res)[[i]] <- names(services)[[i]]
  }

  return(invisible(res))
}

handle_bboxSR <- function(data, projected) {
  if (!is.na(sf::st_crs(data)$epsg)) {
    sf::st_crs(data)$epsg
  } else if (projected) {
    rlang::warn(c(
      "Assuming CRS of EPSG 5071",
      i = "Set bboxSR explicity to override"
      ))
    5071
  } else {
    rlang::warn(c(
      "Assuming CRS of EPSG 4326",
      i = "Set bboxSR explicity to override"
    ))
    4326
  }
}

#' Split a bounding box into a smaller set of component tiles.
#'
#' This function splits the area contained within a bounding box into a set of
#' tiles.
#'
#' @param bbox The bounding box to split into tiles.
#' @param side_length The length of each side of the output tiles.
#' @param resolution How many meters are represented by each pixel? The default
#' value of 1 means that 1 pixel = 1 meter, while a value of 2 means that
#' 1 pixel = 2 meters, and so on.
#' @param projected Logical: is `data` in a projected coordinate reference
#' system? If `NULL`, the default, inferred from [sf::st_is_longlat].
#'
#' @keywords internal
#'
#' @return A data frame containing x and y tile indices, bounding boxes, and
#' side lengths.
#'
#' @noRd
split_bbox <- function(data, side_length, resolution = 1, projected) {
  if (!methods::is(data, "terrainr_bounding_box")) {
    bbox <- terrainr_bounding_box(
      c("lng" = data[["xmin"]], "lat" = data[["ymin"]]),
      c("lng" = data[["xmax"]], "lat" = data[["ymax"]])
    )
  } else {
    bbox <- data
  }
  tl <- terrainr_coordinate_pair(c(bbox@tr@lat, bbox@bl@lng))

  if (projected) {
    img_width <- round(
      (bbox@tr@lng - tl@lng) / resolution,
      digits = 0
    )
    img_height <- round(
      (tl@lat - bbox@bl@lat) / resolution,
      digits = 0
    )
  } else {
    img_width <- round(
      calc_haversine_distance(tl, bbox@tr) / resolution,
      digits = 0
    )
    img_height <- round(
      calc_haversine_distance(tl, bbox@bl) / resolution,
      digits = 0
    )
  }

  x_tiles <- ceiling(img_width / side_length)
  y_tiles <- ceiling(img_height / side_length)

  tile_boxes <- expand.grid(
    x_tiles = 1:x_tiles,
    y_tiles = 1:y_tiles,
    min_x = NA,
    max_x = NA,
    min_y = NA,
    max_y = NA,
    img_width = NA,
    img_height = NA
  )

  if (projected) {
    for (i in 1:x_tiles) {
      tile_boxes[tile_boxes$x_tiles == i, ]$min_x <-
        bbox@bl@lng + (side_length * (i - 1) * resolution)
      if (i == x_tiles) {
        tile_boxes[tile_boxes$x_tiles == i, ]$max_x <- bbox@tr@lng
      } else {
        tile_boxes[tile_boxes$x_tiles == i, ]$max_x <-
          bbox@bl@lng + (side_length * i * resolution)
      }
    }
    for (j in 1:y_tiles) {
      tile_boxes[tile_boxes$y_tiles == j, ]$max_y <-
        bbox@tr@lat - (side_length * (j - 1) * resolution)
      if (j == y_tiles) {
        tile_boxes[tile_boxes$y_tiles == j, ]$min_y <- bbox@bl@lat
      } else {
        tile_boxes[tile_boxes$y_tiles == j, ]$min_y <-
          bbox@tr@lat - (side_length * j * resolution)
      }
    }
  } else {
    for (i in 1:x_tiles) {
      tile_boxes[tile_boxes$x_tiles == i, ]$min_x <- point_from_distance(
        bbox@bl,
        side_length * (i - 1) * resolution,
        90
      )@lng
      if (i == x_tiles) {
        tile_boxes[tile_boxes$x_tiles == i, ]$max_x <- bbox@tr@lng
      } else {
        tile_boxes[tile_boxes$x_tiles == i, ]$max_x <- point_from_distance(
          bbox@bl,
          side_length * i * resolution,
          90
        )@lng
      }
    }
    for (j in 1:y_tiles) {
      tile_boxes[tile_boxes$y_tiles == j, ]$max_y <- point_from_distance(
        bbox@tr,
        side_length * (j - 1) * resolution,
        180
      )@lat
      if (j == y_tiles) {
        tile_boxes[tile_boxes$y_tiles == j, ]$min_y <- bbox@bl@lat
      } else {
        tile_boxes[tile_boxes$y_tiles == j, ]$min_y <- point_from_distance(
          bbox@tr,
          side_length * j * resolution,
          180
        )@lat
      }
    }
  }

  for (i in seq_len(nrow(tile_boxes))) {
    x_idx <- tile_boxes[i, ]$x_tiles
    y_idx <- tile_boxes[i, ]$y_tiles
    tile_boxes[i, ]$img_width <- ifelse(
      (img_width - (x_idx * side_length)) < 0,
      img_width - ((x_idx - 1) * side_length),
      side_length
    )
    tile_boxes[i, ]$img_height <- ifelse(
      img_height - (y_idx * side_length) < 0,
      img_height - ((y_idx - 1) * side_length),
      side_length
    )
  }

  return(list(
    tile_boxes,
    x_tiles,
    y_tiles
  ))
}