Voronoi point order

DESCRIPTION

@@ -112,7 +112,7 @@ LinkingTo:
 VignetteBuilder: 
     knitr
 Encoding: UTF-8
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 Roxygen: list(markdown = TRUE)
 Config/testthat/edition: 2
 Config/needs/coverage: XML
@@ -162,6 +162,7 @@ Collate:
     'gdal_utils.R'
     'nearest.R'
     'normalize.R'
+    'sf_aware.R'
     'sf-package.R'
     'defunct.R'
     'z_range.R'

NAMESPACE

@@ -401,6 +401,7 @@ export(gdal_utils)
 export(gdal_write)
 export(gdal_write_mdim)
 export(get_key_pos)
+export(make_sf_aware)
 export(plot_sf)
 export(rawToHex)
 export(read_sf)

NEWS.md

@@ -1,5 +1,7 @@
 # version 1.0-17
 
+* add flag for preservation of point order in `st_voronoi` #1371 for GEOS >= 3.12
+
 * add `st_transform()` method for `bbox` objects; this uses OGRCoordinateTransformation::TransformBounds(), densifying first and antemeridian proof; #2415
 
 * `st_filter.sf()` correctly scopes `x` and `y` arguments using !! operator; #2416

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,33 +471,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

R/sf.R

@@ -339,41 +339,45 @@ st_sf = function(..., agr = NA_agr_, row.names,
 #' h[g,]
 #' @export
 "[.sf" = function(x, i, j, ..., drop = FALSE, op = st_intersects) {
-	nargs = nargs()
-	agr = st_agr(x)
-	if (!missing(i) && (inherits(i, "sf") || inherits(i, "sfc") || inherits(i, "sfg")))
-		i = lengths(op(x, i, ...)) != 0
-	sf_column = attr(x, "sf_column")
-	geom = st_geometry(x)
-	if (!missing(i) && nargs > 2) { # e.g. a[3:4,] not a[3:4]
-		if (is.character(i))
-			i = match(i, row.names(x))
-		geom = geom[i]
-	}
 
-	# x = as.data.frame(x)
-	class(x) = setdiff(class(x), "sf") # one step down
-	x = if (missing(j)) {
-		if (nargs == 2) # `[`(x,i)
-			x[i] # do sth else for tbl?
-		else
-			x[i, , drop = drop]
-	} else
-		x[i, j, drop = drop]
+	if (is_sf_aware()) {
+		nargs = nargs()
+		agr = st_agr(x)
+		if (!missing(i) && (inherits(i, "sf") || inherits(i, "sfc") || inherits(i, "sfg")))
+			i = lengths(op(x, i, ...)) != 0
+		sf_column = attr(x, "sf_column")
+		geom = st_geometry(x)
+		if (!missing(i) && nargs > 2) { # e.g. a[3:4,] not a[3:4]
+			if (is.character(i))
+				i = match(i, row.names(x))
+			geom = geom[i]
+		}
 
-	if (!missing(j))
-		agr = agr[j]
-	else if (!missing(i) && nargs <= 2)
-		agr = agr[i] # e.g., obj["name"]
+		# x = as.data.frame(x)
+		class(x) = setdiff(class(x), "sf") # one step down
+		x = if (missing(j)) {
+			if (nargs == 2) # `[`(x,i)
+				x[i] # do sth else for tbl?
+			else
+				x[i, , drop = drop]
+		} else
+			x[i, j, drop = drop]
 
-	if (inherits(x, "sfc")) # drop was TRUE, and we selected geom column only
-		x
-	else if (! drop) {
-		x[[ sf_column ]] = geom
-		x = st_sf(x, sf_column_name = sf_column, sfc_last = FALSE)
-		st_set_agr(x, agr[match(setdiff(names(x), sf_column), names(agr))])
+		if (!missing(j))
+			agr = agr[j]
+		else if (!missing(i) && nargs <= 2)
+			agr = agr[i] # e.g., obj["name"]
+
+		if (inherits(x, "sfc")) # drop was TRUE, and we selected geom column only
+			x
+		else if (! drop) {
+			x[[ sf_column ]] = geom
+			x = st_sf(x, sf_column_name = sf_column, sfc_last = FALSE)
+			st_set_agr(x, agr[match(setdiff(names(x), sf_column), names(agr))])
+		} else
+			structure(x, class = setdiff(class(x), "sf"))
 	} else
-		structure(x, class = setdiff(class(x), "sf"))
+		NextMethod()
 }
 
 #' @export

R/sf_aware.R

@@ -0,0 +1,70 @@
+# see https://github.com/r-spatial/sf/issues/2131
+
+# Known namespaces that need sticky columns. This is set to the
+# packages that are already depending on sf-awareness at the time
+# the feature was introduced (Jul 23, 2024). 
+# New packages wanting to opt-in should use `make_sf_aware()` in their .onLoad().
+known_sf_aware <- c(".globalenv", "sf",
+					"alcyon", "animalEKF", "arcpullr", "basf", "bcmaps", "bispdep", "canadamaps", "CCAMLRGIS", "cffdrs", "chilemapas", "dggridR", "ecochange", "geotopbricks", "ggmapinset", "habCluster", "hgwrr", "LabourMarketAreas", "M3", "MazamaSpatialUtils", "micromap", "nngeo", "onmaRg", "PointedSDMs", "ppgm", "raptr", "RCzechia", "rgplates", "rKIN", "rLFT", "rts2", "segmetric", "sftime", "siland", "spatialreg", "spdep", "spMaps", "spsurvey", "stars", "starsExtra", "stcos", "surveyvoi", "tilegramsR", "transfR", "waver", "wdnr.gis", "wdpar", "windAC", "abmR", "abstr", "adw", "agricolaeplotr", "alarmdata", "amt", "anipaths", "aopdata", "appeears", "arcgisgeocode", "arcgislayers", "arcgisutils", "areal", "arealDB", "ARPALData", "ARUtools", "ascotraceR", "atakrig", "atpolR", "automap", "bangladesh", "basemaps", "bayesmove", "BayesX", "bcdata", "bcputility", "bdc", "bdl", "bdvis", "BeeBDC", "BFS", "bfsMaps", "BIEN", "bigDM", "BIOMASS", "bioregion", "blackmarbler", "blockCV", "bluebike", "BoundaryStats", "bRacatus", "btb", "camtrapR", "canadianmaps", "capm", "CARBayes", "CARBayesdata", "CARBayesST", "card", "cartograflow", "cartogram", "cartogramR", "cartographer", "cartographr", "cartography", "CAST", "CatastRo", "CDCPLACES", "cdrcR", "CDSE", "censable", "centr", "cft", "changeRangeR", "chessboard", "chirps", "cleangeo", "clhs", "cmsafvis", "cnmap", "CoastlineFD", "comorosmaps", "concaveman", "conleyreg", "constrainedKriging", "CopernicusDEM", "CopernicusMarine", "covid19br", "covid19sf", "covidcast", "crawl", "crimedata", "cropDemand", "CropScapeR", "cropZoning", "crsuggest", "CRTspat", "cshapes", "CSHShydRology", "csodata", "csquares", "ctmm", "cubble", "CvmortalityMult", "cyclestreets", "dafishr", "damAOI", "datazoom.amazonia", "daymetr", "densityarea", "DEPONS2R", "Directional", "disaggregation", "dispeRse", "distanceto", "divseg", "divvy", "dots", "downscale", "dsims", "dsmSearch", "dssd", "dwp", "dynamicSDM", "ebirdst", "echor", "edbuildmapr", "ediblecity", "EEAaq", "eiExpand", "eixport", "eks", "elevatr", "EmissV", "emstreeR", "enmSdmX", "envi", "epikit", "epiR", "epm", "eSDM", "evolMap", "exactextractr", "expowo", "extRatum", "FedData", "feltr", "fgdr", "FIESTA", "FIESTAutils", "fisheye", "fishRman", "fitbitViz", "flexpolyline", "flightplot", "FLightR", "fmesher", "forestecology", "ForestTools", "FORTLS", "fsr", "fude", "galah", "gbm.auto", "gdalUtilities", "geneHapR", "GeNetIt", "GeoAdjust", "geoAr", "geobr", "geocausal", "geocmeans", "geodimension", "geodiv", "geodrawr", "geofi", "geogenr", "geogrid", "geojsonio", "geomander", "GEOmap", "geomaroc", "geomerge", "geomultistar", "geonetwork", "geoperu", "geostan", "geouy", "gfcanalysis", "ggautomap", "ggfields", "ggOceanMaps", "GGoutlieR", "ggseg", "ggspatial", "GIFT", "giscoR", "GISINTEGRATION", "GISSB", "glottospace", "googletraffic", "gps.track", "GPSeqClus", "grainscape", "graph4lg", "GREENeR", "gridpattern", "gstat", "gtfs2emis", "gtfs2gps", "gtfstools", "gwavr", "GWmodel", "GWnnegPCA", "GWpcor", "gwpcormapper", "GWSDAT", "h3jsr", "happign", "HDSpatialScan", "helsinki", "hemispheR", "hereR", "hero", "himach", "hosm", "hwsdr", "hydroloom", "hyfo", "hypsoLoop", "IceSat2R", "icosa", "idbr", "imcRtools", "importinegi", "inlabru", "INLAspacetime", "inldata", "intamap", "intamapInteractive", "ipdw", "IRexamples", "ispdata", "ISRaD", "itsdm", "jmastats", "jpgrid", "jpmesh", "klexdatr", "KMLtoSHAPE", "kokudosuuchi", "LAGOSNE", "lakemorpho", "LandComp", "landsepi", "lazysf", "lconnect", "leafem", "leafgl", "leafpm", "leafpop", "leastcostpath", "letsR", "lgcp", "lidaRtRee", "lidR", "link2GI", "LMMsolver", "LMoFit", "lwgeom", "macleish", "macroBiome", "MainExistingDatasets", "malariaAtlas", "mapboxapi", "mapchina", "mapedit", "MapGAM", "mapi", "mapiso", "mapme.biodiversity", "mapmixture", "mapping", "mapsapi", "mapscanner", "mapsf", "mapSpain", "mapStats", "maptiles", "mapview", "MassWateR", "MazamaSpatialPlots", "medfateland", "meteo", "meteoland", "meteospain", "metR", "MetricGraph", "mgwrhw", "mgwrsar", "micromapST", "MigConnectivity", "misuvi", "mkde", "mlr3spatial", "modisfast", "MODISTools", "Momocs", "monographaR", "Morpho", "motif", "move2", "movecost", "movegroup", "MTA", "naijR", "naturaList", "ncdfgeom", "ndi", "neonPlantEcology", "netmap", "nhdplusTools", "nhdR", "NipponMap", "njgeo", "nlrx", "nominatimlite", "nswgeo", "oceanexplorer", "oceanic", "oceanis", "oceanmap", "odbr", "ofpetrial", "ohsome", "ohun", "openairmaps", "opendatatoronto", "openeo", "opentripplanner", "Orcs", "osmextract", "OSMscale", "osrm", "otpr", "overturemapsr", "ows4R", "ozmaps", "palaeoverse", "PAMscapes", "pargasite", "pastclim", "patternize", "pavo", "pct", "pgirmess", "PL94171", "PlanetNICFI", "planscorer", "plantTracker", "pliman", "plotdap", "plusCode2", "populR", "potential", "povmap", "pRecipe", "PReMiuM", "pressuRe", "prevR", "prioritizr", "prioritizrdata", "prisonbrief", "pseudohouseholds", "pspatreg", "ptools", "PublicWorksFinanceIT", "pycno", "qlcVisualize", "qualmap", "r5r", "raceland", "rangeBuilder", "rangeMapper", "rasterbc", "rasterDT", "rasterpic", "raybevel", "raytracing", "rcarbon", "RchivalTag", "rdhs", "rdwplus", "readwritesqlite", "redist", "redistmetrics", "redistverse", "redlistr", "ref.ICAR", "Relectoral", "remap", "rerddapXtracto", "rfishnet2", "rflexscan", "rgeoboundaries", "rgeoda", "rgeopat2", "rgugik", "ripc", "RiskMap", "riverdist", "rivnet", "rlandfire", "Rlibkdv", "rmapshaper", "rmapzen", "rnaturalearth", "rnrfa", "roads", "robis", "rolap", "rosmium", "roughsf", "rpostgis", "RPyGeo", "Rsagacmd", "rsat", "rsi", "rsocialwatcher", "rSPARCS", "rstac", "RStoolbox", "rTLsDeep", "rtop", "rWCVP", "RWmisc", "sabre", "sampbias", "satres", "SchoolDataIT", "scider", "SDLfilter", "SDPDmod", "secr", "secrdesign", "secrlinear", "seedreg", "sfarrow", "sfcentral", "sfdct", "sfdep", "sfhotspot", "sfislands", "sfnetworks", "sftrack", "sgapi", "sgsR", "shoredate", "simodels", "SimSurvey", "sits", "smartmap", "smile", "SMITIDstruct", "smoothr", "soilassessment", "sorvi", "spacejamr", "SpatFD", "SpatGC", "spatgeom", "SpatGRID", "spatialEco", "SpatialFeatureExperiment", "SpatialGraph", "SpatialKDE", "SpatialPosition", "SpatialRDD", "spatialrisk", "spatialsample", "SpaTopic", "spatsoc", "spatsurv", "spbal", "spectator", "spectralR", "sphet", "spmodel", "spmoran", "spnaf", "spNetwork", "sqlhelper", "SRTsim", "SSDM", "SSIMmap", "SSN2", "SSNbayes", "stampr", "starsTileServer", "stats19", "stelfi", "StormR", "stplanr", "stppSim", "streamDepletr", "streetscape", "sugarbag", "SUMMER", "SUNGEO", "suntools", "supercells", "SurfaceTortoise", "surveyPrev", "swfscDAS", "swfscMisc", "SWMPrExtension", "SWTools", "tanaka", "TDLM", "tectonicr", "tenm", "terrainr", "tidycensus", "tidygeoRSS", "tidyrgee", "tidysdm", "tidyterra", "tidytransit", "tidyUSDA", "tigris", "tilemaps", "tinytiger", "tmap", "tmaptools", "tongfen", "track2KBA", "trackdf", "trackeRapp", "transformr", "treePlotArea", "TreeRingShape", "treesliceR", "trigpoints", "TUFLOWR", "uavRmp", "uci", "ulex", "UnalR", "upstartr", "ursa", "usmapdata", "valhallr", "valuemap", "VancouvR", "vein", "velociraptr", "versioning", "VicmapR", "vietnameseConverter", "viewscape", "voluModel", "Voyager", "walkboutr", "wallace", "waywiser", "weed", "WEGE", "WeightedTreemaps", "WES", "wflo", "wildlifeDI", "wingen", "WorldMapR", "zippeR", "zonebuilder")
+
+# Register environment as sf-aware. We use a lexical flag instead of
+# inspecting the the name of the calling namespace and add it to
+# `the$sf_aware` because the flag is more flexible. It makes it
+# possible for arbitrary environments to opt into (or disable) the sf
+# API like this:
+#
+# ```
+# local({ make_sf_aware(); sf[...] })
+# ```
+#
+# Packages should call it like this in their `.onLoad()` hook (and
+# thus before their namespace is sealed):
+#
+# ```
+# make_sf_aware(topenv(environment()))
+# ```
+
+#' for packages: use the sticky geometry [ behaviour of sf in package code
+#'
+#' for packages: use the sticky geometry [ behaviour of sf in package code
+#' @param env environment
+#' @param value logical; default `TRUE`
+#' @export
+make_sf_aware <- function(env = parent.frame(), value = TRUE) {
+	env$.__sf_aware__. <- value
+}
+
+is_sf_aware <- function(env = parent.frame(2)) {
+	top <- topenv(env)
+
+	# Check for overrides
+	top_name <- env_name(top)
+	if (!is.null(top_name) && top_name %in% known_sf_aware)
+		TRUE
+	else {
+		if (!requireNamespace("rlang", quietly = TRUE))
+			stop("rlang is not installed: install first?")
+		# Now check for the lexical flag. This could be rewritten
+		# without rlang as a loop over parents that stops at `topenv()`.
+		flag <- rlang::env_get(
+			env,
+			".__sf_aware__.",
+			default = FALSE,
+			inherit = TRUE,
+			last = top
+		)
+		stopifnot(is.logical(flag) && length(flag) == 1 && !is.na(flag))
+		flag
+	}
+}
+
+env_name <- function(env) {
+	ns <- topenv(env)
+
+	if (isNamespace(ns))
+		getNamespaceName(ns)
+	else if (identical(ns, globalenv()))
+		".globalenv"
+	else
+		NULL
+}

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))