modis with gdalcubes

gdalcubes.R

#remotes::install_github("OldLipe/rstac@dev")
#remotes::install_github("appelmar/gdalcubes_R")

source("new_sign_planetary_computer.R")

sites <- readr::read_csv(paste0("https://github.com/eco4cast/neon4cast-noaa-download/",
                                "raw/master/noaa_download_site_list.csv"))

z <- sites |> dplyr::filter(site_id == "HARV")
y <- z$latitude
x <- z$longitude
buffer <- 1
box <- c(x - buffer, y- buffer, x+buffer, y+ buffer)


library(rstac)
library(gdalcubes)
library(stars)
matches <-
  stac("https://planetarycomputer.microsoft.com/api/stac/v1") |>
  stac_search(collections = "modis-15A2H-061",
              datetime = "2022-01-01/2022-11-08",
              bbox = c(box)) |>
  get_request() |>
  items_fetch() |>
  items_sign(sign_fn = new_sign_planetary_computer())

# test we can read a single image as expected
url1 <- matches$features[[1]]$assets$Lai_500m$href
remote <- stars::read_stars(paste0("/vsicurl/", url1))
plot(remote)

# Time for some gdalcubes
cube <- gdalcubes::stac_image_collection(matches$features, asset_names = "Lai_500m")

# make a cube_view for use on this image.
# nx, ny, and dt define a new grid size
v <- cube_view(srs = "EPSG:4326",
               extent = list(t0 = "2022-01-01", t1 = "2022-12-31",
                             left = box[1], right = box[3],
                             top = box[4], bottom = box[2]),
               nx = 2400, ny = 2400, dt= "P28D",
               aggregation = "median", resampling = "near"
               )

Q <- raster_cube(cube,v)
plot(Q, zlim=c(0,25))

patch-sign-planetarycomputer.R

ms_token_endpoint <- "https://planetarycomputer.microsoft.com/api/sas/v1/token"
ms_max_timeleft <- 300
ms_blob_name <- ".blob.core.windows.net"
ms_public_assets <- "ai4edatasetspublicassets.blob.core.windows.net"

get_ms_info <- function(asset) {
  parsed_url <- httr::parse_url(asset[["href"]])
  host_spplited <- strsplit(
    x = parsed_url[["hostname"]], split = ".", fixed = TRUE
  )
  path_spplited <- strsplit(parsed_url[["path"]], split = "/", fixed = TRUE)

  list(
    account = host_spplited[[1]][[1]],
    container = path_spplited[[1]][[1]]
  )
}

get_ms_account <- function(ms_info) {
  ms_info[["account"]]
}

get_ms_container <- function(ms_info) {
  ms_info[["container"]]
}

is_public_assets <- function(parsed_url) {
  !endsWith(parsed_url[["hostname"]], ms_blob_name) ||
    parsed_url[["hostname"]] == ms_public_assets
}


new_sign_planetary_computer <- function(..., token_url = NULL) {
  default_endpoint <- ms_token_endpoint
  if (!is.null(token_url)) {
    default_endpoint <- token_url
  }

  token <- list()

  # parse href to separate each query element, this will be used to dont
  # append the same token for an asset
  parse <- function(obj_req) {
    # transform to a datetime object
    obj_req[["msft:expiry"]] <- strptime(
      obj_req[["msft:expiry"]], "%Y-%m-%dT%H:%M:%SZ"
    )

    token_str <- paste0("?", obj_req[["token"]])
    obj_req[["token_value"]] <- httr::parse_url(token_str)[["query"]]

    obj_req
  }

  new_token <- function(asset) {
    ms_info <- get_ms_info(asset)
    account <- get_ms_account(ms_info)
    container <- get_ms_container(ms_info)

    url <- paste(
      default_endpoint, account, container, sep = "/"
    )

    tryCatch({
      res_content <- httr::content(httr::GET(url, ...), encoding = "UTF-8")
    },
    error = function(e) {
      .error("Request error. %s", e$message)
    })

    if (!"token" %in% names(res_content))
      .error("No collection found with id '%s'", asset$collection)

    token[[account]][[container]] <<- parse(res_content)
  }

  exists_token <- function(asset) {
    ms_info <- get_ms_info(asset)
    account <- get_ms_account(ms_info)
    container <- get_ms_container(ms_info)
    account %in% names(token) && container %in% names(token[[account]])
  }

  is_token_expired <- function(asset) {

    get_token_expiry <- function(asset) {
      ms_info <- get_ms_info(asset)
      account <- get_ms_account(ms_info)
      container <- get_ms_container(ms_info)
      token[[account]][[container]][["msft:expiry"]]
    }

    difftime_token <- difftime(
      time1 = get_token_expiry(asset),
      time2 = as.POSIXlt(Sys.time(), tz = "UTC"),
      units = "secs"
    )

    difftime_token < ms_max_timeleft
  }

  get_token_value <- function(asset) {
    ms_info <- get_ms_info(asset)
    account <- get_ms_account(ms_info)
    container <- get_ms_container(ms_info)
    token[[account]][[container]][["token_value"]]
  }

  sign_asset <- function(asset) {
    parsed_url <- httr::parse_url(asset[["href"]])
    if (is_public_assets(parsed_url)) {
      return(asset)
    }
    if (!exists_token(asset) || is_token_expired(asset)) {
      new_token(asset)
    }
    token_value <- get_token_value(asset)

    # if the href is already sign it will not be modified
    parsed_url$query <- rstac:::.modify_list(parsed_url[["query"]], token_value)

    asset[["href"]] <- httr::build_url(parsed_url)
    asset
  }

  sign_item <- function(item) {
    item[["assets"]] <- lapply(item[["assets"]], sign_asset)

    return(item)
  }

  return(sign_item)
}