diff --git a/xrspatial/geotiff/_backends/dask.py b/xrspatial/geotiff/_backends/dask.py
index 0cc4ab6bf..59a085cbc 100644
--- a/xrspatial/geotiff/_backends/dask.py
+++ b/xrspatial/geotiff/_backends/dask.py
@@ -40,7 +40,7 @@ def _read_geotiff_dask(source: str, *,
                        allow_experimental_codecs: bool = False,
                        allow_internal_only_jpeg: bool = False,
                        band_nodata: str | None = None,
-                       mask_nodata: bool = True,
+                       mask_nodata: bool = False,
                        mask_and_scale: bool = False,
                        parse_coordinates: bool = True) -> xr.DataArray:
     """Read a GeoTIFF as a dask-backed DataArray for out-of-core processing.
@@ -99,14 +99,15 @@ def _read_geotiff_dask(source: str, *,
         mixed-band-metadata check. Forwarded
         verbatim to ``_read_vrt`` when the source is a ``.vrt`` file.
         Passing it with a non-VRT GeoTIFF source raises ``ValueError``.
-    mask_nodata : bool, default True
+    mask_nodata : bool, default False
         [stable] If True, replace the nodata sentinel with NaN per
         chunk (integer rasters get promoted to ``float64``). If False,
         skip the sentinel-to-NaN step so the source dtype survives.
         The raw sentinel is still carried on ``attrs['nodata']``
-        either way. Pass ``mask_nodata=False`` together with
-        ``dtype=<integer>`` to keep an integer source dtype; the
-        default promotes to ``float64`` and the cast then raises.
+        either way. The default matches ``open_geotiff`` (unmasked,
+        rioxarray-compatible): a bare ``_read_geotiff_dask(path)`` keeps
+        the source dtype and leaves the sentinel pixels untouched. Pass
+        ``mask_nodata=True`` to restore the promote-to-NaN behaviour.
     mask_and_scale : bool, default False
         [advanced] If True, apply the source's GDAL ``SCALE`` / ``OFFSET``
         (``data * scale + offset``) lazily on the assembled dask array and
diff --git a/xrspatial/geotiff/_backends/gpu.py b/xrspatial/geotiff/_backends/gpu.py
index d229271c2..9283fdec6 100644
--- a/xrspatial/geotiff/_backends/gpu.py
+++ b/xrspatial/geotiff/_backends/gpu.py
@@ -72,7 +72,7 @@ def _read_geotiff_gpu(source: str, *,
                       allow_experimental_codecs: bool = False,
                       allow_internal_only_jpeg: bool = False,
                       band_nodata: str | None = None,
-                      mask_nodata: bool = True,
+                      mask_nodata: bool = False,
                       gpu: str = _GPU_DEPRECATED_SENTINEL,
                       ) -> xr.DataArray:
     """Read a GeoTIFF with GPU-accelerated decompression via Numba CUDA.
@@ -180,14 +180,15 @@ def _read_geotiff_gpu(source: str, *,
         with values ``'auto'`` / ``'strict'`` and was easy to confuse
         with the boolean ``gpu=`` kwarg on ``open_geotiff`` /
         ``to_geotiff`` / ``_read_vrt``.
-    mask_nodata : bool, default True
+    mask_nodata : bool, default False
         [experimental] If True, replace the nodata sentinel with NaN
         (integer rasters get promoted to ``float64`` first). If False,
         keep the source dtype and leave the raw sentinel in the data.
-        ``attrs['nodata']`` carries the sentinel either way. Pass
-        ``mask_nodata=False`` together with ``dtype=<integer>`` to
-        preserve an integer source dtype on a file with a matching
-        sentinel.
+        ``attrs['nodata']`` carries the sentinel either way. The default
+        matches ``open_geotiff`` (unmasked, rioxarray-compatible): a bare
+        ``_read_geotiff_gpu(path)`` keeps the source dtype and the raw
+        sentinel. Pass ``mask_nodata=True`` to restore the
+        promote-to-NaN behaviour.
     allow_rotated : bool, default False
         [experimental] Read-side opt-in for rotated / sheared
         ``ModelTransformationTag`` files. Forwarded through both GPU
diff --git a/xrspatial/geotiff/_backends/vrt.py b/xrspatial/geotiff/_backends/vrt.py
index 704e4eb67..07b24ce40 100644
--- a/xrspatial/geotiff/_backends/vrt.py
+++ b/xrspatial/geotiff/_backends/vrt.py
@@ -131,7 +131,7 @@ def _read_vrt(source: str, *,
               allow_experimental_codecs: bool = False,
               allow_internal_only_jpeg: bool = False,
               band_nodata: str | None = None,
-              mask_nodata: bool = True) -> xr.DataArray:
+              mask_nodata: bool = False) -> xr.DataArray:
     """Read a GDAL Virtual Raster Table (.vrt) into an xarray.DataArray.
 
     Release-contract tier (see
@@ -233,14 +233,15 @@ def _read_vrt(source: str, *,
         to keep the legacy flatten-to-band-0 semantics explicitly. Any
         other value raises ``ValueError`` at the boundary so typos
         surface up front instead of degrading silently into strict mode.
-    mask_nodata : bool, default True
+    mask_nodata : bool, default False
         [advanced] If True, run the integer-sentinel-to-NaN promotion
         on the assembled mosaic. If False, skip it and keep the source
         dtype with the raw sentinel still in the data. ``attrs['nodata']``
-        carries the sentinel either way. Pass ``mask_nodata=False``
-        together with ``dtype=<integer>`` when you need to preserve an
-        integer source dtype on a VRT whose declared sentinel matches
-        actual pixels. Float source bands are NaN-aware
+        carries the sentinel either way. The default matches
+        ``open_geotiff`` (unmasked, rioxarray-compatible): a bare
+        ``_read_vrt(path)`` keeps the source dtype and the raw sentinel.
+        Pass ``mask_nodata=True`` to restore the promote-to-NaN
+        behaviour on integer mosaics. Float source bands are NaN-aware
         by virtue of how the internal reader handles them, so this kwarg
         is most useful for integer-dtype mosaics.
     allow_rotated : bool, default False
diff --git a/xrspatial/geotiff/tests/gpu/test_codec.py b/xrspatial/geotiff/tests/gpu/test_codec.py
index 9a33649d9..7284da679 100644
--- a/xrspatial/geotiff/tests/gpu/test_codec.py
+++ b/xrspatial/geotiff/tests/gpu/test_codec.py
@@ -1172,8 +1172,11 @@ def _read_cpu_gpu_lerc(path):
     from xrspatial.geotiff._reader import read_to_array
 
     cpu, _geo = read_to_array(path, allow_experimental_codecs=True)
+    # The backend default is unmasked (#2976); request masking so the
+    # GPU read promotes the sentinel to NaN to match the CPU mask.
     gpu_da = _read_geotiff_gpu(
         path, gpu='strict', allow_experimental_codecs=True,
+        mask_nodata=True,
     )
     gpu_host = gpu_da.data.get()
     return cpu, gpu_host
diff --git a/xrspatial/geotiff/tests/gpu/test_kernels_and_kwargs.py b/xrspatial/geotiff/tests/gpu/test_kernels_and_kwargs.py
index ee26d3810..d45430f07 100644
--- a/xrspatial/geotiff/tests/gpu/test_kernels_and_kwargs.py
+++ b/xrspatial/geotiff/tests/gpu/test_kernels_and_kwargs.py
@@ -1643,16 +1643,18 @@ def invalid_pred(a):
     write(arr, path, compression="lerc", tiled=True, tile_size=8,
           nodata=-9999.0)
 
+    # The backend default is unmasked (#2976); request masking on both
+    # the eager and chunked GPU reads so the sentinel promotes to NaN.
     eager = _read_geotiff_gpu(
         path, on_gpu_failure='strict',
-        allow_experimental_codecs=True,
+        allow_experimental_codecs=True, mask_nodata=True,
     ).data.get()
 
     ifd, geo_info, header = _parse_for_gds_1896(path)
     chunked_da = _read_geotiff_gpu_chunked_gds(
         path, ifd, geo_info, header,
         dtype=None, chunks=4, window=None, band=None,
-        name=None, max_pixels=None,
+        name=None, max_pixels=None, mask_nodata=True,
     )
     chunked = chunked_da.data.compute().get()
 
@@ -2459,19 +2461,20 @@ def test_read_geotiff_gpu_mask_nodata_false_preserves_uint16_2052(
 
 
 @_gpu_only
-def test_read_geotiff_gpu_default_mask_nodata_true_still_promotes_2052(
+def test_read_geotiff_gpu_default_unmasked_keeps_uint16_2976(
         uint16_with_matching_sentinel_2052):
-    """The GPU default is unchanged: ``mask_nodata=True`` promotes."""
+    """The GPU default matches ``open_geotiff`` (unmasked): a bare
+    ``_read_geotiff_gpu(path)`` keeps the source dtype and sentinel."""
     import cupy
 
     from xrspatial.geotiff import _read_geotiff_gpu
 
-    path, _ = uint16_with_matching_sentinel_2052
+    path, arr = uint16_with_matching_sentinel_2052
     da = _read_geotiff_gpu(path)
 
-    assert da.dtype == np.float64
-    nan_count = int(cupy.isnan(da.data).sum().get())
-    assert nan_count == len(_SENTINEL_POS_2052)
+    assert da.dtype == np.uint16
+    assert int(cupy.isnan(da.data.astype(cupy.float64)).sum().get()) == 0
+    np.testing.assert_array_equal(da.data.get(), arr)
 
 
 @_gpu_only
@@ -2519,20 +2522,23 @@ def test_read_geotiff_gpu_dask_mask_nodata_false_preserves_uint16_2052(
 
 
 @_gpu_only
-def test_read_geotiff_gpu_dask_default_mask_nodata_true_still_promotes_2052(
+def test_read_geotiff_gpu_dask_default_unmasked_keeps_uint16_2976(
         uint16_with_matching_sentinel_2052):
-    """The dask+GPU default still promotes the graph dtype to float64."""
+    """The dask+GPU default matches ``open_geotiff`` (unmasked): the
+    graph dtype stays uint16 and the sentinel survives."""
     import cupy
 
     from xrspatial.geotiff import _read_geotiff_gpu
 
-    path, _ = uint16_with_matching_sentinel_2052
+    path, arr = uint16_with_matching_sentinel_2052
     da = _read_geotiff_gpu(path, chunks=2)
 
-    assert da.dtype == np.float64
+    assert da.dtype == np.uint16
     computed = da.compute()
-    nan_count = int(cupy.isnan(computed.data).sum().get())
-    assert nan_count == len(_SENTINEL_POS_2052)
+    assert computed.dtype == np.uint16
+    assert int(cupy.isnan(
+        computed.data.astype(cupy.float64)).sum().get()) == 0
+    np.testing.assert_array_equal(computed.data.get(), arr)
 
 
 @_gpu_only
@@ -2565,16 +2571,18 @@ def test_read_vrt_mask_nodata_false_preserves_uint16_2052(
     assert da.attrs["nodata"] == 0
 
 
-def test_read_vrt_default_mask_nodata_true_still_promotes_2052(
+def test_read_vrt_default_unmasked_keeps_uint16_2976(
         uint16_vrt_with_matching_sentinel_2052):
-    """The VRT default unchanged: ``mask_nodata=True`` promotes."""
+    """The VRT default matches ``open_geotiff`` (unmasked): a bare
+    ``_read_vrt(path)`` keeps the source dtype and sentinel."""
     from xrspatial.geotiff import _read_vrt
 
-    vrt_path, _ = uint16_vrt_with_matching_sentinel_2052
+    vrt_path, arr = uint16_vrt_with_matching_sentinel_2052
     da = _read_vrt(vrt_path)
 
-    assert da.dtype == np.float64
-    assert int(np.isnan(da.values).sum()) == len(_SENTINEL_POS_2052)
+    assert da.dtype == np.uint16
+    assert int(np.isnan(np.asarray(da.values, dtype=float)).sum()) == 0
+    np.testing.assert_array_equal(np.asarray(da.values), arr)
 
 
 def test_open_geotiff_vrt_mask_nodata_false_threads_through_2052(
@@ -2617,17 +2625,20 @@ def test_read_vrt_chunked_mask_nodata_false_preserves_uint16_2052(
     np.testing.assert_array_equal(np.asarray(computed.values), arr)
 
 
-def test_read_vrt_chunked_default_mask_nodata_true_still_promotes_2052(
+def test_read_vrt_chunked_default_unmasked_keeps_uint16_2976(
         uint16_vrt_with_matching_sentinel_2052):
-    """The chunked-VRT default still promotes to float64."""
+    """The chunked-VRT default matches ``open_geotiff`` (unmasked): the
+    graph dtype stays uint16 and the sentinel survives."""
     from xrspatial.geotiff import _read_vrt
 
-    vrt_path, _ = uint16_vrt_with_matching_sentinel_2052
+    vrt_path, arr = uint16_vrt_with_matching_sentinel_2052
     da = _read_vrt(vrt_path, chunks=2)
 
-    assert da.dtype == np.float64
+    assert da.dtype == np.uint16
     computed = da.compute()
-    assert int(np.isnan(computed.values).sum()) == len(_SENTINEL_POS_2052)
+    assert computed.dtype == np.uint16
+    assert int(np.isnan(np.asarray(computed.values, dtype=float)).sum()) == 0
+    np.testing.assert_array_equal(np.asarray(computed.values), arr)
 
 
 def test_open_geotiff_vrt_chunked_mask_nodata_false_threads_through_2052(
diff --git a/xrspatial/geotiff/tests/gpu/test_reader.py b/xrspatial/geotiff/tests/gpu/test_reader.py
index 506a1c18f..b75875aaf 100644
--- a/xrspatial/geotiff/tests/gpu/test_reader.py
+++ b/xrspatial/geotiff/tests/gpu/test_reader.py
@@ -1472,12 +1472,13 @@ def test_chunked_gpu_dtype_matches_cpu_dask_1909(
     from xrspatial.geotiff import _read_geotiff_dask
     from xrspatial.geotiff._backends.gpu import _read_geotiff_gpu_chunked_gds
 
-    cpu = _read_geotiff_dask(uint16_no_sentinel_path_1909, chunks=4)
+    cpu = _read_geotiff_dask(
+        uint16_no_sentinel_path_1909, chunks=4, mask_nodata=True)
     ifd, geo_info, header = _parse_for_gds_1909(uint16_no_sentinel_path_1909)
     gpu = _read_geotiff_gpu_chunked_gds(
         uint16_no_sentinel_path_1909, ifd, geo_info, header,
         dtype=None, chunks=4, window=None, band=None,
-        name=None, max_pixels=None,
+        name=None, max_pixels=None, mask_nodata=True,
     )
     assert cpu.data.dtype == gpu.data.dtype, (
         f"CPU dask declared {cpu.data.dtype} but GPU dask declared "
diff --git a/xrspatial/geotiff/tests/parity/test_backend_matrix.py b/xrspatial/geotiff/tests/parity/test_backend_matrix.py
index 4a19ed2d5..3bdb52590 100644
--- a/xrspatial/geotiff/tests/parity/test_backend_matrix.py
+++ b/xrspatial/geotiff/tests/parity/test_backend_matrix.py
@@ -1849,7 +1849,7 @@ def _ap_open_vrt(path, meta):
     )
     with open(vrt_path, 'w') as f:
         f.write(xml)
-    return _read_vrt(vrt_path)
+    return _read_vrt(vrt_path, mask_nodata=True)
 
 
 _AP_BACKENDS = (
diff --git a/xrspatial/geotiff/tests/parity/test_finalization.py b/xrspatial/geotiff/tests/parity/test_finalization.py
index fad11ba4f..95599b043 100644
--- a/xrspatial/geotiff/tests/parity/test_finalization.py
+++ b/xrspatial/geotiff/tests/parity/test_finalization.py
@@ -1030,7 +1030,7 @@ def test_dtype_cast_absent_without_caller_dtype(tmp_path, opener):
     when masking auto-promotes the graph dtype to float64."""
     path = str(tmp_path / "tmp_2178_no_cast.tif")
     _make_int_with_nodata_tiff(path)
-    out = opener(path)
+    out = opener(path, mask_nodata=True)
     # Masking promoted the int source to float64 on the graph dtype,
     # but the caller did not ask for a cast.
     assert out.dtype == np.float64
diff --git a/xrspatial/geotiff/tests/read/test_nodata.py b/xrspatial/geotiff/tests/read/test_nodata.py
index d56bc6c0a..ee38b902b 100644
--- a/xrspatial/geotiff/tests/read/test_nodata.py
+++ b/xrspatial/geotiff/tests/read/test_nodata.py
@@ -484,7 +484,7 @@ def test_dask_mask_nodata_true_reports_true(tmp_path):
     path = str(tmp_path / "tmp_2092_dask_masked.tif")
     _make_float_raster_with_nodata_2092(path)
 
-    out = _read_geotiff_dask(path, chunks=2)
+    out = _read_geotiff_dask(path, chunks=2, mask_nodata=True)
     assert out.attrs.get('masked_nodata') is True
     computed = out.values
     assert np.isnan(computed).sum() == 2
@@ -622,7 +622,7 @@ def test_gpu_mask_nodata_true_reports_true(tmp_path):
     path = str(tmp_path / "tmp_2092_gpu_masked.tif")
     _make_float_raster_with_nodata_2092(path)
 
-    out = _read_geotiff_gpu(path)
+    out = _read_geotiff_gpu(path, mask_nodata=True)
     assert out.attrs.get('masked_nodata') is True
 
 
@@ -782,7 +782,7 @@ def test_dask_leaves_pixels_present_unset(tmp_path):
     ``nodata_pixels_present`` stays unset by design."""
     path = str(tmp_path / "tmp_2135_dask_present.tif")
     _make_float_raster_2135(path)
-    out = _read_geotiff_dask(path, chunks=2)
+    out = _read_geotiff_dask(path, chunks=2, mask_nodata=True)
     assert out.attrs.get('masked_nodata') is True
     assert 'nodata_pixels_present' not in out.attrs
 
@@ -892,7 +892,7 @@ def test_gpu_float_sentinel_present_masked(tmp_path):
 
     path = str(tmp_path / "tmp_2135_gpu_float_present.tif")
     _make_float_raster_2135(path)
-    out = _read_geotiff_gpu(path)
+    out = _read_geotiff_gpu(path, mask_nodata=True)
     assert out.attrs.get('masked_nodata') is True
     assert out.attrs.get('nodata_pixels_present') is True
 
@@ -1913,7 +1913,7 @@ def test_eager_masks_int_sentinel_to_nan(self, int_tif):
 
     def test_dask_matches_eager(self, int_tif):
         eager = open_geotiff(int_tif, masked=True)
-        lazy = _read_geotiff_dask(int_tif, chunks=2)
+        lazy = _read_geotiff_dask(int_tif, chunks=2, mask_nodata=True)
         # Same on-disk sentinel propagated.
         assert lazy.attrs["nodata"] == eager.attrs["nodata"]
         # Same lifecycle decision: dask graph promoted to float64.
@@ -1927,7 +1927,7 @@ def test_gpu_matches_eager(self, int_tif):
         from xrspatial.geotiff import _read_geotiff_gpu
 
         eager = open_geotiff(int_tif, masked=True)
-        gpu = _read_geotiff_gpu(int_tif)
+        gpu = _read_geotiff_gpu(int_tif, mask_nodata=True)
         assert gpu.attrs["nodata"] == eager.attrs["nodata"]
         np.testing.assert_array_equal(
             np.isnan(eager.data), np.isnan(gpu.data.get()),
@@ -1944,7 +1944,7 @@ def test_eager(self, float_tif):
 
     def test_dask(self, float_tif):
         eager = open_geotiff(float_tif, masked=True)
-        lazy = _read_geotiff_dask(float_tif, chunks=2)
+        lazy = _read_geotiff_dask(float_tif, chunks=2, mask_nodata=True)
         np.testing.assert_array_equal(
             np.isnan(eager.data), np.isnan(lazy.compute().data),
         )
@@ -1955,7 +1955,7 @@ def test_gpu(self, float_tif):
         from xrspatial.geotiff import _read_geotiff_gpu
 
         eager = open_geotiff(float_tif, masked=True)
-        gpu = _read_geotiff_gpu(float_tif)
+        gpu = _read_geotiff_gpu(float_tif, mask_nodata=True)
         np.testing.assert_array_equal(
             np.isnan(eager.data), np.isnan(gpu.data.get()),
         )
@@ -2043,7 +2043,7 @@ def test_dask_matches_eager(self, tmp_path):
         path = str(tmp_path / "miw_dask_2226.tif")
         self._build(path)
         eager = open_geotiff(path, masked=True)
-        lazy = _read_geotiff_dask(path, chunks=2)
+        lazy = _read_geotiff_dask(path, chunks=2, mask_nodata=True)
         np.testing.assert_array_equal(
             np.isnan(eager.data), np.isnan(lazy.compute().data),
         )
@@ -2158,7 +2158,7 @@ def simple_vrt(tmp_path):
 
 class TestVRTEagerParity:
     def test_vrt_masks_sentinel_to_nan(self, simple_vrt):
-        da = _read_vrt(simple_vrt)
+        da = _read_vrt(simple_vrt, mask_nodata=True)
         assert da.dtype.kind == "f"
         assert np.isnan(da.data[0, 2])
 
@@ -2472,7 +2472,7 @@ class TestDaskNumpy:
     def test_float_source_with_sentinel(self, tmp_path):
         path = str(tmp_path / "tnss1988_dask_float_sentinel.tif")
         _write_float_tiff_1988(path, with_sentinel=True)
-        da = _read_geotiff_dask(path, chunks=2)
+        da = _read_geotiff_dask(path, chunks=2, mask_nodata=True)
         assert da.attrs["nodata"] == _SENTINEL_1988
         assert da.attrs["masked_nodata"] is True
 
@@ -2494,7 +2494,7 @@ def test_int_source_with_in_range_sentinel(self, tmp_path):
         """
         path = str(tmp_path / "tnss1988_dask_int_in_range.tif")
         _write_int_tiff_1988(path, with_sentinel_hit=False)
-        da = _read_geotiff_dask(path, chunks=2)
+        da = _read_geotiff_dask(path, chunks=2, mask_nodata=True)
         assert da.attrs["nodata"] == 65535
         assert da.dtype.kind == "f"
         assert da.attrs["masked_nodata"] is True
@@ -2561,7 +2561,7 @@ def test_float32_vrt_int_source_with_hit(self, tmp_path):
         )
         vrt = _build_vrt_1988(tmp_path, src, "Float32", 65535,
                               filename="tnss1988_vrt_hit.vrt")
-        r = _read_vrt(vrt)
+        r = _read_vrt(vrt, mask_nodata=True)
         assert r.attrs["nodata"] == 65535.0
         assert r.dtype.kind == "f"
         assert r.attrs["masked_nodata"] is True
@@ -2622,7 +2622,7 @@ def test_chunked_int_source_in_range_sentinel(self, tmp_path):
         )
         vrt = _build_vrt_1988(tmp_path, src, "UInt16", 65535,
                               filename="tnss1988_vrt_chunked.vrt")
-        r = _read_vrt(vrt, chunks=2)
+        r = _read_vrt(vrt, chunks=2, mask_nodata=True)
         assert r.attrs["nodata"] == 65535.0
         # Chunked path promotes to float64 declared dtype.
         assert r.dtype == np.float64
@@ -2638,7 +2638,7 @@ def test_int_source_with_hit(self, tmp_path):
         from xrspatial.geotiff import _read_geotiff_gpu
         path = str(tmp_path / "tnss1988_gpu_int_hit.tif")
         _write_int_tiff_1988(path, with_sentinel_hit=True)
-        da = _read_geotiff_gpu(path)
+        da = _read_geotiff_gpu(path, mask_nodata=True)
         assert da.attrs["nodata"] == 65535
         assert np.dtype(str(da.dtype)).kind == "f"
         assert da.attrs["masked_nodata"] is True
@@ -2662,7 +2662,7 @@ def test_dask_gpu_in_range_sentinel(self, tmp_path):
         from xrspatial.geotiff import _read_geotiff_gpu
         path = str(tmp_path / "tnss1988_gpu_dask_int.tif")
         _write_int_tiff_1988(path, with_sentinel_hit=False)
-        da = _read_geotiff_gpu(path, chunks=2)
+        da = _read_geotiff_gpu(path, chunks=2, mask_nodata=True)
         assert da.attrs["nodata"] == 65535
         assert np.dtype(str(da.dtype)).kind == "f"
         assert da.attrs["masked_nodata"] is True
diff --git a/xrspatial/geotiff/tests/read/test_rioxarray_compat_2961.py b/xrspatial/geotiff/tests/read/test_rioxarray_compat_2961.py
index b2551ac01..9402ef197 100644
--- a/xrspatial/geotiff/tests/read/test_rioxarray_compat_2961.py
+++ b/xrspatial/geotiff/tests/read/test_rioxarray_compat_2961.py
@@ -14,8 +14,11 @@
 import pytest
 import xarray as xr
 
-from xrspatial.geotiff import _build_vrt, open_geotiff, to_geotiff
+from xrspatial.geotiff import (
+    _build_vrt, _read_geotiff_dask, _read_geotiff_gpu, _read_vrt,
+    open_geotiff, to_geotiff)
 from xrspatial.geotiff._runtime import GeoTIFFFallbackWarning
+from xrspatial.geotiff.tests._helpers.markers import requires_gpu
 
 
 def _int_sentinel_tiff(path, sentinel=255):
@@ -64,6 +67,66 @@ def test_default_does_not_mask(tmp_path):
     assert out.attrs.get("nodata") == 255
 
 
+# ---------------------------------------------------------------------------
+# direct backend defaults match open_geotiff's unmasked default (#2976)
+#
+# The three direct backend entry points (_read_geotiff_dask,
+# _read_geotiff_gpu, _read_vrt) used to default to mask_nodata=True while
+# open_geotiff defaults to masked=False. A bare backend call therefore
+# returned a different dtype + NaN-substituted values than the public path.
+# These tests pin the backends to the public unmasked default.
+# ---------------------------------------------------------------------------
+
+def test_read_geotiff_dask_default_matches_open_geotiff_2976(tmp_path):
+    """Bare ``_read_geotiff_dask`` keeps the source dtype and sentinel,
+    matching ``open_geotiff(path, chunks=...)``."""
+    path = _int_sentinel_tiff(str(tmp_path / "t2976_dask.tif"))
+    public = open_geotiff(path, chunks=2).compute()
+    direct = _read_geotiff_dask(path, chunks=2).compute()
+
+    assert direct.dtype == public.dtype == np.uint8
+    assert (direct.data == 255).any()
+    assert not np.isnan(direct.data.astype(float)).any()
+    assert direct.attrs.get("masked_nodata") is False
+    assert direct.attrs.get("nodata") == public.attrs.get("nodata") == 255
+    np.testing.assert_array_equal(direct.data, public.data)
+
+
+def test_read_vrt_default_matches_open_geotiff_2976(tmp_path):
+    """Bare ``_read_vrt`` keeps the source dtype and sentinel, matching
+    ``open_geotiff(<vrt>)``."""
+    src = _int_sentinel_tiff(str(tmp_path / "t2976_vrt_src.tif"))
+    vrt = _build_vrt(str(tmp_path / "t2976.vrt"), source_files=[src])
+    public = open_geotiff(vrt)
+    direct = _read_vrt(vrt)
+
+    assert direct.dtype == public.dtype == np.uint8
+    assert (np.asarray(direct.values) == 255).any()
+    assert not np.isnan(np.asarray(direct.values, dtype=float)).any()
+    assert direct.attrs.get("masked_nodata") is False
+    assert direct.attrs.get("nodata") == public.attrs.get("nodata") == 255
+    np.testing.assert_array_equal(
+        np.asarray(direct.values), np.asarray(public.values))
+
+
+@requires_gpu
+def test_read_geotiff_gpu_default_matches_open_geotiff_2976(tmp_path):
+    """Bare ``_read_geotiff_gpu`` keeps the source dtype and sentinel,
+    matching ``open_geotiff(path, gpu=True)``."""
+    path = _int_sentinel_tiff(str(tmp_path / "t2976_gpu.tif"))
+    public = open_geotiff(path, gpu=True)
+    direct = _read_geotiff_gpu(path)
+
+    assert direct.dtype == public.dtype == np.uint8
+    direct_np = direct.data.get()
+    public_np = public.data.get()
+    assert (direct_np == 255).any()
+    assert not np.isnan(direct_np.astype(float)).any()
+    assert direct.attrs.get("masked_nodata") is False
+    assert direct.attrs.get("nodata") == public.attrs.get("nodata") == 255
+    np.testing.assert_array_equal(direct_np, public_np)
+
+
 def test_masked_true_promotes_and_masks(tmp_path):
     path = _int_sentinel_tiff(str(tmp_path / "t2961_masked.tif"))
     out = open_geotiff(path, masked=True)
diff --git a/xrspatial/geotiff/tests/release_gates/test_features.py b/xrspatial/geotiff/tests/release_gates/test_features.py
index d0e0ae94c..239915b55 100644
--- a/xrspatial/geotiff/tests/release_gates/test_features.py
+++ b/xrspatial/geotiff/tests/release_gates/test_features.py
@@ -2783,7 +2783,7 @@ def test_dask_nodata(self, tmp_path):
         path = str(tmp_path / 'dask_nodata.tif')
         write(arr, path, compression='none', tiled=False, nodata=-9999.0)
 
-        result = _read_geotiff_dask(path, chunks=2)
+        result = _read_geotiff_dask(path, chunks=2, mask_nodata=True)
         computed = result.compute()
         assert np.isnan(computed.values[0, 1])
         assert np.isnan(computed.values[1, 0])
diff --git a/xrspatial/geotiff/tests/release_gates/test_stable_features.py b/xrspatial/geotiff/tests/release_gates/test_stable_features.py
index 8a8f0743d..68bdc52d7 100644
--- a/xrspatial/geotiff/tests/release_gates/test_stable_features.py
+++ b/xrspatial/geotiff/tests/release_gates/test_stable_features.py
@@ -1023,13 +1023,13 @@ def test_release_gate_eager_dask_full_parity(
             f"`python -m xrspatial.geotiff.tests.golden_corpus.generate`"
         )
 
-    # ``_read_geotiff_dask`` keeps the legacy masking-on default
-    # (``mask_nodata=True``); ``open_geotiff`` flipped its default to
-    # ``masked=False``. Mirror the dask default on the eager call so the
-    # two backends are compared under the same masking policy, while
-    # still honouring an explicit ``mask_nodata=`` in ``open_kwargs``
-    # (the masked-nodata-lifecycle scenario).
-    eager_kwargs = {"masked": open_kwargs.get("mask_nodata", True),
+    # Both ``_read_geotiff_dask`` and ``open_geotiff`` now default to
+    # unmasked (``mask_nodata=False`` / ``masked=False``, see #2976).
+    # Mirror the dask masking choice on the eager call so the two
+    # backends are compared under the same masking policy, while still
+    # honouring an explicit ``mask_nodata=`` in ``open_kwargs`` (the
+    # masked-nodata-lifecycle scenario).
+    eager_kwargs = {"masked": open_kwargs.get("mask_nodata", False),
                     **{k: v for k, v in open_kwargs.items()
                        if k != "mask_nodata"}}
     eager = open_geotiff(str(path), **eager_kwargs)
@@ -1565,9 +1565,12 @@ def _overview_read_levels_eager(path: str) -> dict:
 
 
 def _overview_read_levels_dask(path: str) -> dict:
-    out = {0: _read_geotiff_dask(path, chunks=8)}
+    # Mirror the eager helper's ``masked=True``; the backend default is
+    # unmasked (#2976) so masking must be requested explicitly.
+    out = {0: _read_geotiff_dask(path, chunks=8, mask_nodata=True)}
     for i, _ in enumerate(_OVERVIEW_FACTORS, start=1):
-        out[i] = _read_geotiff_dask(path, chunks=8, overview_level=i)
+        out[i] = _read_geotiff_dask(
+            path, chunks=8, overview_level=i, mask_nodata=True)
     return out
 
 
diff --git a/xrspatial/geotiff/tests/vrt/test_dtype_conversion.py b/xrspatial/geotiff/tests/vrt/test_dtype_conversion.py
index 1eebc6270..fef2eb151 100644
--- a/xrspatial/geotiff/tests/vrt/test_dtype_conversion.py
+++ b/xrspatial/geotiff/tests/vrt/test_dtype_conversion.py
@@ -329,7 +329,7 @@ def test_uint64_nodata_masks_max_sentinel_in_data(tmp_path):
     src = tmp_path / 'src.tif'
     _dtype_validation_write(b, src)
     vrt = _dtype_validation_build_single_band_vrt(tmp_path, dtype_attr='UInt64', src_path=str(src), nodata=big)  # noqa: E501
-    r = _read_vrt(vrt)
+    r = _read_vrt(vrt, mask_nodata=True)
     assert r.dtype == np.float64, f'sentinel hit must promote to float64, got {r.dtype}'
     assert np.isnan(r.values[1, 1]), f'the 2**64-1 cell must be masked to NaN; got {r.values[1, 1]!r}'  # noqa: E501
     assert r.values[0, 0] == 1.0
@@ -347,7 +347,7 @@ def test_int64_min_nodata_masks_correctly(tmp_path):
     src = tmp_path / 'src.tif'
     _dtype_validation_write(b, src)
     vrt = _dtype_validation_build_single_band_vrt(tmp_path, dtype_attr='Int64', src_path=str(src), nodata=info.min)  # noqa: E501
-    r = _read_vrt(vrt)
+    r = _read_vrt(vrt, mask_nodata=True)
     assert r.dtype == np.float64
     assert np.isnan(r.values[0, 0])
     assert r.values[0, 1] == -1.0
@@ -366,7 +366,7 @@ def test_int32_negative_nodata_still_masks(tmp_path):
     src = tmp_path / 'src.tif'
     _dtype_validation_write(b, src)
     vrt = _dtype_validation_build_single_band_vrt(tmp_path, dtype_attr='Int32', src_path=str(src), nodata=-9999)  # noqa: E501
-    r = _read_vrt(vrt)
+    r = _read_vrt(vrt, mask_nodata=True)
     assert r.dtype == np.float64
     assert np.isnan(r.values[0, 1])
     assert np.isnan(r.values[1, 0])
@@ -509,7 +509,7 @@ def test_float32_vrt_uint16_source_masks_in_range_sentinel(tmp_path):
     """
     src = _int_source_float_dtype_write_uint16_with_sentinel(tmp_path)
     vrt = _int_source_float_dtype_build_vrt(tmp_path, src, 'Float32', 65535)
-    r = _read_vrt(vrt)
+    r = _read_vrt(vrt, mask_nodata=True)
     assert r.dtype == np.float32, f'Float32-declared VRT should return float32, got {r.dtype}'
     assert np.isnan(r.values[1, 1]), f'Sentinel pixel (uint16 65535 -> float32) should be NaN-masked; got values[1, 1]={r.values[1, 1]}'  # noqa: E501
     assert r.attrs.get('nodata') == 65535.0
@@ -520,7 +520,7 @@ def test_float64_vrt_int16_source_masks_negative_sentinel(tmp_path):
     """Float64 VRT, int16 source with negative sentinel: pixel becomes NaN."""
     src = _int_source_float_dtype_write_int16_with_sentinel(tmp_path, sentinel=-1)
     vrt = _int_source_float_dtype_build_vrt(tmp_path, src, 'Float64', -1)
-    r = _read_vrt(vrt)
+    r = _read_vrt(vrt, mask_nodata=True)
     assert r.dtype == np.float64
     assert np.isnan(r.values[1, 1]), f'Sentinel pixel (-1) should be NaN-masked; got values[1, 1]={r.values[1, 1]}'  # noqa: E501
     assert r.attrs.get('nodata') == -1.0
@@ -568,7 +568,7 @@ def test_float_vrt_int_source_dask_path_masks_sentinel(tmp_path):
     """
     src = _int_source_float_dtype_write_uint16_with_sentinel(tmp_path)
     vrt = _int_source_float_dtype_build_vrt(tmp_path, src, 'Float32', 65535)
-    r = _read_vrt(vrt, chunks=2)
+    r = _read_vrt(vrt, chunks=2, mask_nodata=True)
     assert r.dtype == np.float32
     val = r.values
     assert np.isnan(val[1, 1])
@@ -596,11 +596,11 @@ def test_float_vrt_int_source_with_band_select(tmp_path):
     vrt_path = str(tmp_path / 'mb.vrt')
     with open(vrt_path, 'w') as f:
         f.write(vrt_xml)
-    r0 = _read_vrt(vrt_path, band=0, band_nodata='first')
+    r0 = _read_vrt(vrt_path, band=0, band_nodata='first', mask_nodata=True)
     assert r0.dtype == np.float32
     assert np.isnan(r0.values[1, 1])
     assert r0.attrs.get('nodata') == 65535.0
-    r1 = _read_vrt(vrt_path, band=1, band_nodata='first')
+    r1 = _read_vrt(vrt_path, band=1, band_nodata='first', mask_nodata=True)
     assert r1.dtype == np.float32
     assert np.isnan(r1.values[1, 1])
     assert r1.attrs.get('nodata') == 65000.0
@@ -869,7 +869,7 @@ def test_multiband_uint16_per_band_sentinel_each_masked(tmp_path):
     as NaN but band 1's (1,1) cell as the literal 65000.0.
     """
     vrt_path = _multiband_int_nodata_write_two_band_per_band_nodata_vrt(tmp_path)
-    r = _read_vrt(vrt_path, band_nodata='first')
+    r = _read_vrt(vrt_path, band_nodata='first', mask_nodata=True)
     assert r.shape == (2, 2, 2)
     assert r.dtype == np.float64, f'expected float64 promotion, got {r.dtype}'
     assert np.isnan(r.values[1, 1, 0]), "band 0's sentinel pixel was not NaN-masked."
@@ -889,7 +889,7 @@ def test_multiband_int32_negative_per_band_sentinel(tmp_path):
     range guard accepts negatives.
     """
     vrt_path = _multiband_int_nodata_write_two_band_per_band_nodata_vrt(tmp_path, dtype_str='Int32', np_dtype=np.int32, band0_sentinel=-9999, band1_sentinel=-7777, band0_other=(10, 20, 30), band1_other=(40, 50, 60))  # noqa: E501
-    r = _read_vrt(vrt_path, band_nodata='first')
+    r = _read_vrt(vrt_path, band_nodata='first', mask_nodata=True)
     assert r.dtype == np.float64
     assert np.isnan(r.values[1, 1, 0])
     assert np.isnan(r.values[1, 1, 1])
@@ -909,7 +909,7 @@ def test_multiband_only_one_band_has_sentinel_present(tmp_path):
     import os
     p1 = os.path.join(os.path.dirname(vrt_path), 'vrt_b1_1611.tif')
     write(b1_no_sentinel, p1, nodata=65000, compression='none', tiled=False)
-    r = _read_vrt(vrt_path, band_nodata='first')
+    r = _read_vrt(vrt_path, band_nodata='first', mask_nodata=True)
     assert r.dtype == np.float64, "Even when only band 0 has a present sentinel, the array still needs promotion so band 0's NaN can be expressed."  # noqa: E501
     assert np.isnan(r.values[1, 1, 0])
     assert r.values[1, 1, 1] == 99.0
@@ -944,7 +944,7 @@ def test_multiband_per_band_out_of_range_sentinel_is_no_op(tmp_path):
     xml = xml.replace('<NoDataValue>10</NoDataValue>', '<NoDataValue>-9999</NoDataValue>')
     with open(vrt_path, 'w') as f:
         f.write(xml)
-    r = _read_vrt(vrt_path, band_nodata='first')
+    r = _read_vrt(vrt_path, band_nodata='first', mask_nodata=True)
     assert np.isnan(r.values[1, 1, 0])
     assert r.values[1, 1, 1] == 10.0 or r.values[1, 1, 1] == 10
 
@@ -957,8 +957,8 @@ def test_multiband_band_kwarg_still_per_band_post_pr1602(tmp_path):
     sentinel.
     """
     vrt_path = _multiband_int_nodata_write_two_band_per_band_nodata_vrt(tmp_path)
-    r0 = _read_vrt(vrt_path, band=0, band_nodata='first')
-    r1 = _read_vrt(vrt_path, band=1, band_nodata='first')
+    r0 = _read_vrt(vrt_path, band=0, band_nodata='first', mask_nodata=True)
+    r1 = _read_vrt(vrt_path, band=1, band_nodata='first', mask_nodata=True)
     assert r0.dtype == np.float64
     assert r1.dtype == np.float64
     assert r0.attrs.get('nodata') == 65535.0
diff --git a/xrspatial/geotiff/tests/vrt/test_metadata.py b/xrspatial/geotiff/tests/vrt/test_metadata.py
index 940602f6c..0181a4bd6 100644
--- a/xrspatial/geotiff/tests/vrt/test_metadata.py
+++ b/xrspatial/geotiff/tests/vrt/test_metadata.py
@@ -269,7 +269,7 @@ def test_vrt_chunked_float_source_mask_off_reports_false(tmp_path):
 def test_vrt_chunked_float_source_mask_on_reports_true(tmp_path):
     """Canonical direction on the chunked path: masking on, attr True."""
     vrt = _masked_nodata_attr_write_float_vrt(tmp_path, 'tmp_2159_chunked_float_src_masked.tif', 'tmp_2159_chunked_masked.vrt')  # noqa: E501
-    out = _read_geotiff_dask(vrt, chunks=2)
+    out = _read_geotiff_dask(vrt, chunks=2, mask_nodata=True)
     assert out.attrs.get('nodata') == -9999.0
     assert out.attrs.get('masked_nodata') is True
 
@@ -343,7 +343,7 @@ def test_read_vrt_band0_uses_band0_nodata(tmp_path):
     at the call site that the test is exercising the legacy behaviour.
     """
     vrt_path = _band_nodata_write_two_band_per_band_nodata_vrt(tmp_path)
-    r = _read_vrt(vrt_path, band=0, band_nodata='first')
+    r = _read_vrt(vrt_path, band=0, band_nodata='first', mask_nodata=True)
     assert r.dtype == np.float64
     assert r.attrs.get('nodata') == 65535.0
     assert np.isnan(r.values[1, 1])
@@ -357,7 +357,7 @@ def test_read_vrt_band1_uses_band1_nodata(tmp_path):
     [9,65000]] and attrs['nodata']=65535.
     """
     vrt_path = _band_nodata_write_two_band_per_band_nodata_vrt(tmp_path)
-    r = _read_vrt(vrt_path, band=1, band_nodata='first')
+    r = _read_vrt(vrt_path, band=1, band_nodata='first', mask_nodata=True)
     assert r.dtype == np.float64, 'band=1 read kept uint16 dtype; per-band nodata regression.'
     assert r.attrs.get('nodata') == 65000.0, f"attrs['nodata'] was {r.attrs.get('nodata')}, expected 65000 from band 1's <NoDataValue>."  # noqa: E501
     assert np.isnan(r.values[1, 1]), "band 1's sentinel pixel was not NaN-masked; promotion ran against the wrong sentinel."  # noqa: E501
@@ -510,7 +510,7 @@ def test_vrt_uint16_nodata_promotes_to_float64(tmp_path):
     assert np.isnan(eager.values[1, 0])
     vrt_path = str(tmp_path / 'src_1564.vrt')
     _build_vrt(vrt_path, [tif])
-    via_vrt = _read_vrt(vrt_path)
+    via_vrt = _read_vrt(vrt_path, mask_nodata=True)
     assert via_vrt.dtype == np.float64, f'VRT integer-with-nodata should promote to float64; got {via_vrt.dtype}'  # noqa: E501
     assert np.isnan(via_vrt.values[1, 0]), f'VRT sentinel pixel should be NaN; got {via_vrt.values[1, 0]} (literal sentinel survived)'  # noqa: E501
     assert via_vrt.attrs.get('nodata') == 65535.0
@@ -538,7 +538,7 @@ def test_vrt_float_nodata_still_masks(tmp_path):
     to_geotiff(da, tif, compression='none', nodata=-9999.0)
     vrt_path = str(tmp_path / 'srcf_1564.vrt')
     _build_vrt(vrt_path, [tif])
-    via_vrt = _read_vrt(vrt_path)
+    via_vrt = _read_vrt(vrt_path, mask_nodata=True)
     assert via_vrt.dtype == np.float32
     assert np.isnan(via_vrt.values[0, 2])
     assert np.isnan(via_vrt.values[1, 1])
@@ -642,22 +642,25 @@ def _mask_nodata_float_build_vrt(tmp_path, source_path, vrt_dtype, nodata_value,
     return p
 
 
-def test_default_mask_nodata_true_rewrites_float_sentinel(tmp_path):
-    """The default behaviour (mask_nodata=True) still substitutes NaN.
+def test_default_unmasked_keeps_float_sentinel_2976(tmp_path):
+    """The default (``mask_nodata=False``) matches ``open_geotiff``.
 
-    Pins the existing contract so the fix below does not regress the
-    masking happy path.
+    A bare ``_read_vrt`` keeps the float sentinel literal rather than
+    substituting NaN. ``mask_nodata=True`` restores the NaN substitution
+    (see ``test_default_mask_nodata_true_rewrites_float_sentinel``'s
+    masked sibling).
     """
     src, _ = _mask_nodata_float_write_float32_with_sentinel(tmp_path)
     vrt = _mask_nodata_float_build_vrt(tmp_path, src, 'Float32', -9999.0)
     r = _read_vrt(vrt)
     assert r.dtype == np.float32
-    assert np.isnan(r.values[1, 1])
-    assert np.isnan(r.values[2, 1])
+    assert not np.isnan(r.values).any()
+    assert r.values[1, 1] == np.float32(-9999.0)
+    assert r.values[2, 1] == np.float32(-9999.0)
     assert r.values[0, 0] == 1.0
     assert r.values[1, 0] == 4.0
     assert r.attrs.get('nodata') == -9999.0
-    assert r.attrs.get('masked_nodata') is True
+    assert r.attrs.get('masked_nodata') is False
 
 
 def test_eager_mask_nodata_false_preserves_float_sentinel(tmp_path):
@@ -741,7 +744,7 @@ def test_masked_vs_unmasked_differ_only_at_sentinels(tmp_path):
     """
     src, _ = _mask_nodata_float_write_float32_with_sentinel(tmp_path)
     vrt = _mask_nodata_float_build_vrt(tmp_path, src, 'Float32', -9999.0)
-    masked = _read_vrt(vrt).values
+    masked = _read_vrt(vrt, mask_nodata=True).values
     unmasked = _read_vrt(vrt, mask_nodata=False).values
     nan_positions = np.isnan(masked)
     sentinel_positions = unmasked == np.float32(-9999.0)
@@ -785,21 +788,22 @@ def test_int_source_float_vrt_mask_nodata_false_keeps_literal(tmp_path):
     assert r.attrs.get('masked_nodata') is False
 
 
-def test_int_source_float_vrt_default_still_promotes(tmp_path):
-    """Default ``mask_nodata=True`` still NaN-masks the int->float promotion.
+def test_int_source_float_vrt_default_unmasked_2976(tmp_path):
+    """Default ``mask_nodata=False`` matches ``open_geotiff``.
 
-    Baseline that documents the default contract for the integer
-    source path: the int->float NaN-promotion behavior is unchanged
-    when the opt-out is not requested.
+    An integer source feeding a Float32 VRT casts the literal sentinel
+    to ``65535.0`` rather than NaN-masking it. ``mask_nodata=True``
+    restores the int->float NaN promotion.
     """
     src, _ = _mask_nodata_float_write_uint16_with_sentinel(tmp_path)
     vrt = _mask_nodata_float_build_vrt(tmp_path, src, 'Float32', 65535, filename='int_float_default_2158.vrt', shape=(2, 2))  # noqa: E501
     r = _read_vrt(vrt)
     assert r.dtype == np.float32
-    assert np.isnan(r.values[1, 1])
+    assert not np.isnan(r.values).any()
+    assert r.values[1, 1] == np.float32(65535.0)
     assert r.values[0, 0] == 1.0
     assert r.attrs.get('nodata') == 65535.0
-    assert r.attrs.get('masked_nodata') is True
+    assert r.attrs.get('masked_nodata') is False
 
 
 # ---------------------------------------------------------------------------
@@ -1406,9 +1410,11 @@ def _metadata_parity_read_gpu_eager(vrt_path: str):
     ``open_geotiff(..., gpu=True)`` rejects ``.vrt`` sources up front
     (the dispatcher routes ``.vrt`` to ``_read_vrt`` and ``_read_vrt``
     owns the ``gpu`` kwarg, see ``_backends/vrt.py``). Use the direct
-    entry point here so the GPU eager path is exercised.
+    entry point here so the GPU eager path is exercised. ``mask_nodata``
+    is passed explicitly to match the masked numpy / dask baselines
+    above (the backend default is unmasked, see #2976).
     """
-    return _read_vrt(vrt_path, gpu=True)
+    return _read_vrt(vrt_path, gpu=True, mask_nodata=True)
 
 
 _BACKENDS = [pytest.param('numpy', _metadata_parity_read_eager_numpy, id='numpy'), pytest.param('dask', _metadata_parity_read_dask, id='dask'), pytest.param('gpu', _metadata_parity_read_gpu_eager, id='gpu', marks=requires_gpu)]  # noqa: E501
diff --git a/xrspatial/geotiff/tests/vrt/test_parity.py b/xrspatial/geotiff/tests/vrt/test_parity.py
index 8be3cd99b..a76c308e2 100644
--- a/xrspatial/geotiff/tests/vrt/test_parity.py
+++ b/xrspatial/geotiff/tests/vrt/test_parity.py
@@ -845,7 +845,7 @@ def test_band_nodata_first_band_attrs(tmp_path):
     """``band=1`` with ``band_nodata='first'`` surfaces band 1's sentinel
     on attrs and masks against it."""
     vrt_path = _write_two_band_per_band_nodata_vrt(tmp_path)
-    r = _read_vrt(vrt_path, band=1, band_nodata='first')
+    r = _read_vrt(vrt_path, band=1, band_nodata='first', mask_nodata=True)
     assert r.attrs['nodata'] == 65000.0
     assert r.attrs['masked_nodata'] is True
     assert np.isnan(r.values[1, 1])
@@ -855,7 +855,8 @@ def test_band_nodata_first_band_attrs(tmp_path):
 def test_band_nodata_chunked_first_band_attrs(tmp_path):
     """The chunked path threads the same per-band sentinel onto attrs."""
     vrt_path = _write_two_band_per_band_nodata_vrt(tmp_path)
-    r = _read_vrt(vrt_path, band=1, band_nodata='first', chunks=2)
+    r = _read_vrt(vrt_path, band=1, band_nodata='first', chunks=2,
+                  mask_nodata=True)
     assert r.attrs['nodata'] == 65000.0
     assert r.attrs['masked_nodata'] is True
     assert 'nodata_pixels_present' not in r.attrs
diff --git a/xrspatial/geotiff/tests/write/test_nodata.py b/xrspatial/geotiff/tests/write/test_nodata.py
index 26bd4c57a..1e5c54fd3 100644
--- a/xrspatial/geotiff/tests/write/test_nodata.py
+++ b/xrspatial/geotiff/tests/write/test_nodata.py
@@ -413,7 +413,7 @@ def test_uint64_max_masked_via_dask(self, tmp_path):
         da_in = xr.DataArray(arr, dims=("y", "x"))
         path = os.path.join(str(tmp_path), "t.tif")
         to_geotiff(da_in, path, nodata=2**64 - 1)
-        out = _read_geotiff_dask(path, chunks=16).compute()
+        out = _read_geotiff_dask(path, chunks=16, mask_nodata=True).compute()
         assert out.dtype == np.float64
         assert np.isnan(out.values[0, 0])
         assert out.values[1, 1] == 100.0
@@ -424,7 +424,7 @@ def test_int64_max_masked_via_dask(self, tmp_path):
         da_in = xr.DataArray(arr, dims=("y", "x"))
         path = os.path.join(str(tmp_path), "t.tif")
         to_geotiff(da_in, path, nodata=2**63 - 1)
-        out = _read_geotiff_dask(path, chunks=16).compute()
+        out = _read_geotiff_dask(path, chunks=16, mask_nodata=True).compute()
         assert out.dtype == np.float64
         assert np.isnan(out.values[0, 0])
 
@@ -449,7 +449,7 @@ def test_uint64_max_round_trip_via_vrt(self, tmp_path):
             xml = f.read()
         assert "<NoDataValue>18446744073709551615</NoDataValue>" in xml
 
-        out = _read_vrt(vrt_path)
+        out = _read_vrt(vrt_path, mask_nodata=True)
         assert out.dtype == np.float64
         assert np.isnan(out.values[0, 0])
         assert out.values[1, 1] == 100.0
@@ -469,7 +469,7 @@ def test_int64_max_round_trip_via_vrt(self, tmp_path):
             xml = f.read()
         assert "<NoDataValue>9223372036854775807</NoDataValue>" in xml
 
-        out = _read_vrt(vrt_path)
+        out = _read_vrt(vrt_path, mask_nodata=True)
         assert out.dtype == np.float64
         assert np.isnan(out.values[0, 0])
         assert out.values[1, 1] == 100.0
@@ -486,7 +486,7 @@ def test_uint64_max_masked_via_gpu(self, tmp_path):
         path = os.path.join(str(tmp_path), "t.tif")
         to_geotiff(da_in, path, nodata=2**64 - 1)
 
-        gpu_da = _read_geotiff_gpu(path)
+        gpu_da = _read_geotiff_gpu(path, mask_nodata=True)
         host = gpu_da.data.get()
         assert host.dtype == np.float64
         assert np.isnan(host[0, 0])