diff --git a/R/RcppExports.R b/R/RcppExports.R
index f4d9d07d3..57a366b3c 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -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)
 }
 
diff --git a/R/stars.R b/R/stars.R
index 4d171b164..7085c31e7 100644
--- a/R/stars.R
+++ b/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
diff --git a/man/gdal.Rd b/man/gdal.Rd
index 8d682dfa0..6d70c6861 100644
--- a/man/gdal.Rd
+++ b/man/gdal.Rd
@@ -68,7 +68,11 @@ gdal_polygonize(
 
 gdal_rasterize(sf, x, gt, file, driver = "GTiff", options = character())
 
-gdal_extract(f, pts, bilinear = FALSE)
+gdal_extract(
+  f,
+  pts,
+  resampling = c("nearest", "bilinear", "cubic", "cubicspline")
+)
 
 gdal_read_mdim(
   file,
@@ -147,7 +151,8 @@ gdal_create(f, nxy, values, crs, xlim, ylim)
 
 \item{pts}{points matrix}
 
-\item{bilinear}{logical; use bilinear interpolation, rather than nearest neighbor?}
+\item{resampling}{character; resampling method; for method cubic or cubicspline,
+\code{stars_proxy} objects should be used and GDAL should have version >= 3.10.0}
 
 \item{array_name}{array name}
 
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 56586fa7a..039d37e11 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -1339,14 +1339,14 @@ BEGIN_RCPP
 END_RCPP
 }
 // CPL_extract
-NumericMatrix CPL_extract(CharacterVector input, NumericMatrix xy, bool interpolate);
+NumericMatrix CPL_extract(CharacterVector input, NumericMatrix xy, CharacterVector interpolate);
 RcppExport SEXP _sf_CPL_extract(SEXP inputSEXP, SEXP xySEXP, SEXP interpolateSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< CharacterVector >::type input(inputSEXP);
     Rcpp::traits::input_parameter< NumericMatrix >::type xy(xySEXP);
-    Rcpp::traits::input_parameter< bool >::type interpolate(interpolateSEXP);
+    Rcpp::traits::input_parameter< CharacterVector >::type interpolate(interpolateSEXP);
     rcpp_result_gen = Rcpp::wrap(CPL_extract(input, xy, interpolate));
     return rcpp_result_gen;
 END_RCPP
diff --git a/src/stars.cpp b/src/stars.cpp
index 3671154a5..184baae29 100644
--- a/src/stars.cpp
+++ b/src/stars.cpp
@@ -721,8 +721,8 @@ double get_bilinear(GDALRasterBand *poBand, double Pixel, double Line,
 		dY -= 0.5;
 
 	// read:
-	if (GDALRasterIO(poBand, GF_Read, iPixel, iLine, 2, 2,
-			(void *) pixels, 2, 2, GDT_CFloat64, sizeof(double), 0) != CE_None)
+	if (poBand->RasterIO(GF_Read, iPixel, iLine, 2, 2,
+			(void *) pixels, 2, 2, GDT_CFloat64, sizeof(double), 0, NULL) != CE_None)
 		stop("Error reading!");
 	// f(0,0): pixels[0], f(1,0): pixels[1], f(0,1): pixels[2], f(1,1): pixels[3]
 	if (na_set && (pixels[0] == na_value || pixels[1] == na_value ||
@@ -736,7 +736,7 @@ double get_bilinear(GDALRasterBand *poBand, double Pixel, double Line,
 }
 
 // [[Rcpp::export]]
-NumericMatrix CPL_extract(CharacterVector input, NumericMatrix xy, bool interpolate = false) {
+NumericMatrix CPL_extract(CharacterVector input, NumericMatrix xy, CharacterVector interpolate) {
 	// mostly taken from gdal/apps/gdallocationinfo.cpp
 
 	GDALDataset *poDataset = (GDALDataset *) GDALOpenEx(input[0], GA_ReadOnly,
@@ -749,6 +749,17 @@ NumericMatrix CPL_extract(CharacterVector input, NumericMatrix xy, bool interpol
 	NumericMatrix ret(xy.nrow(), poDataset->GetRasterCount());
 	int xsize = poDataset->GetRasterXSize();
 	int ysize = poDataset->GetRasterYSize();
+	GDALRIOResampleAlg RA;
+	if (interpolate[0] == "nearest")
+		RA = GRIORA_NearestNeighbour;
+	else if (interpolate[0] == "bilinear")
+		RA = GRIORA_Bilinear;
+	else if (interpolate[0] == "cubic")
+		RA = GRIORA_Cubic;
+	else if (interpolate[0] == "cubicspline")
+		RA = GRIORA_CubicSpline;
+	else
+		stop("interpolation method not supported"); // #nocov
 
 	double gt[6];
 	poDataset->GetGeoTransform(gt);
@@ -779,16 +790,17 @@ NumericMatrix CPL_extract(CharacterVector input, NumericMatrix xy, bool interpol
 				pixel = NA_REAL;
 			else { // read pixel:
 #if GDAL_VERSION_NUM >= 3100000
-				if (poBand->InterpolateAtPoint(Pixel, Line, interpolate ? GRIORA_NearestNeighbour : GRIORA_Bilinear, &pixel, nullptr) != CE_None)
+				if (poBand->InterpolateAtPoint(Pixel, Line, RA, &pixel, nullptr) != CE_None)
 						// tbd: handle GRIORA_Cubic, GRIORA_CubicSpline
 					stop("Error in InterpolateAtPoint()");
 #else
-				if (interpolate)
-					// stop("interpolate not implemented");
+				if (RA == GRIORA_Cubic || RA == GRIORA_CubicSpline)
+					stop("cubic or cubicspline requires GDAL >= 3.10.0");
+				if (RA == GRIORA_Bilinear)
 					pixel = get_bilinear(poBand, Pixel, Line, iPixel, iLine,
 						xsize, ysize, nodata_set, nodata);
-				else if (GDALRasterIO(poBand, GF_Read, iPixel, iLine, 1, 1,
-						&pixel, 1, 1, GDT_CFloat64, 0, 0) != CE_None)
+				else if (poBand->RasterIO(GF_Read, iPixel, iLine, 1, 1,
+						&pixel, 1, 1, GDT_CFloat64, 0, 0, NULL) != CE_None)
 					stop("Error reading!");
 #endif
 				if (nodata_set && pixel == nodata)