BoundaryDetection: Areal boundary detection in R

#' @useDynLib BoundaryDetection
#' @importFrom Rcpp sourceCpp
NULL

# Gaussian Smoothing (R implementation)
# http://homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm
# 'One method to remove noise is by convolving the original image with a mask that represents a low-pass filter or smoothing operation. For example, the Gaussian mask comprises elements determined by a Gaussian function. This convolution brings the value of each pixel into closer harmony with the values of its neighbors. In general, a smoothing filter sets each pixel to the average value, or a weighted average, of itself and its nearby neighbors; the Gaussian filter is just one possible set of weights.'
gaussian_filter <- function(x, cendroids, stdev = 0.1) {
    if (stdev == 0)
        return(out)
    out <- x
    sigma <- diag(rep(stdev^2, 2))
    # gaussian smoothing filter
    for(i in 1:length(x)) {
        # get euclidean vector
        euclidean_vector <- sweep(cendroids, 2, cendroids[i,], "-")
        dist <- sqrt(rowSums(euclidean_vector^2))
        # select polygons with distance below 3*stdev
        # In theory, the normal distribution is non-zero everywhere but in practice it is close to zero more than three standard deviations from the mean.
        sel <- dist < (3*stdev)
        if(sum(sel) > 1 & !is.na(x[i])) {
            # Gaussian Smoothing
            gauss <- mvtnorm::dmvnorm(euclidean_vector[sel,], c(0,0), sigma)
            out[i] <- weighted.mean(x[sel], gauss)
        }
    }
    return(out)
}

# Derivate of multivariate gaussian (R implementation)
dmvnorm_deriv <- function(X, mean, sigma) {
    if (is.vector(X)) X <- matrix(X, ncol = length(X))
    if (missing(mean)) mean <- rep(0, length = ncol(X))
    if (missing(sigma)) sigma <- diag(ncol(X))
    n <- nrow(X)
    mvnorm <- mvtnorm::dmvnorm(X, mean = mean, sigma = sigma)
    deriv <- array(NA,c(n,ncol(X)))
    for (i in 1:n)
        deriv[i,] <- -mvnorm[i] * solve(sigma,(X[i,]-mean))
    return(deriv)
}

#' Invert shapefile
#'
#' \code{invert_shapefile} inverts a \code{SpatialPolygonsDataFrame}. This function can be used to create a shapefile with natural boundaries/barriers.
#'
#' @param shp \code{SpatialPolygonsDataFrame} object (sf objects as converted).
#' @param extend Factor by which the bounding box is extended.
#' @return    \code{SpatialPolygonsDataFrame}
#'
#' @export
invert_shapefile <- function(shp, extend = 0.01) {
    if (is(shp, "sf")) shp <- as(shp, "Spatial")
    box         <- sp::bbox(shp)
    margin      <- (box[,2]-box[,1]) * extend
    proj4string <- sp::CRS(sp::proj4string(shp))
    coords <- rbind(
        c(box[1,1] - margin[1], box[2,1] - margin[2]),
        c(box[1,2] + margin[1], box[2,1] - margin[2]),
        c(box[1,2] + margin[1], box[2,2] + margin[2]),
        c(box[1,1] - margin[1], box[2,2] + margin[2]))
    p    <- sp::Polygons(list(sp::Polygon(coords)), ID = 'all')
    box  <- sp::SpatialPolygons(list(p), proj4string = proj4string)
    negative <- rgeos::gDifference(box, shp)
    return(negative)
    # plot(box); plot(shp, lwd = 0.1, add = TRUE)
}


#' Convert SpatialPolygonsDataFrame to SpatialLinesDataFrame with border segments
#'
#' \code{border_lines} converts a SpatialPolygonsDataFrame to a SpatialLinesDataFrame with one element for each border between neighbouring areas.
#'
#' @param sp Object of type \code{SpatialPolygonsDataFrame} (sf objects as converted)
#' @param longlat Use Euclidean or Great Circle distance for calculation of line length. If FALSE, Euclidean distance, if TRUE Great Circle distance in kilometers.
#' @return Object of class \code{SpatialLinesDataFrame} with one element for each border between neighbouring areas. 
#'
#' @importFrom magrittr "%<>%"
#' @importFrom magrittr "%>%"
#' @export
border_lines <- function(sp, longlat = TRUE) {
    # Coerce simple feature geometries to corresponding Spatial* objects
    if (is(sp, "sf")) sp <- as(sp, "Spatial")
    P  <- sp::polygons(sp)
    # get adjacency matrix A
    # nb <- spdep::poly2nb(sp, row.names = rownames(sp), queen = FALSE)
    # A  <- nb2mat::nb2mat(nb, style = "B", zero.policy = TRUE)
    nb <- spdep::poly2nb(sp, queen = FALSE)
    # create data.frame with adjacent areas
    greater_than <- function(a, b) a[a > b]
    data <- data.frame(i = 1:length(nb), j = NA) %>%
        group_by(i, j) %>%
        do(expand.grid(i = .$i, j = greater_than(nb[[.$i]], .$i))) %>%
        as.data.frame()
    # area borders as SpatialLines
    lines <- apply(data, 1, function(d) {
        i <- as.numeric(d["i"])
        j <- as.numeric(d["j"])
        # get list of coordinates for polygons
        c1 <- plyr::llply(P@polygons[[i]]@Polygons, sp::coordinates)
        c2 <- plyr::llply(P@polygons[[j]]@Polygons, sp::coordinates)
        # get borders for each combination of polygons
        grid <- expand.grid(s1 = 1:length(c1), s2 = 1:length(c2))
        line <- apply(grid, 1, function(obj) {
            a   <- c1[[obj["s1"]]]
            b   <- c2[[obj["s2"]]]
            # select intersecting rows
            sel <- a[, 1] %in% b[, 1] & a[, 2] %in% b[, 2]
            if(sum(sel) == 0) return(NULL)
            # create Line object for each sequence of matching coordinates
            runs <- rle(sel)
            runs <- data.frame(
                    val = runs$values,
                    i   = c(1, cumsum(runs$length) + 1)[-(length(runs$length) + 1)],
                    len = runs$length) %>%
                dplyr::filter(val)
            # coordinates for each sequence
            pos    <- plyr::alply(as.matrix(runs), 1, . %>% {.[["i"]] : (.[["i"]] + .[["len"]] - 1)})
            coords <- plyr::llply(pos, . %>% a[., , drop = FALSE])
            # remove duplicate line elements
            len_one <- plyr::laply(coords, nrow) == 1
            if (!all(len_one) & any(len_one)) {
                B <- do.call(rbind, coords)
                coords <- plyr::llply(coords, function(A) {
                    if(nrow(A) > 1) return(A)
                    cond <- sum(A[, 1] == B[, 1] & A[, 2] == B[, 2]) > 1
                    if(cond) return(NULL)
                    return(A)
                })
                coords <- coords[!sapply(coords, is.null)]
            }
            # return list of Line objects (one element for each sequence of coordinates)
            plyr::llply(coords, sp::Line)
        })
        segments <- unlist(line[!sapply(line, is.null)], recursive=FALSE)
        sp::Lines(segments, ID = sprintf("i%s_j%s", i, j))
    })
    sl <- sp::SpatialLines(lines[!sapply(lines, is.null)], proj4string = sp::CRS(sp::proj4string(sp)))
    # SpatialLinesDataFrame from SpatialLines and data
    sldf <- sp::SpatialLinesDataFrame(sl, data, match.ID = FALSE)
    # sldf$length <- SpatialLinesLengths(sldf, longlat = longlat)
    return(sldf)
}