using remote (vsicurl) calculations on sentinel-2 data with the {stars} package?

Question

I'm using rstac to access Sentinel-2 data in a desired bounding box and date range and computing NDVI. This is relatively[*] clean and straight forward for me when using the {terra} package, but I'd like to use the {stars} syntax instead (more on why further down).

First, a quick {rstac} query to get URLs to the data:

library(rstac)
library(sf)
library(stars)
library(terra) 

bbox <- st_bbox(c(xmin=-86.94663, ymin=33.43930, 
                  xmax=-86.67684, ymax=33.62239),
                crs=4326) # Birmingham, AL
matches <- 
  stac("https://planetarycomputer.microsoft.com/api/stac/v1/") |>
  stac_search(collections = "sentinel-2-l2a",
              bbox = bbox,
              datetime = "2019-06-01/2019-08-01") |>
  get_request() |>
  items_sign(sign_fn = sign_planetary_computer())

This returns lots of matches and appropriate metadata, for simplicity I'll just pick one (#59) with low eo:cloudcover from the properties metadata:

best_img <- matches$features[[59]] 

Now I'll use the vsicurl mechanism to access the red and near-infrared bands without downloading the whole file. The images are much larger than my search box, so I'll also want to crop out those pixels I won't be using to avoid unnecessary computation.

My first step is ugly. To crop my image using {terra}, I need a SpatVec cookie-cutter to pass to crop(). I already have bbox above as the sf-type bounding box, I do the following to get it in the projection that matches the Sentinel2 asset, but this feels very hacky. I'd love a concise, pure-terra version but this works:

red <- read_stars( paste0("/vsicurl/", best_img$assets$B04$href) )
bbox_proj <- bbox |> st_as_sfc() |> st_transform(st_crs(red)) |> vect()

Anyway, cropping vector in hand, NDVI calculation in terra is quite elegant and fast (on a good network connection using minimal RAM):

red <- rast( paste0("/vsicurl/", best_img$assets$B04$href) ) |> crop(bbox_proj)
nir <- rast( paste0("/vsicurl/", best_img$assets$B08$href) ) |> crop(bbox_proj)
ndvi_fun <- function(x, y) (x - y) / (x + y)
ndvi <- lapp(c(nir, red), fun = ndvi_fun)
ndvi  |> plot()

So my main question is what is the equivalent syntax to the identical calculation using {stars} ? So far I have only come up with the code below, which notably only works when using the local copy I had to first create, (and thus not surprisingly, is also much slower!)

# ugh why can't we combine these in a single read_stars?
red <- read_stars( paste0("/vsicurl/", best_img$assets$B04$href) )
nir <- read_stars( paste0("/vsicurl/", best_img$assets$B08$href) )

bbox_proj <- bbox |> st_as_sfc() |> st_transform(st_crs(red))

# combine 'by hand' and then crop... 
remote <- c(r1,r2, along=3) |> st_crop(bbox_proj)

# ugh! ugh! why do we have to use local copy for this to work??
stars::write_stars(remote, "test.tif")
local <- read_stars("test.tif")

# Um, I think this is correct NDVI method, hope I didn't reverse the bands...
# also surprised this is considerably slower and uses much more RAM
calc_ndvi <- function(x) (x[2] - x[1])/(x[2] + x[1])
ndvi <- st_apply(local, 1:2,  FUN = calc_ndvi)
plot(ndvi, col =  rgb(0, (0:100)/100, 0))

I'm surely missing something in my stars syntax that is resulting in this being slower, somewhat more verbose to express, and to only work when st_apply() operates on the local copy and not the remote object.

style & preferences

it's maybe reasonable to ask why do this in {stars} if it is working in {terra} -- part of this is me learning stars, but I am also an instructor and always find it cumbersome to teach my students both sf and terra syntax. terra also does not warn about miss-matched CRS if you try the above crop without reprojecting the bounding box CRS, which is a common error for students. In both cases I find the re-projection of the bounding box for the crop to be more cumbersome than I like too. In particular it seems awkward to access the file 'twice', once to read the crs and then again to crop, I expect a more elegant syntax is possible but haven't figured it out.

Answer 1

This does not answer your question but here is a more concise approach with terra, using the project<SpatExtent> method that I introduced in terra 1.6.31 (currently the development version).

library(rstac)
library(terra) 
#terra 1.6.31

bbox <- c(xmin=-86.94663, ymin=33.43930, xmax=-86.67684, ymax=33.62239)

matches <- stac("https://planetarycomputer.microsoft.com/api/stac/v1/") |>
        stac_search(collections = "sentinel-2-l2a",
            bbox = bbox, datetime = "2019-06-01/2019-08-01") |>
        get_request() |>  items_sign(sign_fn = sign_planetary_computer())

best_img <- matches$features[[59]] 

Get the first data source, and project the lon/lat search extent to the coordinate reference system of the data. Note the use of (new) argument xy=TRUE to indicate that the bbox numbers are in (xmin,ymin,xmax,ymax ) order whereas "terra" by default expects (xmin,xmax,ymin,ymax).

red <- rast( paste0("/vsicurl/", best_img$assets$B04$href) ) 
e <- ext(bbox, xy=TRUE) |> 
     project("+proj=longlat", crs(red))

Crop the first data source and download and crop others

red <- crop(red, e)
nir <- rast( paste0("/vsicurl/", best_img$assets$B08$href) ) |> crop(e)

And use the data

ndvi_fun <- function(x, y) (x - y) / (x + y)
ndvi <- lapp(c(nir, red), fun = ndvi_fun)

The above use of lapp is great, but, for a simple function like that, the below is faster

ndvi <- (red-nir) / (red+nir)

If you are going with lapp, you might consider doing that like this

rednir <- paste0("/vsicurl/", c(best_img$assets$B04$href, best_img$assets$B08$href)) |> 
          rast() |> crop(e, names=c("red", "nir"))    
ndvi <- lapp(rednir, ndvi_fun)

Answer 2

Just a follow-up here to note that {stars} code works entirely as expected when using the latest GitHub version (0.5.7). The benchmark is slightly faster than {terra} but both are so quick that it is hard to compare. Adjusting the downsampling parameter makes it easy to reduce memory use associated with the final plot, and the remainder has very low resource footprint.

Both the {stars} and {terra} performance and syntax is pretty impressive here.

library(rstac)
library(stars)

## STAC Search
coords <- c(xmin=-86.94663, ymin=33.43930, 
            xmax=-86.67684, ymax=33.62239)
bbox <- st_bbox(coords, crs=4326) # Birmingham, AL
matches <- 
  stac("https://planetarycomputer.microsoft.com/api/stac/v1/") |>
  stac_search(collections = "sentinel-2-l2a",
              bbox = bbox,
              datetime = "2019-06-01/2019-08-01") |>
  get_request() |>
  items_sign(sign_fn = sign_planetary_computer())
best_img <- matches$features[[59]] # I picked one with low cloud clover
B04 <- paste0("/vsicurl/", best_img$assets$B04$href)
B08 <- paste0("/vsicurl/", best_img$assets$B08$href)

## Time to read, transform, crop, and compute NDVI!
bench::bench_time({

img <- read_stars(c(B04,B08), along=3)  
bbox_proj <- bbox |> st_as_sfc() |> st_transform(st_crs(img))
remote <- img |> st_crop(bbox_proj)
calc_ndvi <- function(x) (x[2] - x[1])/(x[2] + x[1])
ndvi <- st_apply(remote, 1:2,  FUN = calc_ndvi)

})

## above is all 'lazy eval' and nearly instantaneous. 
## downsample by higher power of 2 for lower resolution / less ram
bench::bench_time(

plot(ndvi, col =  rgb(0, (0:100)/100, 0), downsample = 2^3)

)

Footnote

The stac API search seems to fail for some users (outside US/Canada), not sure why that would be. The sentinel-2-l2a data should be available from other catalogs as well which have a STAC catalog but not a STAC API, e.g.