Descubrir productos y coberturas

El primer endpoint que vamos a utilizar es product, el cual acepta una solicitud GET() y devuelve la lista de productos³ disponibles.

get_products = GET(paste0(app_url, "product"))

status_code(get_products)

## [1] 200

Toda respuesta de un método HTTP trae un código indicando el tipo de resultado. La convención es que los códigos 2XX son 👍 y los 4XX son 👎. Arriba, el código 200 indica una operación exitosa. El resto de la respuesta representa la información que solicitamos; como content() transforma esa información en list, he definido una función para transformarla en data.frame:

content2DF = function(content) {
   DF = matrix(unlist(content), length(content), byrow = TRUE)
   colnames(DF) = names(content[[1]])
   as.data.frame(DF)
}

products = content(get_products) |> content2DF() |>
   select(ProductAndVersion, Platform:TemporalGranularity)

En este contexto, por “plataforma” generalmente nos referimos a un satélite artificial que orbita la Tierra y porta uno o más sensores; un “producto” es el conjunto de datos capturados por alguno de esos sensores y procesados en un centro especializado. Actualmente AppEEARS tiene 124 productos correspondientes a 14 plataformas:

unique(products$Platform)

##  [1] "GPW"                  "Combined MODIS"       "Terra MODIS"
##  [4] "Aqua MODIS"           "NASADEM"              "SMAP"
##  [7] "SRTM"                 "ASTER GDEM"           "S-NPP NASA VIIRS"
## [10] "Landsat ARD"          "DAYMET"               "SSEBop ET"
## [13] "eMODIS Smoothed NDVI" "ECOSTRESS"

select(filter(products, Platform == "SRTM"), -Platform)

##   ProductAndVersion     Description RasterType Resolution TemporalGranularity
## 1    SRTMGL1_NC.003 Elevation (DEM)       Tile        30m              Static
## 2 SRTMGL1_NUMNC.003    Source (DEM)       Tile        30m              Static
## 3    SRTMGL3_NC.003 Elevation (DEM)       Tile        90m              Static
## 4 SRTMGL3_NUMNC.003    Source (DEM)       Tile        90m              Static

Arriba hemos filtrado los productos correspondientes a la plataforma SRTM. Se trata de una misión que colocó un radar en órbita, a bordo del transbordador espacial Endeavour. El objetivo fue crear un Modelo Digital de Elevación (DEM) casi global, aplicando la técnica de interferometría radar. Como productos de esta plataforma tenemos dos DEMs⁴, con resolución espacial de 1 y 3 segundos de arco (~30 y ~90 metros). Nosotros solicitaremos un extracto de datos de 90 metros, "SRTMGL3_NC.003", observando que siempre se debe incluir la versión del producto (los últimos tres números).

Figure 1: Parte del hardware de la misión SRTM, fotografiado desde la cabina del Endeavour; vía Wikimedia Commons

Cada producto se conforma de una o más coberturas o layers; muchas veces cada layer es una banda espectral del sensor en cuestión. Con el endpoint product también podemos obtener las coberturas de cierto producto; a continuación, solicitamos las coberturas del producto "SRTMGL3_NC.003" y seleccionamos las variables más importantes de la respuesta:

get_product_layers = GET(paste0(app_url, "product/SRTMGL3_NC.003"))

product_layers = lapply(content(get_product_layers), \(x) x[-5] ) |> content2DF()

select(product_layers, Layer, Description:IsQA, Units:YSize)

##         Layer Description FillValue  IsQA  Units ValidMax ValidMin XSize YSize
## 1 SRTMGL3_DEM   Elevation    -32768 FALSE Meters    32767   -32767  1201  1201

La única cobertura se llama "SRTMGL3_DEM" y precisamente es el DEM. Hasta este punto hemos identificado el producto y la cobertura que nos interesa.

Autenticarse, a.k.a. Log in

El endpoint product no requiere que el usuario se identifique, pero para solicitar datos raster esto sí será un requisito. La manera de autenticarse es utilizar un segundo endpoint, login, que recibe el nombre de usuario y la contraseña de Earthdata Login, y devuelve un token. En el siguiente código aprovechamos askpass() para que la experiencia se parezca a hacer log in en un sitio web.

if (!file.exists("appeears_files/user_token.RDS")) {

   user_secret = paste0(
      askpass::askpass("Username:"), ":", askpass::askpass("Password:")
   ) |> base64_enc()

   post_login = POST(
      paste0(app_url, "login"),
      config = add_headers(
         "Authorization" = paste("Basic", user_secret),
         "Content-Type" = "application/x-www-form-urlencoded"
      ),
      body = "grant_type=client_credentials"
   )
   if (floor(status_code(post_login) / 100) != 2) stop(content(post_login)$message)

   saveRDS(user_token <- content(post_login)$token, "appeears_files/user_token.RDS")

} else user_token = readRDS("appeears_files/user_token.RDS")

El token es una cadena de caracteres que autoriza nuestro uso de la API a partir de ahora. Como es válido por 48 horas, lo guardamos en user_token.RDS para poder reusarlo; pasado dicho tiempo, hay que eliminarlo y repetir esta operación.

Crear y enviar una tarea

Una tarea es una lista de parámetros que definen cuáles datos solicitaremos a la API. Para construir una tarea en R podemos crear un elemento de clase list con la siguiente estructura:

task = list(
   task_type = "area",
   task_name = "st_helens_dem",
   params = list(
      dates = data.frame(startDate = "09-04-2021", endDate = "09-05-2021"),
      layers = data.frame(product = "SRTMGL3_NC.003", layer = "SRTMGL3_DEM"),
      geo = fromJSON("st_helens_aoi.json.txt"),
      output = list(format = list(type = "geotiff"), projection = "geographic")
   )
)

El parámetro task_name es el único que puede ser arbitrario, mientras que task_type puede ser area o point; el resto de parámetros que utilizo son obligatorios para una tarea de tipo area:

• dates define el rango temporal para extraer los datos; en este caso —como la misión SRTM inició y finalizó en febrero de 2000— podemos definir cualquier intervalo, entre startDate y endDate, posterior a "02-22-2000"; ambos valores deben estar en formato "MM-DD-YYYY".

• layers se compone de dos vectores, product y layer, con la misma longitud; aquí insertaremos el producto y la cobertura que identificamos anteriormente.

• geo define la región espacial, en formato GeoJSON y EPSG:4326, para extraer los datos; aquí leeremos con fromJSON() el polígono que guardamos en el primer paso.

• output define el tipo de raster que se obtendrá; type puede ser geotiff o netcdf4, y projection permite reproyectar⁵ el raster; en este caso, geographic equivale al EPSG:4326.

Una vez creada la tarea, podemos enviarla a través de una solicitud POST() al endpoint task. Esa solicitud debe ser autorizada con el token que obtuvimos en el paso anterior, y debe ser codificada como JSON (utilizamos toJSON() para transformar la tarea).

post_task = POST(
   paste0(app_url, "task"),
   config = add_headers(
      "Authorization" = paste("Bearer", user_token),
      "Content-Type" = "application/json"
   ),
   body = toJSON(task, auto_unbox = TRUE, digits = NA)
)
if (floor(status_code(post_task) / 100) != 2) stop(content(post_task)$message)

¿Cómo le va a mi tarea?

Este paso es opcional. Dependiendo del tamaño de la región espacial, del rango temporal y número de coberturas solicitadas, el procesamiento de la tarea puede demorar segundos o algunos minutos. El siguiente código consulta el estado del procesamiento cada 30 segundos, usando el id de la tarea. Para tareas más complejas, puede aumentar el intervalo de consulta modificando el Sys.sleep().

task_id = content(post_task)$task_id

task_status = "not done"

while (task_status != "done") {
   
   Sys.sleep(30)
   
   get_task_status = GET(
      paste0(app_url, "task/", task_id),
      config = add_headers("Authorization" = paste("Bearer", user_token))
   )
   task_status = content(get_task_status)$status
   
   message(as.character(Sys.time(), "[%T] "), task$task_name, " status: ", task_status)

   if (task_status == "error") stop(content(get_task_status)$error)
}

## [16:20:36] st_helens_dem status: processing

## [16:21:07] st_helens_dem status: done

¡Está lista! Descargar el contenido

Cuando la tarea tenga status done, los datos solicitados estarán listos para ser descargados. Para cada tarea exitosa, AppEEARS crea un paquete o bundle de archivos que incluye, además del raster que necesitamos, información de calidad y metadata. El endpoint bundle —añadiendo el id de la tarea— recibe una solicitud GET() y devuelve el contenido del bundle.

get_bundle = GET(paste0(app_url, "bundle/", task_id))

bundle = content2DF(content(get_bundle)$files) |>
   mutate(file_name = sub("^.+/(.+\\.\\d{3})_(.+)_doy", "\\1-\\2/\\2-doy", file_name))

select(bundle, file_name, file_size)

##                                                          file_name file_size
## 1    SRTMGL3_NC.003-SRTMGL3_DEM/SRTMGL3_DEM-doy2000042_aid0001.tif    396174
## 2 SRTMGL3_NUMNC.003-SRTMGL3_NUM/SRTMGL3_NUM-doy2000042_aid0001.tif     49808
## 3                         SRTMGL3-NUMNC-003-SRTMGL3-NUM-lookup.csv      1974
## 4                  SRTMGL3-NUMNC-003-SRTMGL3-NUM-Statistics-QA.csv       536
## 5                                    SRTMGL3-NC-003-Statistics.csv       340
## 6                                   st-helens-dem-granule-list.txt       776
## 7                                       st-helens-dem-request.json      1059
## 8                        st-helens-dem-SRTMGL3-NC-003-metadata.xml     22031
## 9                                                        README.md     23442

Se modificó file_name para descargar los archivos en directorios, según producto y cobertura. También notarán que este endpoint no requiere autenticación; la consecuencia es que cualquiera que tenga el id de nuestra tarea podrá descargar el contenido asociado. Eso es precisamente lo que hace el siguiente código, pero solo descargaremos los cinco primeros archivos del bundle:

dir.create(bundle_path <- file.path("appeears_files", task$task_name))

for (file_id in bundle$file_id[1:5]) {
   
   file_name = bundle$file_name[bundle$file_id == file_id]
   
   message(as.character(Sys.time(), "[%T] download: "), file_name)

   if (!file.path(bundle_path, dirname(file_name)) |> dir.exists()) {
      file.path(bundle_path, dirname(file_name)) |> dir.create()
   }
   get_bundle_file = GET(
      paste0(app_url, "bundle/", task_id, "/", file_id),
      # progress(), # opcional
      write_disk(file.path(bundle_path, file_name))
   )
}

## [16:21:07] download: SRTMGL3_NC.003-SRTMGL3_DEM/SRTMGL3_DEM-doy2000042_aid0001.tif

## [16:21:11] download: SRTMGL3_NUMNC.003-SRTMGL3_NUM/SRTMGL3_NUM-doy2000042_aid0001.tif

## [16:21:11] download: SRTMGL3-NUMNC-003-SRTMGL3-NUM-lookup.csv

## [16:21:11] download: SRTMGL3-NUMNC-003-SRTMGL3-NUM-Statistics-QA.csv

## [16:21:12] download: SRTMGL3-NC-003-Statistics.csv