diff --git a/xrspatial/geotiff/_vrt.py b/xrspatial/geotiff/_vrt.py index 534d09b2..fc97251a 100644 --- a/xrspatial/geotiff/_vrt.py +++ b/xrspatial/geotiff/_vrt.py @@ -148,6 +148,12 @@ class _Source: # ComplexSource extras scale: float | None = None offset: float | None = None + # GDAL ```` values are ``NearestNeighbour`` + # (default), ``Bilinear``, ``Cubic``, ``CubicSpline``, ``Lanczos``, + # ``Average``, ``Mode``. We parse the value so we can advertise it, + # but only nearest-neighbour is implemented in the placement step + # (issue #1694). Higher-quality resamplers are tracked for follow-up. + resample_alg: str | None = None @dataclass @@ -330,6 +336,7 @@ def parse_vrt(xml_str: str, vrt_dir: str = '.') -> VRTDataset: # ComplexSource extras scale = None offset = None + resample_alg = None if tag == 'ComplexSource': scale_str = _text(src_elem, 'ScaleOffset') offset_str = _text(src_elem, 'ScaleRatio') @@ -338,6 +345,11 @@ def parse_vrt(xml_str: str, vrt_dir: str = '.') -> VRTDataset: scale = float(offset_str) if scale_str: offset = float(scale_str) + # ```` is informational for the placement + # step: read_vrt currently always nearest-neighbours, so + # we only record the value for later honouring. See + # issue #1694. + resample_alg = _text(src_elem, 'ResampleAlg') sources.append(_Source( filename=filename, @@ -347,6 +359,7 @@ def parse_vrt(xml_str: str, vrt_dir: str = '.') -> VRTDataset: nodata=src_nodata, scale=scale, offset=offset, + resample_alg=resample_alg, )) bands.append(_VRTBand( @@ -367,6 +380,69 @@ def parse_vrt(xml_str: str, vrt_dir: str = '.') -> VRTDataset: ) +def _resample_nearest(src_arr: np.ndarray, + out_h: int, out_w: int) -> np.ndarray: + """Nearest-neighbour resample ``src_arr`` to ``(out_h, out_w)``. + + Mirrors GDAL's ``SimpleSource`` semantics: each output pixel samples + the source pixel whose centre is closest in the source grid. The + centre-of-pixel mapping is ``src_idx = (out_idx + 0.5) * src_size / + out_size``; ``floor`` of that, clamped into range, gives the index. + + Used by :func:`read_vrt` to honour ````/```` + scaling. See issue #1694 -- before this helper, a source array was + pasted directly into the destination with no resample step, which + raised a broadcast error for downsampled sources and left holes for + upsampled ones. + + Integer-ratio cases short-circuit to ``np.repeat`` (integer + upsample) or strided slicing (integer downsample) to avoid building + an index array for the common GDAL ``SrcRect`` shapes. + """ + src_h, src_w = src_arr.shape[:2] + # Reject an empty source up front: a SimpleSource with SrcRect + # ``xSize=0`` / ``ySize=0`` (or a windowed read that clamps to an + # empty slice) would otherwise reach the integer-ratio fast paths + # below and divide / modulo by zero on ``src_h`` or ``src_w``. Raise + # an explicit ValueError so the bad VRT surfaces with a clear cause + # instead of an opaque ZeroDivisionError. + if src_h == 0 or src_w == 0: + raise ValueError( + "_resample_nearest received an empty source array; " + "the VRT SourceFilename probably has SrcRect with zero size" + ) + if src_h == out_h and src_w == out_w: + return src_arr + # Fast paths for integer ratios -- common when a VRT downsamples by + # 2x / 4x or upsamples by an integer factor. ``np.repeat`` is a + # cheaper nearest-neighbour upsample than the gather below. + if (out_h % src_h == 0 and out_w % src_w == 0 + and out_h >= src_h and out_w >= src_w): + ry = out_h // src_h + rx = out_w // src_w + return np.repeat(np.repeat(src_arr, ry, axis=0), rx, axis=1) + if (src_h % out_h == 0 and src_w % out_w == 0 + and src_h >= out_h and src_w >= out_w): + sy = src_h // out_h + sx = src_w // out_w + # Sample the centre of each output pixel's source block (offset + # by half a block, floored) so the integer-ratio downsample + # agrees with the general non-integer path for the same shapes. + return src_arr[sy // 2::sy, sx // 2::sx][:out_h, :out_w] + # General case: build per-axis nearest indices and gather. ``floor`` + # plus a clamp matches GDAL's nearest-neighbour rounding for + # non-integer ratios. + y_idx = np.minimum( + np.floor((np.arange(out_h) + 0.5) * src_h / out_h).astype(np.intp), + src_h - 1, + ) + x_idx = np.minimum( + np.floor((np.arange(out_w) + 0.5) * src_w / out_w).astype(np.intp), + src_w - 1, + ) + return src_arr[y_idx[:, None], x_idx[None, :]] + + def read_vrt(vrt_path: str, *, window=None, band: int | None = None, max_pixels: int | None = None) -> tuple[np.ndarray, VRTDataset]: @@ -467,15 +543,38 @@ def read_vrt(vrt_path: str, *, window=None, if clip_r0 >= clip_r1 or clip_c0 >= clip_c1: continue # no overlap - # Map back to source coordinates - # Scale factor: source pixels per destination pixel - scale_y = sr.y_size / dr.y_size if dr.y_size > 0 else 1.0 - scale_x = sr.x_size / dr.x_size if dr.x_size > 0 else 1.0 - - src_r0 = sr.y_off + int((clip_r0 - dst_r0) * scale_y) - src_c0 = sr.x_off + int((clip_c0 - dst_c0) * scale_x) - src_r1 = sr.y_off + int((clip_r1 - dst_r0) * scale_y) - src_c1 = sr.x_off + int((clip_c1 - dst_c0) * scale_x) + # SrcRect / DstRect scaling. When sizes match (the common + # SimpleSource case GDAL writes for tile mosaics), the source + # rect maps 1:1 to the destination rect and we can issue a + # windowed read for only the overlap region -- cheaper than + # decoding the whole source rect. When the sizes differ the + # source band is being upsampled or downsampled into its + # destination cell, so we read the full source rect, apply + # nodata masking, resample to ``(dr.y_size, dr.x_size)`` with + # nearest-neighbour (matching GDAL's SimpleSource semantics), + # and *then* take the clip subwindow. Reading the whole + # source rect first keeps the resample math local to the + # source/destination shapes; trying to resample a windowed + # source slice independently introduces fence-post errors at + # the clip boundary and breaks the integer-ratio fast paths. + # See issue #1694. + needs_resample = (sr.y_size != dr.y_size + or sr.x_size != dr.x_size) + if needs_resample: + # TODO(#1704): when the caller passes a small ``window=`` + # this still reads the full SrcRect. Computing the + # inverse of the nearest-neighbour mapping to read only + # the SrcRect subset that feeds ``window`` would cut the + # decode/memory cost for large SrcRects. + read_r0 = sr.y_off + read_c0 = sr.x_off + read_r1 = sr.y_off + sr.y_size + read_c1 = sr.x_off + sr.x_size + else: + read_r0 = sr.y_off + (clip_r0 - dst_r0) + read_c0 = sr.x_off + (clip_c0 - dst_c0) + read_r1 = sr.y_off + (clip_r1 - dst_r0) + read_c1 = sr.x_off + (clip_c1 - dst_c0) # Read from source file using windowed read. # @@ -496,7 +595,7 @@ def read_vrt(vrt_path: str, *, window=None, try: src_arr, _ = read_to_array( src.filename, - window=(src_r0, src_c0, src_r1, src_c1), + window=(read_r0, read_c0, read_r1, read_c1), band=src.band - 1, # convert 1-based to 0-based ) except ( @@ -529,6 +628,16 @@ def read_vrt(vrt_path: str, *, window=None, # fallback: the earlier ``src.nodata or nodata`` shortcut treated # ``0.0`` as falsy and silently replaced it with the band-level # sentinel (issue #1655). + # + # Apply nodata masking *before* the resample. Nearest-neighbour + # carries the sentinel value through unchanged, which is what + # we want -- the resampled pixel inherits the NaN (or integer + # sentinel) of whichever source pixel it sampled. Resampling + # first and masking after would still work for the float case + # but would force the integer-into-float promotion to allocate + # the resampled-shape buffer before discovering the mask, + # which is a waste when the source rect is much larger than + # the destination rect. src_nodata = src.nodata if src.nodata is not None else nodata if src_nodata is not None and src_arr.dtype.kind == 'f': src_arr = src_arr.copy() @@ -570,6 +679,24 @@ def read_vrt(vrt_path: str, *, window=None, if src.offset is not None and src.offset != 0.0: src_arr = src_arr.astype(np.float64) + src.offset + if needs_resample: + # Resample the full source rect to the destination rect + # shape, then clip to the window subregion. Doing the + # resample on the full rect keeps the nearest-neighbour + # index math in one place; the post-resample slice is a + # plain numpy view. + src_arr = _resample_nearest(src_arr, dr.y_size, dr.x_size) + # Take the clip subwindow out of the resampled array. + # ``dst_r0`` / ``dst_c0`` anchor the resampled array in + # destination coordinates; the clip rect is in + # destination coordinates too, so the subtract gives the + # right offsets. + sub_r0 = clip_r0 - dst_r0 + sub_c0 = clip_c0 - dst_c0 + sub_r1 = clip_r1 - dst_r0 + sub_c1 = clip_c1 - dst_c0 + src_arr = src_arr[sub_r0:sub_r1, sub_c0:sub_c1] + # Place into output out_r0 = clip_r0 - r0 out_c0 = clip_c0 - c0 diff --git a/xrspatial/geotiff/tests/test_vrt_scaled_rects_1694.py b/xrspatial/geotiff/tests/test_vrt_scaled_rects_1694.py new file mode 100644 index 00000000..d607539b --- /dev/null +++ b/xrspatial/geotiff/tests/test_vrt_scaled_rects_1694.py @@ -0,0 +1,332 @@ +"""Regression tests for issue #1694. + +``read_vrt`` did not resample source pixel data when a source band's +```` size differed from its ```` size. Downsampling +raised ``ValueError: could not broadcast input array from shape (S,S) +into shape (D,D)`` and upsampling silently left holes -- only the +top-left ``sr.x_size``/``sr.y_size`` pixels of each destination cell +were written. + +The fix: + +* when ``sr.size != dr.size``, read the full source rect, apply nodata + masking, resample to ``(dr.y_size, dr.x_size)`` with nearest-neighbour + (matching GDAL's SimpleSource semantics), and then clip; +* the same-size case still uses windowed reads as before. +""" +from __future__ import annotations + +import numpy as np +import pytest + +from xrspatial.geotiff._vrt import _resample_nearest, read_vrt +from xrspatial.geotiff._writer import write + + +def _write_vrt(tmp_path, xml: str, name: str = 'test.vrt') -> str: + p = str(tmp_path / name) + with open(p, 'w') as f: + f.write(xml) + return p + + +def test_downsample_4x4_to_2x2_does_not_raise_and_uses_nearest(tmp_path): + """SrcRect 4x4 -> DstRect 2x2: result is (2,2), nearest-neighbour. + + Before the fix the source (4,4) array was assigned directly into the + (2,2) destination slice, raising the broadcast error documented in + issue #1694. + """ + src = np.arange(16, dtype=np.uint16).reshape(4, 4) + src_path = str(tmp_path / 'src.tif') + write(src, src_path, compression='none', tiled=False) + + vrt_xml = f""" + 0.0, 2.0, 0.0, 0.0, 0.0, -2.0 + + + {src_path} + 1 + + + + +""" + vrt_path = _write_vrt(tmp_path, vrt_xml, 'down.vrt') + + result, _ = read_vrt(vrt_path) + + assert result.shape == (2, 2), ( + f"expected (2,2), got {result.shape}; resample step missing." + ) + # Nearest-neighbour with the centre-of-output-pixel rule samples + # source indices floor((i+0.5)*4/2) = floor((i+0.5)*2) -> 1, 3 for + # i=0, 1. So we expect src[1, 1], src[1, 3], src[3, 1], src[3, 3]. + expected = np.array([[src[1, 1], src[1, 3]], + [src[3, 1], src[3, 3]]], dtype=np.uint16) + np.testing.assert_array_equal(result, expected) + + +def test_upsample_2x2_to_4x4_repeats_each_source_pixel(tmp_path): + """SrcRect 2x2 -> DstRect 4x4: each source pixel repeated 2x2. + + Before the fix only the top-left 2x2 of the destination was written + and the rest stayed at the fill value (0 for integer, NaN for + float). + """ + src = np.array([[1, 2], [3, 4]], dtype=np.uint16) + src_path = str(tmp_path / 'src.tif') + write(src, src_path, compression='none', tiled=False) + + vrt_xml = f""" + 0.0, 1.0, 0.0, 0.0, 0.0, -1.0 + + + {src_path} + 1 + + + + +""" + vrt_path = _write_vrt(tmp_path, vrt_xml, 'up.vrt') + + result, _ = read_vrt(vrt_path) + + assert result.shape == (4, 4) + expected = np.array([ + [1, 1, 2, 2], + [1, 1, 2, 2], + [3, 3, 4, 4], + [3, 3, 4, 4], + ], dtype=np.uint16) + np.testing.assert_array_equal(result, expected) + # No holes -- every cell was written. + assert not (result == 0).any(), ( + "upsample left zero-filled cells; resample not propagated." + ) + + +def test_non_integer_scale_3x3_to_2x2_no_holes(tmp_path): + """Non-integer source / destination ratio: covers index-mapping path. + + With src=(3,3) -> dst=(2,2), neither integer-ratio fast path applies. + Confirms the general nearest-neighbour gather produces the correct + shape, no holes, no out-of-bounds writes. + """ + src = np.arange(9, dtype=np.uint16).reshape(3, 3) + src_path = str(tmp_path / 'src.tif') + write(src, src_path, compression='none', tiled=False) + + vrt_xml = f""" + 0.0, 1.5, 0.0, 0.0, 0.0, -1.5 + + + {src_path} + 1 + + + + +""" + vrt_path = _write_vrt(tmp_path, vrt_xml, 'nonint.vrt') + + result, _ = read_vrt(vrt_path) + + assert result.shape == (2, 2) + # Nearest mapping: floor((i+0.5) * 3/2) = floor((i+0.5)*1.5) + # i=0 -> floor(0.75) = 0 + # i=1 -> floor(2.25) = 2 + # So output samples src[0,0], src[0,2], src[2,0], src[2,2]. + expected = np.array([[src[0, 0], src[0, 2]], + [src[2, 0], src[2, 2]]], dtype=np.uint16) + np.testing.assert_array_equal(result, expected) + + +def test_per_band_scale_mix(tmp_path): + """Mixed: band 1 downsampled, band 2 at native resolution. + + Both bands must land in the right places without a broadcast error + and without bleeding band 1's resampled values into band 2. + """ + # Band 1 source: 4x4 -- will be downsampled to 2x2 destination. + band1_src = (np.arange(16, dtype=np.uint16) * 10).reshape(4, 4) + # Band 2 source: 2x2 -- native resolution. + band2_src = np.array([[100, 200], [300, 400]], dtype=np.uint16) + + p1 = str(tmp_path / 'b1.tif') + p2 = str(tmp_path / 'b2.tif') + write(band1_src, p1, compression='none', tiled=False) + write(band2_src, p2, compression='none', tiled=False) + + vrt_xml = f""" + 0.0, 1.0, 0.0, 0.0, 0.0, -1.0 + + + {p1} + 1 + + + + + + + {p2} + 1 + + + + +""" + vrt_path = _write_vrt(tmp_path, vrt_xml, 'mix.vrt') + + result, _ = read_vrt(vrt_path) + + assert result.shape == (2, 2, 2) + # Band 1 nearest-neighbour from 4x4 -> 2x2: src[1,1], src[1,3], src[3,1], src[3,3] + expected_b1 = np.array([[band1_src[1, 1], band1_src[1, 3]], + [band1_src[3, 1], band1_src[3, 3]]], + dtype=np.uint16) + np.testing.assert_array_equal(result[..., 0], expected_b1) + # Band 2 native: untouched. + np.testing.assert_array_equal(result[..., 1], band2_src) + + +def test_window_on_downsampled_source_returns_correct_subwindow(tmp_path): + """``window=(0,0,1,1)`` on a 4x4 -> 2x2 source returns the (0,0) cell. + + The destination cell maps to the source pixel that the resample + routine would sample for that location. Confirms the clip-after- + resample ordering: clipping in source coordinates first (as the old + code effectively did) would feed the wrong source slice into the + resampler. + """ + src = np.arange(16, dtype=np.uint16).reshape(4, 4) + src_path = str(tmp_path / 'src.tif') + write(src, src_path, compression='none', tiled=False) + + vrt_xml = f""" + 0.0, 2.0, 0.0, 0.0, 0.0, -2.0 + + + {src_path} + 1 + + + + +""" + vrt_path = _write_vrt(tmp_path, vrt_xml, 'win.vrt') + + result, _ = read_vrt(vrt_path, window=(0, 0, 1, 1)) + + assert result.shape == (1, 1) + # Full-resample value at destination (0,0) is src[1,1] == 5. + assert result[0, 0] == src[1, 1] + + +def test_nodata_preserved_across_downsample(tmp_path): + """Source sentinel pixels survive the resample as NaN in the result. + + Source is uint16 with sentinel=65535. Pixels at the sampled-from + positions whose values are 65535 must appear as NaN in the float64 + VRT output. + """ + sentinel = np.uint16(65535) + src = np.array([ + [10, 20, 30, 40], + [50, sentinel, 70, sentinel], + [90, 100, 110, 120], + [130, sentinel, 150, sentinel], + ], dtype=np.uint16) + src_path = str(tmp_path / 'src_nd.tif') + write(src, src_path, nodata=int(sentinel), compression='none', + tiled=False) + + # Float64 VRT so the integer-into-float promotion path runs and + # leaves NaN at the sentinel pixels. Use ```` on the source + # so the masking branch fires regardless of band-level . + vrt_xml = f""" + 0.0, 2.0, 0.0, 0.0, 0.0, -2.0 + + -9999 + + {src_path} + 1 + + + 65535 + + +""" + vrt_path = _write_vrt(tmp_path, vrt_xml, 'nd.vrt') + + result, _ = read_vrt(vrt_path) + + assert result.shape == (2, 2) + assert result.dtype == np.float64 + # Nearest sampler picks src[1,1], src[1,3], src[3,1], src[3,3]. + # Of those, src[1,1]=65535, src[1,3]=65535, src[3,1]=65535, + # src[3,3]=65535 -- so every output pixel is the sentinel and must + # be NaN after masking. + assert np.isnan(result).all(), ( + f"sentinel did not survive resample as NaN; got {result!r}" + ) + + +def test_nodata_with_mixed_sentinel_and_valid_pixels(tmp_path): + """Mixed sentinel / valid source -> mixed NaN / valid destination. + + Confirms the mask resamples *with* the data, not against the + pre-resampled source. + """ + sentinel = np.uint16(65535) + # Build a 4x4 source where sample sites (1,1), (1,3), (3,1), (3,3) + # are valid, sentinel, valid, sentinel respectively. + src = np.zeros((4, 4), dtype=np.uint16) + src[1, 1] = 11 + src[1, 3] = sentinel + src[3, 1] = 31 + src[3, 3] = sentinel + src_path = str(tmp_path / 'src_mixed.tif') + write(src, src_path, nodata=int(sentinel), compression='none', + tiled=False) + + vrt_xml = f""" + 0.0, 2.0, 0.0, 0.0, 0.0, -2.0 + + -9999 + + {src_path} + 1 + + + 65535 + + +""" + vrt_path = _write_vrt(tmp_path, vrt_xml, 'nd_mixed.vrt') + + result, _ = read_vrt(vrt_path) + + assert result.shape == (2, 2) + assert result[0, 0] == 11.0 + assert np.isnan(result[0, 1]) + assert result[1, 0] == 31.0 + assert np.isnan(result[1, 1]) + + +@pytest.mark.parametrize('shape', [(0, 5), (5, 0), (0, 0)]) +def test_resample_nearest_rejects_empty_source(shape): + """``_resample_nearest`` raises ValueError on an empty source array. + + A SimpleSource with ``SrcRect xSize=0`` or ``ySize=0`` -- or a + windowed read that clamps to an empty slice -- would otherwise feed + a zero-dim array to the integer-ratio fast paths, which compute + ``out_h % src_h`` and divide by ``src_h``/``src_w`` and so would + raise an opaque ``ZeroDivisionError``. Surface the bad input with + a clear ``ValueError`` instead. + """ + src_arr = np.zeros(shape, dtype=np.float64) + with pytest.raises(ValueError, match='empty source array'): + _resample_nearest(src_arr, 2, 2)