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

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.

This does not answer your question but here is a more concise approach with terra, using the project<SpatExtent> method.

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)