R/rarify.R

In rLakeHabitat: Interpolate Bathymetry and Quantify Physical Aquatic Habitat

#' Rarify Depth Data
#'
#' Reduce density of mapped depth data to improve accuracy and computation time.
#'
#' @param outline shapefile outline of a waterbody. Accepts a SpatVector, an sf object, or anything terra::vect() can read.
#' @param df dataframe of coordinates and depths for a given waterbody. Coordinates are assumed to be in the same CRS as 'outline'.
#' @param x character giving name of longitude column
#' @param y character giving name of latitude column
#' @param z character giving name of depth column
#' @param res number describing the target spacing between rarified points, in meters, default = 10
#' @details
#' The function automatically detects whether 'outline' (and therefore 'df') is in a geographic (decimal degree) or
#' projected (meters) coordinate system, the same way interpBathy() does. If geographic, points are rarified in the
#' waterbody's best-fit UTM zone so that 'res' is honored as a true physical spacing in meters, then the rarified
#' points are returned in the original CRS of 'outline'.
#' @return dataframe of rarified xyz coordinates (columns named 'x', 'y', 'z' regardless of the input column names),
#' in the same CRS as the original input
#' @author Sean Bertalot & Tristan Blechinger, Department of Zoology & Physiology, University of Wyoming
#' @export
#' @import dplyr
#' @rawNamespace import(terra, except = c(union,intersect, animate))
#' @examples
#' \donttest{
#' outline <- terra::vect(system.file("extdata", "example_outline.shp", package = 'rLakeHabitat'))
#' depths <- read.csv(system.file("extdata", "example_depths.csv", package = 'rLakeHabitat'))
#' rarify(outline = outline, df = depths, x = "x", y = "y", z = "z", res = 100)}

rarify <- function(outline, df, x, y, z, res = 10){

  #transform outline shapefile into vector
  if(!inherits(outline, "SpatVector")){
    if(inherits(outline, "sf") && requireNamespace("sf", quietly = TRUE)){
      outline <- sf::st_zm(outline, drop = TRUE, what = "ZM")
    }
    outline <- terra::vect(outline)
  }
  else{
    outline <- outline
  }

  #store the original CRS so rarified points can be returned in it
  original_crs <- terra::crs(outline)

  #checks
  if(!inherits(df, "data.frame"))
    stop("df must be a dataframe")
  if(!inherits(x, "character"))
    stop("x must be a character giving the longitude column name")
  if(!inherits(y, "character"))
    stop("y must be a character giving the latitude column name")
  if(!inherits(z, "character"))
    stop("z must be a character giving the depth column name")
  if(x %in% names(df) == FALSE)
    stop("The value of x does not appear to be a valid column name")
  if(y %in% names(df) == FALSE)
    stop("The value of y does not appear to be a valid column name")
  if(z %in% names(df) == FALSE)
    stop("The value of z does not appear to be a valid column name")
  if(!inherits(df[, x], "numeric"))
    stop("data in x column is not formatted as numeric")
  if(!inherits(df[, y], "numeric"))
    stop("data in y column is not formatted as numeric")
  if(!inherits(df[, z], "numeric"))
    stop("data in z column is not formatted as numeric")
  if(!inherits(outline, "SpatVector"))
    stop("outline is not a SpatVector or cannot be transformed")
  if(is.na(res) || is.null(res))
    stop("res must be specified as a numeric value")
  if(!is.numeric(res))
    stop("res must be specified as a numeric value")

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

  # Function to determine the best UTM CRS for a given vector
  get_best_utm <- function(outline) {
    centroid <- terra::centroids(outline)
    lon <- terra::crds(centroid)[1]
    lat <- terra::crds(centroid)[2]
    utm_zone <- base::floor((lon + 180) / 6) + 1
    hemisphere <- base::ifelse(lat >= 0, 32600, 32700)
    epsg_code <- hemisphere + utm_zone
    return(terra::crs(paste0("EPSG:", epsg_code)))
  }

  # automatically detect geographic vs. projected CRS
  if(terra::is.lonlat(outline)){
    best_crs <- get_best_utm(outline)
    outline <- terra::project(outline, best_crs)
  }

  # if outline was reprojected above, reproject the point data to match
  if(!identical(terra::crs(outline), original_crs)){
    pts <- terra::vect(df, geom = c(x, y), crs = original_crs)
    pts <- terra::project(pts, terra::crs(outline))
    proj_coords <- terra::crds(pts)
    df[[x]] <- proj_coords[, 1]
    df[[y]] <- proj_coords[, 2]
  }

  ## Helper: figure out how many rows/columns are needed to hit the requested resolution
  get_res <- function(outline, res) {
    ext_o <- terra::ext(outline)
    ext_length <- base::abs(ext_o$xmin - ext_o$xmax)
    ext_height <- base::abs(ext_o$ymax - ext_o$ymin)

    set_ext_x <- ext_length / res
    set_ext_y <- ext_height / res

    xy <- c(set_ext_x, set_ext_y)
    return(xy)
  }

  xy <- get_res(outline, res)
  empty_raster <- terra::rast(ext(outline), ncol = xy[1], nrow = xy[2], crs = terra::crs(outline))

  # take in xyz dataframe and make it spatial dataframe using vect
  points_unrarified <- terra::vect(df, geom = c(x, y), crs = terra::crs(outline))

  # creates a raster of the shape outline in the grid empty_raster
  ras <- terra::rasterize(outline, empty_raster)

  # masks the raster for the shapefile (everything outside the reservoir = NA)
  grid <- terra::mask(ras, outline)

  # rasterize points using the empty grid generated using the get_res coords
  rasterized_points <- terra::rasterize(points_unrarified, grid, field = z, fun = "mean")

  # turn back into xyz dataframe
  points_df <- as.data.frame(rasterized_points, xy = TRUE)
  names(points_df)[ncol(points_df)] <- "z"

  # reproject rarified points back to the original CRS
  if(!identical(terra::crs(outline), original_crs)){
    pts_out <- terra::vect(points_df, geom = c("x", "y"), crs = terra::crs(outline))
    pts_out <- terra::project(pts_out, original_crs)
    proj_coords <- terra::crds(pts_out)
    points_df$x <- proj_coords[, 1]
    points_df$y <- proj_coords[, 2]
  }

  return(points_df)
}