diff --git a/DESCRIPTION b/DESCRIPTION
index 0972d614..a5d2fec7 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -35,7 +35,7 @@ LazyData: true
 Depends:
     R (>= 3.3.0),
     abind,
-    sf (>= 1.0-15)
+    sf (>= 1.0-19)
 Imports:
     methods,
     parallel,
diff --git a/R/extract.R b/R/extract.R
index 476268d0..28eb1e3f 100644
--- a/R/extract.R
+++ b/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
diff --git a/man/st_extract.Rd b/man/st_extract.Rd
index 686f6c04..ef71edb9 100644
--- a/man/st_extract.Rd
+++ b/man/st_extract.Rd
@@ -14,7 +14,8 @@ st_extract(x, ...)
   bilinear = FALSE,
   time_column = attr(at, "time_column") \%||\% attr(at, "time_col"),
   interpolate_time = bilinear,
-  FUN = mean
+  FUN = mean,
+  resampling = c("nearest", "bilinear", "cubic", "cubicspline")
 )
 }
 \arguments{
@@ -31,6 +32,9 @@ st_extract(x, ...)
 \item{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.}
 
 \item{FUN}{function used to aggregate pixel values when geometries of \code{at} intersect with more than one pixel}
+
+\item{resampling}{character; resampling method; for method cubic or cubicspline,
+`stars_proxy` objects should be used and GDAL should have version >= 3.10.0}
 }
 \value{
 if \code{at} is of class \code{matrix}, a matrix with extracted values is returned;