Access to MOD09 data broken? #232

jgomezdans · 2023-07-05T12:36:07Z

jgomezdans
Jul 5, 2023

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.

Answered by TomAugspurger

Jul 5, 2023

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=

View full answer

TomAugspurger · 2023-07-05T13:50:21Z

TomAugspurger
Jul 5, 2023

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

2 replies

jgomezdans Jul 5, 2023
Author

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.

TomAugspurger Jul 5, 2023

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#

Answer selected by jgomezdans

jgomezdans · 2023-07-05T15:10:26Z

jgomezdans
Jul 5, 2023
Author

Thanks for this! I appreciate your fast response!

0 replies

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Access to MOD09 data broken? #232

Uh oh!

{{title}}

Uh oh!

Replies: 2 comments 2 replies

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Select a reply

Uh oh!

Access to MOD09 data broken? #232

Uh oh!

jgomezdans Jul 5, 2023

Replies: 2 comments · 2 replies

Uh oh!

TomAugspurger Jul 5, 2023

Uh oh!

jgomezdans Jul 5, 2023 Author

Uh oh!

TomAugspurger Jul 5, 2023

Uh oh!

jgomezdans Jul 5, 2023 Author

jgomezdans
Jul 5, 2023

Replies: 2 comments 2 replies

TomAugspurger
Jul 5, 2023

jgomezdans Jul 5, 2023
Author

jgomezdans
Jul 5, 2023
Author