R/saveContours.R

In rLakeHabitat: Interpolate Bathymetry and Quantify Physical Aquatic Habitat

#' Generate and Export Bathymetric Contours
#'
#' Generate depth contours from an interpolated DEM and save them to a file format usable outside R - a GPS unit or chartplotter (.gpx), Google Earth or similar (.kml), or standard GIS formats (.shp, .gpkg).
#'
#' @param dem a SpatRaster of interpolated bathymetry (e.g. from interpBathy()). If 'dem' has multiple layers (e.g. the depth/error output of OK or UK), the layer named 'depth' is used if present, otherwise the first layer.
#' @param by numeric value giving a regular contour interval, e.g. by = 5 contours every 5 depth units starting at 5 (the 0 contour is skipped, since it's just the shoreline itself). Exactly one of 'by' or 'levels' must be specified.
#' @param levels numeric vector of specific depth values to contour, e.g. c(1, 5, 10, 20). Exactly one of 'by' or 'levels' must be specified.
#' @param file_type character giving the output format: "gpx", "kml", "shp", or "gpkg". Default = "gpx".
#' @param file_name character giving the output file name, with or without a path. Default = "contours".
#' @details
#' GPX and KML both require geographic (WGS84 lon/lat) coordinates - if 'dem' is in a projected CRS, the contours
#' are automatically reprojected to EPSG:4326 before being written (.shp/.gpkg outputs keep 'dem's original CRS).
#' @return an sf object of the generated contour lines (with a 'depth' column), invisibly saved to 'file_name' as a
#' side effect. Returned regardless of 'file_type' so you can plot or inspect it without re-reading the saved file.
#' @author Tristan Blechinger, Department of Zoology & Physiology, University of Wyoming
#' @export
#' @rawNamespace import(terra, except = c(union,intersect, animate))
#' @examples
#' \donttest{
#' if (requireNamespace("sf", quietly = TRUE)) {
#' outline <- terra::vect(system.file("extdata", "example_outline.shp", package = 'rLakeHabitat'))
#' data <- read.csv(system.file("extdata", "example_depths.csv", package = 'rLakeHabitat'))
#' dem <- interpBathy(outline, data, "x", "y", "z", zeros = FALSE, separation = 10,
#' res = 10, method = "IDW", nmax = 4)
#' saveContours(dem, by = 2, file_type = "gpx", file_name = tempfile("contours"))
#' }}

saveContours <- function(dem, by = NULL, levels = NULL, file_type = "gpx", file_name = "contours"){

  if(!requireNamespace("sf", quietly = TRUE))
    stop("This function requires the 'sf' package. Install it with install.packages('sf').")

  # checks
  if(!inherits(dem, "SpatRaster"))
    stop("dem must be a SpatRaster (e.g., the output of interpBathy())")
  if(!is.character(file_type) || length(file_type) != 1 || !tolower(file_type) %in% c("gpx", "kml", "shp", "gpkg"))
    stop("file_type must be one of 'gpx', 'kml', 'shp', or 'gpkg'")
  if(!is.character(file_name) || length(file_name) != 1)
    stop("file_name must be a single character string")
  if(is.null(by) && is.null(levels))
    stop("either 'by' or 'levels' must be specified")
  if(!is.null(by) && !is.null(levels))
    stop("only one of 'by' or 'levels' should be specified, not both")
  if(!is.null(by)){
    if(!is.numeric(by) || length(by) != 1 || by <= 0)
      stop("by must be a single positive numeric value")
  }
  if(!is.null(levels)){
    if(!is.numeric(levels) || length(levels) == 0)
      stop("levels must be a numeric vector")
  }

  test_crs <- terra::crs(dem)
  if(is.na(test_crs) || test_crs == ""){
    stop("CRS of 'dem' is unable to be defined.")
  }

  file_type <- tolower(file_type)

  # if dem has multiple layers (e.g. OK/UK's depth + error output), use the
  # 'depth' layer if present, otherwise fall back to the first layer
  if(terra::nlyr(dem) > 1){
    if("depth" %in% names(dem)){
      dem <- dem[["depth"]]
    } else {
      dem <- dem[[1]]
      warning("dem has multiple layers; using the first layer ('", names(dem)[1], "') for contouring.")
    }
  }

  # determine contour levels
  if(!is.null(by)){
    max_depth <- as.numeric(terra::global(dem, "max", na.rm = TRUE)[[1]])
    if(by > max_depth){
      levels <- numeric(0)
    } else {
      levels <- seq(by, max_depth, by = by)
    }
    if(length(levels) == 0)
      stop("no contour levels fall within dem's depth range given 'by' - check 'by' against dem's max depth")
  }

  # generate contour lines
  contours <- terra::as.contour(dem, levels = levels)
  if(is.null(contours) || length(contours) == 0)
    stop("no contours could be generated at the requested levels - check 'levels'/'by' against dem's depth range")

  names(contours)[1] <- "depth"

  # GPX and KML both require geographic (WGS84) coordinates
  if(file_type %in% c("gpx", "kml")){
    if(!terra::is.lonlat(contours)){
      contours <- terra::project(contours, "EPSG:4326")
    }
  }

  contours_sf <- sf::st_as_sf(contours)

  out_file <- paste0(file_name, ".", file_type)
  if(file.exists(out_file)) file.remove(out_file)

  if(file_type == "gpx"){

    # Explode to LINESTRING via MULTILINESTRING first
    contours_sf <- suppressWarnings(sf::st_cast(contours_sf, "MULTILINESTRING"))
    contours_sf <- suppressWarnings(sf::st_cast(contours_sf, "LINESTRING"))

    contours_sf$name <- paste0("Contour_", contours_sf$depth, "m")
    # GPX only understands a small fixed schema - keep just what's needed
    contours_sf <- contours_sf["name"]

    sf::st_write(contours_sf, out_file, driver = "GPX", layer_options = "FORCE_GPX_ROUTE=YES")
  }
  else if(file_type == "kml"){
    # GDAL's KML writer looks for a field named 'Name' (capitalized) to label placemarks
    contours_sf$Name <- paste0("Contour_", contours_sf$depth, "m")

    sf::st_write(contours_sf, out_file, driver = "KML")
  }
  else{
    contours_sf$name <- paste0("Contour_", contours_sf$depth, "m")

    sf::st_write(contours_sf, out_file)
  }

  message("Contours saved to ", out_file)

  return(contours_sf)
}