Lslide: Functions for object-oriented image analysis in landslide science

Documented in neighborDirection

#' Calculation of the direction of every object to its neighbors
#'
#' This function neighborDirection calculate the direction in degree to its surrounding neighbors or a selection
#'
#' @param spdf \linkS4class{SpatialPolygonsDataFrame}, or full path of a shapefile
#' @param spdf.bb \linkS4class{SpatialPolygonsDataFrame} from which the bounding boxes are created
#' @param col.name column name for angles (in degree, N = 0, E = 90, S = 180, W = 270)
#' @param modus "bb" for using a boundig box or "nb" for using object neighborhood \emph{nb}. Modus \emph{bb} returns data.frame based on spdf. Modus \emph{nb} returns list based on nb.
#' @param tol tolerance for select neighbor of specific range. Default: 360
#' @param sp.PoS \linkS4class{SpatialPoints} to which the direction is calculated. As standard they are calculated from spdf input. Default: NULL
#' @param nb object neighborhood based on \code{\link[spdep]{poly2nb}}, for modus "nb. Default: Null
#' @param bb bounding box for modus "bb". Default: Null
#' @param cores number of cores for parallel processing. Default: 1 (sequential)
#' @param quiet no outputs in console. Default: TRUE
#'
#' @return
#' In the case of modus \emph{nb} a list is returned. If modus \emph{bb} is set, a data.frame based on the input spdf is returned
#'
#' @note
#' \itemize{
#'   \item \linkS4class{SpatialPolygonsDataFrame} input must have an \emph{ID} field with unique numbers
#' }
#'
#'
#' @keywords neighbor direction, object-oriented image analysis
#'
#'
#' @export
neighborDirection <- function(spdf, col.name, modus = "nb", tol = 360, spdf.bb = NULL, sp.PoS = NULL, nb = NULL, bb = NULL, cores = 1, quiet = TRUE, ...)
{

  # get start time of process
  process.time.start <- proc.time()

  if(is.null(sp.PoS ))
  {
    if(quiet == FALSE){cat("Point on surface...\n")}
    sp.PoS <- rgeos::gPointOnSurface(spgeom = spdf, byid = TRUE)
    spdf.PoS <- SpatialPointsDataFrame(sp.PoS, data.frame(ID = row.names(spdf)))
  }

  if(modus == "bb" & is.null(bb))
  {
    if(is.null(spdf.bb)){stop("Input for spdf.bb is missing!\n")}

    if(quiet == FALSE){cat("Creation of bounding box...\n")}
    bb <- getBoundingBox(shape = spdf.bb, col.name = col.name, ...)
  } else {
    spdf.bb <- bb
  }

  if(modus == "nb" & is.null(nb))
  {
    if(quiet == FALSE){cat("Creation of neighborhood...\n")}
    nb_speed_up <- rgeos::gUnarySTRtreeQuery(spdf) # speed up function for poly2nb
    nb <- spdep::poly2nb(spdf, queen = TRUE, foundInBox = nb_speed_up) # neighborhood based on queen continuity
  }


  # initialize variables for loop
  if(modus == "bb"){obj.inTol.T <- c()} # variable containing information about which neighbors are in tolerance to angle

  # obj.nb.angle <- list()

  # get object that intersect with class object
  if(modus == "bb")
    {
    if(quiet == FALSE){cat("Intersect geometries...\n")}
    obj.inter <- rgeos::gIntersects(spgeom1 = spdf.bb, spgeom2 = spdf, byid = TRUE, returnDense = FALSE)


    # index.class <- as.numeric(row.names(spdf[which(spdf@data$ID == spdf.bb@data$ID),]))
    # xy.class <- sp::coordinates(sp.PoS[index.class,])
    index.class <- match(spdf.bb@data$ID, spdf@data$ID)
    xy.class <- sp::coordinates(sp.PoS[index.class,])

    xy.class.angle <- spdf[index.class,]@data[[col.name]]

    # intersected objects and remove classes
    obj.inter.unlCl <- unique(unlist(obj.inter))
    obj.inter.unlCl <- setdiff(obj.inter.unlCl, index.class)

    # get subset of points intersecting the class bounding box
    obj.inter.pts <- sp.PoS[obj.inter.unlCl,]

    # extract cooridates of points
    xy.obj.inter <- sp::coordinates(obj.inter.pts)

    # get amount of polygons in shape
    run <- length(obj.inter)
  }

  if(modus == "nb")
  {
    # align variable names
    obj.inter <- nb
    index.class <- as.numeric(row.names(spdf))
    xy.class.angle <- spdf@data[[col.name]]
    xy.class <- sp::coordinates(sp.PoS)

    # extract cooridates of points
    xy.obj.inter <- sp::coordinates(spdf.PoS)

    # get amount of polygons in shape
    run <- length(spdf)
  }

  # start loop ------
  if(quiet == FALSE){cat("Start loop...\n")}


  # init parallel
  switch(Sys.info()[[1]],
         Windows = type  <- "PSOCK",
         Linux = type  <- "FORK",
         Mac = type <- "FORK")

  cl <- parallel::makeCluster(cores, type = type)

  if(Sys.info()[[1]] == "Windows")
  {
    parallel::clusterExport(cl = cl, envir = environment(),
                            varlist = c('obj.inter', 'index.class', 'xy.class.angle', 'xy.class',
                                        'tol', 'modus'))
  }


  obj.nb.angle <- parallel::parLapply(cl = cl, X = 1:run, function(i, obj.inter, index.class, xy.class.angle, xy.class, tol, modus){
  # obj.nb.angle <- lapply(X = 1:run, FUN = function(i, obj.inter, index.class, xy.class.angle, xy.class, tol, modus){
  # for(i in 1:run)
  # {
    # print(i)
    if(quiet == FALSE)
    {
      if(i == round(run * 0.25)) {cat("... 25 % done\n")}
      if(i == round(run * 0.5)) {cat("... 50 % done\n")}
      if(i == round(run * 0.75)) {cat("... 75 % done\n")}
    }

    # get object that intersect with class object
    # obj.inter <- unique(unlist(rgeos::gIntersects(spgeom1 = spdf.bb[i,], spgeom2 = spdf, byid = TRUE, returnDense = FALSE)))

    # remove class object itself
    # index.class <- as.numeric(row.names(spdf[which(spdf@data$ID == spdf.bb[i,]@data$ID),]))
    # obj.inter.sub <- setdiff(obj.inter[[i]], index.class)
    obj.inter.i <- obj.inter[[i]]
    obj.inter.i <- setdiff(obj.inter.i, index.class[i])


    if(!is.null(obj.inter.i) & xy.class.angle[i] != -9999 & !is.na(xy.class.angle[i]))
    {
      # get cooridante from actual i-class
      # xy.class <- sp::coordinates(sp.PoS[index.class,])
      # xy.class.angle <- spdf[index.class,]@data[col.name][1,1]
      xy.class.angle.i <- xy.class.angle[i]

      # get subset of points intersecting the class bounding box
      # obj.inter.pts <- sp.PoS[obj.inter,]

      # extract cooridates of points
      # xy.obj.inter <- sp::coordinates(obj.inter.pts)
      xy.obj.inter.i <- xy.obj.inter[match(obj.inter.i, as.numeric(row.names(xy.obj.inter))),]

      if(class(xy.obj.inter.i ) == "matrix")
      {
        # substract class cooridante from matrix
        xy.obj.inter.i[,1] <- xy.obj.inter.i[,1] - xy.class[i,1] # x coordiante
        xy.obj.inter.i[,2] <- xy.obj.inter.i[,2] - xy.class[i,2] # y coordiante

        xy.obj.angle.i <- (atan2(x = xy.obj.inter.i[, 1], y = xy.obj.inter.i[, 2]) * ((-180)/pi) + 90) %% 360


      } else{
        xy.obj.inter.i[1] <- xy.obj.inter.i[1] - xy.class[i,1] # x coordiante
        xy.obj.inter.i[2] <- xy.obj.inter.i[2] - xy.class[i,2] # y coordiante

        xy.obj.angle.i <- (atan2(x = xy.obj.inter.i[1], y = xy.obj.inter.i[2]) * ((-180)/pi) + 90) %% 360

      }




      # # #  calculate mean angle
      # The negative on angle deals with the fact that we are changing from counterclockwise to clockwise.
      # The +90 deals with the offset of ninety degrees. And lastly we need to mod by 360 to keep our angle in the desired range
      # Angles in R are counterclockwise from E:0|360, N 90, W: 180, S: 270
      # get direction
      # xy.obj.angle.i <- (atan2(x = xy.obj.inter.i[, 1], y = xy.obj.inter.i[, 2]) * ((-180)/pi) + 90) %% 360
      # browser()
      # compare with angle direction and selection
      names(xy.obj.angle.i) <- obj.inter.i # changed!!!!

      anglediff <- (xy.class.angle.i - xy.obj.angle.i + 180 + 360) %% 360 - 180


      xy.obj.inTol <- sapply(X = anglediff, FUN = function(x, tol)
                        {
                          ifelse(x <= tol && x >= -tol, 1, 0)
                      }, tol = tol)

      if(length(which(xy.obj.inTol == 1)) >= 1)
      {
        # add results to variables
        if(modus == "bb")
        {
          obj.inTol.T.sub <- as.numeric(unique(names(xy.obj.inTol[which(xy.obj.inTol == 1)])))
          # if(length(obj.inTol.T.sub) == 0){obj.inTol.T.sub <- NA}

          obj.inTol.T <- c(unique(obj.inTol.T), obj.inTol.T.sub) # indices!
        }

        xy.obj.angle.sub <- xy.obj.angle.i[which(xy.obj.inTol == 1)]
        # obj.nb.angle[[i]] <- list(Object = index.class[i], NeighborDirection = xy.obj.angle.sub)
        return(list(Object = index.class[i], NeighborDirection = xy.obj.angle.sub))
      } else {
        # obj.nb.angle[[i]] <- list(Object = index.class[i], NeighborDirection = 0)
        return(list(Object = index.class[i], NeighborDirection = 0))
      }


    } else {
      # obj.nb.angle[[i]] <- list(Object = index.class[i], NeighborDirection = NA)
      return(list(Object = index.class[i], NeighborDirection = NA))
      }# end of if

  # } # end of loop
    }, obj.inter = obj.inter, index.class = index.class, xy.class.angle = xy.class.angle,
       xy.class = xy.class, tol = tol, modus = modus)


  parallel::stopCluster(cl)


  # get time of process
  process.time.run <- proc.time() - process.time.start
  if(quiet == FALSE){cat("------ Run of ClassNeighborFunction: " , round(x = process.time.run["elapsed"][[1]]/60, digits = 4), " Minutes ------ \n")}

  if(modus == "bb")
  {
    df <- data.frame(Index = row.names(spdf), ID = spdf@data$ID, NbInDir = 0, NbDir = NA)
    df[obj.inTol.T, ]$NbInDir <- 1
    df[unlist(lapply(obj.nb.angle, `[[`, 1)), ]$NbDir <- lapply(obj.nb.angle, `[[`, 2)

    return(df)
  }

  if(modus == "nb")
  {
    return(obj.nb.angle)
  }


} # end of function


# # # # # # # #
# EXAMPLE ----------------------


# setwd("D:/Users/rAVer/Desktop/LSlide New")
#
# source("module_BoundingBox.R")
#
# # save.image("NBDirection.RData")
# load("NBDirection.RData")
#
# pacman::p_load("sf", "sp", "rgdal", "rgeos", "spdep", "maptools")
#
# # spdf <- rgdal::readOGR(dsn = getwd(), layer = "L3_seg_sel_out")
# spdf <- sf::st_read(dsn = paste0(getwd(), "/Input"), layer = "L3_seg_sel_out")
# # seg1.sf <- sf::st_transform(seg1.sf, sp::CRS(paste0("+init=", tolower(epsg.code)))@projargs)
#
# # transform to SP
# spdf <- as(spdf, 'Spatial')
# scarp <- subset(spdf, spdf@data$Class == 1) # 640
#
# # spdf.bb <- getBoundingBox(shape = scarp, col.name = "Flow")
#
#
#
# df <- NeighborDirection(spdf = spdf, spdf.bb = scarp, tol = 45, col.name = "Flow", modus = "bb", quiet = FALSE)
# nb.dir <- NeighborDirection(spdf = spdf, col.name = "Flow", modus = "nb", quiet = FALSE)
#
#
#
#
# # writeSpatialShape(x = bb, fn = paste(getwd(), "/Output/bb.shp", sep = "/"))
# rgdal::writeOGR(obj = spdf.PoS, dsn = paste0(getwd(), "/Output"), layer = "spdf_PoS", driver = "ESRI Shapefile")
# rgdal::writeOGR(obj = bb, dsn = paste0(getwd(), "/Output"), layer = "bb", driver = "ESRI Shapefile")