Stacking all bands and saving as COG

[guillemc23]

What I'm trying to do is to stack all of the bands in a single COG, but, for me this is not working:

# Create the reader
reader = Reader()

# Open your product
prod = reader.open(product_path, remove_tmp=True)
ok_bands = prod.get_existing_bands()
prod.stack(ok_bands, stack_path="stack_COG.tiff", driver="COG", **{TO_REFLECTANCE: False})

Not sure what is the proper way of doing this, but I would like as well to resample all the bands to the higher resolution one instead of the lower resolution one. 🤔

I'm wondering as well, if there is a way to set the colorintrp to the proper Red, Green and Blue bands 🤔

[remi-braun]

• It's creating a classif GeoTiff ? Can I have the GDALinfo of your output ? I really don't see why it wouldn't work 🤔
• To change the resolution you can pass pixel_size=... in meters

[guillemc23]

• Adding pixel_size = 10 works, you're right!
• What do you mean a classif GeoTIFF? I'm just trying to stack all of the bands in a Sentinel-2 product, handling the resampling and keeping the values as they are so that I can store it in uint16. I am able to create the TIFF with all of the bands. I believe the problem comes when adding the driver="COG", in this case it just won't work. I'll add the full stacktrace tomorrow.

[guillemc23]

This is the complete stack trace. It works without driver="COG":

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File c:\Users\gmcastano\Desktop\EOreader_testbench\.venv\Lib\site-packages\xarray\backends\file_manager.py:211, in CachingFileManager._acquire_with_cache_info(self, needs_lock)
    210 try:
--> 211     file = self._cache[self._key]
    212 except KeyError:

File c:\Users\gmcastano\Desktop\EOreader_testbench\.venv\Lib\site-packages\xarray\backends\lru_cache.py:56, in LRUCache.__getitem__(self, key)
     55 with self._lock:
---> 56     value = self._cache[key]
     57     self._cache.move_to_end(key)

KeyError: [<function open at 0x00000280A066A8E0>, ('S2A_MSIL2A_20240716T100601_N0510_R022_T33TUM_20240716T160254.SAFE\\GRANULE\\L2A_T33TUM_A047352_20240716T101543\\IMG_DATA\\R10m\\T33TUM_20240716T100601_B04_10m.jp2',), 'r', (('driver', 'COG'), ('sharing', False)), '2f164536-6afd-4198-b0b4-c7cb4f4a09f9']

During handling of the above exception, another exception occurred:

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:359, in rasterio._err.exc_wrap_pointer()

CPLE_OpenFailedError: 'S2A_MSIL2A_20240716T100601_N0510_R022_T33TUM_20240716T160254.SAFE\GRANULE\L2A_T33TUM_A047352_20240716T101543\IMG_DATA\R10m\T33TUM_20240716T100601_B04_10m.jp2' not recognized as being in a supported file format.

During handling of the above exception, another exception occurred:

RasterioIOError                           Traceback (most recent call last)
File c:\Users\gmcastano\Desktop\EOreader_testbench\.venv\Lib\site-packages\eoreader\utils.py:186, in read(raster_path, pixel_size, size, resampling, masked, indexes, **kwargs)
    185 warnings.simplefilter("ignore", category=NotGeoreferencedWarning)
--> 186 arr = rasters.read(
    187     raster_path,
    188     resolution=pixel_size,
    189     resampling=resampling,
    190     masked=masked,
    191     indexes=indexes,
    192     size=size,
    193     window=window,
    194     chunks=chunks,
    195     **_prune_keywords(
    196         additional_keywords=["window", "chunks", "resolution"], **kwargs
    197     ),
    198 )
    200 # In EOReader, we don't care about the band coordinate of a band loaded from a stack.
    201 # Overwrite it (don't keep the number "2" if we loaded the second band of the stack)

File c:\Users\gmcastano\Desktop\EOreader_testbench\.venv\Lib\site-packages\sertit\rasters.py:857, in __read__any_raster_to_rio_ds.<locals>.wrapper(any_raster_type, *args, **kwargs)
    856     with rasterio.open(str(any_raster_type)) as ds:
--> 857         out = function(ds, *args, **kwargs)
    858 # Input is a tuple: we consider it's composed of an output of rasterio.read function, a numpy array and a metadata dict

File c:\Users\gmcastano\Desktop\EOreader_testbench\.venv\Lib\site-packages\sertit\rasters.py:984, in read(ds, resolution, size, window, resampling, masked, indexes, chunks, as_type, **kwargs)
    980 # Read data (and load it to discard lock)
    981 with (
    982     xr.set_options(keep_attrs=True),
    983     rioxarray.set_options(export_grid_mapping=False),
--> 984     rioxarray.open_rasterio(
    985         ds,
    986         default_name=path.get_filename(ds.name),
    987         chunks=chunks,
    988         masked=masked,
    989         **kwargs,
    990     ) as xda,
    991 ):
    992     orig_dtype = xda.encoding.get(
    993         "rasterio_dtype", xda.encoding.get("dtype", xda.dtype)
    994     )

File c:\Users\gmcastano\Desktop\EOreader_testbench\.venv\Lib\site-packages\rioxarray\_io.py:1135, in open_rasterio(filename, parse_coordinates, chunks, cache, lock, masked, mask_and_scale, variable, group, default_name, decode_times, decode_timedelta, band_as_variable, **open_kwargs)
   1134     manager = URIManager(file_opener, filename, mode="r", kwargs=open_kwargs)
-> 1135 riods = manager.acquire()
   1136 captured_warnings = rio_warnings.copy()

File c:\Users\gmcastano\Desktop\EOreader_testbench\.venv\Lib\site-packages\xarray\backends\file_manager.py:193, in CachingFileManager.acquire(self, needs_lock)
    179 """Acquire a file object from the manager.
    180 
    181 A new file is only opened if it has expired from the
   (...)
    191     An open file object, as returned by ``opener(*args, **kwargs)``.
    192 """
--> 193 file, _ = self._acquire_with_cache_info(needs_lock)
    194 return file

File c:\Users\gmcastano\Desktop\EOreader_testbench\.venv\Lib\site-packages\xarray\backends\file_manager.py:217, in CachingFileManager._acquire_with_cache_info(self, needs_lock)
    216     kwargs["mode"] = self._mode
--> 217 file = self._opener(*self._args, **kwargs)
    218 if self._mode == "w":
    219     # ensure file doesn't get overridden when opened again

File c:\Users\gmcastano\Desktop\EOreader_testbench\.venv\Lib\site-packages\rasterio\env.py:463, in ensure_env_with_credentials.<locals>.wrapper(*args, **kwds)
    462 with env_ctor(session=session):
--> 463     return f(*args, **kwds)

File c:\Users\gmcastano\Desktop\EOreader_testbench\.venv\Lib\site-packages\rasterio\__init__.py:368, in open(fp, mode, driver, width, height, count, crs, transform, dtype, nodata, sharing, opener, **kwargs)
    367 if mode == "r":
--> 368     dataset = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)
    369 elif mode == "r+":

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

RasterioIOError: 'S2A_MSIL2A_20240716T100601_N0510_R022_T33TUM_20240716T160254.SAFE\GRANULE\L2A_T33TUM_A047352_20240716T101543\IMG_DATA\R10m\T33TUM_20240716T100601_B04_10m.jp2' not recognized as being in a supported file format.

The above exception was the direct cause of the following exception:

InvalidProductError                       Traceback (most recent call last)
Cell In[27], line 11
      8 ok_bands = prod.get_existing_bands()
      9 ok_bands = [RED, GREEN, BLUE, NIR, ALL_CLOUDS]
---> 11 s2_prod = prod.stack(ok_bands, stack_path="stack_4bands_clouds_COG_v2.tiff", pixel_size=10, driver="COG", **{TO_REFLECTANCE: False})
     12 # prod.stack(ok_bands, stack_path="stack.tiff", **{TO_REFLECTANCE: False})
     13 # prod.stack(ok_bands, stack_path="stack_COG.tiff", driver="COG", **{TO_REFLECTANCE: False})

File c:\Users\gmcastano\Desktop\EOreader_testbench\.venv\Lib\site-packages\eoreader\products\product.py:1718, in Product.stack(self, bands, pixel_size, size, stack_path, save_as_int, **kwargs)
   1715 bands = self.to_band(bands)
   1717 # Create the analysis stack
-> 1718 band_xds = self.load(bands, pixel_size=pixel_size, size=size, **kwargs)
   1720 # Stack bands
   1721 if save_as_int:

File c:\Users\gmcastano\Desktop\EOreader_testbench\.venv\Lib\site-packages\eoreader\products\product.py:935, in Product.load(self, bands, pixel_size, size, **kwargs)
    933 # Load bands (only once ! and convert the bands to be loaded to correct format)
    934 unique_bands = misc.unique(bands)
--> 935 band_xds = self._load(unique_bands, pixel_size, size, **kwargs)
    937 # Rename all bands and add attributes
    938 for key, val in band_xds.items():

File c:\Users\gmcastano\Desktop\EOreader_testbench\.venv\Lib\site-packages\eoreader\products\product.py:1033, in Product._load(self, bands, pixel_size, size, **kwargs)
   1031 if unique_bands:
   1032     LOGGER.debug(f"Loading bands {to_str(unique_bands)}")
-> 1033     loaded_bands = self._load_bands(
   1034         unique_bands, pixel_size=pixel_size, size=size, **kwargs
   1035     )
   1037     # Compute index (they conserve the nodata)
   1038     if index_list:
   1039         # Collocate bands before indices to ensure the same size to perform operations between bands

File c:\Users\gmcastano\Desktop\EOreader_testbench\.venv\Lib\site-packages\eoreader\products\optical\s2_product.py:1215, in S2Product._load_bands(self, bands, pixel_size, size, **kwargs)
   1212 band_paths = self.get_band_paths(bands, pixel_size=pixel_size, **kwargs)
   1214 # Open bands and get array (resampled if needed)
-> 1215 band_arrays = self._open_bands(
   1216     band_paths, pixel_size=pixel_size, size=size, **kwargs
   1217 )
   1219 return band_arrays

File c:\Users\gmcastano\Desktop\EOreader_testbench\.venv\Lib\site-packages\eoreader\products\optical\optical_product.py:314, in OpticalProduct._open_bands(self, band_paths, pixel_size, size, **kwargs)
    311 for band, band_path in band_paths.items():
    312     # Read band
    313     LOGGER.debug(f"Read {band.name}")
--> 314     band_arr = self._read_band(
    315         band_path, band=band, pixel_size=pixel_size, size=size, **kwargs
    316     )
    318     if not pixel_size:
    319         pixel_size = band_arr.rio.resolution()[0]

File c:\Users\gmcastano\Desktop\EOreader_testbench\.venv\Lib\site-packages\eoreader\products\optical\s2_product.py:732, in S2Product._read_band(self, band_path, band, pixel_size, size, **kwargs)
    729         raise ex
    731 # Read band
--> 732 return utils.read(
    733     geocoded_path,
    734     pixel_size=pixel_size,
    735     size=size,
    736     resampling=Resampling.bilinear,
    737     **kwargs,
    738 )

File c:\Users\gmcastano\Desktop\EOreader_testbench\.venv\Lib\site-packages\eoreader\utils.py:216, in read(raster_path, pixel_size, size, resampling, masked, indexes, **kwargs)
    210 except errors.RasterioIOError as ex:
    211     if not isinstance(raster_path, str) and (
    212         str(raster_path).endswith("jp2")
    213         or str(raster_path).endswith("tif")
    214         and raster_path.exists()
    215     ):
--> 216         raise InvalidProductError(f"Corrupted file: {raster_path}") from ex
    217     else:
    218         raise

InvalidProductError: Corrupted file: S2A_MSIL2A_20240716T100601_N0510_R022_T33TUM_20240716T160254.SAFE\GRANULE\L2A_T33TUM_A047352_20240716T101543\IMG_DATA\R10m\T33TUM_20240716T100601_B04_10m.jp2

[remi-braun]

I didn't know that driver was a used keyword when reading 🤔
I've purged it. It works for me. Could you try again (you may need to install tifffile library) ?

PS : you could also do:

# Create the reader
reader = Reader()

# Open your product
prod = reader.open(product_path, remove_tmp=True)
ok_bands = prod.get_existing_bands()
stack = prod.stack(ok_bands, **{TO_REFLECTANCE: False})

# Sertit (you may need to install `tifffile` library)
rasters.write(stack, "stack_COG.tiff", driver="COG", dtype="float32")
# or using directly 
stack.rio.to_raster("stack_COG.tiff", driver="COG", dtype="float32")

[guillemc23]

there is definitely something weird with these methods. I've tested the three following methods to save a stack with the following execution times:

stack = prod.stack([RED, GREEN, BLUE, NIR, ALL_CLOUDS], stack_path="s2_stack.tiff", clean_optical="raw", save_as_int=True, **{TO_REFLECTANCE: False}) #1m 13.4s
rasters.write(stack, "s2_stack_COG.tiff", write_cogs_with_dask=True, driver="COG", dtype="int16") #3m 7.8s
stack.rio.to_raster("s2_stack_COG_rio.tiff", driver="COG", dtype="int16") #1m 16.4s

When opening the results in QGIS, the metadata is kept but the band names get messed up when trying to save as COG with these methods, I might be missing something 🤔

On the other hand, I'm not sure why the two COG methods yield different sized files:
s2_stack_COG.tiff ▶️ 1010772 KB
s2_stack_COG_rio.tiff ▶️ 1356501 KB

[remi-braun]

• To rename correctly the bands, you have to amend the long_name attribute in your xarray
• Writing with dask is usually slower but allows out-of-memory computation
• The size of the output file classically uses GDAL and hence depends on its internal configuration. Writing the COG with dask doesn't use GDAL, so the behind-the-scene parameters differ.

[guillemc23]

Thank you, I think it's much more clear now 😄 I think it would be nice to add the option to write directly as a COG but preserving the band metadata like rio_cogeo deoes and maybe adding the ColorInterp metadata, since it is useful for visualization purposes!

[remi-braun]

I will just pass the keywords to what the writer can handle 😉
The good way will always to pass driver="COG" to the stack function.
For the ColorInterp metadata:

• stacking with dask doesn't seem to handle it: https://odc-geo.readthedocs.io/en/latest/_api/odc.geo.cog.save_cog_with_dask.html
• rasterio seems to handle that : https://rasterio.readthedocs.io/en/stable/topics/color.html#writing-colormaps
• rioxarray sadly not Add ability to describe color interpretation corteva/rioxarray#191 / Read colormap and colorinterp from GeoTIFF as an attribute corteva/rioxarray#414