Harmonized Landsat-8 Sentinel-2 Ideas

[h-a-graham]

Hey @mikemahoney218! So, the HLS data (https://hls.gsfc.nasa.gov/) are a game changer for getting a high return rate of images and building lovely mosaics - the provided cloud masks are also excellent.

However, it is a little different compared to working with other STAC collections as there are two distinct collections and these have different band numbers and mappings.

The following reprex works great but I wanted to get your thoughts on a few things:

1. Should we add more high-level get_stac_data functions to the package?
2. And if so, how would we tackle this?
   • just offer the two collections as get_hsl_ss_imagery and get_hsl_sl_imagery
   • Three versions? the two above and get_hsl_imagery to wrap both of the others?
3. Or does adding all of this stuff become burdensome, and are better of just chucking something like the following reprex in a vignette?

I'm happy to write a PR but wanted to get your thoughts first.

So yeah let me know your thoughts and as always, thanks for all your fab work on this package!

# Description: This script is used to download HLS data from NASA's Earthdata
# STAC API and build a cloud-free composite image.

aoi <- sf::st_bbox(
  c(
    xmin = 35.401471, ymin = -3.352265,
    xmax = 35.492570, ymax = -3.261740
  ),
  crs = "EPSG:4326"
) |>
  sf::st_as_sfc()


# band mapping for HLS SL data
hlssl_band_mapping <- c(
  B01 = "A", B02 = "B", B03 = "G", B04 = "R",
  B05 = "N2", B06 = "S1", B07 = "S2",
  B09 = "C", B10 = "T1",
  B11 = "T2"
)

# band mapping for HLS SS data
hlsss_band_mapping <- c(
  B01 = "A", B02 = "B", B03 = "G", B04 = "R",
  B05 = "RE1", B06 = "RE2", B07 = "RE3",
  B08 = "N", B8A = "N2", B09 = "WV",
  B10 = "C", B11 = "S1", B12 = "S2"
)

Mask function for HLS data
raster a SpatRaster object. The mask asset of a given HLS item.
bits a numeric vector. The bits to be used for masking.
a SpatRaster with logical values.
This function is used to mask HLS data using the Fmask band.
the bit values relate to the folloeing feautres in the Fmask band:
0: Cirrus, 1: Cloud, 2: Adjacent to cloud/shadow, 3: Cloud shadow,
4: Snow/ice, 5: Water, 6-7: aerosol level
The defaults are quite agressive but do a good job of removing clouds and
shadows.

hls_mask_fun <- function(raster, bits = c(0, 1, 2, 3, 4, 5, 7)) {
  bm <- sum(
    terra::rast(lapply(
      bits,
      \(b) {
        terra::app(
          raster,
          function(x) bitwAnd(x, bitwShiftL(1, b)) > 0
        )
      }
    )),
    na.rm = TRUE
  )
  bm == 0
}

Cloud cover filter for HLS data
A function factory creating a rsi query function that filters HLS data
based on cloud cover proportion.
eo_cloud_cover a numeric value. The maximum cloud cover proportion to
filter HLS data by. must be between 0 and 100.
a function that can be used as the query_function argument in
rsi::get_stac_data.

hls_cloud_cover_filter <- function(eo_cloud_cover = 25) {
  if (!rlang::is_bare_numeric(eo_cloud_cover)) {
    rlang::abort("`eo_cloud_cover`` must be a numeric value",
      class = "eo_cloud_cover_not_numeric"
    )
  }
  if (rlang::is_false(eo_cloud_cover >= 0 && eo_cloud_cover <= 100)) {
    rlang::abort("`eo_cloud_cover` must be between 0 and 100",
      class = "eo_cloud_cover_out_of_range"
    )
  }
  function(bbox, stac_source, collection, start_date, end_date, ...) {
    request <- rstac::stac(stac_source) |>
      rstac::stac_search(
        collections = collection,
        bbox = bbox[1:4],
        datetime = paste(start_date, end_date, sep = "/")
      )

    its <- rstac::items_fetch(rstac::post_request(request)) |>
      rstac::items_filter(filter_fn = function(x) {
        x$properties$`eo:cloud_cover` <= eo_cloud_cover
      })

    if (rstac::items_length(its) == 0) {
      cli::cli_abort("No items found with cloud cover <= {eo_cloud_cover}")
    } else {
      cli::cli_alert_info(
        "Found {rstac::items_length(its)} items
      with cloud cover < {eo_cloud_cover}"
      )
    }

    return(its)
  }
}

# GDAL options for downloading HLS data - you must set the EARTHDATA_TOKEN
# environment variable to your Earthdata token
nasa_ed_gdal_config_options <- c(rsi::rsi_gdal_config_options(),
  GDAL_HTTP_HEADERS = glue::glue(
    "Authorization: Bearer {Sys.getenv('EARTHDATA_TOKEN')}"
  )
)

# shared args
target_dir <- "hls-tz-24"
if (!dir.exists(target_dir)) dir.create(target_dir)
start_date <- "2024-01-01"
end_date <- "2024-02-29"
mask_band <- "Fmask"


# get Landsat data
hlssl_paths <- rsi::get_stac_data(
  aoi = aoi,
  start_date = start_date,
  end_date = end_date,
  asset_names = hlssl_band_mapping,
  stac_source = "https://cmr.earthdata.nasa.gov/cloudstac/LPCLOUD/",
  collection = "HLSL30.v2.0",
  query_function = hls_cloud_cover_filter(eo_cloud_cover = 25),
  output_filename = file.path(target_dir, "TZ-HLS-SL-cc25.tif"),
  mask_band = mask_band,
  mask_function = hls_mask_fun,
  composite_function = NULL,
  gdal_config_options = nasa_ed_gdal_config_options
)
#> ℹ Found 2 items
#> with cloud cover < 25

# get Sentinel-2 data
hlsss_paths <- rsi::get_stac_data(
  aoi = aoi,
  start_date = start_date,
  end_date = end_date,
  asset_names = hlsss_band_mapping,
  stac_source = "https://cmr.earthdata.nasa.gov/cloudstac/LPCLOUD/",
  collection = "HLSS30.v2.0",
  query_function = hls_cloud_cover_filter(eo_cloud_cover = 25),
  output_filename = file.path(target_dir, "TZ-HLS-SS-cc25.tif"),
  mask_band = mask_band,
  mask_function = hls_mask_fun,
  composite_function = NULL,
  gdal_config_options = nasa_ed_gdal_config_options
)
#>   |                                                                              |======================================================================| 100%
#> ℹ Found 4 items
#> with cloud cover < 25
#> Warning: Failed to download HLS.S30.T36MYB.2024036T075121.v2.0 from
#> 2024-02-05T08:10:50.306Z


# get names for unique and shared bands between sentinel and landasat
hlsss_unique <- setdiff(hlsss_band_mapping, hlssl_band_mapping)
hlssl_unique <- setdiff(hlssl_band_mapping, hlsss_band_mapping)
hls_shared <- intersect(hlssl_band_mapping, hlsss_band_mapping)

Build a mosaic from a list of rasters and specified band names
param rlist a character vector of raster file paths
param bands a character vector of band names to extract from the rasters
return a SpatRaster object

structured_mosaic <- function(rlist, bands) {
  if (!is.character(rlist)) rlang::abort("`rlist` must be a character vector")
  if (!is.character(bands)) rlang::abort("`bands` must be a character vector")
  r <- terra::mosaic(
    x = terra::sprc(lapply(rlist, \(x) terra::rast(x)[[bands]])),
    fun = "median",
  )
  names(r) <- bands
  return(r)
}

# construct the final mosaic and order it.
tz_mos <- c(
  structured_mosaic(c(hlssl_paths, hlsss_paths), hls_shared),
  structured_mosaic(hlssl_paths, hlssl_unique),
  structured_mosaic(hlsss_paths, hlsss_unique)
) |>
  terra::subset(
    c(
      "A", "B", "G", "R", "RE1", "RE2", "RE3", "N", "N2", "S1", "S2",
      "WV", "C", "T1", "T2"
    ),
    filename = file.path(target_dir, "TZ-HLS-mosaic.tif"),
    overwrite = TRUE
  )
#> ext 0: 35.4015 35.4926 -3.35228 -3.26174 
#> ext 1: 35.4015 35.4926 -3.35228 -3.26174 
#> ext 2: 35.4015 35.4926 -3.35228 -3.26174 
#> ext 3: 35.4015 35.4926 -3.35228 -3.26174 
#> ext 4: 35.4015 35.4926 -3.35228 -3.26174 
#> ext 0: 35.4015 35.4926 -3.35228 -3.26174 
#> ext 1: 35.4015 35.4926 -3.35228 -3.26174 
#> ext 0: 35.4015 35.4926 -3.35228 -3.26174 
#> ext 1: 35.4015 35.4926 -3.35228 -3.26174 
#> ext 2: 35.4015 35.4926 -3.35228 -3.26174

# plotting
terra::plotRGB(tz_mos,
  r = 4, g = 3, b = 2,
  stretch = "lin", zlim = c(0.003, 0.1)
)

terra::plot(tz_mos,
  col = hcl.colors(100, "Mako")
)

^{Created on 2024-10-16 with reprex v2.1.1}

[h-a-graham]

Another (probably better, but longer) example here: https://gist.github.com/h-a-graham/35b6ee6b5fd4fac75eb87629265fb8ee
using more appropriate CRS.