Merge branch 'main' into docs-geo

.Rbuildignore

@@ -28,3 +28,4 @@ changes1.txt
 changes2.txt
 changes3.txt
 changes.txt
+aware.patch

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: sf
-Version: 1.0-18
+Version: 1.0-20
 Title: Simple Features for R
 Authors@R:
     c(person(given = "Edzer",
@@ -47,11 +47,14 @@ Authors@R:
              role = "ctb",
              comment = c(ORCID = "0000-0002-9415-4582"))
 			 )
-Description: Support for simple features, a standardized way to
-    encode spatial vector data. Binds to 'GDAL' for reading and writing
-    data, to 'GEOS' for geometrical operations, and to 'PROJ' for
-    projection conversions and datum transformations. Uses by default the 's2'
-    package for spherical geometry operations on ellipsoidal (long/lat) coordinates.
+Description: Support for simple feature access, a standardized way to
+    encode and analyze spatial vector data. Binds to 'GDAL' <doi:
+    10.5281/zenodo.5884351> for reading and writing data, to 'GEOS'
+    <doi: 10.5281/zenodo.11396894> for geometrical operations,
+    and to 'PROJ' <doi: 10.5281/zenodo.5884394> for projection
+    conversions and datum transformations. Uses by default the 's2'
+    package for geometry operations on geodetic (long/lat degree)
+	coordinates.
 License: GPL-2 | MIT + file LICENSE
 URL: https://r-spatial.github.io/sf/, https://github.com/r-spatial/sf
 BugReports: https://github.com/r-spatial/sf/issues

NAMESPACE

@@ -270,6 +270,9 @@ S3method(st_m_range,sfc)
 S3method(st_make_valid,sf)
 S3method(st_make_valid,sfc)
 S3method(st_make_valid,sfg)
+S3method(st_minimum_bounding_circle,sf)
+S3method(st_minimum_bounding_circle,sfc)
+S3method(st_minimum_bounding_circle,sfg)
 S3method(st_minimum_rotated_rectangle,sf)
 S3method(st_minimum_rotated_rectangle,sfc)
 S3method(st_minimum_rotated_rectangle,sfg)
@@ -485,6 +488,7 @@ export(st_linestring)
 export(st_m_range)
 export(st_make_grid)
 export(st_make_valid)
+export(st_minimum_bounding_circle)
 export(st_minimum_rotated_rectangle)
 export(st_multilinestring)
 export(st_multipoint)

NEWS.md

@@ -1,6 +1,27 @@
+# version 1.0-20
+
+* `st_buffer()` for geodetic coordinates allows `max_dist` and `min_level` 
+   to be specified by feature; #2488 and https://github.com/r-spatial/s2/pull/264
+
+* `distinct.sf()` allows for comparing exact equality of geometries when `exact = TRUE`; #2484
+
+* `st_minimum_bounding_circle()` returns geometries representing the smallest circle that contains the input; #2473
+
+# version 1.0-19
+
+* fix type checks in C++ GDAL area and length computation functions, anticipating GDAL 3.10.0; #2466, #2468, #2469 by @rsbivand and @rouault
+
+* improve test on empty geometries, which changed in 1.0-18; #2463
+
+* `gdal_utils()` `ogrinfo` has an argument `read_only` which, when `TRUE` (or `options` includes `"-ro"`), opens a datasource in read-only mode (#2460; `sf` did this before 1.0-17); by default a datasource is opened in update (read-write) mode (since sf 1.0-17; #2420)
+
+* 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 the entire Earth surface, as used by `s2geometry`; #2441
+* 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
 
 * `st_sfc()` has an argument `oriented` which, when set to `TRUE`, adds an attribute `oriented=TRUE` to an `sfc` object, indicating that this object should not be reoriented in conversion to `s2_geography` (avoiding using the global option `s2_oriented`); `st_as_sfc.bbox()` sets this to `TRUE`; #2441
 

R/RcppExports.R

@@ -117,8 +117,8 @@ CPL_gdalinfo <- function(obj, options, oo, co) {
     .Call(`_sf_CPL_gdalinfo`, obj, options, oo, co)
 }
 
-CPL_ogrinfo <- function(obj, options, oo, co) {
-    .Call(`_sf_CPL_ogrinfo`, obj, options, oo, co)
+CPL_ogrinfo <- function(obj, options, oo, co, read_only = FALSE) {
+    .Call(`_sf_CPL_ogrinfo`, obj, options, oo, co, read_only)
 }
 
 CPL_gdaladdo <- function(obj, method, overviews, bands, oo, co, clean = FALSE, read_only = FALSE) {
@@ -393,7 +393,7 @@ CPL_write_gdal <- function(x, fname, driver, options, Type, dims, from, gt, p4s,
     invisible(.Call(`_sf_CPL_write_gdal`, x, fname, driver, options, Type, dims, from, gt, p4s, na_val, scale_offset, create, only_create))
 }
 
-CPL_extract <- function(input, xy, interpolate = FALSE) {
+CPL_extract <- function(input, xy, interpolate) {
     .Call(`_sf_CPL_extract`, input, xy, interpolate)
 }
 

R/bbox.R

@@ -71,7 +71,7 @@ st_bbox.MULTIPOINT = bbox.Mtrx
 st_bbox.LINESTRING = bbox.Mtrx
 #' @export
 #' @name st_bbox
-st_bbox.POLYGON = function(obj, ...) if (sf_use_s2() && length(obj) == 1 && nrow(obj[[1]]) == 2) FULL_bbox_ else bbox.MtrxSet(obj)
+st_bbox.POLYGON = function(obj, ...) if (st_is_full(obj)) FULL_bbox_ else bbox.MtrxSet(obj)
 #' @export
 #' @name st_bbox
 st_bbox.MULTILINESTRING = bbox.MtrxSet

R/crs.R

@@ -218,6 +218,9 @@ make_crs = function(x) {
 	if (!is.na(start_crs) && !is.na(end_crs) && start_crs != end_crs)
 		warning("st_crs<- : replacing crs does not reproject data; use st_transform for that", call. = FALSE)
 
+	if (is.na(end_crs) && !is.na(start_crs) && isTRUE(st_is_longlat(start_crs)) && any(st_is_full(x)))
+		stop("To set the crs to NA, first remove the full polygons; see: st_is_full()")
+
 	structure(x, crs = end_crs)
 }
 

R/gdal_utils.R

@@ -31,6 +31,7 @@ resampling_method = function(option = "near") {
 #' @param quiet logical; if \code{TRUE}, suppress printing the output for \code{info} and \code{mdiminfo}, and suppress printing progress
 #' @param processing character; processing options for \code{demprocessing}
 #' @param colorfilename character; name of color file for \code{demprocessing} (mandatory if \code{processing="color-relief"})
+#' @param read_only logical; only for `ogrinfo`: if `TRUE`, source is opened in read-only mode
 #' @return \code{info} returns a character vector with the raster metadata; all other utils return (invisibly) a logical indicating success (i.e., \code{TRUE}); in case of failure, an error is raised.
 #' @export
 #' @seealso \link{gdal_addo} for adding overlays to a raster file; \link{st_layers} to query geometry type(s) and crs from layers in a (vector) data source
@@ -76,7 +77,7 @@ gdal_utils = function(util = "info", source, destination, options = character(0)
 		quiet = !(util %in% c("info", "gdalinfo", "ogrinfo", "vectorinfo", 
 							  "mdiminfo")) || ("-multi" %in% options),
 		processing = character(0), colorfilename = character(0),
-		config_options = character(0)) {
+		config_options = character(0), read_only = FALSE) {
 
 	stopifnot(is.character(options), is.character(config_options))
 	if (!quiet && "-multi" %in% options)
@@ -103,7 +104,8 @@ gdal_utils = function(util = "info", source, destination, options = character(0)
 
 	ret = switch(util,
 			gdalinfo =, info = CPL_gdalinfo(if (missing(source)) character(0) else source, options, oo, config_options),
-			vectorinfo =, ogrinfo = CPL_ogrinfo(if (missing(source)) character(0) else source, options, oo, config_options),
+			vectorinfo =, ogrinfo = CPL_ogrinfo(if (missing(source)) character(0) else source, options, oo, config_options, 
+												isTRUE(read_only) || "-ro" %in% options),
 			warp = CPL_gdalwarp(source, destination, options, oo, doo, config_options, quiet, "-overwrite" %in% options),
 			warper = CPL_gdal_warper(source, destination, as.integer(resampling_method(options)),
 				oo, doo, config_options, quiet), # nocov

R/geom-measures.R

@@ -23,7 +23,7 @@ st_dimension = function(x, NA_if_empty = TRUE)
 #' @export
 #' @return If the coordinate reference system of \code{x} was set, these functions return values with unit of measurement; see \link[units]{set_units}.
 #'
-#' st_area returns the area of a geometry, in the coordinate reference system used; in case \code{x} is in degrees longitude/latitude, \link[lwgeom:geod]{st_geod_area} is used for area calculation.
+#' st_area returns the area of each feature geometry, computed in the coordinate reference system used. In case \code{x} has geodetic coordinates (unprojected), then if `sf_use_s2()` is `FALSE` \link[lwgeom:geod]{st_geod_area} is used for area calculation, if it is `TRUE` then \link[s2:s2_is_collection]{s2_area} is used: the former assumes an ellipsoidal shape, the latter a spherical shape of the Earth. In case of projected data, areas are computed in flat space. The argument `...` can be used to specify `radius` to \link[s2:s2_is_collection]{s2_area}, to modify the Earth radius.
 #' @examples
 #' b0 = st_polygon(list(rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,-1))))
 #' b1 = b0 + 2

R/geom-transformers.R

@@ -19,18 +19,18 @@
 #' @param mitreLimit numeric; limit of extension for a join if \code{joinStyle} 'MITRE' is used (default 1.0, minimum 0.0); see details
 #' @param singleSide logical; if \code{TRUE}, single-sided buffers are returned for linear geometries,
 #' in which case negative \code{dist} values give buffers on the right-hand side, positive on the left; see details
-#' @param ... passed on to  [s2::s2_buffer_cells()]
+#' @param ... in `st_buffer` passed on to  [s2::s2_buffer_cells()], otherwise ignored
 #' @return an object of the same class of \code{x}, with manipulated geometry.
 #' @export
-#' @details \code{st_buffer} computes a buffer around this geometry/each geometry. If any of \code{endCapStyle},
-#' \code{joinStyle}, or \code{mitreLimit} are set to non-default values ('ROUND', 'ROUND', 1.0 respectively) then
-#' the underlying 'buffer with style' GEOS function is used.
-#' If a negative buffer returns empty polygons instead of shrinking, set sf_use_s2() to FALSE
-#' See \href{https://postgis.net/docs/ST_Buffer.html}{postgis.net/docs/ST_Buffer.html} for details.
+#' @details \code{st_buffer} computes a buffer around this geometry/each geometry. Depending on the spatial
+#' coordinate system, a different engine (GEOS or S2) can be used, which have different function
+#' arguments. The \code{nQuadSegs}, \code{endCapsStyle}, \code{joinStyle}, \code{mitreLimit} and
+#' \code{singleSide} parameters only work if the GEOS engine is used (i.e. projected coordinates or
+#' when \code{sf_use_s2()} is set to \code{FALSE}). See \href{https://postgis.net/docs/ST_Buffer.html}{postgis.net/docs/ST_Buffer.html}
+#' for details. The \code{max_cells} and \code{min_level} parameters ([s2::s2_buffer_cells()]) work with the S2
+#' engine (i.e. geographic coordinates) and can be used to change the buffer shape (e.g. smoothing). If
+#' a negative buffer returns empty polygons instead of shrinking, set \code{sf_use_s2()} to \code{FALSE}.
 #' 
-#' \code{nQuadSegs}, \code{endCapsStyle}, \code{joinStyle}, \code{mitreLimit} and \code{singleSide} only
-#' work when the GEOS back-end is used: for projected coordinates or when \code{sf_use_s2()} is set
-#' to \code{FALSE}.
 #' @examples
 #'
 #' ## st_buffer, style options (taken from rgeos gBuffer)
@@ -313,7 +313,7 @@ st_simplify.sf = function(x, preserveTopology, dTolerance = 0.0) {
 
 #' @name geos_unary
 #' @export
-#' @param bOnlyEdges logical; if TRUE, return lines, else return polygons
+#' @param bOnlyEdges logical; if \code{TRUE}, return lines, else return polygons
 #' @details \code{st_triangulate} triangulates set of points (not constrained). \code{st_triangulate} requires GEOS version 3.4 or above
 st_triangulate = function(x, dTolerance = 0.0, bOnlyEdges = FALSE)
 	UseMethod("st_triangulate")
@@ -451,11 +451,40 @@ st_minimum_rotated_rectangle.sf = function(x, dTolerance, ...) {
 	st_set_geometry(x, st_minimum_rotated_rectangle(st_geometry(x)), ...)
 }
 
+#' @name geos_unary
+#' @details \code{st_minimum_bounding_circle} 
+#' returns a geometry which represents the "minimum bounding circle",
+#' the smallest circle that contains the input.
+#' @export
+st_minimum_bounding_circle = function(x, ...)
+	UseMethod("st_minimum_bounding_circle")
+
+#' @export
+st_minimum_bounding_circle.sfg = function(x, ...) {
+	get_first_sfg(st_minimum_bounding_circle(st_sfc(x), ...))
+}
+
+#' @export
+st_minimum_bounding_circle.sfc = function(x, ...) {
+	if (compareVersion(CPL_geos_version(), "3.8.0") > -1) { # >=
+		if (isTRUE(st_is_longlat(x)))
+			warning("st_minimum_rotated_rectangle does not work correctly for longitude/latitude data")
+		st_sfc(CPL_geos_op("bounding_circle", x, 0L, integer(0),
+			dTolerance = 0., logical(0), bOnlyEdges = as.integer(FALSE)))
+	} else
+		stop("for st_minimum_bounding_circle, GEOS version 3.8.0 or higher is required")
+}
+
+#' @export
+st_minimum_bounding_circle.sf = function(x, ...) {
+	st_set_geometry(x, st_minimum_bounding_circle(st_geometry(x)), ...)
+}
 
 #' @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
-#' @details \code{st_voronoi} creates voronoi tesselation. \code{st_voronoi} requires GEOS version 3.5 or above
+#' @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)
 #' x = st_multipoint(matrix(runif(10),,2))
@@ -475,37 +504,52 @@ 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
-#' @details \code{st_polygonize} creates polygon from lines that form a closed ring. In case of \code{st_polygonize}, \code{x} must be an object of class \code{LINESTRING} or \code{MULTILINESTRING}, or an \code{sfc} geometry list-column object containing these
+#' @details \code{st_polygonize} creates a polygon from lines that form a closed ring. In case of \code{st_polygonize}, \code{x} must be an object of class \code{LINESTRING} or \code{MULTILINESTRING}, or an \code{sfc} geometry list-column object containing these
 #' @export
 #' @examples
 #' mls = st_multilinestring(list(matrix(c(0,0,0,1,1,1,0,0),,2,byrow=TRUE)))
@@ -698,7 +742,6 @@ st_node.sf = function(x) {
 #' @details \code{st_segmentize} adds points to straight lines
 #' @export
 #' @param 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}.
-#' @param ... ignored
 #' @examples
 #' sf = st_sf(a=1, geom=st_sfc(st_linestring(rbind(c(0,0),c(1,1)))), crs = 4326)
 #' if (require(lwgeom, quietly = TRUE)) {

R/init.R

@@ -40,6 +40,7 @@ pathGrob <- NULL
 	load_gdal() 
 	if ((s2 <- Sys.getenv("_SF_USE_S2")) != "")
 		options(sf_use_s2 = s2 != "false")
+	FULL_bbox_ <<- st_set_crs(FULL_bbox_, "OGC:CRS84")
 }
 
 .onUnload = function(libname, pkgname) {

R/read.R

@@ -61,7 +61,7 @@ set_utf8 = function(x) {
 #'   and if necessary corrected (when seen from above: exterior ring counter
 #'   clockwise, holes clockwise)
 #' @details for \code{geometry_column}, see also
-#' \url{https://trac.osgeo.org/gdal/wiki/rfc41_multiple_geometry_fields}
+#' \url{https://gdal.org/en/latest/development/rfc/rfc41_multiple_geometry_fields.html}
 #'
 #' for values for \code{type} see
 #' \url{https://en.wikipedia.org/wiki/Well-known_text#Well-known_binary}, but

R/sfc.R

@@ -641,14 +641,9 @@ st_is_full.sfg = function(x, ..., is_longlat = NULL) {
 #' @export
 #' @name st_is_full
 st_is_full.sfc = function(x, ...) {
-	if (sf_use_s2() && inherits(x, c("sfc_POLYGON", "sfc_GEOMETRY"))) {
-		is_longlat = if (!is.null(attr(x, "crs")))
-				st_is_longlat(x)
-			else
-				NA
-		#sapply(x, st_is_full.sfg, ..., is_longlat = is_longlat)
+	if (sf_use_s2() && inherits(x, c("sfc_POLYGON", "sfc_GEOMETRY")))
 		sfc_is_full(x)
-	} else
+	else
 		rep_len(FALSE, length(x))
 }
 
@@ -661,5 +656,5 @@ st_is_full.sf = function(x, ...) {
 #' @export
 #' @name st_is_full
 st_is_full.bbox = function(x, ...) {
-	sf_use_s2() && st_is_longlat(x) && all(x == c(-180,-90,180,90))
+	isTRUE(sf_use_s2() && st_is_longlat(x) && all(x == c(-180,-90,180,90)))
 }

R/spatstat.R

@@ -144,13 +144,12 @@ as.ppp.sf = function(X, ...) {
 	if (st_dimension(X[1,]) == 2)
 		X = X[-1,]
 	st_geometry(X) = NULL # remove geometry column
-	if (ncol(X) > 1)
-		warning("only first attribute column is used for marks")
 
-	if (ncol(X) == 0)
+	if (ncol(X) == 0) {
 		pp
-	else 
-		spatstat.geom::setmarks(pp, X[1])
+	} else {
+		spatstat.geom::setmarks(pp, X)
+	}
 }
 
 as.owin.POLYGON = function(W, ..., fatal, check_polygons = TRUE) {

R/stars.R

@@ -340,9 +340,10 @@ gdal_rasterize = function(sf, x, gt, file, driver = "GTiff", options = character
 #' @rdname gdal
 #' @param f gdal raster data source filename
 #' @param pts points matrix
-#' @param bilinear logical; use bilinear interpolation, rather than nearest neighbor?
-gdal_extract = function(f, pts, bilinear = FALSE) {
-	CPL_extract(f, pts, as.logical(bilinear))
+#' @param resampling character; resampling method; for method cubic or cubicspline, 
+#' `stars_proxy` objects should be used and GDAL should have version >= 3.10.0
+gdal_extract = function(f, pts, resampling = c("nearest", "bilinear", "cubic", "cubicspline")) {
+	CPL_extract(f, pts, match.arg(resampling))
 }
 
 #' @rdname gdal