use InterpolateAtPoint if GDAL >= 3.10.0

DESCRIPTION

@@ -35,7 +35,7 @@ LazyData: true
 Depends:
     R (>= 3.3.0),
     abind,
-    sf (>= 1.0-15)
+    sf (>= 1.0-19)
 Imports:
     methods,
     parallel,

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 with time or date values that will be matched to values of the first temporal dimension (matching classes \code{POSIXct}, \code{POSIXt}, \code{Date}, or \code{PCICt}), 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 FALSE, time instances are matched using the coinciding or the last preceding time in the data cube.
 #' @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"))
 		pts = at
 	else {
@@ -54,7 +62,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
@@ -63,7 +71,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