Roadmap for DMR++/Virtual Dataset Support in earthaccess

[PaulH97]

Hi everyone,

I'm currently integrating earthaccess into our tool Terragon and would love to use the open_virtual_mfdataset() function for streamlined cloud data access. However, for several collections (such as HLS or Sentinel-1), I encounter the well-known error that a required DMR++ file is missing.

Is there a general plan or roadmap for expanding DMR++/virtual dataset support to more collections?
Should we expect wider compatibility in the near future, or is this functionality likely to remain experimental and limited to a subset of datasets?

Knowing this will help us decide whether to rely on these experimental features, or to develop our own workflows using lower-level tools like rasterio.

Thanks for your efforts and for any guidance you can provide!

Best regards,
Paul

[betolink]

This is a great question and here is a partial answer: yes, there are plans to extend dmrpp support to more collections and standardize the way they are generated and provided. This falls into OPeNDAP and NASA Earth roadmaps so I don't have a timeline from the earthaccess side. Perhaps @jgallagher59701 @Mikejmnez knows more about the overall plan.

Aside from this, I'd think that for Cloud-Optimized Geo-TIF (HLS) there is no need for a chunk mapping file (dmrpp), since the rasterio/gdal drivers can take advantage of the internal metadata and get to the bytes in an efficient -cloud optimized- way.
There are still questions on how to group many tiles into a virtual data cube, I know there are projects like odc-stac and stackstac that do something like this by reading and aligning them into a cube https://discourse.pangeo.io/t/comparing-odc-stac-load-and-stackstac-for-raster-composite-workflow/4097 This could be the workflow you guys may be after if your sources are only rasters.

If your workflows include data from HDF5 or NetCDF files then the dmrpp path becomes way more attractive than generating the reference ourselves (kerchunk). We're working on improving these access patterns even further with the support of icechunk (thanks @jbusecke!).

[PaulH97]

Thank you very much for your quick response and feedback. Would you have a minimal code example of how we can use odc-stac or stackstac with earthaccess? I am currently trying to implement this without earthaccess by getting the returned STAC items from the HLSS30_2.0 collection at https://cmr.earthdata.nasa.gov/cloudstac/LPCLOUD into a cube with odc.stac.load, but unfortunately it is not working.

[betolink]

Not working because of authentication? auth is kind of an issue with lower level libraries. Generally speaking if you only need one VM to process your data we can make it work with an .netrc file with EDL credentials. Here is a complete example with odc-stac that uses Dask to parallelize work (mostly taken from @ebolch tutorial on HLS data: https://github.com/nasa/HLS-Data-Resources/blob/main/python/tutorials/EVI_timeseries_with_odc_stac.ipynb

from pystac_client import Client
import os
import odc.stac
import xarray as xr
import dask
import dask.distributed


client = dask.distributed.Client()
# display(client)

catalog = Client.open("https://cmr.earthdata.nasa.gov/stac/LPCLOUD/")
# Configure odc.stac rio env - requires a .netrc file, sends info to dask client
odc.stac.configure_rio(cloud_defaults=True,
                       verbose=True,
                       client=client,
                       GDAL_DISABLE_READDIR_ON_OPEN='EMPTY_DIR',
                       GDAL_HTTP_COOKIEFILE=os.path.expanduser('~/cookies.txt'),
                       GDAL_HTTP_COOKIEJAR=os.path.expanduser('~/cookies.txt'))

bbox = [-105.5, 39.4, -105.1, 39.7]  # [minx, miny, maxx, maxy]

search = catalog.search(
    collections=["HLSS30_2.0","HLSL30_2.0"],
    bbox=bbox,
    datetime="2021-05-01/2021-09-30",
    limit=100
)
items = list(search.items())

print(f"Found {len(items)} HLS granules")

for item in items:
    if "HLS.L30" in item.id:
        item.assets["NIR"] = item.assets.pop("B05")
    if "HLS.S30" in item.id:
        item.assets["NIR"] = item.assets.pop("B8A")

cube = odc.stac.load(
    items,
    crs="utm",
    resolution=30,
    resampling="bilinear",
    bands=("B02", "B04","NIR", "Fmask"),
    chunks={"band":1,"x":512,"y":512},
    groupby="solar_day"
)

print(cube)
# if you want to print a band
# print(cube.B02.isel(time=0).plot.imshow())

I've been prototyping a workflow that does the same without odc-stac (rasterio, rioxarray) https://gist.github.com/betolink/0c3cb07a2c0717c4dfb9840e55e08a68 and will modify it to use NASA tokens. that way we won't need to have a netrc file on each VM if we decide to scale out. I hope this helps.