Last active
June 8, 2023 00:44
-
-
Save cboettig/147547c70084c44b7544648a84712981 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) | |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
## 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) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
This could definitely be cleaned up. Note that the
mysign
is merely my hack around a current bug inrstac
, which should go away soon. That may also resolve the rate-limiting issue that currently requires theSys.sleep()
.