Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 31 additions & 10 deletions xrspatial/geotiff/_backends/vrt.py
Original file line number Diff line number Diff line change
Expand Up @@ -571,15 +571,15 @@ def read_vrt(source: str, *,
_vrt_is_rotated = (
gt is not None and (gt[2] != 0.0 or gt[4] != 0.0)
)
height, width = arr.shape[:2]
if window is not None:
r0 = max(0, window[0])
c0 = max(0, window[1])
coord_window = (r0, c0, r0 + height, c0 + width)
else:
coord_window = None
if gt is not None:
origin_x, res_x, _, origin_y, _, res_y = gt
height, width = arr.shape[:2]
if window is not None:
r0 = max(0, window[0])
c0 = max(0, window[1])
coord_window = (r0, c0, r0 + height, c0 + width)
else:
coord_window = None
# Rotated VRTs emit int64 pixel coords to match the eager
# non-VRT rotated path. Without this gate
# the VRT branch handed back float projected coords while
Expand All @@ -594,7 +594,18 @@ def read_vrt(source: str, *,
has_georef=not _vrt_is_rotated,
)
else:
coords = {}
# No ``<GeoTransform>`` at all: synthesise integer pixel coords
# ``0..N-1`` so a no-georef VRT carries the same x/y arrays the
# non-VRT no-georef read path emits (see ``_coords.py``
# ``has_georef=False`` fallback). Previously this branch handed
# back an empty coord dict, so downstream code that assumes x/y
# exist broke on no-georef VRTs but worked on the equivalent
# plain GeoTIFF.
coords = _coords_from_pixel_geometry(
0.0, 0.0, 1.0, 1.0, height, width,
window=coord_window,
has_georef=False,
)

# Select the per-band nodata sentinel. When a specific band is
# selected, source its nodata from that band's ``<NoDataValue>``
Expand Down Expand Up @@ -1099,13 +1110,12 @@ def _read_vrt_chunked(source, *, window, band, name, chunks, gpu, dtype,
# extent. Mirrors the eager branch in ``read_vrt`` so chunked and
# eager reads share the same x/y arrays.
gt = vrt.geo_transform
coords = {}
_vrt_is_rotated = (
gt is not None and (gt[2] != 0.0 or gt[4] != 0.0)
)
coord_window = (win_r0, win_c0, win_r0 + full_h, win_c0 + full_w)
if gt is not None:
origin_x, res_x, _, origin_y, _, res_y = gt
coord_window = (win_r0, win_c0, win_r0 + full_h, win_c0 + full_w)
# Rotated VRTs emit int64 pixel coords to match the eager
# non-VRT rotated path. Without this gate the
# chunked VRT branch handed back float projected coords on a
Expand All @@ -1117,6 +1127,17 @@ def _read_vrt_chunked(source, *, window, band, name, chunks, gpu, dtype,
window=coord_window,
has_georef=not _vrt_is_rotated,
)
else:
# No ``<GeoTransform>`` at all: synthesise integer pixel coords
# so the chunked no-georef VRT read carries the same x/y arrays
# the eager and non-VRT no-georef paths emit (see ``_coords.py``
# ``has_georef=False`` fallback). Previously this branch left
# ``coords = {}``, dropping x/y on no-georef VRTs.
coords = _coords_from_pixel_geometry(
0.0, 0.0, 1.0, 1.0, full_h, full_w,
window=coord_window,
has_georef=False,
)

# Surface the nodata sentinel for the selected band. The chunked
# path declares ``float64`` up front whenever any selected band has
Expand Down
55 changes: 55 additions & 0 deletions xrspatial/geotiff/tests/vrt/test_parity.py
Original file line number Diff line number Diff line change
Expand Up @@ -743,6 +743,61 @@ def test_vrt_chunked_none_matches_dask(tmp_path):
assert tiff_attrs == vrt_attrs


def test_vrt_eager_none_synthesizes_pixel_coords(tmp_path):
"""Issue #2818: a no-``<GeoTransform>`` VRT read eagerly must
synthesise integer x/y pixel coords, matching the non-VRT
no-georef read instead of dropping coords entirely."""
tiff, vrt = _make_none_pair(tmp_path, 'none_coords_2818')
ref = open_geotiff(tiff)
vrt_da = read_vrt(vrt)
assert ref.attrs['georef_status'] == GEOREF_STATUS_NONE
assert vrt_da.attrs['georef_status'] == GEOREF_STATUS_NONE
assert 'x' in vrt_da.coords and 'y' in vrt_da.coords
np.testing.assert_array_equal(
_coord_view(vrt_da, 'x'), _coord_view(ref, 'x'))
np.testing.assert_array_equal(
_coord_view(vrt_da, 'y'), _coord_view(ref, 'y'))
assert vrt_da.coords['x'].dtype == ref.coords['x'].dtype
assert vrt_da.coords['y'].dtype == ref.coords['y'].dtype


def test_vrt_chunked_none_synthesizes_pixel_coords(tmp_path):
"""Issue #2818: a no-``<GeoTransform>`` VRT read chunked must
synthesise the same integer x/y pixel coords as the non-VRT dask
no-georef read rather than dropping coords entirely."""
tiff, vrt = _make_none_pair(tmp_path, 'none_coords_chunked_2818')
ref = read_geotiff_dask(tiff, chunks=2)
vrt_da = read_vrt(vrt, chunks=2)
assert ref.attrs['georef_status'] == GEOREF_STATUS_NONE
assert vrt_da.attrs['georef_status'] == GEOREF_STATUS_NONE
assert 'x' in vrt_da.coords and 'y' in vrt_da.coords
np.testing.assert_array_equal(
_coord_view(vrt_da, 'x'), _coord_view(ref, 'x'))
np.testing.assert_array_equal(
_coord_view(vrt_da, 'y'), _coord_view(ref, 'y'))
assert vrt_da.coords['x'].dtype == ref.coords['x'].dtype
assert vrt_da.coords['y'].dtype == ref.coords['y'].dtype


def test_vrt_none_windowed_synthesizes_offset_pixel_coords(tmp_path):
"""Issue #2818: a windowed no-georef VRT read shifts the synthesised
integer coords to the window offset, matching the non-VRT windowed
read, on both the eager and chunked paths."""
tiff, vrt = _make_none_pair(tmp_path, 'none_coords_win_2818')
window = (1, 2, 4, 4)
ref = open_geotiff(tiff, window=window)
eager = read_vrt(vrt, window=window)
chunked = read_vrt(vrt, window=window, chunks=2)
for label, actual in (('eager', eager), ('chunked', chunked)):
assert 'x' in actual.coords and 'y' in actual.coords, label
np.testing.assert_array_equal(
_coord_view(actual, 'x'), _coord_view(ref, 'x'),
err_msg=label)
np.testing.assert_array_equal(
_coord_view(actual, 'y'), _coord_view(ref, 'y'),
err_msg=label)


def test_vrt_chunked_rotated_dropped(tmp_path):
_, vrt = _make_rotated_pair(tmp_path, 'rot_chunked_2180')
attrs = dict(read_vrt(vrt, allow_rotated=True, chunks=2).attrs)
Expand Down
Loading