# R/dividePolygon.R In ggRaver/Lslide: Functions for object-oriented image analysis in landslide science

#### Documented in dividePolygon

#' Devide Polygon into N Equal Area Parts
#'
#' This function divides a polygon feature into n-party of approximatly the same area.
#' The division is performed along the y-axis (N-S). Therefore, polygons should be rotated
#' appropriate.
#'
#'
#' @param x object of class sf
#' @param angle Can be one single value, or a vector with the same length as the input object x. North is 0 degree, and East 90. Default: 0 degree for North-South. For East-West use 90 degree.
#' @param n number of parts
#' @param trans angle input is transformed, so that N-S cutting is appropriate. Default: TRUE.
#' @param btol tolerance value for lower and upper boundary in optim().Default: 0.0001
#' @param pgtol tolerance value to reduce the sum of residuals in optimizing function
#' @param maxit maximum iteration to convergence in finding optimal points. Default: 100
#' @param cores number of cores for parallel processing. Default: 1 (sequential)
#' @param quiet no outputs in console. Default: TRUE
#' @return
#' Geometry similar to sf input including divisions.
#'
#'
#' @note
#' \itemize{
#'   \item Algorithm does not find optimum if polygon split results in more objects than the given number
#'          n. This can happen when a polygon is very ragged.
#'   \item Algorithm cut polygon on the N-S-Axis (0-180 degree). Therefore, trans is used to adapt the rotation angle
#'         so that cutting is appropriate.
#'
#' }
#'
#'
#' @keywords rainfall tresholds, rainfall event, landslide, automatic appraoch
#'
#'
#' @export
#'
dividePolygon <- function(x, angle = 0, n, trans = TRUE, btol = 0.0001, pgtol = 1e-9, maxit = 100, cores = 1, quiet = TRUE)
{

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

# # # # # # # #
## check input
# # # # # # # #
n.obj <- nrow(x) # number of features

# check angle input
if(length(angle) == 1)
{
angle <- rep(x = angle, n.obj) # create angle vector
} else{
if(length(angle) != n.obj){stop('Length of "angle" is not equal to "x" input features!')}
}

# transform angle input
if(trans == TRUE)
{
# 0 to 90 degree
angle <- ifelse(angle >= 0 & angle <= 90, 90 - angle, angle)

# 90 to 180 degree
angle <- ifelse(angle > 90 & angle < 180, 180 - (angle - 90), angle)

# 180 to 270 degree
angle <- ifelse(angle >= 180 & angle <= 270, 270 - angle, angle)

# 270 to 360 degree
angle <- ifelse(angle > 270 & angle <= 360, 180 - (angle - 270), angle)
}

# # # # # # # #
## create optimizer function
# # # # # # # #

# par - points of first lines guess
# x.all - x extent - min and max values
# n number of lines/parts
# x.i sf object of run i
opt.fun <- function(par, x.extent, n, x.i, splitByLines)
{
# browser()
## create lines
# get points to optimize
y.par <- par

# split sf object by lines
sfc.i <- splitByLines(x = x.extent, y = y.par, x.i = x.i)

# get area of splitted polygons and squared residuals
# sfc.i$A_RES <- ((sfc.i$A_FUNCT/n) - sfc.i$A_OPT)^2 sfc.i$A_RES <- ((sfc.i$A_FUNCT/nrow(sfc.i)) - sfc.i$A_OPT)^2

## sum of squared residuals to optimize
sum(sfc.i$A_RES) } #end of opt.fun # # # # # # # # ## rotate function (radiant) # # # # # # # # # x - radiant input rot <- function(x){matrix(c(cos(x), sin(x), -sin(x), cos(x)), 2, 2)} rot.inv <- function(x){matrix(c(cos(x), -sin(x), sin(x), cos(x)), 2, 2)} # # # # # # # # ## split sf object by lines # # # # # # # # # x - x coordiantes from extent # y - line coordiantes from guess # sf.i - sf-object to be split splitByLines <- function(x, y, x.i, final = FALSE) { # get line matrix m.xy <- expand.grid(x, y) m.lines <- lapply(X = seq(from = 1, to = n, by = 2), FUN = function(x, m.xy ){return(as.matrix(m.xy[c(x, x+1), ]))}, m.xy = m.xy) if(final) { # get lines lines <- sf::st_multilinestring(m.lines) x.i.dif <- list() x.i.run <- x.i for(i in 1:length(lines)) { # get line string from multiline ls.i <- sf::st_geometry(sf::st_linestring(lines[[i]])) split.i <- sf::st_collection_extract(x = lwgeom::st_split(x = x.i.run, y = ls.i), type = "POLYGON") if((i+1) <= length(lines)) { # get next line for intersection ls.i.next <- lines[[(i+1)]] # get interecting polygon that must be excluded from split inter <- unlist(sf::st_intersects(x = sf::st_linestring(ls.i.next), y = split.i)) n.split <- 1:nrow(split.i) # write result of run x.i.dif[[i]] <- split.i[!n.split %in% inter,] # set polygon to be cut in next run x.i.run <- split.i[n.split %in% inter,] } else{ # write result of run x.i.dif[[i]] <- split.i } } # end of for loop sfc.i <- do.call(rbind, x.i.dif) # combine final result of split } else { # get lines lines <- sf::st_sfc(sf::st_multilinestring(m.lines), crs = sf::st_crs(x.i)) ## buffer line to get thin polygons lines.buf <- sf::st_buffer(x = lines, dist = 0.000001) # create a very thin polygon ## splitting # calculate the differences of polygons and casting x.i.dif <- suppressWarnings(sf::st_difference(x = x.i, y = lines.buf)) if(class(x.i.dif$geometry)[1] == "sfc_GEOMETRYCOLLECTION"){
sfc.i <- sf::st_collection_extract(x = x.i.dif, "POLYGON")
} else {
sfc.i <- suppressWarnings(sf::st_cast(x = x.i.dif, to = "POLYGON"))
}

}

## calculate area
sfc.i$A_OPT <- as.numeric(sf::st_area(x = sfc.i)) return(sfc.i) } # end of function splitByLines # # # # # # # # ## divide features # # # # # # # # # calculate area x$A_FUNCT <- as.numeric(sf::st_area(x))
x.cntrd <- sf::st_centroid(x)

# get crs
x.crs <- sf::st_crs(x)

# init parallel
# 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('x', 'x.crs', 'angle', 'n', 'btol','pgtol', 'maxit', 'x.cntrd', 'rot', 'rot.inv', 'opt.fun', 'splitByLines'))

parallel::clusterEvalQ(cl, library("sf", "lwgeom"))
}

# x.div <- parallel::mclapply(mc.cores = mc.cores, X = 1:n.obj, fun = function(i, x, angle, n, btol, pgtol, maxit, x.cntrd, rot, opt.fun, splitByLines){
x.div <- parallel::parLapply(cl = cl, X = 1:n.obj, fun = function(i, x.sf, crs, angle, n, btol, pgtol, maxit, x.cntrd, rot, rot.inv, opt.fun, splitByLines){
# x.div <- lapply(X = 1:n.obj, FUN = function(i, x.sf, crs, angle, n, btol, pgtol, maxit, x.cntrd, rot, rot.inv, opt.fun, splitByLines){
# browser()
# get objects
x.i <- x.sf[i,]

# check if object is valid, and make it valuid
if(suppressWarnings(!all(sf::st_is_valid(x.i)))){
x.i <- lwgeom::st_make_valid(x.i)

if(class(x.i$geometry)[1] == "sfc_GEOMETRYCOLLECTION") { x.i <- sf::st_collection_extract(x = x.i, type = "POLYGON") } } angle.i <- angle[i] x.cntrd.i <- sf::st_geometry(x.cntrd[i,]) # rotate object x.i.rot <- (sf::st_geometry(x.i) - x.cntrd.i) * rot(angle.i*pi/180) + x.cntrd.i x.i.rot <- sf::st_sf(sf::st_set_geometry(x.i, NULL) , geometry = x.i.rot, crs = sf::st_crs(x.i)) if(suppressWarnings(!all(sf::st_is_valid(x.i.rot)))){ x.i.rot <- lwgeom::st_make_valid(x.i.rot) if(class(x.i.rot$geometry)[1] == "sfc_GEOMETRYCOLLECTION")
{
x.i.rot <- sf::st_collection_extract(x = x.i.rot, type = "POLYGON")
}
}

# test
# plot(sf::st_geometry(x.i))

# get bounding box and start values
bb.i <- sf::st_bbox(x.i.rot)

y.start <- bb.i[2]
y.end <- bb.i[4]
y.par <- seq(from = y.start, to = bb.i[4], length.out = n + 1)[2:n]
x.extent <- c(bb.i[1], bb.i[3])

# # #
# search for optimized splitting lines
# # #
y.par.opt <- optim(par = y.par, fn = opt.fun, method = "L-BFGS-B",
control = list(pgtol = pgtol, maxit = maxit),
lower = (y.start + btol), upper = (y.end - btol),
x.extent = x.extent, n = n, x.i = x.i.rot,
splitByLines = splitByLines)$par # # # # split polygon # # # x.i.rot.div <- splitByLines(x = x.extent, y = y.par.opt, x.i = x.i.rot, final = TRUE) # x.i.rot.div <- lwgeom::st_make_valid(sf::st_buffer(x.i.rot.div, dist = 0.000001)) if(suppressWarnings(!all(sf::st_is_valid(x.i.rot.div)))) { x.i.rot.div <- lwgeom::st_make_valid(x.i.rot.div) } # sf::st_intersects(x.i.rot.div) # plot(st_geometry(sf::st_intersection(x.i.rot.div))) # rotate back x.i.div <- (sf::st_geometry(x.i.rot.div) - x.cntrd.i) * rot.inv(angle.i*pi/180) + x.cntrd.i x.i.div <- sf::st_sf(sf::st_set_geometry(x.i.rot.div, NULL) , geometry = x.i.div, crs = sf::st_crs(x.i.rot.div)) # set crs x.i.div <- sf::st_set_crs(x.i.div, crs) # set unique ID column x.i.div$ID_Split <- paste0(i, "_", c(1:nrow(x.i.div)))

# check last time geometry
if(suppressWarnings(!all(sf::st_is_valid(x.i.div))))
{
x.i.div <- lwgeom::st_make_valid(x.i.div)

if(class(x.i.div\$geometry)[1] == "sfc_GEOMETRYCOLLECTION")
{
x.i.div <- sf::st_collection_extract(x = x.i.div, type = "POLYGON")
}
}

# test
# plot(sf::st_geometry(x.i.div))

return(x.i.div)

}, x.sf = x, crs = x.crs, angle = angle, n = n, btol = btol, pgtol = pgtol, maxit = maxit, x.cntrd = x.cntrd, rot = rot, rot.inv = rot.inv, opt.fun = opt.fun, splitByLines = splitByLines)

parallel::stopCluster(cl)

# rbind simple features
x.div <- do.call(rbind, x.div)
# plot(sf::st_geometry(x.div))

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

return(x.div)
} # end of function dividePolygon

# setwd("d:/Users/qo23hel/Desktop/Test")