From 2baaec75525dea97b7f0bda5beeed52036469608 Mon Sep 17 00:00:00 2001 From: Roger Bivand Date: Thu, 12 Sep 2024 10:54:24 +0200 Subject: [PATCH] create clean voronoi point order PR #1371 #2426 --- NEWS.md | 2 ++ R/geom-transformers.R | 30 +++++++++++++++++++++++------- man/geos_unary.Rd | 14 ++++++++++++-- src/geos.cpp | 3 +++ tests/geos.R | 14 +++++++++++++- tests/geos.Rout.save | 22 +++++++++++++++++----- 6 files changed, 70 insertions(+), 15 deletions(-) diff --git a/NEWS.md b/NEWS.md index 12add0430..c0d276a6f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,7 @@ # version 1.0-18 +* add flag for preservation of point order in `st_voronoi` #1371 for GEOS >= 3.12 + * fix build failure with GDAL < 3.4.0 #2436 # version 1.0-17 diff --git a/R/geom-transformers.R b/R/geom-transformers.R index af13a68e3..f29e67696 100644 --- a/R/geom-transformers.R +++ b/R/geom-transformers.R @@ -450,6 +450,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 tesselation. \code{st_voronoi} requires GEOS version 3.5 or above #' @examples #' set.seed(1) @@ -470,24 +471,39 @@ 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 @@ -495,8 +511,8 @@ st_voronoi.sfc = function(x, envelope = st_polygon(), dTolerance = 0.0, bOnlyEdg } #' @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 diff --git a/man/geos_unary.Rd b/man/geos_unary.Rd index 0041b1510..e9c4ca8c1 100644 --- a/man/geos_unary.Rd +++ b/man/geos_unary.Rd @@ -49,7 +49,7 @@ st_inscribed_circle(x, dTolerance, ...) st_minimum_rotated_rectangle(x, ...) -st_voronoi(x, envelope, dTolerance = 0, bOnlyEdges = FALSE) +st_voronoi(x, envelope, dTolerance = 0, bOnlyEdges = FALSE, point_order = FALSE) st_polygonize(x) @@ -120,6 +120,8 @@ meters.} \item{of_largest_polygon}{logical; for \code{st_centroid}: if \code{TRUE}, return centroid of the largest (sub)polygon of a \code{MULTIPOLYGON} rather than of the whole \code{MULTIPOLYGON}} \item{dfMaxLength}{maximum length of a line segment. If \code{x} has geographical coordinates (long/lat), \code{dfMaxLength} is either a numeric expressed in meter, or an object of class \code{units} with length units \code{rad} or \code{degree}; segmentation in the long/lat case takes place along the great circle, using \link[lwgeom:geod]{st_geod_segmentize}.} + +\item{point_order}{logical; preserve point order if TRUE and GEOS version >= 3.12; overrides bOnlyEdges} } \value{ an object of the same class of \code{x}, with manipulated geometry. @@ -278,7 +280,15 @@ if (compareVersion(sf_extSoftVersion()[["GEOS"]], "3.5.0") > -1) { # 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) diff --git a/src/geos.cpp b/src/geos.cpp index 703cf1d6a..5846d6952 100644 --- a/src/geos.cpp +++ b/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 diff --git a/tests/geos.R b/tests/geos.R index ccce312be..376f14abb 100644 --- a/tests/geos.R +++ b/tests/geos.R @@ -73,7 +73,8 @@ 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) @@ -81,6 +82,17 @@ if (isTRUE(try(compareVersion(sf_extSoftVersion()["GEOS"], "3.5.0") > -1, silent 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)) diff --git a/tests/geos.Rout.save b/tests/geos.Rout.save index 5d0092256..303390bf9 100644 --- a/tests/geos.Rout.save +++ b/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,7 +162,8 @@ 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) @@ -170,6 +171,17 @@ LINESTRING (0 0, 1 1, 2 0) + 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