Loading data using virtual products

Sign up to the DEA Sandbox to run this notebook interactively from a browser
Compatibility: Notebook currently compatible with both the NCI and DEA Sandbox environments
Products used: ga_s2am_ard_3, ga_s2bm_ard_3

Background

In addition to the Datacube.load function, Open Data Cube also provide a mechanism to load data from virtual products. These products are created by combining existing products in different ways. They are designed to make it easier to assemble data from different bands and different sensors together, but they can be useful even when dealing with a single datacube product.

Description

This notebook introduces virtual products and provides examples of how to achieve simple data transformations with them. This includes:

Loading a simple product
Applying a cloud mask
Applying a dilated cloud mask
Calculating band math and indices
Combining multiple sensors
Using virtual products with Dask
Calculating statistics and custom transforms
Masking nodata
Loading data at native resolution
Reprojecting on the fly

The reference documentation on virtual products can be found here.

Getting started

First we import the required Python packages, then we connect to the database, and load the catalog of virtual products.

[1]:

%matplotlib inline

import sys
import datacube
from datacube.virtual import catalog_from_file

import sys
sys.path.insert(1, '../Tools/')
from dea_tools.plotting import rgb

Connect to the datacube

[2]:

dc = datacube.Datacube(app='Virtual_products')

Load virtual products catalogue

For the examples below, we will load some example data using some pre-defined virtual products from a YAML catalog file, virtual-product-catalog.yaml in the Supplementary_data/Virtual_products directory (view on GitHub). First we load the catalogue of virtual products:

[3]:

catalog = catalog_from_file('../Supplementary_data/Virtual_products/virtual-product-catalog.yaml')

We can list the products and transforms that the catalog defines:

[4]:

list(catalog)

[4]:

['ga_s2am_ard_3_simple',
 's2a_nbart_rgb',
 's2a_nbart_rgb_renamed',
 's2a_nbart_rgb_nodata_masked',
 's2a_nbart_rgb_and_fmask',
 's2a_nbart_rgb_and_fmask_native',
 's2a_nbart_rgb_and_fmask_albers',
 's2a_fmask',
 's2a_cloud_free_nbart_rgb_albers',
 's2a_cloud_free_nbart_rgb_albers_reprojected',
 's2a_cloud_free_nbart_rgb_albers_dilated',
 's2a_NDVI',
 's2a_optical_albers',
 's2a_optical_and_fmask_albers',
 's2a_tasseled_cap',
 's2a_cloud_free_tasseled_cap',
 's2a_mean',
 's2_nbart_rgb',
 'nbart_ndvi_transform',
 'nbart_tasseled_cap_transform']

We will go through them one by one. We also specify our example area of interest throughout this notebook.

[5]:

query = {'x': (153.38, 153.47),
         'y': (-28.83, -28.92),
         'time': ('2018-01', '2018-02')}

Combining sensors

To combine multiple sensors into one virtual product, we can use the collate combinator.

[27]:

product = catalog.s2_nbart_rgb
product

[27]:

collate:
- measurements:
  - nbart_blue
  - nbart_green
  - nbart_red
  product: ga_s2am_ard_3
- measurements:
  - nbart_blue
  - nbart_green
  - nbart_red
  product: ga_s2bm_ard_3

[28]:

data = product.load(dc, **query)
rgb(data, col='time', col_wrap=4)

../../../_images/notebooks_How_to_guides_Virtual_products_56_0.png

You can collate any number of virtual products as long as they share the same set of measurements.

More advanced usage

Dask

Dask is a useful tool when working with large analyses (either in space or time) as it breaks data into manageable chunks that can be easily stored in memory. It can also use multiple computing cores to speed up computations.

To load virtual products using Dask, pass the dask_chunks specification to the load method.

Note: For more information about using Dask, refer to the Parallel processing with Dask notebook.

[29]:

data = catalog.s2a_NDVI.load(dc, dask_chunks={'time': 1, 'x': 2000, 'y': 2000}, **query)
data

[29]:

<xarray.Dataset>
Dimensions:      (time: 6, y: 1001, x: 882)
Coordinates:
  * time         (time) datetime64[ns] 2018-01-06T23:52:37.461000 ... 2018-02...
  * y            (y) float64 6.811e+06 6.811e+06 ... 6.801e+06 6.801e+06
  * x            (x) float64 5.37e+05 5.371e+05 ... 5.458e+05 5.459e+05
    spatial_ref  int32 32756
Data variables:
    NDVI         (time, y, x) float32 dask.array<chunksize=(1, 1001, 882), meta=np.ndarray>
Attributes:
    crs:           epsg:32756
    grid_mapping:  spatial_ref

[30]:

data.compute().NDVI.plot.imshow(col='time', col_wrap=4)

[30]:

<xarray.plot.facetgrid.FacetGrid at 0x7f22f7e19730>

../../../_images/notebooks_How_to_guides_Virtual_products_62_1.png

Statistics

Virtual products support some statistics out-of-the-box (those supported by xarray methods). Here we calculate the mean image of our area of interest. Grouping by year here groups all the data together since our specified time range falls within a year.

[31]:

product = catalog.s2a_mean
product

[31]:

aggregate: datacube.virtual.transformations.XarrayReduction
group_by: datacube.virtual.transformations.year
input:
  input:
    flags:
      fmask: valid
    input:
      measurements:
      - nbart_blue
      - nbart_green
      - nbart_red
      - fmask
      output_crs: EPSG:3577
      product: ga_s2am_ard_3
      resampling:
        '*': average
        fmask: nearest
      resolution:
      - -20
      - 20
    mask_measurement_name: fmask
    transform: datacube.virtual.transformations.MakeMask
  mask_measurement_name: fmask
  preserve_dtype: false
  transform: datacube.virtual.transformations.ApplyMask
method: mean

[32]:

data = product.load(dc, **query)
rgb(data, col='time', col_wrap=4)

../../../_images/notebooks_How_to_guides_Virtual_products_66_0.png

User-defined transformations

We can define transformations ourselves too. For illustrative purposes, we define the NDVI transform here again, in Python. However, this example is going to be a little bit different since the catalog cannot possibly know about the transformation we are going to define here in the notebook.

Our transformation needs to inherit from datacube.virtual.Transformation and to implement the methods defined there. The construct method is a lower-level way of constructing a virtual product. We will need it since we are going to refer to a Python class directly.

[33]:

from datacube.virtual import Transformation
from datacube.virtual import Measurement
from datacube.virtual import construct


class NDVI(Transformation):

    def compute(self, data):
        return (((data.nbart_nir_1 - data.nbart_red) /
                 (data.nbart_nir_1 + data.nbart_red))
                .astype('float32')
                .to_dataset(name='NDVI'))

    def measurements(self, input_measurements):
        return {'NDVI': Measurement(name='NDVI', dtype='float32', units='1')}


product = construct(transform=NDVI,
                    input=construct(product='ga_s2am_ard_3',
                                    measurements=['nbart_red', 'nbart_nir_1']))

[34]:

data = product.load(dc, **query)
data.NDVI.plot.imshow(col='time', col_wrap=4)

[34]:

<xarray.plot.facetgrid.FacetGrid at 0x7f22f45333d0>

../../../_images/notebooks_How_to_guides_Virtual_products_70_1.png

Note that this version does not mask nodata like the expressions version from before automatically does.

Reprojecting on-the-fly

We have so far reprojected our images to Albers first, and then applied our cloud mask to it. For (hopefully) higher spectral fidelity, we may want to apply the cloud mask in the native projection, and subsequently reproject to a common grid to obtain a dataset ready for temporal analysis. The reproject combinator in virtual products allow us to do this.

[41]:

product = catalog.s2a_cloud_free_nbart_rgb_albers_reprojected
product

[41]:

input:
  input:
    flags:
      fmask: valid
    input:
      like: fmask
      measurements:
      - nbart_blue
      - nbart_green
      - nbart_red
      - fmask
      product: ga_s2am_ard_3
      resampling: average
    mask_measurement_name: fmask
    transform: datacube.virtual.transformations.MakeMask
  mask_measurement_name: fmask
  preserve_dtype: false
  transform: datacube.virtual.transformations.ApplyMask
reproject:
  output_crs: EPSG:3577
  resolution:
  - -20
  - 20
resampling: average

[42]:

data = product.load(dc, **query)
rgb(data, col='time', col_wrap=4)

../../../_images/notebooks_How_to_guides_Virtual_products_88_0.png

Transforms defined in the catalog

We can also store unapplied transforms in the catalog:

[43]:

catalog.nbart_tasseled_cap_transform

[43]:

output:
  brightness:
    dtype: float32
    formula: 0.2043 * nbart_blue + 0.4158 * nbart_green + 0.5524 * nbart_red + 0.5741
      * nbart_nir_1 + 0.3124 * nbart_swir_2 + -0.2303 * nbart_swir_3
  greenness:
    dtype: float32
    formula: -0.1603 * nbart_blue + -0.2819 * nbart_green + -0.4934 * nbart_red +
      0.7940 * nbart_nir_1 + -0.0002 * nbart_swir_2 + -0.1446 * nbart_swir_3
  wetness:
    dtype: float32
    formula: 0.0315 * nbart_blue + 0.2021 * nbart_green + 0.3102 * nbart_red + 0.1594
      * nbart_nir_1 + -0.6806 * nbart_swir_2 + -0.6109 * nbart_swir_3
transform: expressions

The syntax for these transforms is the same as those of the products, except the input section is ommitted (compare the recipe of the product s2a_tasseled_cap here), so that we can provide the input that we want:

[44]:

product = catalog.nbart_tasseled_cap_transform(catalog.s2a_optical_albers)
product

[44]:

input:
  measurements:
  - nbart_blue
  - nbart_green
  - nbart_red
  - nbart_nir_1
  - nbart_swir_2
  - nbart_swir_3
  output_crs: EPSG:3577
  product: ga_s2am_ard_3
  resampling: average
  resolution:
  - -20
  - 20
output:
  brightness:
    dtype: float32
    formula: 0.2043 * nbart_blue + 0.4158 * nbart_green + 0.5524 * nbart_red + 0.5741
      * nbart_nir_1 + 0.3124 * nbart_swir_2 + -0.2303 * nbart_swir_3
  greenness:
    dtype: float32
    formula: -0.1603 * nbart_blue + -0.2819 * nbart_green + -0.4934 * nbart_red +
      0.7940 * nbart_nir_1 + -0.0002 * nbart_swir_2 + -0.1446 * nbart_swir_3
  wetness:
    dtype: float32
    formula: 0.0315 * nbart_blue + 0.2021 * nbart_green + 0.3102 * nbart_red + 0.1594
      * nbart_nir_1 + -0.6806 * nbart_swir_2 + -0.6109 * nbart_swir_3
transform: datacube.virtual.transformations.Expressions

This resultant product is equivalent to the s2a_tasseled_cap product, but what we have gained here is that the nbart_tasseled_cap_transform is more re-usable and can be applied to other products as long as they provide the nbart_* bands.

Additional information

License: The code in this notebook is licensed under the Apache License, Version 2.0. Digital Earth Australia data is licensed under the Creative Commons by Attribution 4.0 license.

Contact: If you need assistance, please post a question on the Open Data Cube Discord chat or on the GIS Stack Exchange using the open-data-cube tag (you can view previously asked questions here). If you would like to report an issue with this notebook, you can file one on GitHub.

Last modified: September 2021

Compatible datacube version:

[45]:

print(datacube.__version__)

1.8.6

Loading data using virtual products

Background

Description

Getting started

Connect to the datacube

Load virtual products catalogue

Single datacube product

Data transformations

Cloud-free NBART

Dilation of masks

Band math

Combining sensors

More advanced usage

Dask

Statistics

User-defined transformations

Bonus features

Renaming bands

Masking invalid data

Native load

Reprojecting on-the-fly

Transforms defined in the catalog

Additional information

Tags