create clean voronoi point order PR #1371 #2426

NEWS.md

@@ -8,6 +8,8 @@
 
 * the `sf` -> `ppp` conversion `as.ppp.sf()` accepts a data.frame of marks instead of just 1 column #2450, by @agila5
 
+* add flag for preservation of point order in `st_voronoi` #1371 for GEOS >= 3.12
+
 # version 1.0-18
 
 * support  `POLYGON FULL` simple feature geometry, representing the entire Earth surface, as used by `s2geometry`; see also https://r-spatial.org/r/2024/10/11/polygonfull.html for a longer introduction; #2441

R/geom-transformers.R

@@ -455,6 +455,7 @@ st_minimum_rotated_rectangle.sf = function(x, dTolerance, ...) {
 #' @name geos_unary
 #' @export
 #' @param envelope object of class \code{sfc} or \code{sfg} containing a \code{POLYGON} with the envelope for a voronoi diagram; this only takes effect when it is larger than the default envelope, chosen when \code{envelope} is an empty polygon
+#' @param point_order logical; preserve point order if TRUE and GEOS version >= 3.12; overrides bOnlyEdges
 #' @details \code{st_voronoi} creates voronoi tessellation. \code{st_voronoi} requires GEOS version 3.5 or above
 #' @examples
 #' set.seed(1)
@@ -475,33 +476,48 @@ st_minimum_rotated_rectangle.sf = function(x, dTolerance, ...) {
 #'  # compute Voronoi polygons:
 #'  pols = st_collection_extract(st_voronoi(do.call(c, st_geometry(pts))))
 #'  # match them to points:
-#'  pts$pols = pols[unlist(st_intersects(pts, pols))]
+#'  pts_pol = st_intersects(pts, pols)
+#'  pts$pols = pols[unlist(pts_pol)] # re-order
+#'  if (isTRUE(try(compareVersion(sf_extSoftVersion()["GEOS"], "3.12.0") > -1,
+#'    silent = TRUE))) {
+#'    pols_po = st_collection_extract(st_voronoi(do.call(c, st_geometry(pts)),
+#'      point_order = TRUE)) # GEOS >= 3.12 can preserve order of inputs
+#'    pts_pol_po = st_intersects(pts, pols_po)
+#'    print(all(unlist(pts_pol_po) == 1:(n/2)))
+#'  }
 #'  plot(pts["id"], pch = 16) # ID is color
 #'  plot(st_set_geometry(pts, "pols")["id"], xlim = c(0,1), ylim = c(0,1), reset = FALSE)
 #'  plot(st_geometry(pts), add = TRUE)
 #'  layout(matrix(1)) # reset plot layout
 #' }
-st_voronoi = function(x, envelope, dTolerance = 0.0, bOnlyEdges = FALSE)
+st_voronoi = function(x, envelope, dTolerance = 0.0, bOnlyEdges = FALSE, point_order = FALSE)
 	UseMethod("st_voronoi")
 
 #' @export
-st_voronoi.sfg = function(x, envelope = st_polygon(), dTolerance = 0.0, bOnlyEdges = FALSE)
-	get_first_sfg(st_voronoi(st_sfc(x), st_sfc(envelope), dTolerance, bOnlyEdges = bOnlyEdges))
+st_voronoi.sfg = function(x, envelope = st_polygon(), dTolerance = 0.0, bOnlyEdges = FALSE, point_order = FALSE)
+	get_first_sfg(st_voronoi(st_sfc(x), st_sfc(envelope), dTolerance, bOnlyEdges = bOnlyEdges, point_order = point_order))
 
 #' @export
-st_voronoi.sfc = function(x, envelope = st_polygon(), dTolerance = 0.0, bOnlyEdges = FALSE) {
+st_voronoi.sfc = function(x, envelope = st_polygon(), dTolerance = 0.0, bOnlyEdges = FALSE, point_order = FALSE) {
 	if (compareVersion(CPL_geos_version(), "3.5.0") > -1) {
 		if (isTRUE(st_is_longlat(x)))
 			warning("st_voronoi does not correctly triangulate longitude/latitude data")
+		if (compareVersion(CPL_geos_version(), "3.12.0") > -1) {
+			bOnlyEdges = as.integer(bOnlyEdges)
+			if (point_order) bOnlyEdges = 2L # GEOS enum GEOS_VORONOI_PRESERVE_ORDER
+		} else {
+			if (point_order) 
+				warning("Point order retention not supported for GEOS ", CPL_geos_version())
+		}
 		st_sfc(CPL_geos_voronoi(x, st_sfc(envelope), dTolerance = dTolerance,
 			bOnlyEdges = as.integer(bOnlyEdges)))
 	} else
 		stop("for voronoi, GEOS version 3.5.0 or higher is required")
 }
 
 #' @export
-st_voronoi.sf = function(x, envelope = st_polygon(), dTolerance = 0.0, bOnlyEdges = FALSE) {
-	st_set_geometry(x, st_voronoi(st_geometry(x), st_sfc(envelope), dTolerance, bOnlyEdges))
+st_voronoi.sf = function(x, envelope = st_polygon(), dTolerance = 0.0, bOnlyEdges = FALSE, point_order = FALSE) {
+	st_set_geometry(x, st_voronoi(st_geometry(x), st_sfc(envelope), dTolerance, bOnlyEdges = bOnlyEdges, point_order = point_order))
 }
 
 #' @name geos_unary

src/geos.cpp

@@ -35,6 +35,9 @@
 # if GEOS_VERSION_MINOR >= 11
 #  define HAVE311
 # endif
+# if GEOS_VERSION_MINOR >= 12
+#  define HAVE312
+# endif
 #else
 # if GEOS_VERSION_MAJOR > 3
 #  define HAVE340

tests/geos.R

@@ -73,14 +73,26 @@ st_line_merge(mls)
 if (isTRUE(try(compareVersion(sf_extSoftVersion()["GEOS"], "3.5.0") > -1, silent = TRUE))) {
  # voronoi:
  set.seed(1)
- x = st_multipoint(matrix(runif(10),,2))
+ m = matrix(runif(10),,2)
+ x = st_multipoint(m)
  box = st_polygon(list(rbind(c(0,0),c(1,0),c(1,1),c(0,1),c(0,0))))
  v = st_sfc(st_voronoi(x, st_sfc(box)))
  plot(v, col = 0, border = 1, axes = TRUE)
  plot(box, add = TRUE, col = 0, border = 1) # a larger box is returned, as documented
  plot(x, add = TRUE, col = 'red', cex=2, pch=16)
  plot(st_intersection(st_cast(v), box)) # clip to smaller box
  plot(x, add = TRUE, col = 'red', cex=2, pch=16)
+ v0 = st_sfc(st_voronoi(st_sfc(x), st_sfc(box)))
+ pal <- c("black", "red", "green", "blue", "orange")
+ opar = par(mfrow=c(1,2))
+ plot(st_collection_extract(v0, "POLYGON"), col=pal)
+ text(m[,1], m[,2], label=1:5, col="white")
+ if (isTRUE(try(compareVersion(sf_extSoftVersion()["GEOS"], "3.12.0") > -1, silent = TRUE))) {
+  v2 = st_sfc(st_voronoi(st_sfc(x), st_sfc(box), point_order=TRUE))
+  plot(st_collection_extract(v2, "POLYGON"), col=pal)
+  text(m[,1], m[,2], label=1:5, col="white")
+ }
+ par(opar)
 
  v = st_voronoi(x)
  print(class(v))

tests/geos.Rout.save

@@ -1,7 +1,7 @@
 
-R version 4.2.2 Patched (2022-11-10 r83330) -- "Innocent and Trusting"
-Copyright (C) 2022 The R Foundation for Statistical Computing
-Platform: x86_64-pc-linux-gnu (64-bit)
+R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+Copyright (C) 2024 The R Foundation for Statistical Computing
+Platform: x86_64-pc-linux-gnu
 
 R is free software and comes with ABSOLUTELY NO WARRANTY.
 You are welcome to redistribute it under certain conditions.
@@ -162,14 +162,26 @@ LINESTRING (0 0, 1 1, 2 0)
 > if (isTRUE(try(compareVersion(sf_extSoftVersion()["GEOS"], "3.5.0") > -1, silent = TRUE))) {
 +  # voronoi:
 +  set.seed(1)
-+  x = st_multipoint(matrix(runif(10),,2))
++  m = matrix(runif(10),,2)
++  x = st_multipoint(m)
 +  box = st_polygon(list(rbind(c(0,0),c(1,0),c(1,1),c(0,1),c(0,0))))
 +  v = st_sfc(st_voronoi(x, st_sfc(box)))
 +  plot(v, col = 0, border = 1, axes = TRUE)
 +  plot(box, add = TRUE, col = 0, border = 1) # a larger box is returned, as documented
 +  plot(x, add = TRUE, col = 'red', cex=2, pch=16)
 +  plot(st_intersection(st_cast(v), box)) # clip to smaller box
 +  plot(x, add = TRUE, col = 'red', cex=2, pch=16)
++  v0 = st_sfc(st_voronoi(st_sfc(x), st_sfc(box)))
++  pal <- c("black", "red", "green", "blue", "orange")
++  opar = par(mfrow=c(1,2))
++  plot(st_collection_extract(v0, "POLYGON"), col=pal)
++  text(m[,1], m[,2], label=1:5, col="white")
++  if (isTRUE(try(compareVersion(sf_extSoftVersion()["GEOS"], "3.12.0") > -1, silent = TRUE))) {
++   v2 = st_sfc(st_voronoi(st_sfc(x), st_sfc(box), point_order=TRUE))
++   plot(st_collection_extract(v2, "POLYGON"), col=pal)
++   text(m[,1], m[,2], label=1:5, col="white")
++  }
++  par(opar)
 + 
 +  v = st_voronoi(x)
 +  print(class(v))
@@ -666,4 +678,4 @@ CRS:           NA
 > 
 > proc.time()
    user  system elapsed 
- 19.725   0.860  19.747 
+ 11.890   0.081  12.018