Loading data using virtual products

Background

In addition to the Datacube.load function, recent versions of 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:

  1. Loading a simple product

  2. Applying a cloud mask

  3. Applying a dilated cloud mask

  4. Calculating band math and indices

  5. Combining multiple sensors

  6. Using virtual products with Dask

  7. Calculating statistics and custom transforms

  8. Masking nodata

  9. Loading data at native resolution

  10. 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

sys.path.append('../Scripts')
from dea_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 or in Jupyter). 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 that the catalog defines:

[4]:
list(catalog)
[4]:
['s2a_ard_granule_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_tasseled_cap',
 's2a_cloud_free_tasseled_cap',
 's2a_mean',
 's2_nbart_rgb']

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')}

Single datacube product

A virtual product is defined by what we call a “recipe” for it. If we get a virtual product from the catalog and try to print it, it displays the recipe.

[6]:
product = catalog['s2a_ard_granule_simple']
product
[6]:
product: s2a_ard_granule

This very simple recipe says that this virtual product, called s2a_ard_granule_simple, is really just the datacube product s2_ard_granule. We can use this to load some data.

[7]:
data = product.load(dc, measurements=['fmask'], **query)
data
[7]:
Show/Hide data repr Show/Hide attributes
xarray.Dataset
    • time: 6
    • x: 441
    • y: 501
    • time
      (time)
      datetime64[ns]
      2018-01-06T23:52:41.026000 ... 2018-02-25T23:52:41.026000
      units :
      seconds since 1970-01-01 00:00:00
      array(['2018-01-06T23:52:41.026000000', '2018-01-16T23:52:41.026000000',
             '2018-01-26T23:52:41.026000000', '2018-02-05T23:52:41.026000000',
             '2018-02-15T23:52:41.026000000', '2018-02-25T23:52:41.026000000'],
            dtype='datetime64[ns]')
    • y
      (y)
      float64
      6.811e+06 6.811e+06 ... 6.801e+06
      units :
      metre
      resolution :
      -20.0
      crs :
      PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
      array([6810790., 6810770., 6810750., ..., 6800830., 6800810., 6800790.])
    • x
      (x)
      float64
      5.370e+05 5.371e+05 ... 5.458e+05
      units :
      metre
      resolution :
      20.0
      crs :
      PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
      array([537050., 537070., 537090., ..., 545810., 545830., 545850.])
    • fmask
      (time, y, x)
      uint8
      1 1 1 1 1 1 1 1 ... 2 2 2 2 2 2 2 2
      units :
      1
      nodata :
      0
      flags_definition :
      {'fmask': {'bits': [0, 1, 2, 3, 4, 5, 6, 7], 'values': {'0': 'nodata', '1': 'valid', '2': 'cloud', '3': 'shadow', '4': 'snow', '5': 'water'}, 'description': 'Fmask'}}
      crs :
      PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
      array([[[1, 1, 1, ..., 1, 1, 1],
              [1, 1, 1, ..., 1, 1, 1],
              [1, 1, 1, ..., 1, 1, 1],
              ...,
              [1, 1, 1, ..., 2, 2, 2],
              [1, 1, 1, ..., 2, 2, 2],
              [1, 1, 1, ..., 2, 2, 2]],
      
             [[1, 1, 1, ..., 1, 1, 1],
              [1, 1, 1, ..., 1, 1, 1],
              [1, 1, 1, ..., 1, 1, 1],
              ...,
              [1, 1, 1, ..., 1, 1, 1],
              [1, 1, 1, ..., 1, 1, 1],
              [1, 1, 1, ..., 1, 1, 1]],
      
             [[1, 1, 1, ..., 2, 2, 2],
              [1, 1, 1, ..., 2, 2, 2],
              [1, 1, 1, ..., 2, 2, 2],
              ...,
              [1, 1, 1, ..., 3, 3, 3],
              [1, 1, 1, ..., 3, 3, 3],
              [1, 1, 1, ..., 3, 3, 3]],
      
             [[2, 2, 2, ..., 3, 3, 3],
              [2, 2, 2, ..., 3, 3, 3],
              [2, 2, 2, ..., 3, 3, 3],
              ...,
              [2, 2, 2, ..., 1, 1, 1],
              [2, 2, 2, ..., 1, 1, 1],
              [2, 2, 2, ..., 1, 1, 1]],
      
             [[2, 2, 2, ..., 3, 3, 3],
              [2, 2, 2, ..., 3, 3, 3],
              [2, 2, 2, ..., 2, 2, 2],
              ...,
              [1, 1, 1, ..., 3, 3, 3],
              [1, 1, 1, ..., 3, 3, 3],
              [1, 1, 1, ..., 3, 3, 3]],
      
             [[2, 2, 2, ..., 2, 2, 2],
              [2, 2, 2, ..., 2, 2, 2],
              [2, 2, 2, ..., 2, 2, 2],
              ...,
              [3, 3, 3, ..., 2, 2, 2],
              [3, 3, 3, ..., 2, 2, 2],
              [3, 3, 3, ..., 2, 2, 2]]], dtype=uint8)
  • crs :
    PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]

We specified the measurements to load because the complete list of measurements of this product is quite long!

[8]:
list(dc.index.products.get_by_name('s2a_ard_granule').measurements)
[8]:
['azimuthal_exiting',
 'azimuthal_incident',
 'exiting',
 'incident',
 'relative_azimuth',
 'relative_slope',
 'satellite_azimuth',
 'satellite_view',
 'solar_azimuth',
 'solar_zenith',
 'terrain_shadow',
 'fmask',
 'nbar_contiguity',
 'nbar_coastal_aerosol',
 'nbar_blue',
 'nbar_green',
 'nbar_red',
 'nbar_red_edge_1',
 'nbar_red_edge_2',
 'nbar_red_edge_3',
 'nbar_nir_1',
 'nbar_nir_2',
 'nbar_swir_2',
 'nbar_swir_3',
 'nbart_contiguity',
 'nbart_coastal_aerosol',
 'nbart_blue',
 'nbart_green',
 'nbart_red',
 'nbart_red_edge_1',
 'nbart_red_edge_2',
 'nbart_red_edge_3',
 'nbart_nir_1',
 'nbart_nir_2',
 'nbart_swir_2',
 'nbart_swir_3']

For the purposes of this notebook we will mostly restrict ourselves to the visible NBART bands.

[9]:
product = catalog['s2a_nbart_rgb']
product
[9]:
measurements:
- nbart_blue
- nbart_green
- nbart_red
product: s2a_ard_granule
[10]:
data = product.load(dc, **query)
rgb(data, col='time')
../../_images/notebooks_Frequently_used_code_Virtual_products_23_0.png

We will need the fmask band too if we want cloud-free images. However, fmask has a different resolution so the following does not work:

[11]:
product = catalog['s2a_nbart_rgb_and_fmask']
product
[11]:
measurements:
- nbart_blue
- nbart_green
- nbart_red
- fmask
product: s2a_ard_granule
[12]:
try:
    data = product.load(dc, **query)
except ValueError as e:
    print(e, file=sys.stderr)
Not all bands share the same pixel grid

In general, to load datasets from a product that is not “ingested”, such as those in s2a_ard_granule, datacube needs to be provided with the output_crs and resolution (and possibly align for sub-pixel alignment).

Note: In this notebook so far we have got away without specifying these. This is because the area of interest we chose happens to contain datasets from only one CRS (coordinate reference system), and the virtual product engine could figure out what grid specification to use. See sections on reprojecting and native load later on in this notebook for more on this.

We can set these in the virtual product so not to repeat ourselves:

[13]:
product = catalog['s2a_nbart_rgb_and_fmask_albers']
product
[13]:
measurements:
- nbart_blue
- nbart_green
- nbart_red
- fmask
output_crs: EPSG:3577
product: s2a_ard_granule
resampling:
  '*': average
  fmask: nearest
resolution:
- -20
- 20
[14]:
# we do not need to specify output crs or resolution here because the product does it for us
data = product.load(dc, **query)
rgb(data, col='time', col_wrap=4)
../../_images/notebooks_Frequently_used_code_Virtual_products_29_0.png

Data transformations

Cloud-free NBAR-T

So now that we have got the fmask measurement (i.e. band), we can mask out the cloud. First, let’s create the mask that only keeps the valid pixels.

[15]:
product = catalog['s2a_fmask']
product
[15]:
flags:
  fmask: valid
input:
  measurements:
  - fmask
  product: s2a_ard_granule
mask_measurement_name: fmask
transform: datacube.virtual.transformations.MakeMask
[16]:
data = product.load(dc, **query)
data
[16]:
Show/Hide data repr Show/Hide attributes
xarray.Dataset
    • time: 6
    • x: 441
    • y: 501
    • time
      (time)
      datetime64[ns]
      2018-01-06T23:52:41.026000 ... 2018-02-25T23:52:41.026000
      units :
      seconds since 1970-01-01 00:00:00
      array(['2018-01-06T23:52:41.026000000', '2018-01-16T23:52:41.026000000',
             '2018-01-26T23:52:41.026000000', '2018-02-05T23:52:41.026000000',
             '2018-02-15T23:52:41.026000000', '2018-02-25T23:52:41.026000000'],
            dtype='datetime64[ns]')
    • y
      (y)
      float64
      6.811e+06 6.811e+06 ... 6.801e+06
      units :
      metre
      resolution :
      -20.0
      crs :
      PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
      array([6810790., 6810770., 6810750., ..., 6800830., 6800810., 6800790.])
    • x
      (x)
      float64
      5.370e+05 5.371e+05 ... 5.458e+05
      units :
      metre
      resolution :
      20.0
      crs :
      PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
      array([537050., 537070., 537090., ..., 545810., 545830., 545850.])
    • fmask
      (time, y, x)
      bool
      True True True ... False False
      crs :
      PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
      array([[[ True,  True,  True, ...,  True,  True,  True],
              [ True,  True,  True, ...,  True,  True,  True],
              [ True,  True,  True, ...,  True,  True,  True],
              ...,
              [ True,  True,  True, ..., False, False, False],
              [ True,  True,  True, ..., False, False, False],
              [ True,  True,  True, ..., False, False, False]],
      
             [[ True,  True,  True, ...,  True,  True,  True],
              [ True,  True,  True, ...,  True,  True,  True],
              [ True,  True,  True, ...,  True,  True,  True],
              ...,
              [ True,  True,  True, ...,  True,  True,  True],
              [ True,  True,  True, ...,  True,  True,  True],
              [ True,  True,  True, ...,  True,  True,  True]],
      
             [[ True,  True,  True, ..., False, False, False],
              [ True,  True,  True, ..., False, False, False],
              [ True,  True,  True, ..., False, False, False],
              ...,
              [ True,  True,  True, ..., False, False, False],
              [ True,  True,  True, ..., False, False, False],
              [ True,  True,  True, ..., False, False, False]],
      
             [[False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False],
              ...,
              [False, False, False, ...,  True,  True,  True],
              [False, False, False, ...,  True,  True,  True],
              [False, False, False, ...,  True,  True,  True]],
      
             [[False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False],
              ...,
              [ True,  True,  True, ..., False, False, False],
              [ True,  True,  True, ..., False, False, False],
              [ True,  True,  True, ..., False, False, False]],
      
             [[False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False],
              ...,
              [False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False],
              [False, False, False, ..., False, False, False]]])
  • crs :
    PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]

We now have an fmask band as a boolean mask we want to apply.

Notice that the recipe this time specifies a transform to be applied (datacube.virtual.transformations.MakeMask, which is referred to by its short name make_mask in the file virtual-product-catalog.yaml). That transformation takes certain parameters (flags and mask_measurement_name in this case). It also needs the input, which is this case is the fmask band from the s2a_ard_granule product. You can have other virtual products as input too.

The next step is to apply it to the rest of the bands.

[17]:
product = catalog['s2a_cloud_free_nbart_rgb_albers']
product
[17]:
input:
  flags:
    fmask: valid
  input:
    measurements:
    - nbart_blue
    - nbart_green
    - nbart_red
    - fmask
    output_crs: EPSG:3577
    product: s2a_ard_granule
    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
[18]:
data = product.load(dc, **query)
rgb(data, col='time', col_wrap=4)
../../_images/notebooks_Frequently_used_code_Virtual_products_37_0.png

This time there were two (nested) applications of transforms. The outer transform apply_mask (i.e. datacube.virtual.transformations.ApplyMask) took as input the result of the inner transform make_mask (i.e. datacube.virtual.transformations.MakeMask) that was applied, in turn, to the Albers projected NBART bands. Just like mathematical function applications, this recipe too can be read inside out (or right to left) as a data transformation pipeline: we started from our customized s2a_ard_granule product, then we applied the make_mask transform to it, and then we applied the apply_mask transform to its output.

You could give apply_mask any boolean band (specified by mask_measurement_name) to apply, even if it is not created using make_mask.

Note: The recipe this time was pretty verbose. One advantage of using one catalog for many related virtual products is that parts of their recipes can be shared using some YAML features called anchors and references. See the recipe for the s2a_cloud_free_nbart_rgb_albers product that we just used in the virtual-product-catalog.yaml file (view on GitHub or in Jupyter). In this case, we recognize that the input to the transform is just the product s2a_nbart_rgb_and_fmask_albers that we had defined earlier.

Dilation of masks

The make_mask transform optionally can take an argument called dilation that specifies additional buffer for the mask in pixel units. The following recipe is essentially the same as that for s2a_cloud_free_nbart_rgb_albers except for the dilation set to 10 pixels.

[19]:
product = catalog['s2a_cloud_free_nbart_rgb_albers_dilated']
product
[19]:
dilation: 10
input:
  flags:
    fmask: valid
  input:
    measurements:
    - nbart_blue
    - nbart_green
    - nbart_red
    - fmask
    output_crs: EPSG:3577
    product: s2a_ard_granule
    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
[20]:
data = product.load(dc, **query)
rgb(data, col='time', col_wrap=4)
../../_images/notebooks_Frequently_used_code_Virtual_products_42_0.png

Band math

One particularly useful transformation is called expressions (i.e. datacube.virtual.transformations.Expressions) which allows us to use Python arithmetic on bands.

[21]:
product = catalog['s2a_NDVI']
product
[21]:
input:
  measurements:
  - nbart_red
  - nbart_nir_1
  product: s2a_ard_granule
output:
  NDVI:
    dtype: float32
    formula: (nbart_nir_1 - nbart_red) / (nbart_nir_1 + nbart_red)
transform: datacube.virtual.transformations.Expressions
[22]:
data = product.load(dc, **query)
data.NDVI.plot.imshow(col='time', col_wrap=4)
[22]:
<xarray.plot.facetgrid.FacetGrid at 0x7fb1c89569b0>
../../_images/notebooks_Frequently_used_code_Virtual_products_46_1.png

In a similar fashion, we can define a virtual product that calculates the Tasseled Cap indices too.

[23]:
product = catalog['s2a_tasseled_cap']
product
[23]:
input:
  measurements:
  - nbart_blue
  - nbart_green
  - nbart_red
  - nbart_nir_1
  - nbart_swir_2
  - nbart_swir_3
  output_crs: EPSG:3577
  product: s2a_ard_granule
  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
[24]:
data = product.load(dc, **query)
rgb(data, bands=['brightness', 'greenness', 'wetness'], col='time', col_wrap=4)
../../_images/notebooks_Frequently_used_code_Virtual_products_49_0.png

The cloud-free version would probably be more useful for analysis.

[25]:
product = catalog['s2a_cloud_free_tasseled_cap']
product
[25]:
input:
  input:
    flags:
      fmask: valid
    input:
      measurements:
      - nbart_blue
      - nbart_green
      - nbart_red
      - nbart_nir_1
      - nbart_swir_2
      - nbart_swir_3
      - fmask
      output_crs: EPSG:3577
      product: s2a_ard_granule
      resampling:
        '*': average
        fmask: nearest
      resolution:
      - -20
      - 20
    mask_measurement_name: fmask
    transform: datacube.virtual.transformations.MakeMask
  mask_measurement_name: fmask
  transform: datacube.virtual.transformations.ApplyMask
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
[26]:
data = product.load(dc, **query)

rgb(data, bands=['brightness', 'greenness', 'wetness'], col='time', col_wrap=4)
../../_images/notebooks_Frequently_used_code_Virtual_products_52_0.png

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: s2a_ard_granule
- measurements:
  - nbart_blue
  - nbart_green
  - nbart_red
  product: s2b_ard_granule
[28]:
data = product.load(dc, **query)
rgb(data, col='time', col_wrap=4)
../../_images/notebooks_Frequently_used_code_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]:
product = catalog['s2a_NDVI']
data = product.load(dc, dask_chunks={'time': 1, 'x': 2000, 'y': 2000}, **query)
data
[29]:
Show/Hide data repr Show/Hide attributes
xarray.Dataset
    • time: 6
    • x: 882
    • y: 1001
    • time
      (time)
      datetime64[ns]
      2018-01-06T23:52:41.026000 ... 2018-02-25T23:52:41.026000
      units :
      seconds since 1970-01-01 00:00:00
      array(['2018-01-06T23:52:41.026000000', '2018-01-16T23:52:41.026000000',
             '2018-01-26T23:52:41.026000000', '2018-02-05T23:52:41.026000000',
             '2018-02-15T23:52:41.026000000', '2018-02-25T23:52:41.026000000'],
            dtype='datetime64[ns]')
    • y
      (y)
      float64
      6.811e+06 6.811e+06 ... 6.801e+06
      units :
      metre
      resolution :
      -10.0
      crs :
      PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
      array([6810785., 6810775., 6810765., ..., 6800805., 6800795., 6800785.])
    • x
      (x)
      float64
      5.37e+05 5.371e+05 ... 5.459e+05
      units :
      metre
      resolution :
      10.0
      crs :
      PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
      array([537045., 537055., 537065., ..., 545835., 545845., 545855.])
    • NDVI
      (time, y, x)
      float32
      dask.array<chunksize=(1, 1001, 882), meta=np.ndarray>
      nodata :
      nan
      crs :
      PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
      Array Chunk
      Bytes 21.19 MB 3.53 MB
      Shape (6, 1001, 882) (1, 1001, 882)
      Count 90 Tasks 6 Chunks
      Type float32 numpy.ndarray
      882 1001 6
  • crs :
    PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
[30]:
data.compute().NDVI.plot.imshow(col='time', col_wrap=4)
[30]:
<xarray.plot.facetgrid.FacetGrid at 0x7fb1cc3c4400>
../../_images/notebooks_Frequently_used_code_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: s2a_ard_granule
      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)
/env/lib/python3.6/site-packages/xarray/core/nanops.py:142: RuntimeWarning: Mean of empty slice
  return np.nanmean(a, axis=axis, dtype=dtype)
../../_images/notebooks_Frequently_used_code_Virtual_products_66_1.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='s2a_ard_granule',
                                    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 0x7fb1cc29f898>
../../_images/notebooks_Frequently_used_code_Virtual_products_70_1.png

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

Bonus features

Renaming bands

The recipe for the expressions transform, besides the transform and input entries common to all transfoms, contains the output specification, which is a dictionary mapping output band names as keys to recipes for calculating those bands as values. These recipes can themselves be a dictionary as before, containing a formula and optionally specifying a dtype, a nodata value, and a units annotation. They can also be the name of an input band, in which case the input band is copied unaltered. We can therefore use the expressions transform to copy bands, rename bands, or even drop bands simply by not including it in the output section.

[35]:
product = catalog['s2a_nbart_rgb_renamed']
product
[35]:
input:
  measurements:
  - nbart_blue
  - nbart_green
  - nbart_red
  product: s2a_ard_granule
output:
  blue: nbart_blue
  green: nbart_green
  red: nbart_red
transform: datacube.virtual.transformations.Expressions
[36]:
data = product.load(dc, **query)
data
[36]:
Show/Hide data repr Show/Hide attributes
xarray.Dataset
    • time: 6
    • x: 882
    • y: 1001
    • time
      (time)
      datetime64[ns]
      2018-01-06T23:52:41.026000 ... 2018-02-25T23:52:41.026000
      units :
      seconds since 1970-01-01 00:00:00
      array(['2018-01-06T23:52:41.026000000', '2018-01-16T23:52:41.026000000',
             '2018-01-26T23:52:41.026000000', '2018-02-05T23:52:41.026000000',
             '2018-02-15T23:52:41.026000000', '2018-02-25T23:52:41.026000000'],
            dtype='datetime64[ns]')
    • y
      (y)
      float64
      6.811e+06 6.811e+06 ... 6.801e+06
      units :
      metre
      resolution :
      -10.0
      crs :
      PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
      array([6810785., 6810775., 6810765., ..., 6800805., 6800795., 6800785.])
    • x
      (x)
      float64
      5.37e+05 5.371e+05 ... 5.459e+05
      units :
      metre
      resolution :
      10.0
      crs :
      PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
      array([537045., 537055., 537065., ..., 545835., 545845., 545855.])
    • blue
      (time, y, x)
      int16
      210 202 216 294 ... 810 809 834 933
      units :
      1
      nodata :
      -999
      spectral_definition :
      {'response': [0.001206, 0.00204, 0.002623, 0.002738, 0.002746, 0.002076, 0.003246, 0.002224, 0.002789, 0.003127, 0.002377, 0.003776, 0.002856, 0.003063, 0.00607, 0.009028, 0.019547, 0.038955, 0.084088, 0.176255, 0.292197, 0.364612, 0.382418, 0.385789, 0.393447, 0.400158, 0.410291, 0.424686, 0.449286, 0.481594, 0.505323, 0.523406, 0.529543, 0.534688, 0.533786, 0.534656, 0.5381, 0.543691, 0.557717, 0.578585, 0.601967, 0.616037, 0.621092, 0.613597, 0.596062, 0.575863, 0.558063, 0.546131, 0.542099, 0.553602, 0.571684, 0.598269, 0.633236, 0.67337, 0.711752, 0.738396, 0.758249, 0.768325, 0.773367, 0.780468, 0.788363, 0.795449, 0.809151, 0.824011, 0.837709, 0.844983, 0.847328, 0.840111, 0.825797, 0.804778, 0.78694, 0.771578, 0.761923, 0.765487, 0.781682, 0.810031, 0.850586, 0.901671, 0.95467, 0.987257, 1.0, 0.986389, 0.908308, 0.724, 0.478913, 0.286992, 0.169976, 0.102833, 0.065163, 0.04116, 0.02508, 0.013112, 0.002585, 0.001095, 0.000308, 0.000441, 0.0, 0.0, 0.000443], 'wavelength': [440, 441, 442, 443, 444, 445, 446, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 495, 496, 497, 498, 499, 500, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510, 511, 512, 513, 514, 515, 516, 517, 518, 519, 520, 521, 522, 523, 524, 525, 526, 527, 528, 529, 530, 531, 532, 533, 534, 535, 536, 537, 538]}
      crs :
      PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
      array([[[ 210,  202,  216, ...,  414,  433,  436],
              [ 214,  195,  198, ...,  408,  406,  412],
              [ 209,  190,  207, ...,  423,  421,  427],
              ...,
              [ 277,  293,  265, ...,  217,  213,  233],
              [ 262,  305,  275, ...,  203,  217,  218],
              [ 256,  281,  298, ...,  189,  232,  221]],
      
             [[ 257,  225,  241, ...,  498,  489,  504],
              [ 229,  215,  218, ...,  507,  508,  537],
              [ 235,  234,  224, ...,  503,  504,  572],
              ...,
              [ 325,  320,  322, ...,  227,  248,  268],
              [ 309,  330,  316, ...,  244,  245,  252],
              [ 295,  308,  308, ...,  256,  261,  276]],
      
             [[ 408,  392,  447, ..., 4362, 4740, 4969],
              [ 404,  396,  418, ..., 4331, 4762, 5033],
              [ 411,  427,  418, ..., 4272, 4699, 5011],
              ...,
              [ 556,  513,  497, ...,  272,  289,  329],
              [ 545,  555,  511, ...,  293,  288,  316],
              [ 525,  560,  532, ...,  284,  292,  335]],
      
             [[ 621, 1050,  938, ...,  360,  373,  383],
              [1295, 1469,  896, ...,  354,  348,  353],
              [1546, 1126,  649, ...,  345,  348,  351],
              ...,
              [4595, 4624, 4607, ...,  258,  273,  312],
              [4589, 4602, 4597, ...,  278,  292,  300],
              [4625, 4640, 4618, ...,  278,  298,  310]],
      
             [[3825, 4217, 4435, ...,  373,  366,  365],
              [3886, 4184, 4399, ...,  367,  363,  358],
              [3809, 4110, 4311, ...,  366,  356,  354],
              ...,
              [ 485,  489,  497, ...,  270,  276,  325],
              [ 455,  473,  474, ...,  269,  264,  304],
              [ 447,  481,  488, ...,  280,  275,  291]],
      
             [[ 951,  948,  917, ..., 1543, 1543, 1542],
              [ 923,  924,  898, ..., 1615, 1567, 1559],
              [ 869,  847,  866, ..., 1644, 1569, 1583],
              ...,
              [ 943,  899,  850, ...,  916,  962, 1062],
              [ 903,  898,  869, ...,  857,  900, 1003],
              [ 863,  848,  830, ...,  809,  834,  933]]], dtype=int16)
    • green
      (time, y, x)
      int16
      422 380 416 511 ... 1141 1219 1262
      units :
      1
      nodata :
      -999
      spectral_definition :
      {'response': [0.00084, 0.016372, 0.037688, 0.080665, 0.169509, 0.341374, 0.573031, 0.746711, 0.828036, 0.868228, 0.888565, 0.891572, 0.87815, 0.860271, 0.843698, 0.834035, 0.832617, 0.844968, 0.867734, 0.897361, 0.933938, 0.966923, 0.990762, 1.0, 0.997812, 0.981107, 0.947088, 0.907183, 0.868656, 0.837622, 0.81291, 0.792434, 0.7851, 0.789606, 0.807011, 0.830458, 0.846433, 0.858402, 0.85799, 0.799824, 0.62498, 0.386994, 0.200179, 0.098293, 0.042844, 0.016512], 'wavelength': [537, 538, 539, 540, 541, 542, 543, 544, 545, 546, 547, 548, 549, 550, 551, 552, 553, 554, 555, 556, 557, 558, 559, 560, 561, 562, 563, 564, 565, 566, 567, 568, 569, 570, 571, 572, 573, 574, 575, 576, 577, 578, 579, 580, 581, 582]}
      crs :
      PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
      array([[[ 422,  380,  416, ...,  828,  828,  862],
              [ 414,  362,  381, ...,  822,  817,  867],
              [ 406,  365,  389, ...,  824,  833,  856],
              ...,
              [ 457,  436,  402, ...,  459,  422,  424],
              [ 446,  475,  475, ...,  389,  421,  456],
              [ 405,  442,  508, ...,  370,  441,  446]],
      
             [[ 474,  416,  445, ...,  907,  923,  940],
              [ 423,  390,  415, ...,  923,  948,  952],
              [ 454,  430,  410, ...,  932,  967,  950],
              ...,
              [ 505,  528,  515, ...,  428,  452,  480],
              [ 480,  528,  530, ...,  445,  459,  487],
              [ 477,  488,  523, ...,  455,  504,  503]],
      
             [[ 589,  541,  573, ..., 4304, 4636, 4905],
              [ 564,  547,  553, ..., 4234, 4647, 4989],
              [ 568,  567,  573, ..., 4229, 4651, 4994],
              ...,
              [ 760,  719,  708, ...,  360,  386,  505],
              [ 737,  798,  746, ...,  357,  396,  479],
              [ 685,  775,  778, ...,  383,  452,  447]],
      
             [[ 639,  705,  880, ...,  447,  456,  502],
              [ 662,  827,  960, ...,  427,  432,  451],
              [1092, 1368, 1108, ...,  421,  418,  409],
              ...,
              [4654, 4711, 4727, ...,  461,  462,  555],
              [4637, 4659, 4657, ...,  471,  494,  540],
              [4605, 4610, 4601, ...,  457,  519,  502]],
      
             [[3928, 4304, 4476, ...,  415,  424,  436],
              [3937, 4286, 4524, ...,  416,  414,  425],
              [3954, 4210, 4443, ...,  417,  421,  416],
              ...,
              [ 709,  707,  707, ...,  337,  352,  487],
              [ 678,  729,  731, ...,  356,  361,  434],
              [ 658,  709,  718, ...,  345,  346,  359]],
      
             [[1041, 1011, 1034, ..., 1767, 1771, 1754],
              [ 999,  974, 1015, ..., 1754, 1772, 1781],
              [ 957,  950,  996, ..., 1790, 1799, 1794],
              ...,
              [ 989,  937,  922, ..., 1279, 1370, 1430],
              [1013,  986,  940, ..., 1189, 1293, 1362],
              [ 992, 1001,  966, ..., 1141, 1219, 1262]]], dtype=int16)
    • red
      (time, y, x)
      int16
      230 205 225 340 ... 1299 1281 1333
      units :
      1
      nodata :
      -999
      spectral_definition :
      {'response': [0.002584, 0.034529, 0.14997, 0.464826, 0.817746, 0.965324, 0.983869, 0.9969, 1.0, 0.995449, 0.991334, 0.977215, 0.936802, 0.873776, 0.814166, 0.776669, 0.764864, 0.775091, 0.801359, 0.830828, 0.857112, 0.883581, 0.90895, 0.934759, 0.955931, 0.96811, 0.973219, 0.971572, 0.969003, 0.965712, 0.960481, 0.944811, 0.884152, 0.706167, 0.422967, 0.189853, 0.063172, 0.020615, 0.002034], 'wavelength': [646, 647, 648, 649, 650, 651, 652, 653, 654, 655, 656, 657, 658, 659, 660, 661, 662, 663, 664, 665, 666, 667, 668, 669, 670, 671, 672, 673, 674, 675, 676, 677, 678, 679, 680, 681, 682, 683, 684]}
      crs :
      PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
      array([[[ 230,  205,  225, ...,  511,  540,  565],
              [ 215,  197,  192, ...,  514,  526,  565],
              [ 209,  186,  197, ...,  532,  525,  522],
              ...,
              [ 344,  354,  334, ...,  213,  214,  234],
              [ 321,  368,  357, ...,  186,  213,  209],
              [ 320,  351,  358, ...,  188,  239,  209]],
      
             [[ 306,  241,  317, ...,  633,  632,  630],
              [ 259,  229,  238, ...,  633,  629,  629],
              [ 294,  267,  245, ...,  654,  622,  671],
              ...,
              [ 391,  397,  412, ...,  222,  253,  284],
              [ 381,  400,  407, ...,  233,  250,  277],
              [ 377,  405,  399, ...,  271,  307,  320]],
      
             [[ 386,  336,  395, ..., 4322, 4639, 4908],
              [ 356,  309,  328, ..., 4226, 4673, 4997],
              [ 360,  333,  341, ..., 4181, 4663, 5047],
              ...,
              [ 627,  534,  513, ...,  203,  234,  353],
              [ 606,  601,  541, ...,  214,  243,  330],
              [ 536,  635,  597, ...,  224,  286,  350]],
      
             [[ 378,  403,  600, ...,  318,  332,  379],
              [ 365,  499,  660, ...,  301,  311,  326],
              [ 380,  542,  613, ...,  297,  289,  276],
              ...,
              [4760, 4850, 4847, ...,  209,  235,  311],
              [4766, 4815, 4816, ...,  239,  269,  280],
              [4723, 4740, 4757, ...,  233,  280,  318]],
      
             [[4077, 4437, 4591, ...,  292,  300,  311],
              [4085, 4462, 4657, ...,  290,  298,  304],
              [4112, 4416, 4600, ...,  291,  299,  293],
              ...,
              [ 565,  579,  579, ...,  198,  210,  308],
              [ 540,  585,  587, ...,  206,  218,  277],
              [ 521,  566,  569, ...,  215,  216,  250]],
      
             [[ 979,  953, 1028, ..., 1724, 1724, 1740],
              [ 970,  964,  988, ..., 1695, 1716, 1727],
              [ 960,  965,  990, ..., 1684, 1728, 1714],
              ...,
              [ 769,  733,  796, ..., 1311, 1288, 1303],
              [ 859,  811,  844, ..., 1287, 1283, 1313],
              [ 913,  839,  851, ..., 1299, 1281, 1333]]], dtype=int16)
  • crs :
    PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]

Masking invalid data

The data loaded by datacube usually have a nodata value specified. Often, we want it to be replaced with the NaN value. Since the expressions transform is nodata-aware, we get the same effect just by asking the expressions transform to convert the band to float.

[37]:
product = catalog['s2a_nbart_rgb_nodata_masked']
product
[37]:
input:
  measurements:
  - nbart_blue
  - nbart_green
  - nbart_red
  product: s2a_ard_granule
output:
  blue:
    dtype: float32
    formula: nbart_blue
  green:
    dtype: float32
    formula: nbart_green
  red:
    dtype: float32
    formula: nbart_red
transform: datacube.virtual.transformations.Expressions
[38]:
data = product.load(dc, **query)
data
[38]:
Show/Hide data repr Show/Hide attributes
xarray.Dataset
    • time: 6
    • x: 882
    • y: 1001
    • time
      (time)
      datetime64[ns]
      2018-01-06T23:52:41.026000 ... 2018-02-25T23:52:41.026000
      units :
      seconds since 1970-01-01 00:00:00
      array(['2018-01-06T23:52:41.026000000', '2018-01-16T23:52:41.026000000',
             '2018-01-26T23:52:41.026000000', '2018-02-05T23:52:41.026000000',
             '2018-02-15T23:52:41.026000000', '2018-02-25T23:52:41.026000000'],
            dtype='datetime64[ns]')
    • y
      (y)
      float64
      6.811e+06 6.811e+06 ... 6.801e+06
      units :
      metre
      resolution :
      -10.0
      crs :
      PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
      array([6810785., 6810775., 6810765., ..., 6800805., 6800795., 6800785.])
    • x
      (x)
      float64
      5.37e+05 5.371e+05 ... 5.459e+05
      units :
      metre
      resolution :
      10.0
      crs :
      PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
      array([537045., 537055., 537065., ..., 545835., 545845., 545855.])
    • blue
      (time, y, x)
      float32
      210.0 202.0 216.0 ... 834.0 933.0
      nodata :
      nan
      crs :
      PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
      array([[[ 210.,  202.,  216., ...,  414.,  433.,  436.],
              [ 214.,  195.,  198., ...,  408.,  406.,  412.],
              [ 209.,  190.,  207., ...,  423.,  421.,  427.],
              ...,
              [ 277.,  293.,  265., ...,  217.,  213.,  233.],
              [ 262.,  305.,  275., ...,  203.,  217.,  218.],
              [ 256.,  281.,  298., ...,  189.,  232.,  221.]],
      
             [[ 257.,  225.,  241., ...,  498.,  489.,  504.],
              [ 229.,  215.,  218., ...,  507.,  508.,  537.],
              [ 235.,  234.,  224., ...,  503.,  504.,  572.],
              ...,
              [ 325.,  320.,  322., ...,  227.,  248.,  268.],
              [ 309.,  330.,  316., ...,  244.,  245.,  252.],
              [ 295.,  308.,  308., ...,  256.,  261.,  276.]],
      
             [[ 408.,  392.,  447., ..., 4362., 4740., 4969.],
              [ 404.,  396.,  418., ..., 4331., 4762., 5033.],
              [ 411.,  427.,  418., ..., 4272., 4699., 5011.],
              ...,
              [ 556.,  513.,  497., ...,  272.,  289.,  329.],
              [ 545.,  555.,  511., ...,  293.,  288.,  316.],
              [ 525.,  560.,  532., ...,  284.,  292.,  335.]],
      
             [[ 621., 1050.,  938., ...,  360.,  373.,  383.],
              [1295., 1469.,  896., ...,  354.,  348.,  353.],
              [1546., 1126.,  649., ...,  345.,  348.,  351.],
              ...,
              [4595., 4624., 4607., ...,  258.,  273.,  312.],
              [4589., 4602., 4597., ...,  278.,  292.,  300.],
              [4625., 4640., 4618., ...,  278.,  298.,  310.]],
      
             [[3825., 4217., 4435., ...,  373.,  366.,  365.],
              [3886., 4184., 4399., ...,  367.,  363.,  358.],
              [3809., 4110., 4311., ...,  366.,  356.,  354.],
              ...,
              [ 485.,  489.,  497., ...,  270.,  276.,  325.],
              [ 455.,  473.,  474., ...,  269.,  264.,  304.],
              [ 447.,  481.,  488., ...,  280.,  275.,  291.]],
      
             [[ 951.,  948.,  917., ..., 1543., 1543., 1542.],
              [ 923.,  924.,  898., ..., 1615., 1567., 1559.],
              [ 869.,  847.,  866., ..., 1644., 1569., 1583.],
              ...,
              [ 943.,  899.,  850., ...,  916.,  962., 1062.],
              [ 903.,  898.,  869., ...,  857.,  900., 1003.],
              [ 863.,  848.,  830., ...,  809.,  834.,  933.]]], dtype=float32)
    • green
      (time, y, x)
      float32
      422.0 380.0 416.0 ... 1219.0 1262.0
      nodata :
      nan
      crs :
      PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
      array([[[ 422.,  380.,  416., ...,  828.,  828.,  862.],
              [ 414.,  362.,  381., ...,  822.,  817.,  867.],
              [ 406.,  365.,  389., ...,  824.,  833.,  856.],
              ...,
              [ 457.,  436.,  402., ...,  459.,  422.,  424.],
              [ 446.,  475.,  475., ...,  389.,  421.,  456.],
              [ 405.,  442.,  508., ...,  370.,  441.,  446.]],
      
             [[ 474.,  416.,  445., ...,  907.,  923.,  940.],
              [ 423.,  390.,  415., ...,  923.,  948.,  952.],
              [ 454.,  430.,  410., ...,  932.,  967.,  950.],
              ...,
              [ 505.,  528.,  515., ...,  428.,  452.,  480.],
              [ 480.,  528.,  530., ...,  445.,  459.,  487.],
              [ 477.,  488.,  523., ...,  455.,  504.,  503.]],
      
             [[ 589.,  541.,  573., ..., 4304., 4636., 4905.],
              [ 564.,  547.,  553., ..., 4234., 4647., 4989.],
              [ 568.,  567.,  573., ..., 4229., 4651., 4994.],
              ...,
              [ 760.,  719.,  708., ...,  360.,  386.,  505.],
              [ 737.,  798.,  746., ...,  357.,  396.,  479.],
              [ 685.,  775.,  778., ...,  383.,  452.,  447.]],
      
             [[ 639.,  705.,  880., ...,  447.,  456.,  502.],
              [ 662.,  827.,  960., ...,  427.,  432.,  451.],
              [1092., 1368., 1108., ...,  421.,  418.,  409.],
              ...,
              [4654., 4711., 4727., ...,  461.,  462.,  555.],
              [4637., 4659., 4657., ...,  471.,  494.,  540.],
              [4605., 4610., 4601., ...,  457.,  519.,  502.]],
      
             [[3928., 4304., 4476., ...,  415.,  424.,  436.],
              [3937., 4286., 4524., ...,  416.,  414.,  425.],
              [3954., 4210., 4443., ...,  417.,  421.,  416.],
              ...,
              [ 709.,  707.,  707., ...,  337.,  352.,  487.],
              [ 678.,  729.,  731., ...,  356.,  361.,  434.],
              [ 658.,  709.,  718., ...,  345.,  346.,  359.]],
      
             [[1041., 1011., 1034., ..., 1767., 1771., 1754.],
              [ 999.,  974., 1015., ..., 1754., 1772., 1781.],
              [ 957.,  950.,  996., ..., 1790., 1799., 1794.],
              ...,
              [ 989.,  937.,  922., ..., 1279., 1370., 1430.],
              [1013.,  986.,  940., ..., 1189., 1293., 1362.],
              [ 992., 1001.,  966., ..., 1141., 1219., 1262.]]], dtype=float32)
    • red
      (time, y, x)
      float32
      230.0 205.0 225.0 ... 1281.0 1333.0
      nodata :
      nan
      crs :
      PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]
      array([[[ 230.,  205.,  225., ...,  511.,  540.,  565.],
              [ 215.,  197.,  192., ...,  514.,  526.,  565.],
              [ 209.,  186.,  197., ...,  532.,  525.,  522.],
              ...,
              [ 344.,  354.,  334., ...,  213.,  214.,  234.],
              [ 321.,  368.,  357., ...,  186.,  213.,  209.],
              [ 320.,  351.,  358., ...,  188.,  239.,  209.]],
      
             [[ 306.,  241.,  317., ...,  633.,  632.,  630.],
              [ 259.,  229.,  238., ...,  633.,  629.,  629.],
              [ 294.,  267.,  245., ...,  654.,  622.,  671.],
              ...,
              [ 391.,  397.,  412., ...,  222.,  253.,  284.],
              [ 381.,  400.,  407., ...,  233.,  250.,  277.],
              [ 377.,  405.,  399., ...,  271.,  307.,  320.]],
      
             [[ 386.,  336.,  395., ..., 4322., 4639., 4908.],
              [ 356.,  309.,  328., ..., 4226., 4673., 4997.],
              [ 360.,  333.,  341., ..., 4181., 4663., 5047.],
              ...,
              [ 627.,  534.,  513., ...,  203.,  234.,  353.],
              [ 606.,  601.,  541., ...,  214.,  243.,  330.],
              [ 536.,  635.,  597., ...,  224.,  286.,  350.]],
      
             [[ 378.,  403.,  600., ...,  318.,  332.,  379.],
              [ 365.,  499.,  660., ...,  301.,  311.,  326.],
              [ 380.,  542.,  613., ...,  297.,  289.,  276.],
              ...,
              [4760., 4850., 4847., ...,  209.,  235.,  311.],
              [4766., 4815., 4816., ...,  239.,  269.,  280.],
              [4723., 4740., 4757., ...,  233.,  280.,  318.]],
      
             [[4077., 4437., 4591., ...,  292.,  300.,  311.],
              [4085., 4462., 4657., ...,  290.,  298.,  304.],
              [4112., 4416., 4600., ...,  291.,  299.,  293.],
              ...,
              [ 565.,  579.,  579., ...,  198.,  210.,  308.],
              [ 540.,  585.,  587., ...,  206.,  218.,  277.],
              [ 521.,  566.,  569., ...,  215.,  216.,  250.]],
      
             [[ 979.,  953., 1028., ..., 1724., 1724., 1740.],
              [ 970.,  964.,  988., ..., 1695., 1716., 1727.],
              [ 960.,  965.,  990., ..., 1684., 1728., 1714.],
              ...,
              [ 769.,  733.,  796., ..., 1311., 1288., 1303.],
              [ 859.,  811.,  844., ..., 1287., 1283., 1313.],
              [ 913.,  839.,  851., ..., 1299., 1281., 1333.]]], dtype=float32)
  • crs :
    PROJCS["WGS 84 / UTM zone 56S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",153],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32756"]]

Native load

This feature is in virtual product only at the moment but not in Datacube.load. We plan to integrate this there in the future.

If all the datasets to be loaded are in the same CRS and all the bands to be loaded have the same resolution and sub-pixel alignment, virtual products will load them “natively”, that is, in their common grid specification, as we have seen before in this notebook.

Sometimes the datasets to be loaded are all in the same CRS but the bands do not have the same resolution. We can specify which resolution to use by setting the like keyword argument to the name of a band.

[39]:
product = catalog['s2a_nbart_rgb_and_fmask_native']
product
[39]:
like: fmask
measurements:
- nbart_blue
- nbart_green
- nbart_red
- fmask
product: s2a_ard_granule
resampling:
  '*': average
  fmask: nearest
[40]:
data = product.load(dc, **query)
rgb(data, col='time', col_wrap=4)
../../_images/notebooks_Frequently_used_code_Virtual_products_84_0.png

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: s2a_ard_granule
      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_Frequently_used_code_Virtual_products_88_0.png

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 Slack channel 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: April 2020

Compatible datacube version:

[43]:
print(datacube.__version__)
1.7+262.g1cf3cea8

Tags

Browse all available tags on the DEA User Guide’s Tags Index

Tags: NCI compatible, sandbox compatible, sentinel 2, dea_plotting, rgb, virtual products, NDVI, tasseled cap, cloud masking, Dask, image compositing, statistics, pixel quality, combining data, native load, reprojecting