diff --git a/DESCRIPTION b/DESCRIPTION
index b54ad6b..a4f5f09 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -14,6 +14,8 @@ Suggests:
     mockery,
     knitr,
     kableExtra,
+    leafem,
+    leaflet,
     rmarkdown,
     mapview,
     rstac,
diff --git a/R/ngr_spk_stac_calc.R b/R/ngr_spk_stac_calc.R
index b8f1fcc..38d73eb 100644
--- a/R/ngr_spk_stac_calc.R
+++ b/R/ngr_spk_stac_calc.R
@@ -1,10 +1,10 @@
-#' Spectral index calculation from a STAC item
+#' Retrieve and optionally calculate spectral indices from a STAC item
 #'
-#' Compute a spectral index for a single STAC item. By default, the function reads
-#' red and near-infrared (NIR) assets to compute NDVI, but this is only the
-#' default configuration. The function is designed to read any two or three
-#' assets and apply different spectral index calculations as they are added in
-#' the future.
+#' Retrieve raster assets from a single STAC item with optional clipping to an
+#' AOI and optional spectral index calculation. When `calc = NULL`, the function
+#' simply returns the first asset (`asset_a`), optionally clipped. When `calc`
+#' specifies an index (e.g., `"ndvi"`), the function reads the required assets
+#' and computes the index.
 #'
 #' The default asset names (`"red"` and `"nir08"`) align with common conventions
 #' in Landsat Collection 2 Level-2 STAC items (e.g. the `landsat-c2-l2` collection
@@ -16,15 +16,25 @@
 #' @param aoi [sf] or [NULL] Optional. An AOI object used to restrict reads (by
 #' bbox) and to crop and mask the output. If `NULL`, the full assets are read and
 #' no crop/mask is applied. Default is `NULL`.
-#' @param asset_a [character] Optional. A single string giving the asset name for
-#' the first input band. For NDVI this corresponds to the red band. Default is
-#' `"red"`.
-#' @param asset_b [character] Optional. A single string giving the asset name for
-#' the second input band. For NDVI this corresponds to the NIR band. Default is
-#' `"nir08"`.
+#' @param asset_a [character] A single string giving the asset name for
+#' the first (or only) input band. For NDVI this corresponds to the red band.
+#' Default is `"red"`.
+#' @param asset_b [character] or [NULL] Optional. A single string giving the asset
+#' name for the second input band. For NDVI this corresponds to the NIR band.
+#' For RGB this is the green band. Required when `calc` is not `NULL`. Default
+#' is `"nir08"`.
 #' @param asset_c [character] or [NULL] Optional. A single string giving the asset
-#' name for an optional third input band, used by some calculations. Default is
-#' `NULL`.
+#' name for the third input band. Required when `calc = "rgb"` (blue band).
+#' Default is `NULL`.
+#' @param calc [character] or [NULL] Optional. The calculation to perform:
+#' \itemize{
+#'   \item `NULL`: No calculation; return `asset_a` only (optionally clipped).
+#'   \item `"ndvi"`: Normalized Difference Vegetation Index using
+#'     `(asset_b - asset_a) / (asset_b + asset_a)`.
+#'   \item `"rgb"`: Stack three bands into an RGB composite. Requires `asset_a`
+#'     (red), `asset_b` (green), and `asset_c` (blue).
+#' }
+#' Default is `"ndvi"`.
 #' @param vsi_prefix [character] Optional. A single string giving the GDAL VSI
 #' prefix used by [terra::rast()] to enable streaming, windowed reads over HTTP.
 #' Exposed to allow alternative VSI backends (e.g. `/vsis3/`). Default is
@@ -36,21 +46,22 @@
 #'
 #' @details
 #' This function expects `feature` to provide the required `assets` referenced by
-#' `asset_a`, `asset_b`, and (if used) `asset_c`. When `aoi` is not `NULL`,
-#' `feature` must also have `properties$"proj:epsg"` so the AOI can be transformed
-#' for windowed reads.
+#' `asset_a` and, when `calc` is not `NULL`, also `asset_b` (and `asset_c` for
+#' RGB). When `aoi` is not `NULL`, `feature` must also have
+#' `properties$"proj:epsg"` so the AOI can be transformed for windowed reads.
 #'
-#' The calculation performed is controlled by `calc`. By default, `calc = "ndvi"`
-#' computes the Normalized Difference Vegetation Index using `(b - a) / (b + a)`.
+#' When `calc = NULL`, only `asset_a` is read and returned (optionally clipped
+#' to the AOI). When `calc = "ndvi"`, the function computes the Normalized
+#' Difference Vegetation Index using `(b - a) / (b + a)`. When `calc = "rgb"`,
+#' the three assets are stacked into an RGB composite at native resolution.
 #' This structure allows additional spectral indices (e.g. NDWI, EVI) to be added
 #' in the future without changing the data access or cropping logic.
 #'
 #' When `aoi` is provided, reading is limited to the AOI bounding box in the
 #' feature's projected CRS, and then cropped/masked to the AOI. When `aoi = NULL`,
-#' the full assets are read and the calculated raster is returned for the full
-#' raster extent.
+#' the full assets are read and returned for the full raster extent.
 #'
-#' @return A [terra::SpatRaster] with calculated index values.
+#' @return A [terra::SpatRaster] with asset values or calculated index values.
 #'
 #' @section Production-ready alternative:
 #' This function is a proof of concept demonstrating how to perform raster
@@ -118,6 +129,36 @@
 #'
 #' ndvi_list <- ndvi_list |>
 #'   purrr::set_names(purrr::map_chr(items$features, "id"))
+#'
+#' # Retrieve a single asset (no calculation) from Sentinel-2
+#' stac_query_s2 <- rstac::stac(stac_url) |>
+#'   rstac::stac_search(
+#'     collections = "sentinel-2-l2a",
+#'     datetime = "2023-07-01/2023-07-31",
+#'     intersects = sf::st_geometry(aoi)[[1]],
+#'     limit = 10
+#'   ) |>
+#'   rstac::ext_filter(`eo:cloud_cover` <= 10)
+#'
+#' items_s2 <- stac_query_s2 |>
+#'   rstac::post_request() |>
+#'   rstac::items_fetch() |>
+#'   rstac::items_sign_planetary_computer()
+#'
+#' # Get just the visual (RGB) asset clipped to AOI, no calculation
+#' visual_list <- items_s2$features |>
+#'   purrr::map(ngr_spk_stac_calc, aoi = aoi, asset_a = "visual", calc = NULL)
+#'
+#' # Build RGB composite from native 10m bands (higher resolution than visual)
+#' rgb_list <- items_s2$features |>
+#'   purrr::map(
+#'     ngr_spk_stac_calc,
+#'     aoi = aoi,
+#'     asset_a = "B04",
+#'     asset_b = "B03",
+#'     asset_c = "B02",
+#'     calc = "rgb"
+#'   )
 #' }
 
 ngr_spk_stac_calc <- function(
@@ -132,13 +173,24 @@ ngr_spk_stac_calc <- function(
     timing = FALSE
 ) {
   chk::chk_string(asset_a)
-  chk::chk_string(asset_b)
+  if (!is.null(asset_b)) chk::chk_string(asset_b)
   if (!is.null(asset_c)) chk::chk_string(asset_c)
-  chk::chk_string(calc)
+  if (!is.null(calc)) chk::chk_string(calc)
   chk::chk_string(vsi_prefix)
   chk::chk_flag(quiet)
   chk::chk_flag(timing)
 
+
+  # asset_b is required when calc is specified
+  if (!is.null(calc) && is.null(asset_b)) {
+    cli::cli_abort("`asset_b` is required when `calc` is not NULL.")
+  }
+
+  # asset_c is required when calc = "rgb"
+  if (!is.null(calc) && calc == "rgb" && is.null(asset_c)) {
+    cli::cli_abort("`asset_c` is required when `calc = \"rgb\"`.")
+  }
+
   if (!is.list(feature) || is.null(feature$assets) || !is.list(feature$assets)) {
     cli::cli_abort("`feature` must be a single STAC item (a list with an `assets` element).")
   }
@@ -156,31 +208,59 @@ ngr_spk_stac_calc <- function(
   }
 
   a_href <- feature$assets[[asset_a]]$href
-  b_href <- feature$assets[[asset_b]]$href
   if (is.null(a_href) || !nzchar(a_href)) {
     cli::cli_abort("Missing asset href for `asset_a` `{asset_a}` in item: {id}.")
   }
-  if (is.null(b_href) || !nzchar(b_href)) {
-    cli::cli_abort("Missing asset href for `asset_b` `{asset_b}` in item: {id}.")
+  a_url <- paste0(vsi_prefix, a_href)
+
+  # Only get asset_b when calc is specified
+  b_url <- NULL
+  if (!is.null(calc)) {
+    b_href <- feature$assets[[asset_b]]$href
+    if (is.null(b_href) || !nzchar(b_href)) {
+      cli::cli_abort("Missing asset href for `asset_b` `{asset_b}` in item: {id}.")
+    }
+    b_url <- paste0(vsi_prefix, b_href)
   }
 
-  a_url <- paste0(vsi_prefix, a_href)
-  b_url <- paste0(vsi_prefix, b_href)
+  # Only get asset_c when calc = "rgb"
+  c_url <- NULL
+  if (!is.null(calc) && calc == "rgb") {
+    c_href <- feature$assets[[asset_c]]$href
+    if (is.null(c_href) || !nzchar(c_href)) {
+      cli::cli_abort("Missing asset href for `asset_c` `{asset_c}` in item: {id}.")
+    }
+    c_url <- paste0(vsi_prefix, c_href)
+  }
 
   if (is.null(aoi)) {
-    if (!quiet) cli::cli_alert_info("aoi is NULL: reading full assets: {id}")
+    if (!quiet) cli::cli_alert_info("aoi is NULL: reading full asset(s): {id}")
 
     if (!quiet) cli::cli_alert_info("read asset_a: {id}")
     t0 <- if (timing) base::proc.time()[[3]] else NA_real_
     a <- terra::rast(a_url)
     if (timing && !quiet) cli::cli_alert_info("read asset_a elapsed (s): {format(base::proc.time()[[3]] - t0, digits = 3)}")
 
+    # If no calculation, return asset_a directly
+    if (is.null(calc)) {
+      return(a)
+    }
+
     if (!quiet) cli::cli_alert_info("read asset_b: {id}")
     t1 <- if (timing) base::proc.time()[[3]] else NA_real_
     b <- terra::rast(b_url)
     if (timing && !quiet) cli::cli_alert_info("read asset_b elapsed (s): {format(base::proc.time()[[3]] - t1, digits = 3)}")
 
-    out <- .ngr_spk_calc(calc, a = a, b = b, c = NULL)
+    # Read asset_c if needed (for rgb)
+    c_rast <- NULL
+    if (!is.null(c_url)) {
+      if (!quiet) cli::cli_alert_info("read asset_c: {id}")
+      t2 <- if (timing) base::proc.time()[[3]] else NA_real_
+      c_rast <- terra::rast(c_url)
+      if (timing && !quiet) cli::cli_alert_info("read asset_c elapsed (s): {format(base::proc.time()[[3]] - t2, digits = 3)}")
+    }
+
+    out <- .ngr_spk_calc(calc, a = a, b = b, c = c_rast)
     return(out)
   }
 
@@ -192,12 +272,29 @@ ngr_spk_stac_calc <- function(
   a <- terra::rast(a_url, win = win)
   if (timing && !quiet) cli::cli_alert_info("read asset_a elapsed (s): {format(base::proc.time()[[3]] - t0, digits = 3)}")
 
+  # If no calculation, crop/mask asset_a and return
+  if (is.null(calc)) {
+    aoi_v <- aoi |>
+      sf::st_transform(sf::st_crs(a)) |>
+      terra::vect()
+    return(terra::crop(a, aoi_v) |> terra::mask(aoi_v))
+  }
+
   if (!quiet) cli::cli_alert_info("read asset_b: {id}")
   t1 <- if (timing) base::proc.time()[[3]] else NA_real_
   b <- terra::rast(b_url, win = win)
   if (timing && !quiet) cli::cli_alert_info("read asset_b elapsed (s): {format(base::proc.time()[[3]] - t1, digits = 3)}")
 
-  out <- .ngr_spk_calc(calc, a = a, b = b, c = NULL)
+  # Read asset_c if needed (for rgb)
+  c_rast <- NULL
+  if (!is.null(c_url)) {
+    if (!quiet) cli::cli_alert_info("read asset_c: {id}")
+    t2 <- if (timing) base::proc.time()[[3]] else NA_real_
+    c_rast <- terra::rast(c_url, win = win)
+    if (timing && !quiet) cli::cli_alert_info("read asset_c elapsed (s): {format(base::proc.time()[[3]] - t2, digits = 3)}")
+  }
+
+  out <- .ngr_spk_calc(calc, a = a, b = b, c = c_rast)
 
   aoi_v <- aoi |>
     sf::st_transform(sf::st_crs(out)) |>
@@ -224,6 +321,10 @@ ngr_spk_stac_calc <- function(
     return((b - a) / (b + a))
   }
 
+  if (calc == "rgb") {
+    return(c(a, b, c))
+  }
+
   cli::cli_abort("Unknown calc: {calc}.")
 }
 
diff --git a/man/ngr_spk_stac_calc.Rd b/man/ngr_spk_stac_calc.Rd
index 1185f35..4426b54 100644
--- a/man/ngr_spk_stac_calc.Rd
+++ b/man/ngr_spk_stac_calc.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/ngr_spk_stac_calc.R
 \name{ngr_spk_stac_calc}
 \alias{ngr_spk_stac_calc}
-\title{Spectral index calculation from a STAC item}
+\title{Retrieve and optionally calculate spectral indices from a STAC item}
 \usage{
 ngr_spk_stac_calc(
   feature,
@@ -24,17 +24,28 @@ list).}
 bbox) and to crop and mask the output. If \code{NULL}, the full assets are read and
 no crop/mask is applied. Default is \code{NULL}.}
 
-\item{asset_a}{\link{character} Optional. A single string giving the asset name for
-the first input band. For NDVI this corresponds to the red band. Default is
-\code{"red"}.}
+\item{asset_a}{\link{character} A single string giving the asset name for
+the first (or only) input band. For NDVI this corresponds to the red band.
+Default is \code{"red"}.}
 
-\item{asset_b}{\link{character} Optional. A single string giving the asset name for
-the second input band. For NDVI this corresponds to the NIR band. Default is
-\code{"nir08"}.}
+\item{asset_b}{\link{character} or \link{NULL} Optional. A single string giving the asset
+name for the second input band. For NDVI this corresponds to the NIR band.
+For RGB this is the green band. Required when \code{calc} is not \code{NULL}. Default
+is \code{"nir08"}.}
 
 \item{asset_c}{\link{character} or \link{NULL} Optional. A single string giving the asset
-name for an optional third input band, used by some calculations. Default is
-\code{NULL}.}
+name for the third input band. Required when \code{calc = "rgb"} (blue band).
+Default is \code{NULL}.}
+
+\item{calc}{\link{character} or \link{NULL} Optional. The calculation to perform:
+\itemize{
+\item \code{NULL}: No calculation; return \code{asset_a} only (optionally clipped).
+\item \code{"ndvi"}: Normalized Difference Vegetation Index using
+\code{(asset_b - asset_a) / (asset_b + asset_a)}.
+\item \code{"rgb"}: Stack three bands into an RGB composite. Requires \code{asset_a}
+(red), \code{asset_b} (green), and \code{asset_c} (blue).
+}
+Default is \code{"ndvi"}.}
 
 \item{vsi_prefix}{\link{character} Optional. A single string giving the GDAL VSI
 prefix used by \code{\link[terra:rast]{terra::rast()}} to enable streaming, windowed reads over HTTP.
@@ -48,14 +59,14 @@ Exposed to allow alternative VSI backends (e.g. \verb{/vsis3/}). Default is
 for reads. Default is \code{FALSE}.}
 }
 \value{
-A \link[terra:SpatRaster-class]{terra::SpatRaster} with calculated index values.
+A \link[terra:SpatRaster-class]{terra::SpatRaster} with asset values or calculated index values.
 }
 \description{
-Compute a spectral index for a single STAC item. By default, the function reads
-red and near-infrared (NIR) assets to compute NDVI, but this is only the
-default configuration. The function is designed to read any two or three
-assets and apply different spectral index calculations as they are added in
-the future.
+Retrieve raster assets from a single STAC item with optional clipping to an
+AOI and optional spectral index calculation. When \code{calc = NULL}, the function
+simply returns the first asset (\code{asset_a}), optionally clipped. When \code{calc}
+specifies an index (e.g., \code{"ndvi"}), the function reads the required assets
+and computes the index.
 }
 \details{
 The default asset names (\code{"red"} and \code{"nir08"}) align with common conventions
@@ -64,19 +75,20 @@ on Planetary Computer), but are fully parameterized to support other
 collections with different band naming schemes.
 
 This function expects \code{feature} to provide the required \code{assets} referenced by
-\code{asset_a}, \code{asset_b}, and (if used) \code{asset_c}. When \code{aoi} is not \code{NULL},
-\code{feature} must also have \code{properties$"proj:epsg"} so the AOI can be transformed
-for windowed reads.
-
-The calculation performed is controlled by \code{calc}. By default, \code{calc = "ndvi"}
-computes the Normalized Difference Vegetation Index using \code{(b - a) / (b + a)}.
+\code{asset_a} and, when \code{calc} is not \code{NULL}, also \code{asset_b} (and \code{asset_c} for
+RGB). When \code{aoi} is not \code{NULL}, \code{feature} must also have
+\code{properties$"proj:epsg"} so the AOI can be transformed for windowed reads.
+
+When \code{calc = NULL}, only \code{asset_a} is read and returned (optionally clipped
+to the AOI). When \code{calc = "ndvi"}, the function computes the Normalized
+Difference Vegetation Index using \code{(b - a) / (b + a)}. When \code{calc = "rgb"},
+the three assets are stacked into an RGB composite at native resolution.
 This structure allows additional spectral indices (e.g. NDWI, EVI) to be added
 in the future without changing the data access or cropping logic.
 
 When \code{aoi} is provided, reading is limited to the AOI bounding box in the
 feature's projected CRS, and then cropped/masked to the AOI. When \code{aoi = NULL},
-the full assets are read and the calculated raster is returned for the full
-raster extent.
+the full assets are read and returned for the full raster extent.
 }
 \section{Production-ready alternative}{
 
@@ -137,6 +149,36 @@ ndvi_list <- items$features |>
 
 ndvi_list <- ndvi_list |>
   purrr::set_names(purrr::map_chr(items$features, "id"))
+
+# Retrieve a single asset (no calculation) from Sentinel-2
+stac_query_s2 <- rstac::stac(stac_url) |>
+  rstac::stac_search(
+    collections = "sentinel-2-l2a",
+    datetime = "2023-07-01/2023-07-31",
+    intersects = sf::st_geometry(aoi)[[1]],
+    limit = 10
+  ) |>
+  rstac::ext_filter(`eo:cloud_cover` <= 10)
+
+items_s2 <- stac_query_s2 |>
+  rstac::post_request() |>
+  rstac::items_fetch() |>
+  rstac::items_sign_planetary_computer()
+
+# Get just the visual (RGB) asset clipped to AOI, no calculation
+visual_list <- items_s2$features |>
+  purrr::map(ngr_spk_stac_calc, aoi = aoi, asset_a = "visual", calc = NULL)
+
+# Build RGB composite from native 10m bands (higher resolution than visual)
+rgb_list <- items_s2$features |>
+  purrr::map(
+    ngr_spk_stac_calc,
+    aoi = aoi,
+    asset_a = "B04",
+    asset_b = "B03",
+    asset_c = "B02",
+    calc = "rgb"
+  )
 }
 }
 \seealso{
diff --git a/tests/testthat/test-ngr_spk_stac_calc.R b/tests/testthat/test-ngr_spk_stac_calc.R
index 864562a..0514112 100644
--- a/tests/testthat/test-ngr_spk_stac_calc.R
+++ b/tests/testthat/test-ngr_spk_stac_calc.R
@@ -174,6 +174,116 @@ testthat::test_that(".ngr_spk_calc errors on unknown calc type", {
   )
 })
 
+testthat::test_that(".ngr_spk_calc stacks RGB correctly", {
+  testthat::skip_if_not_installed("terra")
+
+  # Create synthetic rasters for R, G, B bands
+
+  r <- terra::rast(nrows = 2, ncols = 2, vals = c(100, 200, 150, 50))
+  g <- terra::rast(nrows = 2, ncols = 2, vals = c(110, 210, 160, 60))
+  b <- terra::rast(nrows = 2, ncols = 2, vals = c(120, 220, 170, 70))
+
+  result <- ngr:::.ngr_spk_calc("rgb", a = r, b = g, c = b)
+
+  # Should have 3 layers
+
+testthat::expect_equal(terra::nlyr(result), 3)
+
+  # Check values are preserved
+  testthat::expect_equal(as.vector(terra::values(result[[1]])), c(100, 200, 150, 50))
+  testthat::expect_equal(as.vector(terra::values(result[[2]])), c(110, 210, 160, 60))
+  testthat::expect_equal(as.vector(terra::values(result[[3]])), c(120, 220, 170, 70))
+})
+
+# Test calc = NULL (single asset retrieval) ------------------------------------
+
+testthat::test_that("ngr_spk_stac_calc allows asset_b = NULL when calc = NULL", {
+  mock_feature <- list(
+    id = "test_item",
+    assets = list(
+      visual = list(href = "https://example.com/visual.tif")
+    )
+  )
+
+  # Should not error on validation - will only fail when trying to read the URL
+  testthat::expect_error(
+    ngr::ngr_spk_stac_calc(
+      feature = mock_feature,
+      asset_a = "visual",
+      asset_b = NULL,
+      calc = NULL,
+      quiet = TRUE
+    ),
+    "cannot open|HTTP|curl"
+  )
+})
+
+testthat::test_that("ngr_spk_stac_calc errors when asset_b is NULL but calc is specified", {
+  mock_feature <- list(
+    id = "test_item",
+    assets = list(
+      red = list(href = "https://example.com/red.tif")
+    )
+  )
+
+  testthat::expect_error(
+    ngr::ngr_spk_stac_calc(
+      feature = mock_feature,
+      asset_b = NULL,
+      calc = "ndvi",
+      quiet = TRUE
+    ),
+    "asset_b.*required when.*calc.*is not NULL"
+  )
+})
+
+# Test calc = "rgb" validation -------------------------------------------------
+
+testthat::test_that("ngr_spk_stac_calc errors when calc = 'rgb' but asset_c is NULL", {
+  mock_feature <- list(
+    id = "test_item",
+    assets = list(
+      B04 = list(href = "https://example.com/B04.tif"),
+      B03 = list(href = "https://example.com/B03.tif")
+    )
+  )
+
+  testthat::expect_error(
+    ngr::ngr_spk_stac_calc(
+      feature = mock_feature,
+      asset_a = "B04",
+      asset_b = "B03",
+      asset_c = NULL,
+      calc = "rgb",
+      quiet = TRUE
+    ),
+    "asset_c.*required when.*calc.*rgb"
+  )
+})
+
+testthat::test_that("ngr_spk_stac_calc errors when calc = 'rgb' and asset_c href is missing", {
+  mock_feature <- list(
+    id = "test_item",
+    assets = list(
+      B04 = list(href = "https://example.com/B04.tif"),
+      B03 = list(href = "https://example.com/B03.tif"),
+      B02 = list(href = "")
+    )
+  )
+
+  testthat::expect_error(
+    ngr::ngr_spk_stac_calc(
+      feature = mock_feature,
+      asset_a = "B04",
+      asset_b = "B03",
+      asset_c = "B02",
+      calc = "rgb",
+      quiet = TRUE
+    ),
+    "Missing asset href.*asset_c.*B02"
+  )
+})
+
 # Test null coalescing helper --------------------------------------------------
 
 testthat::test_that("%||% returns first value when not NULL", {
diff --git a/vignettes/stac-sentinel2-ortho-timelapse.Rmd b/vignettes/stac-sentinel2-ortho-timelapse.Rmd
new file mode 100644
index 0000000..783c2e4
--- /dev/null
+++ b/vignettes/stac-sentinel2-ortho-timelapse.Rmd
@@ -0,0 +1,241 @@
+---
+title: "STAC ortho timelapse"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{STAC ortho timelapse}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  # echo = FALSE,
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+```{r setup}
+library(ngr)
+```
+
+## Planetary Computer Collections
+
+The [Microsoft Planetary Computer](https://planetarycomputer.microsoft.com/) hosts
+a STAC catalog with numerous Earth observation datasets. We can use `rstac` to
+query the available collections directly from the API:
+
+```{r pc-collections}
+stac_url <- "https://planetarycomputer.microsoft.com/api/stac/v1"
+
+# Fetch all collections from the Planetary Computer STAC catalog
+collections <- rstac::stac(stac_url) |>
+  rstac::collections() |>
+  rstac::get_request()
+
+# determin which dtasets overlap bc
+# bc_bbox <- c(-139, 48, -114, 60)
+
+ #100km buffer from the us border to not overlap all the us dtasets 
+bc_bbox <- c(-139, 49.5, -114, 60)
+
+bc_overlap <-  purrr::keep(collections$collections, ~ {
+
+  bbox <- .x$extent$spatial$bbox[[1]]
+  bbox[1] <= bc_bbox[3] && bbox[3] >= bc_bbox[1] &&
+   bbox[2] <= bc_bbox[4] && bbox[4] >= bc_bbox[2]
+  }) |>
+  purrr::map_chr("id")
+
+# Extract collection IDs and descriptions into a data frame
+collections_df <- purrr::map_df(collections$collections, ~ {
+    interval <- .x$extent$temporal$interval[[1]]
+    data.frame(
+      id = .x$id,
+      start = interval[[1]],
+      end = interval[[2]],
+      description = paste0(.x$title, ": ", paste(.x$keywords %||% NA_character_, collapse = ", ")),
+      check.names = FALSE
+    )
+  }) |>
+    dplyr::mutate(
+      start = as.Date(start),
+      end = as.Date(end),
+      # rowid = dplyr::row_number(),
+      bc_overlap = dplyr::if_else(id %in% bc_overlap, "yes", "no"), .after = id
+    ) |> 
+  dplyr::arrange(id)
+
+
+collections_df |> 
+  kableExtra::kbl() |> 
+  kableExtra::scroll_box(width = "100%", height = "500px") |> 
+  kableExtra::kable_styling(c("condensed", 
+        "responsive"), full_width = T, font_size = 10)
+```
+
+For more details on each collection, visit the
+[Planetary Computer Data Catalog](https://planetarycomputer.microsoft.com/catalog).
+
+This vignette demonstrates computing NDVI from Landsat imagery using STAC APIs. The
+`ngr_spk_stac_calc()` function provides a simple proof-of-concept approach using
+`terra`. For production workflows involving time series, composites, or data cubes,
+see the [gdalcubes](https://github.com/appelmar/gdalcubes) package.
+
+## Single Year Example
+
+Define an area of interest and query the Planetary Computer STAC catalog for
+Landsat Collection 2 Level-2 imagery with low cloud cover.
+
+```{r rstac-query-2016}
+# Define an AOI from a bounding box (WGS84) - Maxxam floodplain near bulkley lake
+bbox <- c(
+  xmin = -126.17545256019142,
+  ymin =  54.36161045287439,
+  xmax = -126.12615394008702,
+  ymax =  54.38908432381547
+)
+
+
+aoi <- sf::st_as_sfc(sf::st_bbox(bbox, crs = 4326)) |>
+  sf::st_as_sf()
+
+stac_url <- "https://planetarycomputer.microsoft.com/api/stac/v1"
+y <- 2016
+date_time <- paste0(y, "-05-01/", y, "-09-15")
+
+stac_query <- rstac::stac(stac_url) |>
+  rstac::stac_search(
+    collections = "sentinel-2-l2a",
+    datetime = date_time,
+    intersects = sf::st_geometry(aoi)[[1]],
+    limit = 200
+  ) |>
+  rstac::ext_filter(`eo:cloud_cover` <= 20)
+
+items <- stac_query |>
+  rstac::post_request() |>
+  rstac::items_fetch() |>
+  rstac::items_sign_planetary_computer()
+
+```
+
+Clip "visual" item to the bbox aoi
+
+```{r clip-2016}
+# items$features is a list of items, but the function expects a single item so we either
+# Use purrr::map() to process each feature or access a single feature directly
+r <- items$features[[1]] |> 
+  ngr_spk_stac_calc(asset_a = "visual", asset_b = NULL, calc = NULL, aoi = aoi)
+
+```
+
+Create a mapview object for each NDVI raster with a red-yellow-green color scale.
+
+```{r mapview1}
+leaflet::leaflet() |>
+  leaflet::addTiles() |>
+  leafem::addRasterRGB(r, quantiles = c(0.02, 0.98))
+
+# static view
+# terra::plotRGB(r, stretch = "lin")
+
+```
+
+## Native Resolution RGB
+
+The "visual" asset is a pre-rendered RGB that may be compressed. For higher
+resolution, we can build an RGB composite from the native 10m bands (B04=Red,
+B03=Green, B02=Blue).
+
+```{r rgb-native}
+r_rgb <- items$features[[1]] |>
+  ngr_spk_stac_calc(
+    asset_a = "B04",
+    asset_b = "B03",
+    asset_c = "B02",
+    calc = "rgb",
+    aoi = aoi
+  )
+```
+
+Compare the pre-rendered "visual" asset (left) with native 10m bands (right):
+
+```{r rgb-plot, fig.width=10, fig.height=5}
+par(mfrow = c(1, 2))
+terra::plotRGB(r, stretch = "lin", main = "Visual Asset")
+terra::plotRGB(r_rgb, stretch = "lin", main = "Native 10m RGB")
+
+```
+
+
+They look basically identical and when we check the resolution we can see that they are equal...
+
+```{r res-compare}
+# compare the resolution
+terra::res(r)
+terra::res(r_rgb)
+```
+
+## Multi-Year Comparison
+
+Query multiple years and retrieve the latest available ortho (visual asset) from
+peak growing-season for each year.
+
+```{r multi-year, eval=T}
+years <- seq(2016, 2025, by = 1)
+
+orthos_by_year <- years |>
+  purrr::map(function(y) {
+    date_time <- paste0(y, "-06-01/", y, "-07-31")
+
+    items <- rstac::stac(stac_url) |>
+      rstac::stac_search(
+        collections = "sentinel-2-l2a",
+        datetime = date_time,
+        intersects = sf::st_geometry(aoi)[[1]],
+        limit = 10
+      ) |>
+      rstac::ext_filter(`eo:cloud_cover` <= 20) |>
+      rstac::post_request() |>
+      rstac::items_fetch() |>
+      rstac::items_sign_planetary_computer()
+
+    # Return NULL if no features found for this year
+    if (length(items$features) == 0) return(NULL)
+
+    # Sort by datetime descending to get the latest scene
+    datetimes <- purrr::map_chr(items$features, ~ .x$properties$datetime)
+    latest_idx <- which.max(as.POSIXct(datetimes))
+
+    # Extract the acquisition date for the layer name
+    acq_date <- as.Date(datetimes[latest_idx])
+
+    # Get the latest feature's visual asset clipped to AOI
+    ortho <- items$features[[latest_idx]] |>
+      ngr_spk_stac_calc(asset_a = "visual", asset_b = NULL, calc = NULL, aoi = aoi)
+
+    list(date = acq_date, raster = ortho)
+  }) |>
+  purrr::compact() |>  # Remove NULLs (years with no data)
+  (\(x) purrr::set_names(x, purrr::map_chr(x, ~ as.character(.x$date))))()
+```
+
+Display the multi-year comparison as toggleable layers on a leaflet map.
+
+```{r map2, out.width="100%", eval=T}
+map <- leaflet::leaflet() |>
+  leaflet::addTiles()
+
+for (dt in names(orthos_by_year)) {
+  map <- map |>
+    leafem::addRasterRGB(orthos_by_year[[dt]]$raster, group = dt, quantiles = c(0.02, 0.98))
+}
+
+map |>
+  leaflet::addLayersControl(
+    baseGroups = names(orthos_by_year),
+    options = leaflet::layersControlOptions(collapsed = FALSE)
+  )
+```
+