|
| 1 | +"""Regression tests for issue #1597. |
| 2 | +
|
| 3 | +``read_geotiff_dask`` on an integer raster with an in-range nodata |
| 4 | +sentinel used to silently lose the mask when the sentinel only appeared |
| 5 | +in non-first chunks. Per-chunk dtype divergence (uint16 vs float64) |
| 6 | +caused dask concatenation to preallocate from the first chunk's actual |
| 7 | +dtype, casting float64 chunks back to int and converting NaN to 0. |
| 8 | +
|
| 9 | +The fix threads the resolved ``target_dtype`` (the dask graph's |
| 10 | +declared dtype) unconditionally through ``_delayed_read_window`` so |
| 11 | +every chunk lands as float64 regardless of whether its mask hit. |
| 12 | +""" |
| 13 | +from __future__ import annotations |
| 14 | + |
| 15 | +import numpy as np |
| 16 | +import pytest |
| 17 | + |
| 18 | +from xrspatial.geotiff import open_geotiff |
| 19 | +from xrspatial.geotiff._writer import write |
| 20 | + |
| 21 | + |
| 22 | +@pytest.fixture |
| 23 | +def uint16_with_sentinel_only_in_corner(tmp_path): |
| 24 | + """Write a uint16 8x8 TIFF whose nodata sentinel is in the |
| 25 | + bottom-right 2x2 quadrant. With ``chunks=4`` the top-left chunk |
| 26 | + never sees a sentinel and used to keep its uint16 dtype. |
| 27 | + """ |
| 28 | + arr = np.arange(64, dtype=np.uint16).reshape(8, 8) + 1 |
| 29 | + arr[6:8, 6:8] = 65535 |
| 30 | + path = str(tmp_path / 'uint16_corner_sentinel_1597.tif') |
| 31 | + write(arr, path, nodata=65535, compression='none', tiled=False) |
| 32 | + return path, arr |
| 33 | + |
| 34 | + |
| 35 | +def test_eager_promotes_to_float64_and_masks(uint16_with_sentinel_only_in_corner): |
| 36 | + """Baseline: the eager path produces float64 with 4 NaNs.""" |
| 37 | + path, _ = uint16_with_sentinel_only_in_corner |
| 38 | + eager = open_geotiff(path) |
| 39 | + assert eager.dtype == np.float64 |
| 40 | + assert np.isnan(eager.values).sum() == 4 |
| 41 | + assert np.isnan(eager.values[6:8, 6:8]).all() |
| 42 | + |
| 43 | + |
| 44 | +def test_dask_chunks_4_matches_eager(uint16_with_sentinel_only_in_corner): |
| 45 | + """The dask compute result matches the eager path bit-for-bit. |
| 46 | +
|
| 47 | + Before the fix this returned a uint16 array with 0s where the |
| 48 | + sentinel had been, because dask coerced the late-arriving float64 |
| 49 | + chunk back to uint16 at concat time. |
| 50 | + """ |
| 51 | + path, _ = uint16_with_sentinel_only_in_corner |
| 52 | + eager = open_geotiff(path) |
| 53 | + dk = open_geotiff(path, chunks=4) |
| 54 | + assert dk.dtype == np.float64 |
| 55 | + computed = dk.compute() |
| 56 | + assert computed.dtype == np.float64 |
| 57 | + np.testing.assert_array_equal(np.isnan(computed.values), |
| 58 | + np.isnan(eager.values)) |
| 59 | + finite = ~np.isnan(eager.values) |
| 60 | + np.testing.assert_array_equal(computed.values[finite], |
| 61 | + eager.values[finite]) |
| 62 | + |
| 63 | + |
| 64 | +def test_dask_chunks_2_per_chunk_dtype_uniform( |
| 65 | + uint16_with_sentinel_only_in_corner): |
| 66 | + """Every dask chunk returns float64 regardless of mask hit. |
| 67 | +
|
| 68 | + Iterates the delayed blocks and asserts each one computes to |
| 69 | + float64; the regression had the first chunk's actual data come back |
| 70 | + as uint16 because the mask never matched there. |
| 71 | + """ |
| 72 | + path, _ = uint16_with_sentinel_only_in_corner |
| 73 | + dk = open_geotiff(path, chunks=2) |
| 74 | + blocks = dk.data.to_delayed().flatten() |
| 75 | + for i, block in enumerate(blocks): |
| 76 | + chunk = block.compute() |
| 77 | + assert chunk.dtype == np.float64, ( |
| 78 | + f"chunk {i} computed as {chunk.dtype}, expected float64; " |
| 79 | + f"per-chunk dtype divergence is the #1597 regression." |
| 80 | + ) |
| 81 | + |
| 82 | + |
| 83 | +def test_dask_keeps_dtype_for_out_of_range_sentinel(tmp_path): |
| 84 | + """Out-of-range sentinels (uint16 + nodata=-9999) stay uint16. |
| 85 | +
|
| 86 | + The fix should not regress #1581: when the sentinel cannot match |
| 87 | + any pixel, no float64 promotion is needed and the dask path keeps |
| 88 | + the file's native dtype. |
| 89 | + """ |
| 90 | + arr = np.array([[1, 2, 3, 4]] * 4, dtype=np.uint16) |
| 91 | + path = str(tmp_path / 'uint16_out_of_range_1597.tif') |
| 92 | + write(arr, path, nodata=-9999, compression='none', tiled=False) |
| 93 | + |
| 94 | + dk = open_geotiff(path, chunks=2) |
| 95 | + assert dk.dtype == np.uint16 |
| 96 | + result = dk.compute() |
| 97 | + assert result.dtype == np.uint16 |
| 98 | + np.testing.assert_array_equal(result.values, arr) |
| 99 | + |
| 100 | + |
| 101 | +def test_dask_float_input_with_sentinel_in_one_chunk(tmp_path): |
| 102 | + """Float rasters with sentinel in non-first chunk also stay float. |
| 103 | +
|
| 104 | + The float path doesn't promote dtype, but it does in-place NaN |
| 105 | + substitution. Verify the substitution holds for chunks with and |
| 106 | + without the sentinel. |
| 107 | + """ |
| 108 | + arr = np.arange(64, dtype=np.float32).reshape(8, 8) + 1 |
| 109 | + arr[6:8, 6:8] = -9999.0 |
| 110 | + path = str(tmp_path / 'float_corner_sentinel_1597.tif') |
| 111 | + write(arr, path, nodata=-9999, compression='none', tiled=False) |
| 112 | + |
| 113 | + eager = open_geotiff(path) |
| 114 | + dk = open_geotiff(path, chunks=4).compute() |
| 115 | + np.testing.assert_array_equal(np.isnan(dk.values), |
| 116 | + np.isnan(eager.values)) |
0 commit comments