Loading metadata for Landsat images loaded with `dc.load` using Open Data Cube

Question

How can I pull in metadata from Open Data Cube for multiple Landsat observations imported using dc.load?

Answer 1

I've previously used loaded metadata for each observation in a xarray loaded using Open Data Cube like this:

import numpy as np
import datacube
from datacube.api.query import query_group_by
from datacube.model.utils import xr_apply

# Connect to datacube
dc = datacube.Datacube(app='Metadata')

# Spatiotemporal query
query = dict(product='ga_ls8c_ard_3', 
             lat=(-35.0, -35.1), 
             lon=(148.0, 148.1), 
             time='2018', 
             group_by='solar_day'
            )

# Load the metadata from each dataset into an xarray after grouping by solar day
gb = query_group_by(**query)
datasets = dc.find_datasets(**query)
dataset_array = dc.group_datasets(datasets, gb)
cloud_cover = xr_apply(dataset_array, lambda t, dd: np.mean([d.metadata.cloud_cover for d in dd]), dtype=float)

# Load the satellite data data
ds = dc.load(datasets=datasets, 
             measurements=['nbart_red'],
             output_crs='EPSG:3577', 
             resolution=(-30, 30), 
             dask_chunks={},
             group_by='solar_day'
            )

# Add metadata to data
ds['cloud_cover'] = cloud_cover
ds

The np.mean([d.metadata.cloud_cover for d in dd]) bit above is because there might be multiple datasets for a single solar day, each with their own cloud cover percent, and we only want a single number out (in this case, mean cloud cover across all those datasets).

(Just keep in mind that "cloud_cover" in the example above is metadata that is calculated across the entire satellite scene, and doesn't necessarily match the clouds in the area you've requested. The load_ard function from dea-tools uses a much more accurate approach by calculating the number of actually cloudy pixels in the area you request).