Access to MOD09 data broken?

[José Gómez-Dans (jgomezdans)]

Hi,
I am trying to extract surface reflectance from the MODIS MOD09 product. I have basically tried to replicate example here, but I get a weird error when building up the xarray Dataset. I did call the planetary sign-in before running this code, and 'items' has a bunch of STAC objects.

data = odc.stac.load(
    items.values(),
    crs="EPSG:3857",
    bbox=bbox,
    resolution=473,
)

yields the following error:

Aborting load due to failure while reading: https://modiseuwest.blob.core.windows.net/modis-061/MYD09A1/18/03/2021337/MYD09A1.A2021337.h18v03.061.2021346055103.hdf?st=2023-07-04T11%3A36%3A54Z&se=2023-07-12T11%3A36%3A54Z&sp=rl&sv=2021-06-08&sr=c&skoid=c85c15d6-d1ae-42d4-af60-e2ca0f81359b&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2023-07-05T11%3A36%3A53Z&ske=2023-07-12T11%3A36%3A53Z&sks=b&skv=2021-06-08&sig=EzKhiwa7ShxbWdtFBkz3PuVYCoKMHDfNZpEDYhtBvsE%3D:1
---------------------------------------------------------------------------
CPLE_OpenFailedError                      Traceback (most recent call last)
File rasterio/_base.pyx:310, in rasterio._base.DatasetBase.__init__()

File rasterio/_base.pyx:221, in rasterio._base.open_dataset()

File rasterio/_err.pyx:221, in rasterio._err.exc_wrap_pointer()

CPLE_OpenFailedError: '/vsicurl/https://modiseuwest.blob.core.windows.net/modis-061/MYD09A1/18/03/2021337/MYD09A1.A2021337.h18v03.061.2021346055103.hdf?st=2023-07-04T11%3A36%3A54Z&se=2023-07-12T11%3A36%3A54Z&sp=rl&sv=2021-06-08&sr=c&skoid=c85c15d6-d1ae-42d4-af60-e2ca0f81359b&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2023-07-05T11%3A36%3A53Z&ske=2023-07-12T11%3A36%3A53Z&sks=b&skv=2021-06-08&sig=EzKhiwa7ShxbWdtFBkz3PuVYCoKMHDfNZpEDYhtBvsE%3D' not recognized as a supported file format.

During handling of the above exception, another exception occurred:

RasterioIOError                           Traceback (most recent call last)
Cell In[18], line 1
----> 1 data = odc.stac.load(
      2     items.values(),
      3     crs="EPSG:3857",
      4     bbox=bbox,
      5     resolution=473,
      6 )

File /srv/conda/envs/notebook/lib/python3.11/site-packages/odc/stac/_load.py:610, in load(items, bands, groupby, resampling, dtype, chunks, pool, crs, resolution, anchor, geobox, bbox, lon, lat, x, y, like, geopolygon, progress, fail_on_error, stac_cfg, patch_url, preserve_original_order, **kw)
    607 if progress is not None:
    608     _work = progress(SizedIterable(_work, total_tasks))
--> 610 for _ in _work:
    611     pass
    613 return _with_debug_info(ds, tasks=_tasks)

File /srv/conda/envs/notebook/lib/python3.11/site-packages/odc/stac/_utils.py:38, in pmap(func, inputs, pool)
     34 """
     35 Wrapper for ThreadPoolExecutor.map
     36 """
     37 if pool is None:
---> 38     yield from map(func, inputs)
     39     return
     41 if isinstance(pool, int):

File /srv/conda/envs/notebook/lib/python3.11/site-packages/odc/stac/_load.py:601, in load.<locals>._do_one(task)
    595 srcs = [
    596     src
    597     for src in (_parsed[idx].get(band, None) for idx, band in task.srcs)
    598     if src is not None
    599 ]
    600 with rio_env(**_rio_env):
--> 601     _ = _fill_2d_slice(srcs, task.dst_gbox, task.cfg, dst_slice)
    602 t, y, x = task.idx_tyx
    603 return (task.band, t, y, x)

File /srv/conda/envs/notebook/lib/python3.11/site-packages/odc/stac/_load.py:698, in _fill_2d_slice(srcs, dst_gbox, cfg, dst)
    695     return dst
    697 src, *rest = srcs
--> 698 _roi, pix = rio_read(src, cfg, dst_gbox, dst=dst)
    700 for src in rest:
    701     # first valid pixel takes precedence over others
    702     _roi, pix = rio_read(src, cfg, dst_gbox)

File /srv/conda/envs/notebook/lib/python3.11/site-packages/odc/stac/_reader.py:194, in rio_read(src, cfg, dst_geobox, dst)
    188     if cfg.fail_on_error:
    189         log.error(
    190             "Aborting load due to failure while reading: %s:%d",
    191             src.uri,
    192             src.band,
    193         )
--> 194         raise e
    196 # Failed to read, but asked to continue
    197 log.warning("Ignoring read failure while reading: %s:%d", src.uri, src.band)

File /srv/conda/envs/notebook/lib/python3.11/site-packages/odc/stac/_reader.py:186, in rio_read(src, cfg, dst_geobox, dst)
    163 """
    164 Internal read method.
    165 
   (...)
    182 
    183 """
    185 try:
--> 186     return _rio_read(src, cfg, dst_geobox, dst)
    187 except rasterio.errors.RasterioIOError as e:
    188     if cfg.fail_on_error:

File /srv/conda/envs/notebook/lib/python3.11/site-packages/odc/stac/_reader.py:219, in _rio_read(src, cfg, dst_geobox, dst)
    209 def _rio_read(
    210     src: RasterSource,
    211     cfg: RasterLoadParams,
   (...)
    215     # if resampling is `nearest` then ignore sub-pixel translation when deciding
    216     # whether we can just paste source into destination
    217     ttol = 0.9 if cfg.nearest else 0.05
--> 219     with rasterio.open(src.uri, "r", sharing=False) as rdr:
    220         assert isinstance(rdr, rasterio.DatasetReader)
    221         ovr_idx: Optional[int] = None

File /srv/conda/envs/notebook/lib/python3.11/site-packages/rasterio/env.py:451, in ensure_env_with_credentials.<locals>.wrapper(*args, **kwds)
    448     session = DummySession()
    450 with env_ctor(session=session):
--> 451     return f(*args, **kwds)

File /srv/conda/envs/notebook/lib/python3.11/site-packages/rasterio/__init__.py:304, in open(fp, mode, driver, width, height, count, crs, transform, dtype, nodata, sharing, **kwargs)
    301 path = _parse_path(raw_dataset_path)
    303 if mode == "r":
--> 304     dataset = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)
    305 elif mode == "r+":
    306     dataset = get_writer_for_path(path, driver=driver)(
    307         path, mode, driver=driver, sharing=sharing, **kwargs
    308     )

File rasterio/_base.pyx:312, in rasterio._base.DatasetBase.__init__()

RasterioIOError: '/vsicurl/https://modiseuwest.blob.core.windows.net/modis-061/MYD09A1/18/03/2021337/MYD09A1.A2021337.h18v03.061.2021346055103.hdf?st=2023-07-04T11%3A36%3A54Z&se=2023-07-12T11%3A36%3A54Z&sp=rl&sv=2021-06-08&sr=c&skoid=c85c15d6-d1ae-42d4-af60-e2ca0f81359b&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2023-07-05T11%3A36%3A53Z&ske=2023-07-12T11%3A36%3A53Z&sks=b&skv=2021-06-08&sig=EzKhiwa7ShxbWdtFBkz3PuVYCoKMHDfNZpEDYhtBvsE%3D' not recognized as a supported file format.

[Tom Augspurger (TomAugspurger)]

Thanks for the report. Are you able to provide the code for getting items so that I can try reproducing this?

[José Gómez-Dans (jgomezdans)]

Thanks for the quick reply. You can just run the linked example notebook and you'll get the same error. I just changed the coordinates, but it also happens if I run this notebook exactly as it is given.

[Tom Augspurger (TomAugspurger)]

Thanks. I can make some progress with this:

import odc.stac
import planetary_computer
import pystac_client
import shapely.geometry

catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)
location = shapely.geometry.Point(-118.2923, 36.5785)
bbox = location.buffer(2).bounds

items = catalog.search(
    collections=["modis-09A1-061"],
    intersects=location,
    datetime="2021-12-03",
).item_collection()


data = odc.stac.load(
    items,
    crs="EPSG:3857",
    bbox=bbox,
    resolution=1000,
    bands=["sur_refl_b01", "sur_refl_b04", "sur_refl_b03"]
)

data

The key change is the inclusion of the bands= argument there to limit which assets were loaded. odc.stac was having trouble loading the (source) HDF5 file. But that example just uses the derived COGs.

After that, I did run into opendatacube/odc-tools#570. If you need to create rgb imagery, you might you might need to do that manually or use something like https://xarray-spatial.org/reference/_autosummary/xrspatial.multispectral.true_color.html#

[José Gómez-Dans (jgomezdans)]

Thanks for this! I appreciate your fast response!