Incorrect handling of pixel Window when loading spectral indices (misalignment after resampling)

[cmayet]

Hi

I’ve run into an issue when loading a spectral index with EOReader using a pixel-based Window.

Summary

When loading a spectral index (e.g. SCoWI), EOReader:

• Automatically loads all required bands
• Resamples them to a common resolution (if needed)
• Aligns them
• Computes the index

This works as expected when using windows :

• bounds (CRS coordinates)
• geometries

However, when using a rasterio.windows.Window (pixel-based window), the behavior seems incorrect.

Minimal example

import logging
import os
import subprocess

from eodag import EODataAccessGateway, setup_logging
from eoreader.bands import *
from eoreader.env_vars import USE_DASK
from eoreader.reader import Reader
from rasterio.windows import Window
from sertit import logs


os.environ[USE_DASK] = "0"


# eoreader_logger
logger = logging.getLogger("eoreader")
logs.init_logger(logger)

# eodag logger
setup_logging(1)
dag = EODataAccessGateway()
search_results = dag.search(
    provider="cop_dataspace",
    collection="S2_MSI_L1C",
    start="2025-08-26",
    end="2025-08-27",
    **{"grid:code": "MGRS-28PCC"},
)

logger.info(f"Working on product {search_results[0].properties['id']}")


path = dag.download(
    search_results[0],
)


def get_gdalinfo(raster_file_path, pattern):
    gdalinfo = subprocess.run(
        ["gdalinfo", raster_file_path],
        capture_output=True,
        text=True,
    )
    info = next((l for l in gdalinfo.stdout.splitlines() if pattern in l), None)
    return info

# Custom user defined read window (is a rasterio Window instance)
window = Window(3208, 2047, 932, 4322)
bounds = (332080.0, 1736310.0, 341400.0, 1779530.0)
logger.info(f"Rasterio reading window (in pixels): {window}")

index_band = SCoWI
bands = get_needed_bands(index_band)
bands.append(SCoWI)
reader = Reader()
with reader.open(
    path, remove_tmp=True, output_path="/tmp/eoreader_pixel_window/"
) as product:
    ds = product.load(bands, window=window,)

    for band in bands:
        print(f"\n\nBand {band}")
        tmp_file = product.get_band_path(band, pixel_size=product.resolution, window=window)
        transform = ds[band].rio.transform()
        shape = ds[band].rio.shape
        print(f"In memory (DataArray) \n\t da.rio.transform() : {transform} \n\t da.rio.shape {shape}")
        print(f"From gdalinfo : ({tmp_file}")
        print(f"\t {get_gdalinfo(tmp_file, "Size is")}")
        print(f"\t {get_gdalinfo(tmp_file, "Origin =")}")
        print(f"\t {get_gdalinfo(tmp_file, "Pixel Size =")}")

Output

Band SpectralBandNames.BLUE
In memory (DataArray) 
	 da.rio.transform() : | 10.00, 0.00, 332080.00|
| 0.00,-10.00, 1779530.00|
| 0.00, 0.00, 1.00| 
	 da.rio.shape (4322, 932)
From gdalinfo : (/tmp/eoreader_pixel_window/tmp_20250826T114651_S2_T28PCC_L1C_150613/20250826T114651_S2_T28PCC_L1C_150613_BLUE_10m_win1e92321750_nodata.tif
	 Size is 932, 4322
	 Origin = (332080.000000000000000,1779530.000000000000000)
	 Pixel Size = (10.000000000000000,-10.000000000000000)


Band SpectralBandNames.GREEN
In memory (DataArray) 
	 da.rio.transform() : | 10.00, 0.00, 332080.00|
| 0.00,-10.00, 1779530.00|
| 0.00, 0.00, 1.00| 
	 da.rio.shape (4322, 932)
From gdalinfo : (/tmp/eoreader_pixel_window/tmp_20250826T114651_S2_T28PCC_L1C_150613/20250826T114651_S2_T28PCC_L1C_150613_GREEN_10m_win1e92321750_nodata.tif
	 Size is 932, 4322
	 Origin = (332080.000000000000000,1779530.000000000000000)
	 Pixel Size = (10.000000000000000,-10.000000000000000)


Band SpectralBandNames.NIR
In memory (DataArray) 
	 da.rio.transform() : | 10.00, 0.00, 332080.00|
| 0.00,-10.00, 1779530.00|
| 0.00, 0.00, 1.00| 
	 da.rio.shape (4322, 932)
From gdalinfo : (/tmp/eoreader_pixel_window/tmp_20250826T114651_S2_T28PCC_L1C_150613/20250826T114651_S2_T28PCC_L1C_150613_NIR_10m_win1e92321750_nodata.tif
	 Size is 932, 4322
	 Origin = (332080.000000000000000,1779530.000000000000000)
	 Pixel Size = (10.000000000000000,-10.000000000000000)


Band SpectralBandNames.SWIR_1
In memory (DataArray) 
	 da.rio.transform() : | 10.00, 0.00, 332080.00|
| 0.00,-10.00, 1779530.00|
| 0.00, 0.00, 1.00| 
	 da.rio.shape (4322, 932)
From gdalinfo : (/tmp/eoreader_pixel_window/tmp_20250826T114651_S2_T28PCC_L1C_150613/20250826T114651_S2_T28PCC_L1C_150613_SWIR_1_10m_win1e92321750_nodata.tif
	 Size is 1864, 8644
	 Origin = (364160.000000000000000,1759060.000000000000000)
	 Pixel Size = (10.000000000000000,-7.966219342896793)


Band SpectralBandNames.SWIR_2
In memory (DataArray) 
	 da.rio.transform() : | 10.00, 0.00, 332080.00|
| 0.00,-10.00, 1779530.00|
| 0.00, 0.00, 1.00| 
	 da.rio.shape (4322, 932)
From gdalinfo : (/tmp/eoreader_pixel_window/tmp_20250826T114651_S2_T28PCC_L1C_150613/20250826T114651_S2_T28PCC_L1C_150613_SWIR_2_10m_win1e92321750_nodata.tif
	 Size is 1864, 8644
	 Origin = (364160.000000000000000,1759060.000000000000000)
	 Pixel Size = (10.000000000000000,-7.966219342896793)


Band SCoWI
In memory (DataArray) 
	 da.rio.transform() : | 10.00, 0.00, 332080.00|
| 0.00,-10.00, 1779530.00|
| 0.00, 0.00, 1.00| 
	 da.rio.shape (4322, 932)
From gdalinfo : (/tmp/eoreader_pixel_window/tmp_20250826T114651_S2_T28PCC_L1C_150613/20250826T114651_S2_T28PCC_L1C_150613_SCoWI_10m_win1e92321750.tif
	 Size is 932, 4322
	 Origin = (332080.000000000000000,1779530.000000000000000)
	 Pixel Size = (10.000000000000000,-10.000000000000000)

Observed discrepancy (example)

• Bands at native 10m resolution (BLUE, GREEN, NIR):
  • ✔ consistent in memory and on disk
• Bands requiring resampling (SWIR):
  • ❌ different size (e.g. ×2)
  • ❌ different origin
  • ❌ inconsistent pixel size

When opening tmp files in qgis, the bands are not aligned.

I suppose the pixel windows are not handled properly, maybe we should consider that a pixel window is valid only for the bands having the product's native resolution, and compute a modified window for the other bands.

Related issue ?

When running the same script (code above) twice, with remove_tmp=False, the first run is OK, the second run raises

Traceback (most recent call last):
  File "scripts/04_eoread_load_scowi.py", line 60, in <module>
    ds = product.load(bands, window=window)
  File ".../site-packages/eoreader/products/product.py", line 1083, in load
    band_xds = self._load(unique_bands, pixel_size, size, **kwargs)
  File ".../site-packages/eoreader/products/product.py", line 1202, in _load
    loaded_bands = self._collocate_bands(loaded_bands)
  File ".../site-packages/eoreader/products/product.py", line 2140, in _collocate_bands
    bands[band_id] = rasters.collocate(reference, band_arr)
  File ".../site-packages/sertit/rasters.py", line 1515, in collocate
    collocated_xds = _collocate_dataarray(
        reference=reference, other=other, resampling=resampling, **kwargs
    )
  File ".../site-packages/sertit/rasters.py", line 1478, in _collocate_dataarray
    collocated_xda = other.rio.reproject_match(reference, resampling=resampling)
  File ".../site-packages/rioxarray/raster_array.py", line 467, in reproject_match
    transform=match_data_array.rio.transform(recalc=True),
  File ".../site-packages/rioxarray/rioxarray.py", line 438, in transform
    src_left, _, _, src_top = self._unordered_bounds(recalc=recalc)
  File ".../site-packages/rioxarray/rioxarray.py", line 810, in _unordered_bounds
    resolution_x, resolution_y = self.resolution(recalc=recalc)
  File ".../site-packages/rioxarray/rioxarray.py", line 777, in resolution
    left, bottom, right, top = self._internal_bounds()
  File ".../site-packages/rioxarray/rioxarray.py", line 745, in _internal_bounds
    raise NoDataInBounds(
rioxarray.exceptions.NoDataInBounds: Unable to determine bounds from coordinates. Data variable: <BLUE_10m_win..._nodata>

[remi-braun]

Pixel-based windows are a bit untested, and I haven't thought of multi-resolution products! Good catch 💯