Error in downloading Sentinel 2 images rstac with sf::gdal_utils

Question

I'm trying to follow the example providede here https://stacspec.org/en/tutorials/1-download-data-using-r/ to download Sentinel 2 data from a defined bounding box determined by a custom shapefile that I have drawn.

I was successfull in using rstac::assets_download() function to download Sentinel 2 data for the whole scene from several months. That's in the first part of the tutorial. In trying to adapt the second part an error is being ouputed that I don't know how to solve.

Here is my code:

require(rstac)

stac_source <- rstac::stac(
  "https://planetarycomputer.microsoft.com/api/stac/v1"
)


collections_query <- stac_source |>
  rstac::collections()


available_collections <- rstac::get_request(collections_query)
available_collections



rstac::stac_search(
  q = stac_source,
  collections = "sentinel-2-l2a",
  datetime = "2023-01-01/2023-01-31",
  limit = 999
)



ashe <- sf::read_sf("E:/temp/test_stac1.shp")

ashe_bbox <- ashe |>
  sf::st_transform(4326) |>
  sf::st_bbox()


stac_query <- rstac::stac(
  "https://planetarycomputer.microsoft.com/api/stac/v1"
) |>
  rstac::stac_search(
    collections = "sentinel-2-l2a",
    bbox = ashe_bbox,
    datetime = "2023-01-01/2023-01-31"
  ) |>
  rstac::get_request()



make_vsicurl_url <- function(base_url) {
  paste0(
    "/vsicurl", 
    "?pc_url_signing=yes",
    "&pc_collection=sentinel-2-l2a",
    "&url=",
    base_url
  )
}


lcpri_url <- rstac::assets_url(executed_stac_query, "B02")




out_file <- tempfile(fileext = ".tif")
sf::gdal_utils(
  "warp",
  source = lcpri_url,
  destination = out_file,
  options = c(
    "-t_srs", sf::st_crs(ashe)$wkt,
    "-te", sf::st_bbox(ashe)
  )
)



terra::rast(out_file) |>
  terra::plot()
ashe |>
  sf::st_geometry() |>
  plot(lwd = 3, add = TRUE)

Error in sf::gdal_utils("warp", source = lcpri_url, destination = out_file,  : 
  gdal_utils warp: an error occured
In addition: Warning messages:
1: In CPL_gdalwarp(source, destination, options, oo, doo, config_options,  :
  GDAL Error 1: HTTP error code : 404
2: In CPL_gdalwarp(source, destination, options, oo, doo, config_options,  :
  GDAL Error 1: HTTP error code : 404
3: In CPL_gdalwarp(source, destination, options, oo, doo, config_options,  :
  GDAL Error 1: HTTP error code : 404
4: In CPL_gdalwarp(source, destination, options, oo, doo, config_options,  :
  GDAL Error 1: HTTP error code : 404
5: In CPL_gdalwarp(source, destination, options, oo, doo, config_options,  :
  GDAL Error 1: HTTP error code : 404
6: In CPL_gdalwarp(source, destination, options, oo, doo, config_options,  :
  GDAL Error 1: HTTP error code : 404

With this code I had sucess download Sentinel2 data:

stac_source <- rstac::stac(
  "https://planetarycomputer.microsoft.com/api/stac/v1"
)
stac_source

str(stac_source)
rstac::get_request(stac_source)


collections_query <- stac_source |>
  rstac::collections()

collections_query


class(stac_source)



class(collections_query)



available_collections <- rstac::get_request(collections_query)
available_collections



rstac::stac_search(
  q = stac_source,
  collections = "sentinel-2-l2a",
  datetime = "2023-02-01/2023-10-31",
  limit = 999
)


ashe <- sf::read_sf("E:/temp/test_stac1.shp")

sf::st_geometry(ashe) |> plot()



ashe_bbox <- ashe |>
  sf::st_transform(4326) |>
  sf::st_bbox()

ashe_bbox


stac_query <- rstac::stac_search(
  q = stac_source,
  collections = "sentinel-2-l2a",
  bbox = ashe_bbox,
  datetime = "2023-10-01/2023-10-07"
)

executed_stac_query <- rstac::get_request(stac_query)

executed_stac_query




signed_stac_query <- rstac::items_sign(
  executed_stac_query,
  rstac::sign_planetary_computer()
)

signed_stac_query


rstac::assets_download(signed_stac_query, c("B02"), output_dir = "E:/temp/")

Here is one resulting image:

The shapefile used is avaiable here : https://drive.google.com/file/d/1NNtnglNJ3b-E9d_Z3fmUoghMYxwRM947/view?usp=sharing

Answer 1

The issue seems to be that the lcpri_url is a vector of length 6, thus resulting in an error for gdalwarp. If you were to subset to a single URL, warping the raster to the extent of ashe is successful. Also, in case you are asking for help here, it is easier to help you if you reduce the amount of code to the minimal amount needed to reproduce your issue.

library(rstac)
library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE

ashe <- sf::read_sf(system.file("shape/nc.shp", package = "sf"))[1, ]
ashe_bbox <- ashe |>
  sf::st_transform(4326) |>
  sf::st_bbox()

stac_query <- rstac::stac(
  "https://planetarycomputer.microsoft.com/api/stac/v1"
) |>
  rstac::stac_search(
    collections = "sentinel-2-l2a",
    bbox = ashe_bbox,
    datetime = "2023-01-01/2023-01-31"
  ) |>
  rstac::get_request()

make_vsicurl_url <- function(base_url) {
  paste0(
    "/vsicurl", 
    "?pc_url_signing=yes",
    "&pc_collection=sentinel-2-l2a",
    "&url=",
    base_url
  )
}

lcpri_urls <- rstac::assets_url(stac_query, "B02")
lcpri_urls <- make_vsicurl_url(lcpri_urls)
length(lcpri_urls)
#> [1] 6

# take only the first file
out_file <- tempfile(fileext = ".tif")
sf::gdal_utils(
  "warp",
  source = lcpri_urls[1],
  destination = out_file,
  options = c(
    "-t_srs", sf::st_crs(ashe)$wkt,
    "-te", sf::st_bbox(ashe)
  )
)

terra::rast(out_file) |>
  terra::plot()

plot(ashe, lwd = 3, add = TRUE, col = NA)
#> Warning in plot.sf(ashe, lwd = 3, add = TRUE, col = NA): ignoring all but the
#> first attribute

^{Created on 2024-04-09 with reprex v2.1.0}

Answer 2

I'm not sure if you must use R for this, but there's definitely a more hassle-free way of downloading the data, and clipping it to your bounding shp file. I see you have QGIS installed; you can install the SCP plugin, select the area of interest, search the Sentinel 2 data, and choose the best capture for you. You can even select which bands to download, if you don't need all twelve.