R/nearest.features.R

In mottensmann/DBChecks: Checking Database

#' Habitat analysis: Nearest features
#'
#' @description
#' Calculates the shortest distance between a point and the edges of a feature.
#' Preprocessing as in habitat.analysis.R
#'
#' @details
#' 1.) If feature is a polygon convert to lines using SAGA
#' 2.) Create points along lines
#' 3.) find nearest match between points of both layers
#'
#' @param features list of ESRI shapes
#' @param overlay overlay ESRI shape file
#' @param point point ESRI shape file
#' @inheritParams meancoordinates
#' @import magrittr
#' @export
#'
#'
nearest.features <- function(features = NULL,
                             overlay =  NULL,
                             point = stringr::str_replace(overlay, "buffer.shp", ".shp"),
                             root =  'C:/OSGeo4W64') {

  ## checks:
  ## ###########################################################################
  if (class(features) != "list") stop("features must be a list of shapes")

  ## set up RQGIS
  ## ###########################################################################
  if ("RQGIS3" %in% rownames(installed.packages())) {
    library(RQGIS3, quietly = T)
    ## Where to find RQGIS
    if (!is.null(root)) {
      RQGIS3::set_env(root = root)
    } else {
      RQGIS3::set_env()
    }
  } else {
    stop("Install RQGIS3 to proceed")
  }
  ## ###########################################################################

  ## define output file
  .RData <- stringr::str_replace(overlay, ".shp", ".edges.RData")


  ## Loop through features
  ## ###########################################################################
  results <- lapply(features, function(one_feature) {
    #cat(one_feature, "\n")

    temp.file <- temp.shp()
    #temp.folder <- stringr::str_split(string = temp.file, pattern = "\\\\", simplify = T)
    temp.folder <- stringr::str_split(string = temp.file, pattern = "/", simplify = T)[1]

    ## keep track of temp files to remove them once not needed anymore
    unlink_file <- file.path(temp.folder, "unlink.txt")
    if (!file.exists(unlink_file)) {
      cat("Create", unlink_file, "\n")
      sink(unlink_file, append = F)
      sink()
    }

    ## Clip to overlay extent
    ## ######################
    x <- RQGIS3::run_qgis(
      alg = "native:clip",
      INPUT = one_feature,
      OVERLAY = overlay,
      OUTPUT = temp.file,
      load_output = T,
      show_output_paths = F)
    ## ######################

    ## Convert to Lines
    ## ####################
    if (nrow(x) > 1) {
    temp.file2 <- temp.shp()
    if (any(names(x) == "Area")) {
      x <- RQGIS3::run_qgis(
        alg = "saga:convertpolygonstolines",
        POLYGONS = temp.file,
        LINES = temp.file2,
        load_output = T,
        show_output_paths = F)
    } else if (any(names(x) %in% c("length", "Length"))) {
      temp.file2 <- temp.file
    } else {
      stop(paste("Unknown geometry in", one_feature))
    }

    ## Points along line
    ## ####################
    temp.file3 <- temp.shp()
    #temp.file3 <- "../../01-PhD/01-Paper/03-Age-of-first-reproduction/data/GIS/nest01-points.shp"
    x <- RQGIS3::run_qgis(
      alg = "native:pointsalonglines",
      INPUT = temp.file2,
      OUTPUT = temp.file3,
      load_output = T,
      show_output_paths = F)

    ## Nearest match
    temp.file4 <- temp.shp()
    #temp.file4 <- "../../01-PhD/01-Paper/03-Age-of-first-reproduction/data/GIS/nest01-knea.shp"
    x <- RQGIS3::run_qgis(
      alg = "qgis:distancematrix",
      TARGET = temp.file3,
      TARGET_FIELD = "ID",
      INPUT = point,
      INPUT_FIELD = "name",
      NEAREST_POINTS = 1,
      OUTPUT = temp.file4,
      show_output_paths = F,
      load_output = T)

    ## discard temp files
    ## ##################
    silent <- lapply(c(temp.file, temp.file2, temp.file3, temp.file4), unlink.shp) %>%
      do.call("rbind",.) %>%
      set_rownames(1:nrow(.))
    ## if some can not be removed retry
    #print(silent)
    to_unlink <- character()
    to_unlink <- dplyr::filter(silent, success == "No")$file
    #####################

    output <- ifelse(length(x$Distance) > 0, min(x$Distance), NA)

    } else {
    ## no geometry detected
      ## discard temp files
      ## ##################
      silent <- unlink.shp(temp.file) %>%
        set_rownames(1:nrow(.))
      ## if some can not be removed retry
      #print(silent)
      to_unlink <- character()
      to_unlink <- dplyr::filter(silent, success == "No")$file
      #####################
      output <- NA
    }

    ## add to file and try to delete
    if (length(to_unlink) > 0) {
      ## write to file
      sink(unlink_file, append = TRUE)
      cat(to_unlink, sep = "\n")
      sink()
      ## read file content
      to_unlink <- readr::read_delim(unlink_file, delim = "\t",
                                     skip_empty_rows = T,
                                     col_names = F,
                                     col_types = "c")[[1]]
      silent <- sapply(to_unlink, unlink)
      to_unlink <- data.frame(file = to_unlink,
                              success = ifelse(silent == 0, "Yes", "No"))
      to_unlink <- dplyr::filter(to_unlink, success == "No")$file
      ## Update
      sink(unlink_file, append = F)
      cat(to_unlink, sep = "\n")
      sink()

    }
    return(output)

  }) %>%
    set_names(names(features)) %>%
    do.call("cbind",.) %>%
    as.data.frame()

  save(results, file = .RData)
  return(results)
}

# rm(list = ls())
# library(DBChecks)
# test <- nearest.features(features = list(
#   Forest = "../../01-PhD/GIS-Projects/04-HabitatImprinting/DLM/Forest.shp",
#   Mainroad = "../../01-PhD/GIS-Projects/04-HabitatImprinting/DLM/MainRoad.shp"),
#   overlay = "../../01-PhD/01-Paper/03-Age-of-first-reproduction/data/GIS/nest1200buffer2000.shp",
#   point = "../../01-PhD/01-Paper/03-Age-of-first-reproduction/data/GIS/nest1200.shp",
#   root = 'C:/OSGeo4W64')
# # # # # # # # #
# folder <- "c:\\Users/MOMO/AppData/Local/Temp/"
# files <- list.files(path = folder, recursive = T, pattern = ".shp$", full.names = T)
# x <- pbapply::pblapply(files, unlink.shp, fileEnding = ".shp")
# #
# unlink.shp(files[2])

#one_feature <- features$Wood