R/areal_edge_detection.R

In jlegewie/BoundaryDetection: Areal boundary detection in R

# scale mask so that positive and negative values have equal weights
scale_mask <- function(mask) {
    pos <- mask > 0
    mask[ pos[,1],1]  <-  1/sum(mask[ pos[,1],1]) * mask[ pos[,1],1]
    mask[!pos[,1],1]  <- -1/sum(mask[!pos[,1],1]) * mask[!pos[,1],1]
    mask[ pos[,2],2]  <-  1/sum(mask[ pos[,2],2]) * mask[ pos[,2],2]
    mask[!pos[,2],2]  <- -1/sum(mask[!pos[,2],2]) * mask[!pos[,2],2]
    mask[is.na(mask)] <- 0
    return(mask)
}

# adapt scale parameter to local area size
stddev_local <- function(i, N, C) {
    sel <- N[[i]]
    if(length(sel) <= 1) return(NA)
    vec <- sweep(C[sel,], 2, C[i,], "-")
    x   <- vec[which.max(abs(vec)[,1] * (1-(abs(vec)[,2]/max(abs(vec)[,2])))), 1]
    y   <- vec[which.max(abs(vec)[,2] * (1-(abs(vec)[,1]/max(abs(vec)[,1])))), 2]
    return(abs(c(x, y)))
}

# crossing of natural barrier
crossing <- function(shp, i, mu, valid, ip, natural_boundaries, nb_fortify) {
    blocked <- rep(TRUE, sum(valid))
    # distance to natural boundary
    vec_nb  <- sweep(nb_fortify[,1:2], 2, mu, "-")
    dist_nb <- sqrt(rowSums(vec_nb^2))
    vec_ip  <- sweep(ip[valid,], 2, mu, "-")
    dist_ip <- sqrt(rowSums(vec_ip^2))
    blocked[dist_ip < min(dist_nb)] = FALSE
    if(all(!blocked)) return(valid)
    # intersection with natural boundary
    width <- -3e-5 * max(bbox(shp[i,])[,2] - bbox(shp[i,])[,1])
    area <- rgeos::gBuffer(shp[i,], width = width, capStyle = "FLAT")
    # plot(shp[i,], border = "red"); plot(area, add=TRUE)
    coords_poly <- ggplot2::fortify(area)
    for (j in 1:nrow(coords_poly)) {
        vertex   <- as.matrix(coords_poly[j, 1:2])
        coords_subset <- ip[valid, , drop = FALSE][blocked, , drop = FALSE]
        if(nrow(coords_subset) == 0) break
        lines <- sapply(1:nrow(coords_subset),
            function(j) Lines(Line(rbind(vertex, coords_subset[j,])), ID = paste0(j)))
        sl    <- sp::SpatialLines(lines, proj4string = CRS(proj4string(shp)))
        sel   <- !rgeos::gIntersects(natural_boundaries, sl, byid = TRUE)[,1]
        blocked[blocked][sel] <- FALSE
    }
    valid[valid][blocked]  <- FALSE
    return(valid)
}


#' Edge detection for areal data
#'
#' \code{areal_edge_detection} ...
#'
#' @param shp     \code{SpatialPolygonsDataFrame} object (sf objects as converted).
#' @param vars    Vector of variable names.
#' @param scale_smoothing  Scale parameter for gaussian smoothing to pre-process data.
#'   \code{scale_smoothing} corresponds to the standard deviation of the normal
#'   distribution. Note that the scale depends on the projection of the shapefile.
#'   0 (the default) uses no smoothing.
#' @param scale_sobel  Scale parameter for the edge detection algorithm.
#'   The parameter corresponds to the standard deviation of the first derivative
#'   of the multivariate normal distribution (the edge detection filter).
#'   Note that the scale depends on the projection of the shapefile. 0 (the default)
#'   automatically adapts to cover the adjacent areas.
#' @param domain  Domain for edge detection algorithms: \code{neighbors} or
#'   \code{unbounded}. \code{neighbors} just uses neighboring areas. \code{unbounded}
#'   uses all areas within 3 times the standard deviation. 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.
#' @param natural_boundaries     natural_boundaries
#' @param N      Adjacency matrix as neighbours list based on \code{\link[spdep]{poly2nb}} (optional, saves time for multiple calls).
#' @param C      SpatialPoints object with coordinates of cendroids (optional, saves time for multiple calls).
#' @param method    Centroid or numerical integration.
#' @param smoothing    Smoothing method. 'gaussian' and 'bilateral' are supported.
#'   'bilateral' requires two scale parameter.
#' @param combine    Function to combine edges based on different variables.
#'   \code{NULL} (the default) returns an array with the edges for each variables.
#' @param res    Resolution for the numerical integration. The number of integration
#'   points is \code{res * res}, which is directly related to speed and precision.
#' @param .progress    Show process bar. 
#' @return       Numerical vector or array with edge intensity variable.
#'
#' @examples
#' areal_edge_detection(...)
#' areal_edge_detection(..., combine = function(x) max(x, na.rm = TRUE))
#' areal_edge_detection(..., combine = function(x) prod(tail(sort(x), 2)))
#'
#' @export
areal_edge_detection <- function(shp, vars, scale_smoothing = 0, scale_sobel = 0,
        domain = "neighbors", natural_boundaries = NULL, N = NULL, C = NULL,
        method = "integration", smoothing = "gaussian", combine = NULL, res = 50, .progress = "none") {
    # check arguments
    if(smoothing == "bilateral" & length(scale_smoothing) != 2)
        stop("Bilateral smoothing requires two scale parameters but length of 'scale_smoothing' is not 2.")
    if(!is.null(N)) if(length(N) != nrow(shp))
        stop("The number of rows and columns in the adjacency matrix 'N' has to equal to the number of units in 'shp'")
    if(!is.null(C)) if(length(C) != nrow(shp))
        stop("The number of units in 'C' has to be equal the number of units in 'shp'")
    # Coerce simple feature geometries to corresponding Spatial* objects
    if (is(shp, "sf")) shp <- as(shp, "Spatial")
    # adjacency matrix, centroids, and fortify natural_boundaries
    if(is.null(N) & (scale_sobel == 0 | domain == "neighbors"))
        N <- spdep::poly2nb(shp)
    if(is.null(C))
        C <- rgeos::gCentroid(shp, byid = TRUE)
    if(!is.null(natural_boundaries))
        nb_fortify <- ggplot2::fortify(natural_boundaries)
    cendroids <- sp::coordinates(C)
    # gaussian smoothing filter (cpp function with r fallback)
    if(scale_smoothing[1] != 0 & smoothing == "gaussian")
        shp@data[,vars] <- apply(shp@data[,vars, drop = FALSE], 2, function(x) gaussian_filter_arma(x, cendroids, scale_smoothing[1]))
    if(scale_smoothing[1] != 0 & smoothing == "bilateral")
        shp@data[,vars] <- apply(shp@data[,vars, drop = FALSE], 2, function(x) bilateral_filter_arma(x, cendroids, scale_smoothing))
    # scale parameter options
    stdev <- rep(scale_sobel, 2)
    # convolution of mask with map
    fail        <- rep(NA, length(vars))
    names(fail) <- vars
    edge <- plyr::aaply(1:nrow(shp), 1, function(i) {
        mu <- cendroids[i,]
        # get standard deviation
        if(scale_sobel == 0) stdev <- stddev_local(i, N, cendroids)
        if(is.na(stdev[1])) return(fail)
        # domain: unbounded
        if(domain == "unbounded") {
            euc_vec   <- sweep(cendroids, 2, mu, "-")
            dist      <- sqrt(rowSums(euc_vec^2))
            sel_areas <- which(dist != 0 & dist < (5 * max(stdev)))
            if (length(sel_areas) == 0) {
                warning('No other areas fall in the kernel domain of area ', i, '.')
                return(rep(0, length(vars)))
            }
        }
        # domain: neighbors
        if(domain == "neighbors") sel_areas <- N[[i]]
        if(length(sel_areas) <= 1) {
            warning('Area ', i, ' has too few neighboring areas.')
            return(fail)
        }
        # integration points
        if(method == "integration") {
            x  <- seq(mu[1] - 3 * stdev[1], mu[1] + 3 * stdev[1], length.out = res)
            y  <- seq(mu[2] - 3 * stdev[2], mu[2] + 3 * stdev[2], length.out = res)
            ip <- as.matrix(expand.grid(x = x, y = y))
            p  <- sp::SpatialPoints(ip, sp::CRS(sp::proj4string(shp)))
            # get data values at sample points (only close areas to improve perf of `over`)
            values <- sp::over(p, shp[sel_areas, vars])
            valid  <- complete.cases(values)
            if(sum(valid) <= 1) return(fail)
            # crossing of natural barrier
            if(!is.null(natural_boundaries))
                valid <- crossing(shp, i, mu, valid, ip, natural_boundaries, nb_fortify)
        }
        if(method == "centroids") {
            ip     <- cendroids[sel_areas,]
            values <- shp@data[sel_areas, vars]
            valid  <- complete.cases(values)
        }
        if(sum(valid) <= 1) return(fail)
        # get weights at sample points
        mask <- dmvnorm_deriv_arma(ip[valid, , drop = FALSE], mu, diag(stdev^2))
        # balance mask (by adding point in focal area)
        mx   <- mask[,1]
        r    <- sum(mx[mx > 0]) / sum(abs(mx[mx < 0]))
        if(!dplyr::between(r, 1 / 3, 3)) {
            w      <- sum(abs(mx[mx < 0])) - sum(mx[mx > 0])
            mask   <- rbind(mask, c(w, 0))
            values <- rbind(values, shp@data[i, vars])
            valid  <- c(valid, TRUE)
        }
        # mx <- mask[,1]; sum(mx[mx>0])/sum(abs(mx[mx<0]))
        my  <- mask[,2]
        r   <- sum(my[my > 0]) / sum(abs(my[my < 0]))
        if(!dplyr::between(r, 1 / 3, 3)) {
            w      <- sum(abs(my[my < 0])) - sum(my[my > 0])
            mask   <- rbind(mask, c(0, w))
            values <- rbind(values, shp@data[i, vars])
            valid  <- c(valid, TRUE)
        }
        # scale mask so that positive and negative values sum to 1 and -1 respectively
        mask <- scale_mask(mask)
        # calculate filter values
        dx <- colSums(apply(values[valid, , drop = FALSE], 2, "*", y = mask[,1]))
        dy <- colSums(apply(values[valid, , drop = FALSE], 2, "*", y = mask[,2]))
        # combine x and y filter
        d <- sqrt(dx^2 + dy^2)
        # return
        return(d)
    }, .progress = .progress)
    # return
    if(!is.null(combine))
        edge <- apply(edge, 1, combine)
    return(edge)
}