How to align X and Y pixel values to a boundary using the Open Data Cube

Question

I am trying to load the following dataset into the DEA Sandbox:

import datacube 
dc = datacube.Datacube(app="Landsat")
import sys
sys.path.append("./Scripts")
from dea_datahandling import load_ard

ds = load_ard(dc = dc,
        products = ['ga_ls8c_ard_3'],
        time = ('2015-05-01', '2015-05-31'),
        x = (691200, 692340),
        y = (6093660, 6094260),
        crs = 'EPSG:28355',
        resampling = 'nearest',
        group_by = 'solar_day',
        measurements = ['nbart_red', 'nbart_nir', 'oa_fmask'],
        output_crs = 'EPSG:28355',
        resolution = (-30, 30))

When I view the x and y values, I notice they don't align up to the boundary I specified:

ds.x
xarray.DataArray'x'x: 38
array([691215., 691245., 691275., 691305., 691335., 691365., 691395., 691425.,
       691455., 691485., 691515., 691545., 691575., 691605., 691635., 691665.,
       691695., 691725., 691755., 691785., 691815., 691845., 691875., 691905.,
       691935., 691965., 691995., 692025., 692055., 692085., 692115., 692145.,
       692175., 692205., 692235., 692265., 692295., 692325.])
ds.y
xarray.DataArray'y'y: 20
array([6094245., 6094215., 6094185., 6094155., 6094125., 6094095., 6094065.,
       6094035., 6094005., 6093975., 6093945., 6093915., 6093885., 6093855.,
       6093825., 6093795., 6093765., 6093735., 6093705., 6093675.])

Is there a way to rectify this? I have a (areal) dataset with veg. structure data on a 30 x 30 m grid, based on these x and y boundaries and would like to calculate and then overlay NDVI from Landsat.

Answer 1

To follow up on @huite-bootsma's comment, I would try plotting your 30 x 30 m gridded vegetation data over the data you've loaded from Open Data Cube. My feeling is that the pixel edges should line up correctly even though the coordinates are offset by 15 m to reflect the pixel midpoints. For example, here's a demonstration of your loaded data overlayed with points starting at your 691200, 6094260:

import datacube 
dc = datacube.Datacube(app="Landsat")
import sys
sys.path.append("./Scripts")
from dea_datahandling import load_ard

ds = load_ard(dc = dc,
        products = ['ga_ls8c_ard_3'],
        time = ('2015-05-01', '2015-05-31'),
        x = (691200, 692340),
        y = (6093660, 6094260),
        crs = 'EPSG:28355',
        resampling = 'nearest',
        group_by = 'solar_day',
        measurements = ['nbart_red', 'nbart_nir', 'oa_fmask'],
        output_crs = 'EPSG:28355',
        resolution = (-30, 30))

# Plot top-left corner of array
ds.nbart_red.isel(time=3, x=[0, 1, 2, 3], y=[0, 1, 2, 3]).plot()

# Overlay points aligned to original coords
plt.plot([691200, 691230, 691260, 691290], 
         [6094260, 6094260, 6094260, 6094260], 'ro', markersize=12)
plt.plot([691200, 691230, 691260, 691290], 
         [6094230, 6094230, 6094230, 6094230], 'ro', markersize=12)

Note that this is only the case when re-projecting the Landsat data like you do above. If you instead wanted to load the data in its native CRS and resolution, there's a special align=(15, 15) param you need to provide to dc.load to account for the fact that Landsat data is stored on file with coordinates that represent pixel centres, not edges. This is the relevant section of our more detailed Introduction to DEA Landsat Surface Reflectance (Geoscience Australia Landsat Collection 3) guide:

In addition to output_crs and resolution, one final parameter is required to construct the spatial grid used for natively loading Landsat data: align. By default, datacube assumes that pixel edges are aligned such that x=0 and y=0 lines fall on pixel edges. However, Landsat data is stored on file with coordinates that define the centre (not the edge) of each pixel. When natively loading data, to ensure the spatial pixel grid created by datacube is exactly the same grid as used by Landsat imagery, we need to shift this grid by 15 m in both directions by specifying align=(15, 15). Otherwise, natively loaded pixels will be offset by half a pixel from their true location.

To determine whether you need to use align when loading DEA Landsat Surface Reflectance data, all the following must apply:

• You are loading Landsat data from the ga_ls5t_ard_3, ga_ls7e_ard_3 and ga_ls8c_ard_3 products that define pixel coordinates by their centers rather than pixel edges

• You want to load data in its native projection without a half pixel offset

• You are supplying a native UTM zone CRS generated by mostcommon_crs or copied from a datacube dataset

• You are supplying a native resolution (e.g. resolution=(-30, 30))