Use harmony to subset TEMPO data from R

harmony_tempo.R

# demo of downloading a single spatially subset granule from TEMPO
# using harmony API

# Need dev version of earthdatalogin for edl_search() to work,
# install with:
# remotes::install_github("boettiger-lab/earthdatalogin")
library(earthdatalogin)
library(httr2)

# space-time bounds (API formatted strings)
# https://harmony.earthdata.nasa.gov/docs
bbox <- c(-100.4, 18.1, -97.4, 20.6)
lat_string <- paste0("lat(", bbox[2], ":", bbox[4], ")")
lon_string <- paste0("lon(", bbox[1], ":", bbox[3], ")")
temporal_string <- "2026-03-09T17:00:00Z,2026-03-09T19:00:00Z"
# example shows Mexico City region 2026-03-09
# https://projects.cosmicds.cfa.harvard.edu/tempo-lite/?lat=19.2541&lon=-99.2473&zoom=8&t=1741525740000&extendedRange=true

# 1. Auth: Get your token
# we use earthdatalogin with NASA EarthData credentials already stored in our environment
# https://boettiger-lab.github.io/earthdatalogin/
edl_set_token(prompt_for_netrc = FALSE)
header_file <- Sys.getenv("GDAL_HTTP_HEADER_FILE")
my_token <- gsub(
  "Authorization: Bearer ",
  "",
  readLines(header_file, warn = FALSE)
)

# check your token by seeing if you can access a file listed in the earthdatalogin README
url <- "https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/HLSL30.020/HLS.L30.T56JKT.2023246T235950.v2.0/HLS.L30.T56JKT.2023246T235950.v2.0.SAA.tif"
terra::rast(url, vsi = TRUE)
# if earthdatalogin isn't set up, instead of giving SpatRaster info this will return
# "Error: [rast] file does not exist"...

# 2. STEP ONE: Find the Granule ID using CMR (NASA's Search Engine)
# This confirms the data exists before we even talk to Harmony

# TEMPO_NO2_L3_V04
shortname <- "TEMPO_NO2_L3"
concept_id <- "C3685896708-LARC_CLOUD"
# https://doi.org/10.5067/IS-40e/TEMPO/NO2_L3.004

# TEMPO_O3PROF_L3_V04
# concept_id <- "C3685896402-LARC_CLOUD"
# https://doi.org/10.5067/IS-40e/TEMPO/O3PROF_L3.004

results <- edl_search(
  short_name = shortname,
  temporal = temporal_string,
  bounding_box = bbox,
  parse_results = FALSE # return as a list
)

length(results)
# Returns three time steps, but we just want the first one for this demo

concept_id <- results[[1]]$collection_concept_id
granule_id <- results[[1]]$title

# 3. STEP TWO: Send that specific Granule to Harmony for subsetting
# Note the URL structure: we use the collection ID + /ogc-api-coverages/...
harmony_url <- paste0(
  "https://harmony.earthdata.nasa.gov/",
  concept_id,
  "/ogc-api-coverages/1.0.0/collections/all/coverage/rangeset"
)

req <- request(harmony_url) %>%
  req_url_query(
    granuleName = granule_id,
    subset = lat_string,
    subset = lon_string,
    format = "application/x-netcdf4"
  ) %>%
  req_auth_bearer_token(my_token) %>%
  req_headers(Accept = "application/json")

# 1. Define your output path
out_file <- "Mexico_City_TEMPO_V04.nc"

# 2. Run the request, but specify the 'path'
# in this case - when there is a single expected output the API returns it directly
# as a netcdf -- so this wouldn't work for a large set of resulting files
resp <- req %>%
  req_perform(path = out_file)

cat("Success! File saved to:", out_file, "\n")

# 3. Verify with terra
# note issue with finding CRS so we use md=T
# https://github.com/rspatial/terra/issues/1795#issuecomment-2918037100
library(terra)
if (file.exists(out_file)) {
  r <- rast(out_file, subds = "/product/vertical_column_troposphere", md = T)
  print(r)
  plot(r, main = "TEMPO V04 - Central Mexico")
}

# end of script