Tooling for converting STAC metadata to ODC data model

~Required because odc.stac.configure_s3_access is part of the public API, which imports from _aws which requires botocore.~

Adds botocore as an optional dependency. Includes:

• Nicer import error message when botocore is not installed
• Skip aws tests when botocore is not installed
• CI with/without botocore
• README instructions on using extra_requires with pip

I'm relatively new to odc-stac and Digital Earth Africa, and trying out the stac client and odc-stac package to access sentinel-1. I started with the digital earth africa sentinel2 example code and that worked fine, but when I tried adapting it for sentinel-1, I'm getting AssertionError: Only axis-aligned geoboxes are currently supported, despite trying both projected and lon/lat specification.

I'm guessing this is user error rather than a bug, so suggestions welcome! :)

Hi there. Thanks a lot for this awesome tool!

I am working with some Sentinel-2 data which I am downsampling for testing purposes.

I wonder whether the following behaviour is expected...

Let's get some STAC items (with a reproducible example from here):

import pystac_client
import planetary_computer

catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)

time_range = "2020-12-01/2020-12-31"
bbox = [-122.2751, 47.5469, -121.9613, 47.7458]

search = catalog.search(collections=["landsat-c2-l2"], bbox=bbox, datetime=time_range)
items = search.get_all_items()

...get the size of one band at full resolution:

out = odc_stac.load(
    items,
    bands=['blue'],
    chunks={'x': 256, 'y': 256}
)
x = out.sizes['x']
y = out.sizes['y']
res = items[0].to_dict()['assets']['blue']['raster:bands'][0]['spatial_resolution']

..and then load the same band with different resolutions:

for factor in [2, 3, 4]:
    print('factor:', factor)
    out = odc_stac.load(
        items,
        bands=['blue'],
        resolution=res * factor,
        chunks={'x': 256, 'y': 256}
    )
    print(f'expected size: {y // factor} x {x // 2}')
    print(f'actual size:   {out.sizes["y"]} x {out.sizes["x"]}')
    print('-' * 30)

This gives me:

factor: 2
expected size: 3986 x 6186
actual size:   3987 x 6187
------------------------------
factor: 3
expected size: 2657 x 6186
actual size:   2658 x 4125
------------------------------
factor: 4
expected size: 1993 x 6186
actual size:   1994 x 3094
------------------------------

Although it is possible to use odc.stac.load to load items whose assets are in requester-pays S3 buckets (e.g. landsat-c2l2-sr), the documentation is not entirely clear on how to accomplish this. Setting the environment variable AWS_REQUEST_PAYER to requester does the trick, but I was unable to use odc.stac.configure_rio(aws={"requester_pays": True}) to load the data. I'm happy to make a docs PR here to address this topic, unless there is a more native way of going about this.

Captures important context about benchmarking tasks and provides some reporting facilities for that data.

I have added docs for odc.stac.bench.* too, but there are some issues with dataclasss and auto-type-hints in sphinx. Would be nice to figure out how to turn off documentation for autogenerated __init__ method.

Testing how well new way of updating notebooks works.

There should be a netlify link posted that allows to review documentation changes before they make it to read the docs.

Notebook changes

Significant improvements in "load data from AWS Sentinel2 collection".

• Don't sure odc.ui.show_datasets instead construct ipyleaflet.Map manually
• Plot STAC footprints and query rectangle
• Demonstrate how to load smaller region
• Fix CSS problems when displaying xarray widget state

To review changes see here: https://618a4f6414a45b24dd03d7fb--odc-stac-docs.netlify.app/notebooks/stac-load-e84-aws.html

• add links to documentation and sample notebooks
• remove crs/resolution parameters from sample
• Update sample notebook for planetary computer, removed no longer needed installation instructions

Closes #38

This throws a ValueError when trying to build a Dataset from a GOES-CMI STAC item from the Planetary Computer:

import pystac
import planetary_computer
import odc.stac


item = planetary_computer.sign(
    pystac.Item.from_file("https://planetarycomputer.microsoft.com/api/stac/v1/collections/goes-cmi/items/OR_ABI-L2-C-M6_G16_s20200012356184")
)

# Drop assets in different projections
# https://github.com/opendatacube/odc-stac/issues/70
bands = {"C01_2km", "C02_2km", "C03_2km"}
item.assets = {k: v for k, v in item.assets.items() if k in bands}

ds = odc.stac.stac_load([item], bands=bands, chunks={})

The error is

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [7], in <cell line: 15>()
     12 bands = {"C01_2km", "C02_2km", "C03_2km"}
     13 item.assets = {k: v for k, v in item.assets.items() if k in bands}
---> 15 ds = odc.stac.stac_load([item], bands=bands, chunks={})

File /srv/conda/envs/notebook/lib/python3.8/site-packages/odc/stac/_load.py:460, in load(items, bands, groupby, resampling, dtype, chunks, pool, crs, resolution, align, geobox, bbox, lon, lat, x, y, like, geopolygon, progress, stac_cfg, patch_url, preserve_original_order, **kw)
    457 items = list(items)
    458 _parsed = list(parse_items(items, cfg=stac_cfg))
--> 460 gbox = output_geobox(
    461     _parsed,
    462     bands=bands,
    463     crs=crs,
    464     resolution=resolution,
    465     align=align,
    466     geobox=geobox,
    467     like=like,
    468     geopolygon=geopolygon,
    469     bbox=bbox,
    470     lon=lon,
    471     lat=lat,
    472     x=x,
    473     y=y,
    474 )
    476 if gbox is None:
    477     raise ValueError("Failed to auto-guess CRS/resolution.")

File /srv/conda/envs/notebook/lib/python3.8/site-packages/odc/stac/_mdtools.py:805, in output_geobox(items, bands, crs, resolution, align, geobox, like, geopolygon, bbox, lon, lat, x, y)
    803 x0, y0, x1, y1 = _compute_bbox(items, crs)
    804 geopolygon = geom.box(x0, y0, x1, y1, crs)
--> 805 return GeoBox.from_geopolygon(geopolygon, resolution=resolution, align=align)

File /srv/conda/envs/notebook/lib/python3.8/site-packages/odc/geo/geobox.py:227, in GeoBox.from_geopolygon(geopolygon, resolution, crs, align, shape, tight, anchor, tol)
    224 else:
    225     geopolygon = geopolygon.to_crs(crs)
--> 227 return GeoBox.from_bbox(
    228     geopolygon.boundingbox,
    229     crs,
    230     shape=shape,
    231     resolution=resolution,
    232     anchor=anchor,
    233     tol=tol,
    234     tight=tight,
    235 )

File /srv/conda/envs/notebook/lib/python3.8/site-packages/odc/geo/geobox.py:151, in GeoBox.from_bbox(bbox, crs, tight, shape, resolution, anchor, tol)
    149     offy, ny = snap_grid(bbox.bottom, bbox.top, ry, None, tol=tol)
    150 else:
--> 151     offx, nx = snap_grid(bbox.left, bbox.right, rx, _snap.x, tol=tol)
    152     offy, ny = snap_grid(bbox.bottom, bbox.top, ry, _snap.y, tol=tol)
    154 affine = Affine.translation(offx, offy) * Affine.scale(rx, ry)

File /srv/conda/envs/notebook/lib/python3.8/site-packages/odc/geo/math.py:172, in snap_grid(x0, x1, res, off_pix, tol)
    169     return x1, max(nx, 1)
    171 off = off_pix * abs(res)
--> 172 _tx, nx = _snap_edge(x0 - off, x1 - off, res, tol)
    173 return _tx + off, nx

File /srv/conda/envs/notebook/lib/python3.8/site-packages/odc/geo/math.py:140, in _snap_edge(x0, x1, res, tol)
    138 assert x1 >= x0
    139 if res > 0:
--> 140     return _snap_edge_pos(x0, x1, res, tol)
    141 _tx, nx = _snap_edge_pos(x0, x1, -res, tol)
    142 tx = _tx + nx * (-res)

File /srv/conda/envs/notebook/lib/python3.8/site-packages/odc/geo/math.py:132, in _snap_edge_pos(x0, x1, res, tol)
    130 assert x1 >= x0
    131 _x0 = floor(maybe_int(x0 / res, tol))
--> 132 _x1 = ceil(maybe_int(x1 / res, tol))
    133 nx = max(1, _x1 - _x0)
    134 return _x0 * res, nx

File /srv/conda/envs/notebook/lib/python3.8/site-packages/odc/geo/math.py:57, in maybe_int(x, tol)
     50 def maybe_int(x: float, tol: float) -> Union[int, float]:
     51     """
     52     Turn almost ints to actual ints.
     53 
     54     pass through other values unmodified.
     55     """
---> 57     x_whole, x_part = split_float(x)
     59     if abs(x_part) < tol:  # almost int
     60         return int(x_whole)

File /srv/conda/envs/notebook/lib/python3.8/site-packages/odc/geo/math.py:39, in split_float(x)
     28 def split_float(x: float) -> Tuple[float, float]:
     29     """
     30     Split float number into whole and fractional parts.
     31 
   (...)
     37     :return: ``whole, fraction``
     38     """
---> 39     x_part = fmod(x, 1.0)
     40     x_whole = x - x_part
     41     if x_part > 0.5:

ValueError: math domain error

I haven't dug too deeply yet, but x there at the end, which I think is something like bbox.right is inf.

https://planetarycomputer.microsoft.com/api/stac/v1/collections/goes-cmi/items/OR_ABI-L2-C-M6_G16_s20200012356184 is a link to the STAC metadata.

probably a bunch of things in odc.stac.testing.{_xr,_cog} should be generalized and moved to odc-geo and made available when rasterio and xarray are installed.

_cog.py is basically copy from datacube with dask parts removed, new things are write_cog_layers which allows arbitrary overviews to be written, this will be useful when testing read from overviews properly, so far that was only tested in the notebook.

also includes a fix for #46

Hi all, first of all thanks for this great tool!

I'm currently trying to aggregate raster values of a xarray.DataSet created with odc-stac based on several hundred STAC items, similar to this reproducible example:

import pystac_client
import planetary_computer
from odc.stac import load

catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1", 
    modifier=planetary_computer.sign_inplace
) 

time_range = "2020-01-01/2020-12-31"
bbox = [-122.2751, 47.5469, -121.9613, 47.7458]

search = catalog.search(collections=["landsat-c2-l2"], bbox=bbox, datetime=time_range)
items = search.get_all_items()

data = load(items, bands=["green"], chunks={"x": 256, "y": 256})
res = data.sum(dim="time").compute()

If, however, one of the many items is corrupt, it is very hard to indentify the faulty data source just from the RasterioIOError that gets returned:

# lets break a random href of an Landsat item
items[0].assets['green'].href = "https://faulty_url"

Excecuting the code above returns:

Exception has occurred: RasterioIOError
HTTP response code: 404

During handling of the above exception, another exception occurred:

  File "/path/to/file", line 16, in <module>
    res = data.sum(dim="time").compute()

Is there an option to extend the logging in odc.stac.load in order to identify which Item/Asset rasterio wasn't able to open?

Overview

This PR fixes how odc-stac uses the Electro-Optical STAC extension to generate band aliases. Specifically:

• Asset keys are used as the canonical name for bands, instead of eo:bands -> name
• eo:bands -> common_name and eo:bands -> name are both used as aliases back to asset keys

This enables support for the common STAC pattern of "one band per asset".

Background

Initially raised in https://github.com/Element84/geo-notebooks/issues/8. USGS has changed the structure of their STAC items such that eo:band -> name does not match the asset key (example item here). Instead, eo:band -> common_name is the asset key.

More generally, odc-stac was incorrectly reading eo:bands information from the Item instead of from the Assets. From the EO extension:

Each Asset should specify its own band object. If the individual bands are repeated in different assets they should all use the same values and include the optional name field to enable clients to easily combine and summarize the bands.

The eo:bands array may optionally be used in the Item Properties to summarize the available bands in the assets. This should be a 'union' of all the possible bands represented in assets. It should be considered merely informative - clients should rely on the eo:bands of each asset.

This PR fixes odc-stac to read band information from each asset.

Details

• If name or common_name matches the asset key, it is not added to the aliases.
• Multi-band assets are not supported via the alias mechanism. It's not clear if this was ever correctly supported previously, but now it is explicitly unsupported with a warning.
• Test case included with a new Landsat item demonstrating its layout.
• Since this does not change the public API, IMO this could be incorporated in a MINOR release, e.g. v0.3.1.

This is related to #99 and #97.

Now that "continue in the face of errors" feature is implemented see (#99). When IO errors happen they are logged and forgotten, however that means there is no easy way of knowing

1. Whether any failures have happened
2. Which inputs were corrupt

Quote from issue #99

For local load case errors can be reported via an attribute, something like xx.load_failures, or via thread-local global variable that reports failures for the last load that happened: odc.stac.load_failures(). With Dask, computation is delayed and distributed, so certainly can not be done as an attribute. In case of Dask we could maintain a global dictionary mapping load-task-id -> [Failure], so would be something like odc.stac.load_failures(xx) that would lookup failures that might have happened when processing xx, complexity with this approach is mainly to do with cache life-cycle management: when and how to purge failures that happened in the past.

Basically having an interface for error reporting that works "the same" for Dask and local load case can be tricky. Maybe we could instead have two different interfaces, one for local load and a separate one for Dask. We can start with local only, as this one will be used by Dask implementation anyway.

Problem Description

odc.stac.load output pixel grid can be confusing to users. Unexpected result often look like an error in implementation: an extra pixel column/row, but in fact it is an expected feature of the loading process and has to do with "pixel grid snapping". There is a common expectation that output will start exactly at min(x) of the input data. This is not always the case. Most commonly this manifests when loading data at reduced resolution, example #94.

Another source of confusion is to do with coordinates of the returned xarray. Coordinates are for pixel centers, not for the edge. Common mistake is to assume that image spans [img.x[0], img.x[-1]]. It does not.

• ❌ [img.x[0], img.x[-1]], pixel centers are not pixel edges
• ✅ [img.x[0]-res/2, img.x[-1]+res/2], better but need to know pixel size
• ✅✅ img.odc.geobox.boundingbox.range_x, best

Worth noting that users of stackstac experience similar confusion, and there is a pinned issue https://github.com/gjoseph92/stackstac/issues/152 helping users.

Pixel Grid Snapping

Just because two images have the same resolution, the same projection and overlap, does not mean that individual pixels overlap exactly with each other:

When requested to load two images above how the output pixel grid should look like? Should we pick pixel grid that aligns with orange or with blue image, or something else? We choose something else.

Say images above have pixel size of 100 units, orange one starts at X=130 and blue one at X=690. Image that will be produced in this case will default to starting at X=100 and so both images will be resampled into that common grid. We pick X=100 because this aligns pixel edges to integer multiples of pixel size.

A common assumption users have, is that instead, output image should start at X=130 as this is the left edge of available input data. It's a sensible assumption. Image like that can be produced by odc.stac, but requires specifying anchor= parameter, value of odc.geo.geobox.AnchorEnum.FLOATING, this turns off pixel snapping. This is not well reflected in the docs. Maybe worth adding snap=False or tight=True argument as a more obvious alias for that behaviour.

The advantage of snapping becomes apparent when combining data from different sources: just use the same resolution to load data from different sources and resulting images are guaranteed to have exactly overlapping pixels, so you can combine those arrays without needing further resampling.

References

• https://github.com/opendatacube/odc-geo/issues/59
• https://odc-stac.readthedocs.io/en/latest/_api/odc.stac.load.html?highlight=anchor
• https://odc-geo.readthedocs.io/en/latest/intro-geobox.html
• https://odc-geo.readthedocs.io/en/latest/_api/odc.geo.geobox.GeoBox.from_bbox.html
• #94
• https://github.com/gjoseph92/stackstac/issues/152

Hi all, as far as I can see it is currently not possible to control the chunking of non-spatial dimensions when loading a dataset, right? So in order to have chunks that, for example, extend through the time dimension, I'd need an extra step like this:

ds = odc.stac.load(items=items, bands=bands, chunks={'x': 2048, 'y': 2048})
ds = ds.chunk(chunks={'time':-1})

I've used ODC in the past and expected chunks={'x': 2048, 'y': 2048, 'time':-1} to work. If this is only going to be added sometime in the future, I think there should be a note in the documentation that currently only chunking of spatial dimensions is supported.

odc.stac.stac_load is an alias of odc.stac.load: https://github.com/opendatacube/odc-stac/blob/064c5f04573f8a873ae52ab76d047510ed033caf/odc/stac/init.py#L13. IMO one should be deprecated in favor of the other, and odc.stac.load feels like a better name to me.

odc-stac should accept STAC items that have local paths. The underlying error comes from datacube, so maybe the fix should go there? If so we can transfer this issue over.

Example script

This script uses a MODIS file from the Planetary Computer (staging), but any tif will do. Requires stactools.

import datetime
import odc.stac
import stactools.core.create

bbox = [-107, 39, -106, 40]
file_name = "MYD10A1.A2020335.h10v04.061.2020357064115_NDSI_Snow_Cover.tif"
item = stactools.core.create.item(file_name)  # asset href is set to "MYD10A1.A2020335.h10v04.061.2020357064115_NDSI_Snow_Cover.tif"
item.common_metadata.start_datetime = datetime.datetime(2020, 11, 30)
item.common_metadata.end_datetime = datetime.datetime(2020, 11, 30, 23, 59, 59)
xarray = odc.stac.load([item], bands="data", resolution=500, bbox=bbox)

Expected behavior

It works

Actual behavior

Traceback (most recent call last):
  File "/Users/gadomski/Code/odc-stac/scripts/local_path.py", line 10, in <module>
    xarray = odc.stac.load([item], bands="data", resolution=500, bbox=bbox)
  File "/Users/gadomski/Code/odc-stac/odc/stac/_load.py", line 462, in load
    return dc_load(
  File "/Users/gadomski/Code/odc-stac/odc/stac/_dcload.py", line 191, in dc_load
    return Datacube.load_data(
  File "/Users/gadomski/.virtualenvs/odc-stac/lib/python3.9/site-packages/datacube/api/core.py", line 768, in load_data
    return Datacube._xr_load(sources, geobox, measurements,
  File "/Users/gadomski/.virtualenvs/odc-stac/lib/python3.9/site-packages/datacube/api/core.py", line 701, in _xr_load
    _fuse_measurement(data_slice, datasets, geobox, m,
  File "/Users/gadomski/.virtualenvs/odc-stac/lib/python3.9/site-packages/datacube/api/core.py", line 902, in _fuse_measurement
    src = new_datasource(BandInfo(ds, measurement.name, extra_dim_index=extra_dim_index))
  File "/Users/gadomski/.virtualenvs/odc-stac/lib/python3.9/site-packages/datacube/drivers/readers.py", line 110, in new_datasource
    return source_type(band)
  File "/Users/gadomski/.virtualenvs/odc-stac/lib/python3.9/site-packages/datacube/storage/_rio.py", line 223, in __init__
    filename = _url2rasterio(band.uri, band.format, band.layer)
  File "/Users/gadomski/.virtualenvs/odc-stac/lib/python3.9/site-packages/datacube/storage/_rio.py", line 297, in _url2rasterio
    raise ValueError("Expect either URL or /vsi path")
ValueError: Expect either URL or /vsi path

The raster extension provides scale and offset information. These scales and offsets should be applied to loaded data, if available. Right now, we're having to work around it like this:

data = odc.stac.load(
    items,
    crs="EPSG:3857",
    bands="500m_16_days_NDVI",
    resolution=500,
    bbox=bbox,
)
raster = RasterExtension.ext(items[0].assets["500m_16_days_NDVI"])
data = data["500m_16_days_NDVI"] * raster.bands[0].scale

Full notebook example that includes the workaround in Cell 4 here.