[Bug] Some windows that do overlap an ingested raster in world coordinates are silently dropped during rslearn raster materialization because the pixel-bounds intersection logic can treat them as non-overlapping.

[robmarkcole]

Some windows that do overlap an ingested raster in world coordinates are silently dropped during rslearn raster materialization because the pixel-bounds intersection logic can treat them as non-overlapping.

What happened

We were materializing an ERA5LandMonthlyMeans raster layer into field windows.

The source ERA5 tile was valid and overlapped the target window in world coordinates. A direct reprojection from the tile-store GeoTIFF into the target window grid produced valid temperature values.

However, rslearn materialized the window as all zeros and still marked the layer completed.

Relevant code

• rslearn/dataset/materialize.py, around line 88
• rslearn/tile_stores/default.py, around line 189

The current flow is:

1. get_raster_bounds() warps the source raster to the target CRS and returns normalized pixel bounds as (minx, miny, maxx, maxy).
2. read_raster_window_from_tiles() intersects those source bounds with the window bounds.
3. If the computed intersection has non-positive width or height, the item is skipped.

Suspected bug

get_raster_bounds() normalizes the warped source bounds to (minx, miny, maxx, maxy), but the window bounds are effectively being carried as (left, top, right, bottom) in pixel space.

That mismatch can make a real overlap look empty.

Concrete example

For one failing window:

• Window bounds: (39098, -479739, 39154, -479639)
• Source ERA5 bounds in the same projection: (37833, -480633, 40370, -479787)

The current intersection logic computes:

intersection = (
    max(bounds[0], src_bounds[0]),
    max(bounds[1], src_bounds[1]),
    min(bounds[2], src_bounds[2]),
    min(bounds[3], src_bounds[3]),
)
# -> (39098, -479739, 39154, -479787)

This yields a negative height (-48), so rslearn skips the raster read.

But a direct reprojection from the tile-store GeoTIFF into the same target window grid produced valid data:

valid pixels: 5600
temperature range: about 282.18 K

So the source tile and window do overlap in world coordinates, but the pixel-bounds intersection drops them.

This causes silent data loss:

• materialized rasters can be all zeros even though valid source data exist
• the layer may still be marked completed
• downstream analysis can treat the zero raster as real data

Suggested fix

Normalize window bounds and source bounds to the same pixel-bounds convention before intersecting them.

At minimum, it would help to:

• make the bounds convention explicit and consistent across materialization
• raise an error or strong warning when a source item overlaps in world coordinates but is skipped due to pixel-bounds intersection
• avoid marking the layer completed if every source item was skipped for this reason

Additional context

This was easier to trigger because the ERA5 request bbox was tight enough that CDS returned a very small raster (1 x 3 cells), so edge cases near the tile boundary became more obvious. But the underlying issue appears to be in pixel-bounds handling during materialization, not in the data source itself.

[favyen2]

Sounds like this might be due to distortion when re-projecting from the original projection of the item to the window's projection. Would #541 resolve it?

[robmarkcole]

Test to reproduce the issue:

from __future__ import annotations

import numpy as np
import shapely
from rasterio.crs import CRS
from upath import UPath

from rslearn.data_sources.data_source import Item
from rslearn.dataset.materialize import read_raster_window_from_tiles
from rslearn.tile_stores.default import DefaultTileStore
from rslearn.tile_stores.tile_store import TileStoreWithLayer
from rslearn.utils.geometry import Projection, STGeometry
from rslearn.utils.raster_array import RasterArray


def test_read_raster_window_from_tiles_drops_overlapping_window(tmp_path) -> None:
    bands = [
        "2m-temperature",
        "total-precipitation",
        "surface-solar-radiation-downwards",
    ]

    source_projection = Projection(CRS.from_epsg(4326), 0.1, -0.1)
    source_bounds = (-945, -434, -942, -433)
    source_item = Item(
        "era5_test",
        STGeometry(source_projection, shapely.box(*source_bounds), None),
    )
    source_raster = RasterArray(
        chw_array=np.array(
            [
                [[282.18, 282.18, 282.18]],
                [[1.0, 1.0, 1.0]],
                [[2.0, 2.0, 2.0]],
            ],
            dtype=np.float32,
        )
    )

    tile_store = DefaultTileStore(path_suffix="tiles")
    tile_store.set_dataset_path(UPath(tmp_path))
    layer_store = TileStoreWithLayer(tile_store, "era5")
    layer_store.write_raster(
        source_item,
        bands,
        source_projection,
        source_bounds,
        source_raster,
    )

    # These target bounds come from a real failing USA field window. A direct read in
    # the target grid returns valid data, but the materialize path decides there is no
    # overlap after intersecting pixel bounds.
    target_projection = Projection(CRS.from_epsg(32615), 10, -10)
    target_bounds = (39098, -479739, 39154, -479639)

    direct = layer_store.read_raster(
        source_item,
        bands,
        target_projection,
        target_bounds,
    ).array
    assert int((direct[0] != 0).sum()) > 0

    source_bounds_in_target_projection = layer_store.get_raster_bounds(
        source_item,
        bands,
        target_projection,
    )
    intersection = (
        max(target_bounds[0], source_bounds_in_target_projection[0]),
        max(target_bounds[1], source_bounds_in_target_projection[1]),
        min(target_bounds[2], source_bounds_in_target_projection[2]),
        min(target_bounds[3], source_bounds_in_target_projection[3]),
    )
    assert intersection[3] <= intersection[1]

    materialized = read_raster_window_from_tiles(
        layer_store,
        source_item,
        bands,
        target_projection,
        target_bounds,
        nodata_vals=[0.0, 0.0, 0.0],
        band_dtype=np.float32,
    )
    assert materialized is None

@favyen2 this test fails on that branch -> the pixel-bounds intersection is still (..., height = -48)
and read_raster_window_from_tiles(...) still returns None