Merge branch 'main' into voronoi_point_order_new

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: sf
-Version: 1.0-18
+Version: 1.0-19
 Title: Simple Features for R
 Authors@R:
     c(person(given = "Edzer",

NAMESPACE

@@ -237,6 +237,10 @@ S3method(st_intersects,sfg)
 S3method(st_is,sf)
 S3method(st_is,sfc)
 S3method(st_is,sfg)
+S3method(st_is_full,bbox)
+S3method(st_is_full,sf)
+S3method(st_is_full,sfc)
+S3method(st_is_full,sfg)
 S3method(st_is_valid,sf)
 S3method(st_is_valid,sfc)
 S3method(st_is_valid,sfg)
@@ -380,6 +384,7 @@ export(.get_layout)
 export(.image_scale)
 export(.image_scale_factor)
 export(.stop_geos)
+export(FULL_bbox_)
 export(NA_agr_)
 export(NA_bbox_)
 export(NA_crs_)
@@ -463,6 +468,7 @@ export(st_intersection)
 export(st_intersects)
 export(st_is)
 export(st_is_empty)
+export(st_is_full)
 export(st_is_longlat)
 export(st_is_simple)
 export(st_is_valid)

NEWS.md

@@ -1,9 +1,25 @@
-# version 1.0-18
+# 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 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
+
 * fix build failure with GDAL < 3.4.0 #2436
 
+* `st_simplify()` now accepts feature-wise tolerance values when `s2` is switched on #2442
+
 # version 1.0-17
 
 * add `st_transform()` method for `bbox` objects; this uses OGRCoordinateTransformation::TransformBounds(), densifying first and antemeridian proof; #2415

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) {
@@ -361,6 +361,10 @@ sfc_is_empty <- function(sfc) {
     .Call(`_sf_sfc_is_empty`, sfc)
 }
 
+sfc_is_full <- function(sfc) {
+    .Call(`_sf_sfc_is_full`, sfc)
+}
+
 points_cpp <- function(pts, gdim = "XY") {
     .Call(`_sf_points_cpp`, pts, gdim)
 }
@@ -389,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 = bbox.MtrxSet
+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
@@ -138,7 +138,7 @@ compute_bbox = function(obj) {
 		sfc_POINT = bb_wrap(bbox.Set(obj)),
 		sfc_MULTIPOINT = bb_wrap(bbox.MtrxSet(obj)),
 		sfc_LINESTRING = bb_wrap(bbox.MtrxSet(obj)),
-		sfc_POLYGON = bb_wrap(bbox.MtrxSetSet(obj)),
+		sfc_POLYGON = if (any(st_is_full(obj))) FULL_bbox_ else bb_wrap(bbox.MtrxSetSet(obj)),
 		sfc_MULTILINESTRING = bb_wrap(bbox.MtrxSetSet(obj)),
 		sfc_MULTIPOLYGON = bb_wrap(bbox.MtrxSetSetSet(obj)),
 		bbox_list(obj)
@@ -219,6 +219,14 @@ NA_bbox_ = structure(rep(NA_real_, 4),
 	crs = NA_crs_,
 	class = "bbox")
 
+#' @name st_bbox
+#' @details \code{NA_bbox_} represents the missing value for a \code{bbox} object
+#' @export
+FULL_bbox_ = structure(c(-180.,-90.,180.,90.),
+	names = c("xmin", "ymin", "xmax", "ymax"),
+	crs = NA_crs_,
+	class = "bbox")
+
 #' @name st_bbox
 #' @export
 format.bbox = function(x, ...) {

R/crs.R

@@ -249,7 +249,7 @@ st_is_longlat = function(x) {
 		NA
 	else {
 		ret = crs_parameters(crs)$IsGeographic
-		if (ret && inherits(x, c("sf", "sfc", "stars"))) {
+		if (ret && inherits(x, c("sf", "sfc", "stars")) && !is.null(attr(x, "bbox"))) {
 			bb = st_bbox(x)
 			# check for potentially meaningless value range:
 			eps = sqrt(.Machine$double.eps)

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-predicates.R

@@ -4,15 +4,24 @@
 #' @examples
 #' ls = st_linestring(rbind(c(0,0), c(1,1), c(1,0), c(0,1)))
 #' st_is_simple(st_sfc(ls, st_point(c(0,0))))
-st_is_simple = function(x) CPL_geos_is_simple(st_geometry(x))
+st_is_simple = function(x) {
+	x = st_geometry(x)
+	not_full = !sfc_is_full(x)
+	ret = rep(TRUE, length(x))
+	ret[not_full] = CPL_geos_is_simple(x[not_full])
+	ret
+}
 
 #' @name geos_query
 #' @export
 #' @return st_is_empty returns for each geometry whether it is empty
 #' @examples
 #' ls = st_linestring(rbind(c(0,0), c(1,1), c(1,0), c(0,1)))
 #' st_is_empty(st_sfc(ls, st_point(), st_linestring()))
-st_is_empty = function(x) CPL_geos_is_empty(st_geometry(x))
+st_is_empty = function(x) sfc_is_empty(st_geometry(x)) 
+  # used to call
+  # CPL_geos_is_empty(st_geometry(x))
+  # but this avoids a R -> WKB -> GEOS conversion
 
 is_symmetric = function(operation, pattern) {
 	if (!is.na(pattern)) {
@@ -165,9 +174,9 @@ st_intersects.sfg = function(x, y, sparse = TRUE, prepared = TRUE, ...)
 
 #' @name geos_binary_pred
 #' @export
-st_disjoint		= function(x, y = x, sparse = TRUE, prepared = TRUE) {
+st_disjoint		= function(x, y = x, sparse = TRUE, prepared = TRUE, ...) {
 	# st_geos_binop("disjoint", x, y, sparse = sparse, prepared = prepared) -> didn't use STRtree
-	int = st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared)
+	int = st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared, ...)
 	# disjoint = !intersects :
 	if (sparse)
 		sgbp(lapply(int, function(g) setdiff(seq_along(st_geometry(y)), g)),
@@ -254,7 +263,7 @@ st_equals_exact = function(x, y, par, sparse = TRUE, prepared = FALSE, ...) {
 #' @name geos_binary_pred
 #' @export
 #' @param dist distance threshold; geometry indexes with distances smaller or equal to this value are returned; numeric value or units value having distance units.
-st_is_within_distance = function(x, y = x, dist, sparse = TRUE, ...) {
+st_is_within_distance = function(x, y = x, dist, sparse = TRUE, ..., remove_self = FALSE) {
 
 	ret = if (isTRUE(st_is_longlat(x))) {
 			units(dist) = as_units("m") # might convert
@@ -269,11 +278,11 @@ st_is_within_distance = function(x, y = x, dist, sparse = TRUE, ...) {
 					lwgeom::st_geod_distance(x, y, tolerance = dist, sparse = TRUE)
 				}
 			sgbp(r, predicate = "is_within_distance", region.id = seq_along(x), 
-				ncol = length(st_geometry(y)))
+				remove_self = remove_self, ncol = length(st_geometry(y)))
 		} else {
 			if (!is.null(st_crs(x)$ud_unit))
 				units(dist) = st_crs(x)$ud_unit # might convert
-			st_geos_binop("is_within_distance", x, y, par = dist, sparse = sparse, ...)
+			st_geos_binop("is_within_distance", x, y, par = dist, sparse = sparse, remove_self = remove_self, ...)
 		}
 	if (!sparse)
 		as.matrix(ret)

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)
@@ -286,7 +286,12 @@ st_simplify.sfc = function(x, preserveTopology, dTolerance = 0.0) {
 	if (ll && sf_use_s2()) {
 		if (!missing(preserveTopology) && isFALSE(preserveTopology))
 			warning("argument preserveTopology cannot be set to FALSE when working with ellipsoidal coordinates since the algorithm behind st_simplify always preserves topological relationships")
-		st_as_sfc(s2::s2_simplify(x, dTolerance), crs = st_crs(x))
+		if (length(dTolerance) == 1) {
+			st_as_sfc(s2::s2_simplify(x, dTolerance), crs = st_crs(x))
+		} else {
+			simplify <- function(x, dTolerance) st_as_sfc(s2::s2_simplify(x, dTolerance))
+			st_as_sfc(mapply(simplify, x, dTolerance), crs = st_crs(x))
+		}
 	} else {
 		if (missing(preserveTopology)) {
 			preserveTopology = FALSE
@@ -308,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,7 +456,7 @@ st_minimum_rotated_rectangle.sf = function(x, dTolerance, ...) {
 #' @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
+#' @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))
@@ -516,7 +521,7 @@ st_voronoi.sf = function(x, envelope = st_polygon(), dTolerance = 0.0, bOnlyEdge
 }
 
 #' @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)))
@@ -709,7 +714,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)) {
@@ -1189,5 +1193,5 @@ st_exterior_ring.sfc = function(x, ...) {
 		else
 			stop(paste("no exterior_ring method for objects of class", class(x)[1]))
 	}
-    st_as_sfc(lapply(x, exterior_sfg), crs = st_crs(x))
+	st_as_sfc(lapply(x, exterior_sfg), crs = st_crs(x))
 }

R/init.R

@@ -26,6 +26,7 @@ setOldClass(c("sfc_GEOMETRY", "sfc"))
 setOldClass(c("sfc_GEOMETRYCOLLECTION", "sfc"))
 setOldClass("sfg")
 setOldClass("crs")
+setOldClass("bbox")
 
 .sf_cache <- new.env(FALSE, parent=globalenv())
 
@@ -39,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
@@ -295,8 +295,6 @@ st_read.character = function(dsn, layer, ..., query = NA, options = NULL, quiet
 	if (length(promote_to_multi) > 1)
 		stop("`promote_to_multi' should have length one, and applies to all geometry columns")
 
-
-
 	if (use_stream) {
 		stream = nanoarrow::nanoarrow_allocate_array_stream()
 		info = CPL_read_gdal_stream(stream, dsn, layer, query, as.character(options), quiet,

R/s2.R

@@ -44,13 +44,13 @@ as_s2_geography.sfg <- function(x, ..., oriented = getOption("s2_oriented", FALS
 }
 
 # dynamically exported in tidyverse.R
-as_s2_geography.sfc <- function(x, ..., oriented = getOption("s2_oriented", FALSE)) {
+as_s2_geography.sfc <- function(x, ..., oriented = getOption("s2_oriented", FALSE) || isTRUE(attr(x, "oriented"))) {
 	st_as_s2.sfc(x, ..., oriented = oriented)
 }
 
 # dynamically exported in tidyverse.R
-as_s2_geography.sf <- function(x, ..., oriented = getOption("s2_oriented", FALSE)) {
-	st_as_s2.sf(x, ..., oriented = oriented)
+as_s2_geography.sf <- function(x, ...) {
+	st_as_s2.sf(x, ...)
 }
 
 #' functions for spherical geometry, using s2 package
@@ -86,7 +86,8 @@ st_as_s2.sf = function(x, ...) st_as_s2(st_geometry(x), ...)
 #' left of the polygon's path.
 #' @param rebuild logical; call \link[s2]{s2_rebuild} on the geometry (think of this as a \code{st_make_valid} on the sphere)
 #' @export
-st_as_s2.sfc = function(x, ..., oriented = getOption("s2_oriented", FALSE), rebuild = FALSE) {
+st_as_s2.sfc = function(x, ..., oriented = getOption("s2_oriented", FALSE) || isTRUE(attr(x, "oriented")), 
+						rebuild = FALSE) {
 	if (!is.na(st_crs(x)) && !st_is_longlat(x))
 		x = st_transform(x, ifelse(st_axis_order(), "OGC:CRS84", "EPSG:4326"))
 	if (length(x) && nchar(class(x[[1]])[1]) > 2) { # Z, M, ZM: