diff --git a/DESCRIPTION b/DESCRIPTION index a761727f1..3d59d06a7 100644 --- a/DESCRIPTION +++ b/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' diff --git a/NAMESPACE b/NAMESPACE index be8f53fce..b1227fbc2 100644 --- a/NAMESPACE +++ b/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) diff --git a/NEWS.md b/NEWS.md index e4ddcdf7a..185be8c90 100644 --- a/NEWS.md +++ b/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 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/R/sf.R b/R/sf.R index a387ebe53..c4b00ac5a 100644 --- a/R/sf.R +++ b/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 diff --git a/R/sf_aware.R b/R/sf_aware.R new file mode 100644 index 000000000..db1bba126 --- /dev/null +++ b/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 +} diff --git a/man/geos_unary.Rd b/man/geos_unary.Rd index 0041b1510..81f8cfc09 100644 --- a/man/geos_unary.Rd +++ b/man/geos_unary.Rd @@ -49,7 +49,13 @@ 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) @@ -115,6 +121,8 @@ meters.} \item{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} +\item{point_order}{logical; preserve point order if TRUE and GEOS version >= 3.12; overrides bOnlyEdges} + \item{directed}{logical; if \code{TRUE}, lines with opposite directions will not be merged} \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}} @@ -278,7 +286,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/man/make_sf_aware.Rd b/man/make_sf_aware.Rd new file mode 100644 index 000000000..2b10f5a28 --- /dev/null +++ b/man/make_sf_aware.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/sf_aware.R +\name{make_sf_aware} +\alias{make_sf_aware} +\title{for packages: use the sticky geometry [ behaviour of sf in package code} +\usage{ +make_sf_aware(env = parent.frame(), value = TRUE) +} +\arguments{ +\item{env}{environment} + +\item{value}{logical; default \code{TRUE}} +} +\description{ +for packages: use the sticky geometry [ behaviour of sf in package code +} 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))