fix: various improvements to avoid erroneous polygons and catch errors

R/allocate.R

@@ -1,7 +1,7 @@
 #' @importFrom stats rnorm
 #' @importFrom dplyr %>%
 #' @importFrom sf st_area
-#' 
+#'
 cellError <- function(a, target) {
   normA <- a / sum(a)
   diff <- abs(normA - target)
@@ -14,9 +14,9 @@ breaking <- function(
   debug = FALSE,
   error_tol,
   prevError) {
-  
+
   if (max) {
-    
+
     # Stop when largest individual cell error is less than 1%
     # (the default)
     err <- cellError(a, target)
@@ -29,7 +29,7 @@ breaking <- function(
     prevError <- err
 
   } else {
-    
+
     normA <- a / sum(a)
     diff <- abs(normA - target)
     # Stop when *change* in *total* cell error is tiny
@@ -46,7 +46,7 @@ breaking <- function(
       ))
     stopping <- abs(sum(diff) - prevError) < 0.001
     prevError <- sum(diff)
-    
+
   }
   list(
     stopping = stopping,
@@ -59,7 +59,7 @@ breaking <- function(
 # adjust by multiple of average absolute weights
 # This avoids problem of getting stuck at a tiny weight
 # (and stabilizes the algorithm generally)
-# difference to original implementation: adjustment of maximal 
+# difference to original implementation: adjustment of maximal
 # step change of weights to prevent crashing of algorithm
 adjustWeights <- function(w, a, target) {
   # OPTION: avoid extreme scaling values -> squareroot function
@@ -115,34 +115,51 @@ shiftWeights <- function(s, w) {
 # just give up after 'maxIteration's
 allocate <- function(
   names, s, w, outer, target,
-  maxIteration, 
+  maxIteration,
   error_tol,
-  debug = FALSE, 
+  min_target = 0.01,
+  debug = FALSE,
   debugCell = FALSE)
 {
   count <- 1
   prevError <- 1
-
+
+  # check for extremely small cell size compared to theoretical average
+  target_fc <- target * length(target)
+  too_small <- target_fc < min_target
+  if (any(too_small)) {
+    message(paste0(
+      "Found extremely small cell (<", round(min_target * 100, 1), "% of average size);\n",
+      "inflating cell size to prevent failure when calculating polygons."
+    ))
+    correction <- ifelse(
+      too_small,
+      min_target/length(target),
+      -min_target/length(target) * sum(too_small)/sum(!too_small)
+    )
+    target <- target + correction
+  }
+
   repeat {
-    
+
     # if all weights are identical the CGAL algorithm often fails
     # in this case we introduce a bit of random variation
     if (length(unique(w)) == 1) {
       w <- w * rnorm(length(w), mean = 1, sd = 0.01)
     }
-    
+
     # call to awv function, the additively weighted voronoi tesselation,
     # wrapped within a trycatch statement to catch errors and start over
     k <- tryCatch(awv(s, w, outer, debug, debugCell),
-      error = function(e) { print(e); NULL}
+      error = function(e) { message(e); NULL}
     )
     if (is.null(k)) {
       return(NULL)
     }
     areas <- lapply(k, sf::st_area)
-    
+
     # if debug=TRUE, every iteration is drawn to the viewport
-    # this can be very time and resource consuming and should be used 
+    # this can be very time and resource consuming and should be used
     # with care. The result resembles the final treemap but is an overlay of
     # many iterations
     if (debug) {
@@ -156,7 +173,7 @@ allocate <- function(
         lwd = 2, col = grey(0.5),
         fill = grey(1, alpha=0.33)
       )
-      
+
       info <-
         rbind(
           area = round(unlist(areas) / sum(unlist(areas)), 4),
@@ -167,30 +184,30 @@ allocate <- function(
       colnames(info) <- names
       print(info)
     }
-    
+
     stop_cond <- breaking(
-      unlist(areas), 
-      target, 
+      unlist(areas),
+      target,
       debug = debug,
       error_tol = error_tol,
       prevError = prevError)
-    
+
     # if stop condition is fulfilled, return result in form of
     # list of polygons and metadata
     if (count == maxIteration || stop_cond$stopping) {
-      
+
       res <- lapply(1:length(names), function(i) {
         list(
           name = names[i], poly = k[[i]],
-          site = c(s$x[[i]], s$y[[i]]), 
-          weight = w[i], area = unlist(areas)[i], 
+          site = c(s$x[[i]], s$y[[i]]),
+          weight = w[i], area = unlist(areas)[i],
           target = target[i],
           count = count)
       }) %>% setNames(names)
       return(res)
-    
+
     } else {
-      
+
       w <- adjustWeights(w, unlist(areas), target)
       s <- shiftSites(s, k)
       w <- shiftWeights(s, w)

R/tesselation.R

@@ -11,7 +11,7 @@
 # Generate one iteration of the Additively Weighted Voronoi diagram
 awv <- function(
   s, w, region, debug = FALSE,
-  debugCell = FALSE) 
+  debugCell = FALSE)
 {
   # combine X, Y coordinates and weights as input for
   # C++ tesselation function
@@ -30,27 +30,29 @@ awv <- function(
 
 tidyCell <-
   function(cell, tolerance) {
-    
+
     # if cell touches the border at two points, we need to close it
     # this is not necessary if cell touches border at 4 points (like a stripe)
     if (sum(
-      cell$border$x %in% c(4000, -4000), 
+      cell$border$x %in% c(4000, -4000),
       cell$border$y %in% c(4000, -4000)) == 2
     ) {
-
-      closeCell(cell$border, cell$vertex, tol = tolerance)
-
+      poly <- closeCell(cell$border, cell$vertex, tol = tolerance)
     } else {
-
       # return a list of the polygon
-      list(
+      poly <- list(
         x = cell$border$x,
         y = cell$border$y,
         end = "boundary"
       )
-
+      if (sf::st_is_valid(convertCell(poly[1:2]))) {
+        return(poly)
+      } else {
+        message("Found invalid polygon (self-intersection)")
+        return(NULL)
+      }
     }
-    
+
   }
 
 # SIDES
@@ -97,7 +99,7 @@ antiSide <- function(corner) {
 closeClock <- function(x, y, start, end, scale = 2000) {
   cornerX <- c(-2 * scale, 2 * scale, 2 * scale,-2 * scale)
   cornerY <- c(2 * scale, 2 * scale,-2 * scale,-2 * scale)
-  
+
   side <- end
   repeat {
     corner <- clockCorner(side)
@@ -116,7 +118,7 @@ closeClock <- function(x, y, start, end, scale = 2000) {
 closeAnti <- function(x, y, start, end, scale = 2000) {
   cornerX <- c(-2 * scale, 2 * scale, 2 * scale,-2 * scale)
   cornerY <- c(2 * scale, 2 * scale,-2 * scale,-2 * scale)
-  
+
   side <- end
   repeat {
     corner <- antiCorner(side)
@@ -155,40 +157,36 @@ closeCell <- function(cell, vertex, tol, scale = 2000) {
   # If not, do second one (and check that vertex is "inside" that result!)
   x <- cell$x
   y <- cell$y
-  
+
   # ASSUME that both first and last vertices are on boundary!
   N <- length(x)
   startSide <- side(x[1], y[1])
   endSide <- side(x[N], y[N])
-  
+
   # exit if not both end points lie on boundary
   if (length(startSide) != 1 | length(endSide) != 1) {
     return(NULL)
   }
-    
+
   # Start and end on same side
   if (startSide == endSide) {
-
-    cell <- list(x = x, y = y)
+
     if (sp::point.in.polygon(vertex[1], vertex[2],
-                         cell$x, cell$y) == 0) {
+        cell$x, cell$y) == 0) {
       boundRect <- to_sfpoly(list(
         x = c(-2 * scale,-2 * scale, 2 * scale, 2 * scale),
         y = c(-2 * scale, 2 * scale, 2 * scale,-2 * scale)
       ))
       # "Subtract" smallCell from bound rect to get largeCell
       cellPoly <- to_sfpoly(cell)
       cellPoly <- sf::st_difference(boundRect, cellPoly)
-
-      pts <- to_coords(cellPoly)
+      cell <- to_coords(cellPoly)
       if (sp::point.in.polygon(vertex[1], vertex[2],
-                           cell$x, cell$y) == 0) {
+          cell$x, cell$y) == 0) {
         stop("Failed to close cell")
       }
     }
-
   } else {
-
     cell <- closeClock(x, y, startSide, endSide)
     if (sp::point.in.polygon(vertex[1], vertex[2],
                          cell$x, cell$y) == 0) {
@@ -222,6 +220,10 @@ trimCells <- function(cells, region) {
     if (inherits(poly, "MULTIPOLYGON")) {
       poly <- suppressWarnings(sf::st_cast(poly, to = "POLYGON"))
     }
+    if (inherits(poly, "GEOMETRYCOLLECTION")) {
+      valid <- which(sapply(poly, function(x) inherits(x, "POLYGON")))[1]
+      poly <- poly[[valid]]
+    }
     poly
   })
 }
@@ -235,21 +237,21 @@ samplePoints <- function(ParentPoly, n, seed, positioning) {
   if (!is.null(seed)) {
     set.seed(seed)
   }
-  
-  # This loop keeps repeating until the correct number of coordinates 
+
+  # This loop keeps repeating until the correct number of coordinates
   # is sampled. The reason is that sp::spsample() does not always sample
   # the correct number of coordinates, but too few or too many
   repeat {
-    
+
     sampled <- tryCatch({
       points <- sp::spsample(
         sp::Polygon(coords = ParentPoly),
-        n = n, 
+        n = n,
         type = ifelse(positioning == "random", "random", "nonaligned")
       )
       points@coords}, error = function(e) NULL
     )
-    
+
     if (is.null(sampled) || nrow(sampled) != n) {
       next
     } else {