In mikejohnson51/opendap.catalog: R OpenDap Client

knitr::opts_chunk$set(
  collapse = TRUE,
  warning = FALSE,
  message = FALSE,
  comment = "#>"
)

library(opendap.catalog)
library(terra)
library(dplyr)

One challenge that opendap catalogs have for end users is that they are often dissagregated by space and/or time. While this structure may (or may not) make sense from the data storage/service side, it is a burden on people who want the data.

Here we look at two data sets MODIS and MACA. The first is dissagregated across space (XY) while the second across time (T).

XY Aggregates

Within MODIS there are 460 non-fill tiles with an approximate 10 degree by 10 degree size at the equator. The public (and keyless) MODIS OpenDAP server stores each tile - aggregated through time - as a single resource.

mod <- grids[grids$grid.id == "XY_modis", ]

bboxs <- vect(sapply(1:nrow(mod), function(x) {
  make_vect(mod[x, ])
}))

AOI_vect <- project(vect(AOI::aoi_get(state = "conus", union = TRUE)), crs(bboxs))
bbox_int = terra::relate(bboxs, AOI_vect, "intersects")

{
  plot(bboxs[bbox_int[, 1], ], main = "15 tiles over CONUS" )
  plot(AOI_vect, add = TRUE)
}

AOI_vect_fl <- project(vect(AOI::aoi_get(state = "FL", union = TRUE)), crs(bboxs))
bbox_int <- terra::relate(bboxs, AOI_vect_fl, "intersects")

{
  plot(bboxs[bbox_int[, 1], ], main = "3 tiles over Florida")
  plot(AOI_vect_fl, add = TRUE)
}

In cases where the AOI (like the state of Florida) crosses mutiple tiles, the multiple resources must be identified, subset and then stitched together! This XY aggregations is one of the perks avaialble with dap().

Example

Lets find a dataset of interest:

(modis_ex <- search("MOD16A2.006 PET"))

And then query that dataset for a spatial and temporal slice:

system.time({
  dap <- dap(
    catolog = modis_ex,
    AOI = AOI::aoi_get(state = "FL"),
    startDate = "2020-01-01",
    endDate = "2020-01-31"
  )
})

Note that the returned object is a single SpatRaster layer for each time period, but in the summary we see this is the result of compositing 3 unique tiles (e.g. resources).

plot(dap)

T Aggregates

Some data sets also tile by time period. Often this occurs when there is a historic period and multiple periods of future forecasts using different climate scenarios.

An example of this is the MACA dataset. Say we want daily specific humidity and rainfall from the MACA down scaling of the BNU-ESM model.

First we need to find the catalog elements:

tmp <- rbind(search("maca daily huss bnu-esm"),
             search("maca daily pr bnu-esm"))

(select(tmp, id, variable, model, scenario, duration))

Here we see that the data set is split at "2006-01-01" into a historic, and multi scenario future. For out example lets chose the historic and rcp85 future:

(maca_ex <- tmp[tmp$scenario %in% c("historical", "rcp85"), ])

Example

Once defined

system.time({
 (dap <- dap(
    catolog = maca_ex,
    AOI = AOI::aoi_get(state = "NC"),
    startDate = "2005-12-25", 
    endDate   = "2006-01-05"
  ))
})

dap