3

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.

1
  • 3
    They do line up: xarray shows midpoint coordinates. So your leftmost cell has its bounds from 691200 to 691230, with its midpoint at 691215. There's always N + 1 bounds for N cells. Since we'd like to selection on coordinates, having equal sized coordinate arrays is very useful. Note that the CF conventions have a description for bounds arrays, but I don't think xarray has everything ironed out yet: github.com/pydata/xarray/issues/1475 Jul 5, 2021 at 20:32

1 Answer 1

4

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)

Grid overlayed with points

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

2
  • Incredible answer!
    – Alex Leith
    Jul 8, 2021 at 1:58
  • Excellent answer - very clear. Thank you so much. Jul 9, 2021 at 9:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.