Merge branch 'r-spatial:main' into main

.github/workflows/R-CMD-check.yaml

@@ -21,7 +21,7 @@ jobs:
       fail-fast: false
       matrix:
         config:
-          - {os: macos-latest,   r: 'release'}
+          # - {os: macos-latest,   r: 'release'}
           - {os: windows-latest, r: 'release'}
           - {os: ubuntu-latest,   r: 'devel', http-user-agent: 'release'}
           - {os: ubuntu-latest,   r: 'release'}

DESCRIPTION

@@ -35,7 +35,8 @@ LazyData: true
 Depends:
     R (>= 3.3.0),
     abind,
-    sf (>= 1.0-15)
+    sf (>= 1.0-19)
+Remotes: r-spatial/sf
 Imports:
     CFtime (>= 1.3.0),
     methods,

NEWS.md

@@ -1,5 +1,9 @@
 # version 0.6-7
 
+* `st_extract()` if used with GDAL 3.10.0 uses InterpolateAtPoints, allowing for cubic and cubicspline interpolators (requiring sf >= 1.0-19).
+
+* `Ops.stars()` (math ops) now also recycle arrays in the first argument; #718
+
 * `c.stars()` verifies semantic equivalence of objects' CRS; #703
 
 * initial support for `read_mdim()` to work with `proxy = TRUE`; #659 

R/extract.R

@@ -12,6 +12,8 @@ st_extract = function(x, ...) UseMethod("st_extract")
 #' @param time_column character or integer; name or index of a column in \code{at} (which must be of class \code{sf} or \code{sfc}) with time or date values that will be matched to values of the first temporal dimension (matching classes \code{POSIXct}, \code{Date}, or \code{CFtime}) in \code{x}, after which this dimension is reduced. This is useful to extract data cube values along a trajectory; see https://github.com/r-spatial/stars/issues/352 .
 #' @param interpolate_time logical; should time be interpolated? If \code{FALSE}, time instances are matched using the coinciding or the last preceding time in the data cube, unless the dimension uses \code{CFtime} with bounds set, in which case the corresponding time interval in which the time or date value falls is returned.
 #' @param FUN function used to aggregate pixel values when geometries of \code{at} intersect with more than one pixel
+#' @param resampling character; resampling method; for method cubic or cubicspline,
+#' `stars_proxy` objects should be used and GDAL should have version >= 3.10.0
 #' @param ... passed on to \link{aggregate.stars} when geometries are not exclusively POINT geometries
 #' @returns if \code{at} is of class \code{matrix}, a matrix with extracted values is returned; 
 #' otherwise: if \code{x} has more dimensions than only x and y (raster), an 
@@ -31,9 +33,15 @@ st_extract = function(x, ...) UseMethod("st_extract")
 #' st_extract(r, st_coordinates(pnt)) # "at" is a matrix: return a matrix
 st_extract.stars = function(x, at, ..., bilinear = FALSE, time_column = 
 		attr(at, "time_column") %||% attr(at, "time_col"),
-		interpolate_time = bilinear, FUN = mean) {
+		interpolate_time = bilinear, FUN = mean,
+		resampling = c("nearest", "bilinear", "cubic", "cubicspline")) {
 
 	stopifnot(inherits(at, c("sf", "sfc", "matrix")))
+	resampling = match.arg(resampling)
+	if (bilinear) {
+		stopifnot(resampling %in% c("nearest", "bilinear"))
+		resampling = "bilinear"
+	}
 	if (inherits(at, "matrix"))
 		# FIXME: matrix form of at supports only x,y rasters, no z,t or higher dimensions
 		pts = at
@@ -63,7 +71,7 @@ st_extract.stars = function(x, at, ..., bilinear = FALSE, time_column =
 		pts = st_coordinates(at)
 	}
 
-	if (bilinear && !inherits(x, "stars_proxy"))
+	if (resampling != "nearest" && !inherits(x, "stars_proxy"))
 		x = st_as_stars_proxy(x)
 
 	min_dist = NULL
@@ -72,7 +80,7 @@ st_extract.stars = function(x, at, ..., bilinear = FALSE, time_column =
 	m = if (inherits(x, "stars_proxy")) {
 			try_result = try(x0 <- st_as_stars(x, downsample = dim(x)/2), silent = TRUE)
 			lapply(x, function(y) do.call(abind, lapply(get_names(y), 
-				gdal_extract, pts = pts, bilinear = bilinear)))
+				gdal_extract, pts, resampling)))
 		} else {
 			x = st_normalize(st_upfront(x))
 			if (is_curvilinear(x)) { # https://github.com/r-spatial/stars/issues/632

R/ops.R

@@ -34,11 +34,22 @@ first_dimensions_match = function(e1, e2) {
 #' to permutate dimensions first. 
 Ops.stars <- function(e1, e2) {
 	if (!missing(e2)) { 
-		if (inherits(e1, "stars") && inherits(e2, "stars") && !first_dimensions_match(e1, e2))
-			stop("(first) dimensions of e1 and e2 do not match")
-		if (!inherits(e2, "stars"))
+		if (inherits(e1, "stars") && inherits(e2, "stars")) {
+			if (!first_dimensions_match(e1, e2))
+				stop("(first) dimensions of e1 and e2 do not match")
+			dim_final = if (prod(dim(e1)) < prod(dim(e2)))
+							st_dimensions(e2)
+						else
+							st_dimensions(e1)
+		} else if (inherits(e1, "stars")) {
+			dim_final = st_dimensions(e1)
+		} else if (inherits(e2, "stars")) {
+			dim_final = st_dimensions(e2)
+		}
+		if (!inherits(e2, c("stars", "units")))
 			e1 = drop_units(e1)
-	}
+	} else
+		dim_final = st_dimensions(e1)
 	ret = if (missing(e2))
 			lapply(e1, .Generic)
 		else if (!inherits(e2, "stars"))
@@ -48,16 +59,16 @@ Ops.stars <- function(e1, e2) {
 			if (!is.null(dim(e1)) &&
 					!isTRUE(all.equal(dim(e1), dim(e2), check.attributes = FALSE))) {
 				stopifnot(length(e2) == 1)
-				lapply(e1, .Generic, e2 = structure(e2[[1]], dim = NULL))
+				lapply(lapply(e1, structure, dim=NULL), .Generic, e2 = structure(e2[[1]], dim = NULL))
 			} else
 				mapply(.Generic, e1, e2, SIMPLIFY = FALSE)
 		}
 	if (any(sapply(ret, function(x) is.null(dim(x))))) # happens if e1[[1]] is a factor; #304
-		ret = lapply(ret, function(x) { dim(x) = dim(e1); x })
+		ret = lapply(ret, function(x) { dim(x) = dim(dim_final); x })
 	if (! inherits(e1, "stars"))
-		st_as_stars(setNames(ret, names(e2)), dimensions = st_dimensions(e2))
+		st_as_stars(setNames(ret, names(e2)), dimensions = dim_final)
 	else
-		st_as_stars(ret, dimensions = st_dimensions(e1))
+		st_as_stars(ret, dimensions = dim_final)
 }
 
 #' Mathematical operations for stars objects

R/raster.R

@@ -358,7 +358,7 @@ st_as_raster = function(x, class, ...) {
 				if (is.character(z)) {
 					names(b) = z
 				} else if (inherits(z, c("POSIXt", "Date"))) {
-					time(b) = z
+					terra::time(b) = z
 				} else {
 					names(b) = paste0(third, z)
 				}

R/read.R

@@ -83,7 +83,7 @@ geoloc_is_2D = function(geolocation, driver) { # the thing has 2-D x and y array
 #' \code{RasterIO} is a list with zero or more of the following named arguments:
 #' \code{nXOff}, \code{nYOff} (both 1-based: the first row/col has offset value 1),
 #' \code{nXSize}, \code{nYSize}, \code{nBufXSize}, \code{nBufYSize}, \code{bands}, \code{resample}.
-#' See \url{https://gdal.org/doxygen/classGDALDataset.html} for their meaning;
+#' See \url{https://gdal.org/en/latest/doxygen/classGDALDataset.html} for their meaning;
 #' \code{bands} is an integer vector containing the band numbers to be read (1-based: first band is 1).
 #' Note that if \code{nBufXSize} or \code{nBufYSize} are specified for downsampling an image,
 #' resulting in an adjusted geotransform. \code{resample} reflects the resampling method and