modis-LAI.R

library(tidyverse)
library(rstac)
library(httr)
library(stars)
library(spData)
box <- st_bbox(us_states)

# hack around rstac bug in signing
mysign <- function(href) {
  x <- parse_url(href)
  y <- glue::glue("{scheme}://{hostname}/{path}", scheme = x$scheme, hostname = x$hostname, path = x$path)
  resp <- GET(paste0("https://planetarycomputer.microsoft.com/api/sas/v1/sign?href=", y))
  stop_for_status(resp)
  url <- httr::content(resp)
  out <- paste0("/vsicurl/",url$href)
  return(out)
}

# search any time in the past 8 days anywhere in continental US:
end <- Sys.Date()
start <- end - 8

matches <-
  stac("https://planetarycomputer.microsoft.com/api/stac/v1") |>
  stac_search(collections = "modis-15A2H-061",
              datetime = paste(start, end, sep="/"),
              bbox = c(box),
              limit = 100) |>
  post_request() |>
  items_sign(sign_fn = sign_planetary_computer())

v <- purrr::map_int(matches$features, list("properties", "modis:vertical-tile"))
h <- purrr::map_int(matches$features, list("properties", "modis:horizontal-tile"))
date <- purrr::map_chr(matches$features, list("properties", "created"))

urls <- map_chr(matches$features, list("assets", "Lai_500m", "href"))
usa <- tibble(v,h, date, urls) |> group_by(v,h) |> slice_max(date)


rast <- vector("list", nrow(usa))
for(i in 1:nrow(usa)) {
  rast[[i]] <- stars::read_stars(mysign(usa$urls[i]))
  Sys.sleep(10) # avoid rate limiting
}

# Combine all 17 tiles into a mosaic of the continental USA
mosaic <- do.call(st_mosaic, compact(rast))

## Optionally, re-project to a more familiar coordinate system
# mosaic <- mosaic |> st_transform(crs="+proj=longlat")

# keep a local copy if desired (e.g. later analysis etc)
stars::write_stars(mosaic, "mosiac.tif")

# quick peek at data
mosaic |> plot()



## tmap-based plot

## Log scale (proof of principle we can do basic math on every pixel)
remap <- function(x, epsilon=1e-1) {
  log(x[1] + epsilon)
}

mosaic |>
    st_apply(1:2, remap) |>
    tm_shape() +
    tm_raster(n = 100,
              palette = viridisLite::mako(100),
              legend.show = FALSE)

modis-w-gdalcubes.R

## Let's get 1-month averages of Leaf Area Index for the whole US from the
## MODIS 500m LAI data:


library(rstac)
library(gdalcubes)
library(spData)
library(stars)
gdalcubes_options(parallel = 24)
box <- spData::us_states |> sf::st_bbox()
start <- "2022-01-01"
end <- "2022-12-01"

# Query STAC
matches <-
  stac("https://planetarycomputer.microsoft.com/api/stac/v1") |>
  stac_search(collections = "modis-15A2H-061",
              datetime = paste(start, end, sep="/"),
              bbox = c(box)) |>
  get_request() |>
  items_fetch() |>
  items_sign(sign_fn = sign_planetary_computer())

# Define datacube
cube <- gdalcubes::stac_image_collection(matches$features,
                                         asset_names = "Lai_500m",
                                         duration = "start")
# Define view
v <- cube_view(srs = "EPSG:4326",
               extent = list(t0 = start, t1 = end,
                             left = box[1], right = box[3],
                             top = box[4], bottom = box[2]),
               dx = 0.1, dy = 0.1, dt= "P1M",
               aggregation = "median", resampling = "near"
)

# tada, abstract, lazy eval datacube
Q <- raster_cube(cube,v)

# animate over time, automatically runs only at resolution needed
Q |>
  animate(col = viridisLite::mako, zlim=c(0,60),
          key.pos = 1, save_as = "anim.gif", fps = 4)



# We can do most any of simple and rich spatial calculations
# Here's a simple average-over-time, but can apply arbitrary functions by pixel too
Q |> reduce_time(c("mean(Lai_500m)")) |>
  plot(zlim=c(0,50), col = viridisLite::mako)

This could definitely be cleaned up. Note that the mysign is merely my hack around a current bug in rstac, which should go away soon. That may also resolve the rate-limiting issue that currently requires the Sys.sleep().

Looping over one image per month for last year shows the seasons: