From d911bb54c9fe0ffeee74a56de733c8cf80d266ca Mon Sep 17 00:00:00 2001
From: Chris Holden <ceholden@gmail.com>
Date: Wed, 4 Nov 2026 12:53:07 -0500
Subject: [PATCH 01/11] perf: avoid duplicate sort in shadow fill
 reconstruction

Peak memory on l9_clear drops from 12.71 to 10.97 GB and wall time from 328
to 302 s; s2 drops from 5.89 to 5.75 GB and 79 to 72 s. Final masks are
identical on all three test scenes.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/fmask/_reconstruct.py | 103 ++++++++++++++++++
 src/fmask/utils.py        |   7 +-
 tests/test_reconstruct.py | 214 ++++++++++++++++++++++++++++++++++++++
 3 files changed, 320 insertions(+), 4 deletions(-)
 create mode 100644 src/fmask/_reconstruct.py
 create mode 100644 tests/test_reconstruct.py

diff --git a/src/fmask/_reconstruct.py b/src/fmask/_reconstruct.py
new file mode 100644
index 0000000..31d1f5c
--- /dev/null
+++ b/src/fmask/_reconstruct.py
@@ -0,0 +1,103 @@
+"""Grayscale reconstruction by erosion that sorts the image once."""
+
+import warnings
+
+import numpy as np
+import skimage
+from skimage.morphology import reconstruction
+
+MIRRORED_VERSION_PREFIX = "0.26."
+
+
+def _mirror_matches(version: str) -> bool:
+    """Return True for the scikit-image versions whose private internals this module mirrors."""
+    return version.startswith(MIRRORED_VERSION_PREFIX)
+
+
+try:
+    from skimage._shared.utils import _supported_float_type
+    from skimage.morphology._grayreconstruct import reconstruction_loop
+
+    MIRROR_AVAILABLE = _mirror_matches(skimage.__version__)
+except ImportError:
+    MIRROR_AVAILABLE = False
+
+if not MIRROR_AVAILABLE:
+    warnings.warn(
+        f"fmask._reconstruct mirrors scikit-image {MIRRORED_VERSION_PREFIX}x internals, so "
+        f"scikit-image {skimage.__version__} uses its slower reconstruction, which gives the "
+        "same results. Update fmask._reconstruct to regain the speedup.",
+        RuntimeWarning,
+        stacklevel=2,
+    )
+
+
+def _needs_skimage(images: np.ndarray) -> bool:
+    """Return True for values whose result depends on which equal value skimage's sort keeps."""
+    return bool(np.isnan(images).any() or np.signbit(images[images == 0]).any())
+
+
+def reconstruct_erosion(seed: np.ndarray, mask: np.ndarray) -> np.ndarray:
+    """Return skimage.morphology.reconstruction(seed, mask, method="erosion") for 2D images."""
+    if not MIRROR_AVAILABLE:
+        return reconstruction(seed, mask, method="erosion")
+    assert tuple(seed.shape) == tuple(mask.shape)
+    if np.any(seed < mask):
+        raise ValueError(
+            "Intensity of seed image must be greater than that "
+            "of the mask image for reconstruction by erosion."
+        )
+    rows, cols = seed.shape
+    padded_cols = cols + 2
+    image_stride = (rows + 2) * padded_cols
+
+    images = np.full(
+        (2, rows + 2, padded_cols), np.max(seed), dtype=_supported_float_type(mask.dtype)
+    )
+    images[0, 1:-1, 1:-1] = seed
+    images[1, 1:-1, 1:-1] = mask
+    if _needs_skimage(images):
+        return reconstruction(seed, mask, method="erosion")
+    images = images.reshape(-1)
+
+    isize = images.size
+    signed_int_dtype = np.result_type(np.min_scalar_type(-isize), np.int32)
+    unsigned_int_dtype = np.dtype(signed_int_dtype.char.upper())
+    neighbor_strides = np.array(
+        [
+            -padded_cols - 1, -padded_cols, -padded_cols + 1,
+            -1, 1,
+            padded_cols - 1, padded_cols, padded_cols + 1,
+        ],
+        signed_int_dtype,
+    )
+
+    index_sorted = np.argsort(images).astype(signed_int_dtype, copy=False)
+    sorted_values = images[index_sorted]
+    del images
+    is_different = sorted_values[:-1] != sorted_values[1:]
+
+    ascending_rank = np.zeros(isize, unsigned_int_dtype)
+    np.cumsum(is_different, out=ascending_rank[1:], dtype=unsigned_int_dtype)
+    value_map = np.concatenate((sorted_values[:1], sorted_values[1:][is_different]))[::-1]
+    del sorted_values, is_different
+
+    # The reconstruction loop visits ranks in ascending order, so erosion needs the
+    # largest value at rank 0.
+    value_rank = np.empty(isize, unsigned_int_dtype)
+    value_rank[index_sorted] = ascending_rank[-1] - ascending_rank
+    del ascending_rank
+
+    prev = np.full(isize, -1, signed_int_dtype)
+    next_ = np.full(isize, -1, signed_int_dtype)
+    prev[index_sorted[1:]] = index_sorted[:-1]
+    next_[index_sorted[:-1]] = index_sorted[1:]
+    start = index_sorted[0]
+    del index_sorted
+
+    reconstruction_loop(value_rank, prev, next_, neighbor_strides, start, image_stride)
+    del prev, next_
+
+    filled = value_map[value_rank[:image_stride]]
+    filled.shape = (rows + 2, padded_cols)
+    return filled[1:-1, 1:-1]
diff --git a/src/fmask/utils.py b/src/fmask/utils.py
index 38b2002..7fdf1ef 100644
--- a/src/fmask/utils.py
+++ b/src/fmask/utils.py
@@ -26,12 +26,13 @@
 import pandas
 from skimage.measure import block_reduce, label, regionprops
 from skimage.feature import local_binary_pattern
-from skimage.morphology import binary_dilation, binary_erosion, reconstruction
+from skimage.morphology import binary_dilation, binary_erosion
 import matplotlib as mpl
 import matplotlib.pyplot as plt
 from matplotlib.colors import ListedColormap
 from . import constant as C
 from ._paths import get_data_root as _get_data_root
+from ._reconstruct import reconstruct_erosion
 from scipy.ndimage import convolve
 
 # ignore the invalid errors
@@ -539,9 +540,7 @@ def imfill(data, obsmask, fill_value=None):
     # Check the assertion
     # assert np.all(seed <= data), "Error: seed must be <= data for erosion."
     # seed = np.minimum(seed, data)
-    data_filled = reconstruction(
-        seed, data, method="erosion"
-    )  # a 3x3 square, of which return is float64
+    data_filled = reconstruct_erosion(seed, data)  # a 3x3 square
     # data_filled = data_filled.astype(
     #    np.float32
     # )  # covnert to the scalar that we want to use
diff --git a/tests/test_reconstruct.py b/tests/test_reconstruct.py
new file mode 100644
index 0000000..fa3d62d
--- /dev/null
+++ b/tests/test_reconstruct.py
@@ -0,0 +1,214 @@
+import numpy as np
+import pytest
+from skimage.morphology import reconstruction
+
+from fmask import _reconstruct, utils
+from fmask._reconstruct import reconstruct_erosion
+
+
+def _hole_seed(mask):
+    seed = mask.copy()
+    if mask.shape[0] > 2 and mask.shape[1] > 2:
+        seed[1:-1, 1:-1] = np.nanmax(mask)
+    return seed
+
+
+def _random(dtype, seed=0, shape=(40, 57)):
+    return np.random.default_rng(seed).random(shape).astype(dtype)
+
+
+def _ties(dtype):
+    return np.random.default_rng(1).integers(0, 4, (40, 57)).astype(dtype)
+
+
+def _with_nan(dtype):
+    mask = _random(dtype, 2)
+    mask[5, 7] = np.nan
+    mask[20:22, 30:33] = np.nan
+    return mask
+
+
+def _signed_zeros(dtype):
+    rng = np.random.default_rng(3)
+    mask = rng.integers(0, 3, (40, 57)).astype(dtype)
+    mask[mask == 0] = rng.choice([-0.0, 0.0], size=int((mask == 0).sum()))
+    return mask
+
+
+def _only_negative_zero(dtype):
+    return np.full((9, 11), -0.0, dtype=dtype)
+
+
+def _positive_zero_and_values(dtype):
+    mask = _ties(dtype)
+    mask[0, 0] = 0.0
+    return mask
+
+
+def _infinities(dtype):
+    mask = _random(dtype, 4)
+    mask[3, 3] = -np.inf
+    mask[0, 4] = -np.inf
+    return mask
+
+
+def _negative_values(dtype):
+    return (_random(dtype, 5) - 0.5).astype(dtype)
+
+
+CASES = {
+    "random": _random,
+    "ties": _ties,
+    "constant": lambda dtype: np.full((9, 11), 2.5, dtype=dtype),
+    "nan": _with_nan,
+    "signed_zeros": _signed_zeros,
+    "only_negative_zero": _only_negative_zero,
+    "positive_zero": _positive_zero_and_values,
+    "infinities": _infinities,
+    "negative_values": _negative_values,
+    "one_by_one": lambda dtype: np.array([[1.5]], dtype=dtype),
+    "one_by_n": lambda dtype: _random(dtype, 6, (1, 9)),
+    "n_by_one": lambda dtype: _random(dtype, 7, (9, 1)),
+}
+
+
+def _assert_same(actual, expected):
+    assert actual.dtype == expected.dtype
+    assert actual.shape == expected.shape
+    assert actual.tobytes() == expected.tobytes()
+
+
+@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+@pytest.mark.parametrize("case", sorted(CASES))
+def test_matches_skimage(case, dtype):
+    mask = CASES[case](dtype)
+    seed = _hole_seed(mask)
+    _assert_same(
+        reconstruct_erosion(seed, mask),
+        reconstruction(seed, mask, method="erosion"),
+    )
+
+
+@pytest.mark.parametrize("seed_dtype", [np.float32, np.float64])
+@pytest.mark.parametrize("mask_dtype", [np.uint8, np.int16, np.bool_, np.float16])
+def test_matches_skimage_integer_and_small_types(mask_dtype, seed_dtype):
+    mask = np.random.default_rng(8).integers(0, 2, (20, 23)).astype(mask_dtype)
+    seed = _hole_seed(mask.astype(seed_dtype))
+    _assert_same(
+        reconstruct_erosion(seed, mask),
+        reconstruction(seed, mask, method="erosion"),
+    )
+
+
+def test_matches_skimage_on_imfill_input():
+    mask = _random(np.float32, 9, (61, 47)) * 100
+    mask[10:30, 10:30] -= 50
+    seed = mask.copy()
+    seed[1:-1, 1:-1] = mask.max()
+    _assert_same(
+        reconstruct_erosion(seed, mask),
+        reconstruction(seed, mask, method="erosion"),
+    )
+
+
+def test_leaves_inputs_unchanged():
+    mask = _random(np.float32, 10)
+    seed = _hole_seed(mask)
+    mask_before, seed_before = mask.copy(), seed.copy()
+    reconstruct_erosion(seed, mask)
+    _assert_same(mask, mask_before)
+    _assert_same(seed, seed_before)
+
+
+def test_seed_below_mask_raises_same_error():
+    mask = np.ones((5, 5), dtype=np.float32)
+    seed = np.zeros((5, 5), dtype=np.float32)
+    with pytest.raises(ValueError) as expected:
+        reconstruction(seed, mask, method="erosion")
+    with pytest.raises(ValueError) as actual:
+        reconstruct_erosion(seed, mask)
+    assert str(actual.value) == str(expected.value)
+
+
+def test_shape_mismatch_raises_same_error():
+    with pytest.raises(AssertionError):
+        reconstruction(np.ones((5, 5)), np.ones((4, 5)), method="erosion")
+    with pytest.raises(AssertionError):
+        reconstruct_erosion(np.ones((5, 5)), np.ones((4, 5)))
+
+
+def test_private_api_available():
+    assert _reconstruct.MIRROR_AVAILABLE
+    from skimage.morphology._grayreconstruct import reconstruction_loop
+
+    ranks = np.array([1, 1, 1, 1, 1, 1, 1, 1, 0], dtype=np.uint32)
+    prev = np.full(9, -1, np.int32)
+    next_ = np.full(9, -1, np.int32)
+    strides = np.array([1, -1], np.int32)
+    reconstruction_loop(ranks, prev, next_, strides, 8, 3)
+
+
+@pytest.mark.parametrize(
+    "version,expected",
+    [("0.26.0", True), ("0.26.3", True), ("0.27.0", False), ("0.25.2", False), ("0.260.0", False)],
+)
+def test_mirror_is_limited_to_the_mirrored_version(version, expected):
+    assert _reconstruct._mirror_matches(version) is expected
+
+
+def test_uses_skimage_when_mirror_is_unavailable(monkeypatch):
+    calls = []
+
+    def spy(seed, mask, method):
+        calls.append(method)
+        return reconstruction(seed, mask, method=method)
+
+    monkeypatch.setattr(_reconstruct, "MIRROR_AVAILABLE", False)
+    monkeypatch.setattr(_reconstruct, "reconstruction", spy)
+    mask = _random(np.float32, 13)
+    seed = _hole_seed(mask)
+    _assert_same(reconstruct_erosion(seed, mask), reconstruction(seed, mask, method="erosion"))
+    assert calls == ["erosion"]
+
+
+@pytest.mark.parametrize("fill_value", [None, -1.0])
+@pytest.mark.parametrize("obs_shape", ["full", "partial", "single"])
+def test_imfill_matches_skimage(obs_shape, fill_value):
+    data = _random(np.float32, 11, (31, 37)) * 10
+    obsmask = np.ones(data.shape, dtype=bool)
+    if obs_shape == "partial":
+        obsmask[:5] = False
+        obsmask[:, -4:] = False
+    elif obs_shape == "single":
+        obsmask[:] = False
+        obsmask[15, 15] = True
+    max_val = data[obsmask].max()
+    filled_input = data.copy()
+    filled_input[~obsmask] = max_val if fill_value is None else fill_value
+    seed = filled_input.copy()
+    seed[1:-1, 1:-1] = max_val
+
+    result = utils.imfill(data, obsmask, fill_value=fill_value)
+
+    _assert_same(result, reconstruction(seed, filled_input, method="erosion"))
+    _assert_same(data, filled_input)
+
+
+def test_imfill_fills_a_closed_hole():
+    data = np.full((5, 5), 5.0, dtype=np.float32)
+    data[2, 2] = 1.0
+    result = utils.imfill(data, np.ones((5, 5), dtype=bool))
+    assert result[2, 2] == 5.0
+
+
+def test_imfill_keeps_a_depression_open_to_the_border():
+    data = np.full((5, 5), 5.0, dtype=np.float32)
+    data[2, :3] = 1.0
+    result = utils.imfill(data, np.ones((5, 5), dtype=bool))
+    assert (result[2, :3] == 1.0).all()
+
+
+def test_imfill_empty_observation_mask_raises():
+    data = _random(np.float32, 12, (8, 8))
+    with pytest.raises(ValueError):
+        utils.imfill(data, np.zeros(data.shape, dtype=bool))

From b22dc8bd55f7eab48905b261b93f49948d46ff7f Mon Sep 17 00:00:00 2001
From: Chris Holden <ceholden@gmail.com>
Date: Wed, 4 Nov 2026 12:53:07 -0500
Subject: [PATCH 02/11] perf: normalize in blocks without full-image copies

UNet classify no longer copies the predictor bands before normalizing, and
normalization no longer copies the cube or scales whole bands in float64 at
once. Byte-for-byte comparison with v5.0.1-2 on synthetic data passed on
perf/equivalence-tests; final masks and timing are to be evaluated on AWS Batch.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/fmask/unetlib.py        |  11 ++-
 src/fmask/utils.py          |  50 +++++++++--
 tests/test_normalize.py     | 164 ++++++++++++++++++++++++++++++++++++
 tests/test_unet_classify.py |  66 +++++++++++++++
 4 files changed, 280 insertions(+), 11 deletions(-)
 create mode 100644 tests/test_normalize.py
 create mode 100644 tests/test_unet_classify.py

diff --git a/src/fmask/unetlib.py b/src/fmask/unetlib.py
index 50607fb..61c1b0c 100644
--- a/src/fmask/unetlib.py
+++ b/src/fmask/unetlib.py
@@ -12,7 +12,7 @@
 from torch.utils.data import DataLoader, Dataset as BaseDataset
 import segmentation_models_pytorch as smp
 from pathlib import Path
-from .utils import normalize_datacube, init_patch_offanchors
+from .utils import normalize_bands, normalize_datacube, init_patch_offanchors
 import time
 from . import constant as C
 np.seterr(invalid='ignore') # ignore the invalid errors
@@ -547,9 +547,14 @@ def classify(self, probability="none", normalize=True, shift=True) -> tuple:
         )
 
         # rescale the date cube
-        datacube = self.image.data.get(self.predictors)
         if normalize:
-            datacube = normalize_datacube(datacube, obsmask=self.image.obsmask)
+            datacube = normalize_bands(
+                self.image.data.data,
+                [self.image.data.bands.index(predictor) for predictor in self.predictors],
+                obsmask=self.image.obsmask,
+            )
+        else:
+            datacube = self.image.data.get(self.predictors)
     
         if C.MSG_FULL:
             print(">>> classifying image by unet")
diff --git a/src/fmask/utils.py b/src/fmask/utils.py
index 7fdf1ef..e61994a 100644
--- a/src/fmask/utils.py
+++ b/src/fmask/utils.py
@@ -35,6 +35,8 @@
 from ._reconstruct import reconstruct_erosion
 from scipy.ndimage import convolve
 
+NORMALIZE_BLOCK_ROWS = 1024
+
 # ignore the invalid errors
 np.seterr(invalid='ignore') 
 
@@ -1185,6 +1187,43 @@ def normalize_image(image, **kwargs):
     image_scaled[np.where(image_scaled > normal_scale_range[1])] = normal_scale_range[1]
     return image_scaled
 
+def normalize_bands(
+    cube, indices, obsmask=None, percentiles=(1, 99), srange=(-1, 1), block_rows=NORMALIZE_BLOCK_ROWS
+):
+    """
+    Normalize selected bands of a datacube into a new array, a block of rows at a time.
+
+    Args:
+        cube (ndarray): The input datacube (band, row, column), which is not changed.
+        indices (list): The bands to normalize, in output order.
+        obsmask (ndarray, optional): The observation mask. Pixels out of the mask get min(srange).
+        percentiles (sequence, optional): The percentiles used for rescaling.
+        srange (sequence, optional): The scale range used for rescaling.
+        block_rows (int, optional): The number of rows scaled in float64 at once.
+
+    Returns:
+        ndarray: The normalized bands, with the dtype of the input.
+    """
+    rows = cube.shape[1]
+    normalized = np.empty((len(indices), *cube.shape[1:]), dtype=cube.dtype)
+    for output_index, band_index in enumerate(indices):
+        band = cube[band_index]
+        valid = band.astype("float") if obsmask is None else band[obsmask].astype("float")
+        min_val, max_val = np.percentile(valid, percentiles, method="linear", overwrite_input=True)
+        del valid
+        for start in range(0, rows, block_rows):
+            stop = start + block_rows
+            # float64 scalars keep this arithmetic in float64 before the cast on assignment.
+            scaled = ((band[start:stop] - min_val) / (max_val - min_val + C.EPS)) * (
+                srange[1] - srange[0]
+            ) + srange[0]
+            scaled[scaled < srange[0]] = srange[0]
+            scaled[scaled > srange[1]] = srange[1]
+            if obsmask is not None:
+                scaled[~obsmask[start:stop]] = min(srange)
+            normalized[output_index, start:stop] = scaled
+    return normalized
+
 def normalize_datacube(datacube, **kwargs):
     """
     Normalize a datacube by rescaling each layer.
@@ -1203,16 +1242,11 @@ def normalize_datacube(datacube, **kwargs):
     percentiles = kwargs.get("percentiles", [1, 99])
     srange = kwargs.get("srange", [-1, 1])
     obsmask = kwargs.get("obsmask", None)
-    _datacube = datacube.copy() # not to alter the orginal values
     if C.MSG_FULL:
         print(f">>> normalizing the datacube to {srange} with percentiles {percentiles}")
-    # rescale the data cube
-    for i in range(0, _datacube.shape[0]):
-        datalayer = normalize_image(_datacube[i, :, :], obsmask=obsmask, percentiles=percentiles, srange=srange)
-        if obsmask is not None:
-            datalayer[~obsmask] = min(srange) # give min value to the pixels out of the mask, as background
-        _datacube[i, :, :] = datalayer
-    return _datacube
+    return normalize_bands(
+        datacube, list(range(datacube.shape[0])), obsmask=obsmask, percentiles=percentiles, srange=srange
+    )
 
 def select_cloud_seed(
     cloud_mask, cloud_prob, label_cloud, label_noncloud, label_filled, seed_level_cloud, seed_level_noncloud
diff --git a/tests/test_normalize.py b/tests/test_normalize.py
new file mode 100644
index 0000000..7c3de12
--- /dev/null
+++ b/tests/test_normalize.py
@@ -0,0 +1,164 @@
+import numpy as np
+import pytest
+
+from fmask import constant as C
+from fmask import utils
+from fmask.satellite import Data
+
+SHAPE = (4, 23, 31)
+
+
+def _cube(seed=0, dtype=np.float32):
+    return (np.random.default_rng(seed).random(SHAPE) * 3000 - 500).astype(dtype)
+
+
+def _partial_mask():
+    mask = np.ones(SHAPE[1:], dtype=bool)
+    mask[:4] = False
+    mask[:, -3:] = False
+    mask[10, 10] = False
+    return mask
+
+
+def _assert_same(actual, expected):
+    assert actual.dtype == expected.dtype
+    assert actual.shape == expected.shape
+    assert actual.tobytes() == expected.tobytes()
+
+
+def _bands(cube, **kwargs):
+    return utils.normalize_bands(cube, list(range(cube.shape[0])), **kwargs)
+
+
+def _scale_in_float64(band, valid, percentiles=(1, 99), srange=(-1, 1)):
+    low, high = np.percentile(valid.astype("float"), percentiles, method="linear")
+    scaled = ((band.astype("float") - low) / (high - low + C.EPS)) * (srange[1] - srange[0]) + srange[0]
+    return np.clip(scaled, srange[0], srange[1])
+
+
+def test_known_values():
+    cube = np.array([[[0.0, 25.0, 50.0, 75.0, 100.0]]], dtype=np.float32)
+    result = _bands(cube, percentiles=(0, 100))
+    expected = (np.array([0.0, 25.0, 50.0, 75.0, 100.0]) / (100.0 + C.EPS)) * 2 - 1
+    _assert_same(result, expected.astype(np.float32).reshape(1, 1, 5))
+
+
+def test_scales_in_float64_before_casting():
+    cube = _cube(1)
+    expected = np.stack(
+        [_scale_in_float64(band, band).astype(np.float32) for band in cube]
+    )
+    _assert_same(_bands(cube), expected)
+
+    low, high = np.percentile(cube[0].astype("float"), (1, 99), method="linear")
+    in_float32 = ((cube[0] - np.float32(low)) / np.float32(high - low + C.EPS)) * 2 - 1
+    assert in_float32.tobytes() != expected[0].tobytes()
+
+
+def test_percentiles_ignore_pixels_outside_the_mask():
+    cube, mask = _cube(2), _partial_mask()
+    with_outliers = cube.copy()
+    with_outliers[:, ~mask] = 1e9
+    _assert_same(_bands(with_outliers, obsmask=mask), _bands(cube, obsmask=mask))
+
+
+@pytest.mark.parametrize("srange", [(-1, 1), (0, 1), (-3, 7)])
+def test_pixels_outside_the_mask_get_the_range_minimum(srange):
+    result = _bands(_cube(3), obsmask=_partial_mask(), srange=srange)
+    assert (result[:, ~_partial_mask()] == min(srange)).all()
+
+
+@pytest.mark.parametrize("percentiles,srange", [((1, 99), (-1, 1)), ((5, 50), (-3, 7))])
+def test_values_inside_the_mask_stay_in_range(percentiles, srange):
+    mask = _partial_mask()
+    result = _bands(_cube(4), obsmask=mask, percentiles=percentiles, srange=srange)
+    assert result.min() >= srange[0]
+    assert result.max() <= srange[1]
+    assert result[:, mask].min() == srange[0]
+    assert result[:, mask].max() == srange[1]
+
+
+def test_constant_band_gets_the_range_minimum():
+    cube = _cube(5)
+    cube[0] = 17.5
+    assert (_bands(cube)[0] == -1).all()
+
+
+@pytest.mark.filterwarnings("ignore::RuntimeWarning")
+def test_nan_inside_the_mask_only_affects_its_band():
+    cube, mask = _cube(6), _partial_mask()
+    clean = _bands(cube, obsmask=mask)
+    cube[1, 12, 12] = np.nan
+    result = _bands(cube, obsmask=mask)
+    assert np.isnan(result[1, mask]).all()
+    assert (result[1, ~mask] == -1).all()
+    _assert_same(result[[0, 2, 3]], clean[[0, 2, 3]])
+
+
+@pytest.mark.filterwarnings("ignore::RuntimeWarning")
+def test_nan_outside_the_mask_is_replaced():
+    cube, mask = _cube(7), _partial_mask()
+    clean = _bands(cube, obsmask=mask)
+    cube[2, ~mask] = np.nan
+    _assert_same(_bands(cube, obsmask=mask), clean)
+
+
+@pytest.mark.parametrize("use_mask", [True, False])
+@pytest.mark.parametrize("block_rows", [1, 7, 1000])
+def test_block_size_does_not_change_results(block_rows, use_mask):
+    cube = _cube(8)
+    obsmask = _partial_mask() if use_mask else None
+    _assert_same(
+        _bands(cube, obsmask=obsmask, block_rows=block_rows),
+        _bands(cube, obsmask=obsmask, block_rows=SHAPE[1]),
+    )
+
+
+@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+def test_result_keeps_the_input_dtype(dtype):
+    assert _bands(_cube(9, dtype)).dtype == dtype
+
+
+def test_normalize_datacube_matches_normalize_bands():
+    cube, mask = _cube(10), _partial_mask()
+    _assert_same(utils.normalize_datacube(cube, obsmask=mask), _bands(cube, obsmask=mask))
+    _assert_same(utils.normalize_datacube(cube), _bands(cube))
+
+
+def test_empty_obsmask_raises():
+    with pytest.raises(IndexError):
+        _bands(_cube(11), obsmask=np.zeros(SHAPE[1:], dtype=bool))
+
+
+@pytest.mark.filterwarnings("ignore::RuntimeWarning")
+def test_input_not_mutated():
+    mask = _partial_mask()
+    cube = _cube(12)
+    cube[1, 12, 12] = np.nan
+    cube_before, mask_before = cube.copy(), mask.copy()
+    _bands(cube, obsmask=mask)
+    _assert_same(cube, cube_before)
+    _assert_same(mask, mask_before)
+
+
+def test_result_does_not_share_memory_with_input():
+    cube = _cube(13)
+    assert not np.shares_memory(_bands(cube), cube)
+
+
+def test_bands_subset_follows_the_requested_order():
+    cube, mask = _cube(14), _partial_mask()
+    _assert_same(
+        utils.normalize_bands(cube, [2, 0], obsmask=mask),
+        _bands(cube[[2, 0]], obsmask=mask),
+    )
+
+
+def test_predictor_order_matches_data_get():
+    cube, mask = _cube(15), _partial_mask()
+    data = Data(cube, ["blue", "green", "red", "nir"])
+    predictors = ["nir", "blue", "red"]
+    _assert_same(
+        utils.normalize_bands(data.data, [data.bands.index(p) for p in predictors], obsmask=mask),
+        _bands(data.get(predictors), obsmask=mask),
+    )
diff --git a/tests/test_unet_classify.py b/tests/test_unet_classify.py
new file mode 100644
index 0000000..97889bd
--- /dev/null
+++ b/tests/test_unet_classify.py
@@ -0,0 +1,66 @@
+from types import SimpleNamespace
+
+import numpy as np
+import pytest
+import torch
+
+from fmask import utils
+from fmask.satellite import Data
+from fmask.unetlib import UNet
+
+BANDS = ["blue", "green", "red", "nir", "swir1"]
+PREDICTORS = ["nir", "blue", "swir1"]
+CLASSES = ["noncloud", "cloud"]
+SHAPE = (70, 83)
+
+
+def _stub(batch_size=3):
+    torch.manual_seed(0)
+    cube = (np.random.default_rng(0).random((len(BANDS), *SHAPE)) * 3000).astype(np.float32)
+    obsmask = np.ones(SHAPE, dtype=bool)
+    obsmask[:6] = False
+    obsmask[:, -5:] = False
+    return SimpleNamespace(
+        image=SimpleNamespace(data=Data(cube, BANDS), obsmask=obsmask),
+        predictors=PREDICTORS,
+        classes=CLASSES,
+        patch_size=32,
+        patch_stride_classify=24,
+        batch_size=batch_size,
+        device="cpu",
+        # A 1x1 kernel scores each pixel alone, so the patch layout cannot affect the answer.
+        model=torch.nn.Conv2d(len(PREDICTORS), len(CLASSES), 1),
+    )
+
+
+def _whole_image_probs(stub, normalize):
+    data = stub.image.data
+    indices = [data.bands.index(predictor) for predictor in PREDICTORS]
+    if normalize:
+        cube = utils.normalize_bands(data.data, indices, obsmask=stub.image.obsmask)
+    else:
+        cube = data.data[indices]
+    with torch.no_grad():
+        return torch.softmax(stub.model(torch.from_numpy(cube)[None]), dim=1)[0].numpy()
+
+
+@pytest.mark.parametrize("batch_size", [1, 3, 500], ids=["single", "short_last_batch", "one_batch"])
+@pytest.mark.parametrize("normalize", [True, False])
+@pytest.mark.parametrize("probability", ["none", "default", "cloud"])
+def test_classify_matches_per_pixel_scores(probability, normalize, batch_size):
+    stub = _stub(batch_size)
+    cube_before = stub.image.data.data.copy()
+    mask = stub.image.obsmask
+    probs = _whole_image_probs(stub, normalize)
+
+    image_class, image_prob = UNet.classify(stub, probability=probability, normalize=normalize)
+
+    assert image_class.dtype == np.uint8
+    assert np.array_equal(image_class[mask], probs.argmax(axis=0)[mask])
+    if probability == "none":
+        assert image_prob is None
+    else:
+        assert image_prob.dtype == np.float32
+        expected = probs.max(axis=0) if probability == "default" else probs[1]
+        assert np.allclose(image_prob[mask], expected[mask], rtol=1e-5, atol=0)
+    assert stub.image.data.data.tobytes() == cube_before.tobytes()

From 71893faf577a750c96350b53f5bea6387c21b802 Mon Sep 17 00:00:00 2001
From: Chris Holden <ceholden@gmail.com>
Date: Wed, 4 Nov 2026 12:53:07 -0500
Subject: [PATCH 03/11] perf: sample LightGBM pixels without copying the cube

Unit tests confirm the model receives the same matrix as v5.0.1-2; final masks
and timing are to be evaluated on AWS Batch.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/fmask/lightgbmlib.py  |  4 +-
 src/fmask/satellite.py    | 16 +++++++
 tests/test_data_sample.py | 89 +++++++++++++++++++++++++++++++++++++++
 3 files changed, 107 insertions(+), 2 deletions(-)
 create mode 100644 tests/test_data_sample.py

diff --git a/src/fmask/lightgbmlib.py b/src/fmask/lightgbmlib.py
index fbd40fc..e558eaa 100644
--- a/src/fmask/lightgbmlib.py
+++ b/src/fmask/lightgbmlib.py
@@ -554,9 +554,9 @@ def classify(self, probability="none", subsampling_mask = None, base = False) ->
             sample_image_col = sample_image_col_all[idx:idx_end]
             if point_predictors: # image-based predictors + latitude + longitude
                 n_pixels = len(sample_image_row)
-                pro_pred = self.model.predict_proba((np.concatenate((self.image.data.get(image_predictors)[:, sample_image_row, sample_image_col], np.full((1, n_pixels), self.image.lat_center, dtype=np.float32), np.full((1, n_pixels), self.image.lon_center, dtype=np.float32)), axis=0)).T, n_jobs=0)
+                pro_pred = self.model.predict_proba((np.concatenate((self.image.data.sample(image_predictors, sample_image_row, sample_image_col), np.full((1, n_pixels), self.image.lat_center, dtype=np.float32), np.full((1, n_pixels), self.image.lon_center, dtype=np.float32)), axis=0)).T, n_jobs=0)
             else: # only image-based predictors
-                pro_pred = self.model.predict_proba((self.image.data.get(image_predictors)[:,sample_image_row, sample_image_col]).T, n_jobs=0)
+                pro_pred = self.model.predict_proba((self.image.data.sample(image_predictors, sample_image_row, sample_image_col)).T, n_jobs=0)
             
             label_pred = labs_model[np.argmax(pro_pred, axis=1)]
 
diff --git a/src/fmask/satellite.py b/src/fmask/satellite.py
index 8e8033a..ad3f90b 100644
--- a/src/fmask/satellite.py
+++ b/src/fmask/satellite.py
@@ -117,6 +117,22 @@ def copyget(self, band):
             else:
                 return None
 
+    def sample(self, bands, rows, cols):
+        """
+        Returns the values of the specified bands at the given pixels, without copying whole bands.
+
+        Args:
+            bands (list): The names of the bands, in output order.
+            rows (ndarray): The row of each pixel.
+            cols (ndarray): The column of each pixel.
+
+        Returns:
+            ndarray: A new array of shape (len(bands), len(rows)).
+
+        """
+        indices = np.array([self.bands.index(band) for band in bands])
+        return self.data[indices[:, None], rows, cols]
+
     def __init__(self, data, bands):
         """
         Initializes a new instance of the Data class.
diff --git a/tests/test_data_sample.py b/tests/test_data_sample.py
new file mode 100644
index 0000000..81166ac
--- /dev/null
+++ b/tests/test_data_sample.py
@@ -0,0 +1,89 @@
+from types import SimpleNamespace
+
+import numpy as np
+import pytest
+
+from fmask.lightgbmlib import LightGBM
+from fmask.satellite import Data
+
+BANDS = ["blue", "green", "red", "nir"]
+
+
+def _data():
+    cube = np.random.default_rng(0).random((4, 30, 40)).astype(np.float32)
+    return Data(cube, BANDS)
+
+
+def _coords():
+    rng = np.random.default_rng(1)
+    rows, cols = rng.integers(0, 30, 200), rng.integers(0, 40, 200)
+    rows[:20], cols[:20] = rows[20:40], cols[20:40]
+    return rows, cols
+
+
+def test_sample_matches_get():
+    data = _data()
+    bands = ["nir", "blue", "red"]
+    rows, cols = _coords()
+    expected = data.get(bands)[:, rows, cols]
+    actual = data.sample(bands, rows, cols)
+    assert actual.dtype == expected.dtype
+    assert actual.shape == (3, 200)
+    assert np.array_equal(actual, expected)
+    assert actual.T.tobytes() == expected.T.tobytes()
+
+
+def test_sample_is_copy():
+    data = _data()
+    before = data.data.copy()
+    rows, cols = _coords()
+    data.sample(["nir", "blue"], rows, cols)[:] = -1
+    assert np.array_equal(data.data, before)
+
+
+class _RecordingModel:
+    classes_ = np.array([0, 1])
+
+    def __init__(self):
+        self.seen = []
+
+    def predict_proba(self, matrix, n_jobs=0):
+        self.seen.append(np.array(matrix))
+        return np.tile([0.25, 0.75], (matrix.shape[0], 1))
+
+
+@pytest.mark.parametrize("use_lat_lon", [False, True])
+def test_classify_gives_model_the_original_matrix(use_lat_lon):
+    data = _data()
+    predictors = ["nir", "blue", "red"] + (["latitude", "longitude"] if use_lat_lon else [])
+    obsmask = np.ones((30, 40), dtype=bool)
+    obsmask[3:7] = False
+    stub = SimpleNamespace(
+        image=SimpleNamespace(
+            data=data, obsmask=obsmask, shape=(30, 40), lat_center=41.5, lon_center=-72.25
+        ),
+        subsampling_size=2,
+        subsampling_min=0,
+        predictors=predictors,
+        classes=["noncloud", "cloud"],
+        model=_RecordingModel(),
+    )
+    LightGBM.classify(stub, probability="cloud")
+
+    mask = np.zeros((30, 40), dtype=bool)
+    mask[::2, ::2] = True
+    rows, cols = np.where(mask & obsmask)
+    expected = data.get(["nir", "blue", "red"])[:, rows, cols]
+    if use_lat_lon:
+        n_pixels = len(rows)
+        expected = np.concatenate(
+            (
+                expected,
+                np.full((1, n_pixels), 41.5, dtype=np.float32),
+                np.full((1, n_pixels), -72.25, dtype=np.float32),
+            ),
+            axis=0,
+        )
+    seen = stub.model.seen[0]
+    assert seen.dtype == expected.dtype
+    assert seen.tobytes() == expected.T.tobytes()

From 42fe40f251776fabb62ad5bb50122d4fb3d522db Mon Sep 17 00:00:00 2001
From: Chris Holden <ceholden@gmail.com>
Date: Wed, 4 Nov 2026 12:53:07 -0500
Subject: [PATCH 04/11] perf: read Landsat bands concurrently

Spectral bands are read and converted in up to four threads, and the sun angle
cosine is computed once. Observation mask, band and saturation updates stay in
the main thread in the original order. Byte-for-byte comparison with v5.0.1-2 on synthetic data passed on
perf/equivalence-tests; final masks and timing are to be evaluated on AWS Batch.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/fmask/satellite.py     | 198 +++++++++++---------
 tests/test_landsat_read.py | 358 +++++++++++++++++++++++++++++++++++++
 2 files changed, 469 insertions(+), 87 deletions(-)
 create mode 100644 tests/test_landsat_read.py

diff --git a/src/fmask/satellite.py b/src/fmask/satellite.py
index ad3f90b..da3c33f 100644
--- a/src/fmask/satellite.py
+++ b/src/fmask/satellite.py
@@ -26,6 +26,8 @@
 # pylint: disable=invalid-name
 # pylint: too-many-branches
 # pylint: too-many-statements
+from collections import deque
+from concurrent.futures import ThreadPoolExecutor
 import os
 from pathlib import Path
 import xml.etree.ElementTree as ET
@@ -53,6 +55,21 @@
 )
 from . import constant as C
 
+READ_THREADS = 4
+
+LANDSAT_SPECTRAL_READS = (
+    ("coastal", "TOA", False),
+    ("blue", "TOA", True),
+    ("green", "TOA", True),
+    ("red", "TOA", True),
+    ("nir", "TOA", False),
+    ("swir1", "TOA", False),
+    ("swir2", "TOA", False),
+    ("tirs1", "BT", False),
+    ("tirs2", "BT", False),
+    ("cirrus", "TOA", False),
+)
+
 
 class Data:
     """
@@ -149,6 +166,43 @@ def __init__(self, data, bands):
 class Satellite:
     """Base class for satellite data object"""
 
+    def _read_bands(self, requests):
+        """Read bands concurrently, each from its own dataset, without changing this object.
+
+        Args:
+            requests (list): Keyword arguments of _read_band_values for each band.
+
+        Yields:
+            The result of _read_band_values for each request, in request order. A new read
+            starts only when fewer than READ_THREADS are unconsumed, so bands read quickly do
+            not pile up in memory behind a slow one.
+        """
+        if C.MSG_FULL:
+            for request in requests:
+                print(f">>> loading {request['band'].lower()} in {request['level'].lower()}")
+        # NumPy error settings are per thread, so workers would otherwise lose the caller's.
+        error_state = np.geterr()
+
+        def read(request):
+            with np.errstate(**error_state):
+                return self._read_band_values(**request)
+
+        with ThreadPoolExecutor(max_workers=READ_THREADS) as pool:
+            pending = deque()
+            for request in requests:
+                if len(pending) == READ_THREADS:
+                    yield pending.popleft().result()
+                pending.append(pool.submit(read, request))
+            while pending:
+                yield pending.popleft().result()
+
+    def _update_obsmask(self, valid):
+        """Overlap the extent of the observation with the valid pixels of a band."""
+        if self.obsmask is None:
+            self.obsmask = valid
+        else:
+            self.obsmask = np.logical_and(self.obsmask, valid)
+
     def __init__(self, folder):
         self.folder: str = folder
         self.name: str = Path(folder).stem
@@ -502,39 +556,51 @@ def read_band(self, band, level="TOA", sza=None, maskupdate=True, saturation = F
         """
         if C.MSG_FULL:
             print(f">>> loading {band.lower()} in {level.lower()}")
-        # specify the band information
-        band = band.upper()  # uppercase
-        # band_id = self.bands[self.bands['NAME'] == band]['ID'].values[0]
+        cos_sza = None if sza is None else np.cos(sza)
+        image, satu, valid = self._read_band_values(band, level, cos_sza, saturation)
+        # update the extent of the observation when new band is touched
+        # in this way, we do not need to load all the bands once time
+        if maskupdate:
+            self._update_obsmask(valid)
+        # append the profile of the geotiff, which will be used as the reference to auxiliary data warping and mask saving.
+        if profile:
+            with rasterio.open(self._band_path(band.upper())) as src:
+                self.profile = src.meta.copy()
+                self.profilefull = src  # save the full profile for the further use
+        return image, satu
+
+    def _band_path(self, band):
+        """Return the path of the band named in uppercase."""
         band_id = self.get_band_id(band)
-        path_band = os.path.join(
+        return os.path.join(
             self.folder,
             self.metadata["LEVEL1_PROCESSING_RECORD"][f"FILE_NAME_BAND_{band_id}"],
         )
 
+    def _read_band_values(self, band, level="TOA", cos_sza=None, saturation=False):
+        """Read and convert a band without changing this object, so it can run in a thread.
+
+        Args:
+            band (str): Band name
+            level (str, optional): Level of the data. Defaults to 'TOA'.
+            cos_sza (2d array, optional): Cosine of the Solar Zenith Angle. Defaults to None.
+            saturation (bool, optional): Return the saturated layer of this band. Defaults to False.
+
+        Returns:
+            tuple: image of the band, saturated layer or None, and valid pixels before conversion
+        """
+        # specify the band information
+        band = band.upper()  # uppercase
+        band_id = self.get_band_id(band)
+
         # read the raw band
-        with rasterio.open(path_band) as src:
+        with rasterio.open(self._band_path(band)) as src:
             image = src.read(1)  # only one band
+            valid = image != src.nodata  # should be 0 when we test the Landsat Collection 2.
         if saturation:
             satu = image == np.iinfo(image.dtype).max # 255 for Landsat 4-7, and 65535 for Landsat 8, according to the data's type
         else:
             satu = None
-        # update LANDSAT object before converting the pixel value
-        # update the extent of the observation when new band is touched
-        # in this way, we do not need to load all the bands once time
-        if maskupdate:
-            if self.obsmask is None:  # initlize the mask of observation extent
-                self.obsmask = (
-                    image != src.nodata
-                )  # should be 0 when we test the Landsat Collection 2.
-            else:  # update by overlapping it with new band's extent
-                self.obsmask = np.logical_and(
-                    self.obsmask, image != src.nodata
-                )  # should be 0 when we test the Landsat Collection 2.
-
-        # append the profile of the geotiff, which will be used as the reference to auxiliary data warping and mask saving.
-        if profile:
-            self.profile = src.meta.copy()
-            self.profilefull = src  # save the full profile for the further use
 
         ## Converting
         image = image.astype(np.float32)
@@ -550,11 +616,11 @@ def read_band(self, band, level="TOA", sza=None, maskupdate=True, saturation = F
                 ]
             )
             image = image * REFLECTANCE_MULT + REFLECTANCE_ADD
-            if sza is None:
-                sza = self.read_angle(
-                    angle="SOLAR_ZENITH", unit="RADIAN"
+            if cos_sza is None:
+                cos_sza = np.cos(
+                    self.read_angle(angle="SOLAR_ZENITH", unit="RADIAN")
                 )  # when we do not speficfy the angle, we can read it
-            image = image / np.cos(sza)  # with a correction for the sun angle
+            image = image / cos_sza  # with a correction for the sun angle
         elif level == "BT":  # Conversion to Top of Atmosphere Brightness Temperature
             # convert to TOA radiance
             RADIANCE_MULT = float(
@@ -581,7 +647,7 @@ def read_band(self, band, level="TOA", sza=None, maskupdate=True, saturation = F
             image = K2_CONSTANT / np.log(K1_CONSTANT / image + 1)  # np.log is ln
             # convert from Kelvin to Celcius
             image = image - 273.15
-        return image, satu  # return final result
+        return image, satu, valid  # return final result
 
     def read_angle(self, angle, unit="RADIAN", profile=False):
         """read angle bands
@@ -633,68 +699,26 @@ def read_datacube(self, predictors, maskupdate=True):
             (len(predictors), self.profile["height"], self.profile["width"]),
             dtype=np.float32,
         )
+        for band, _, saturation in LANDSAT_SPECTRAL_READS:
+            if saturation and band not in predictors:
+                raise ValueError(f"{band} is required for the saturation layer")
         # Spectral bands
-        predictor = "coastal"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], _ = self.read_band(
-                band=predictor, level="TOA", sza=sza, maskupdate=maskupdate
-            )
-        predictor = "blue"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], satu_blue = self.read_band(
-                band=predictor, level="TOA", sza=sza, maskupdate=maskupdate, saturation = True
-            )
+        cos_sza = np.cos(sza)
+        del sza
+        requests = [
+            {"band": band, "level": level, "cos_sza": cos_sza, "saturation": saturation}
+            for band, level, saturation in LANDSAT_SPECTRAL_READS
+            if band in predictors
+        ]
         # saturation for visible bands after the calculation of whiteness
-        statu = BitLayer(satu_blue.shape)
-        statu.append(satu_blue) # set bit 0 for blue
-        del satu_blue
-
-        predictor = "green"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], satu_green = self.read_band(
-                band=predictor, level="TOA", sza=sza, maskupdate=maskupdate, saturation = True
-            )
-        statu.append(satu_green) # set bit 1 for green
-        del satu_green
-        
-        predictor = "red"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], satu_red  = self.read_band(
-                band=predictor, level="TOA", sza=sza, maskupdate=maskupdate, saturation = True
-            )
-        statu.append(satu_red) # set bit 2 for red
-        del satu_red
-    
-        predictor = "nir"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], _  = self.read_band(
-                band=predictor, level="TOA", sza=sza, maskupdate=maskupdate
-            )
-        predictor = "swir1"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], _ = self.read_band(
-                band=predictor, level="TOA", sza=sza, maskupdate=maskupdate
-            )
-        predictor = "swir2"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], _ = self.read_band(
-                band=predictor, level="TOA", sza=sza, maskupdate=maskupdate
-            )
-        predictor = "tirs1"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], _ = self.read_band(
-                band=predictor, level="BT", maskupdate=maskupdate
-            )
-        predictor = "tirs2"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], _ = self.read_band(
-                band=predictor, level="BT", maskupdate=maskupdate
-            )
-        predictor = "cirrus"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], _ = self.read_band(
-                band=predictor, level="TOA", sza=sza, maskupdate=maskupdate
-            )
+        statu = BitLayer(bands.shape[1:])
+        for request, (image, satu, valid) in zip(requests, self._read_bands(requests)):
+            if maskupdate:
+                self._update_obsmask(valid)
+            bands[predictors.index(request["band"]), :, :] = image
+            if satu is not None:
+                statu.append(satu) # set bits 0, 1 and 2 for blue, green and red, in request order
+        del image, satu, valid, cos_sza, requests, request
 
         # spectral index
         predictor = "hot"
diff --git a/tests/test_landsat_read.py b/tests/test_landsat_read.py
new file mode 100644
index 0000000..ea6392b
--- /dev/null
+++ b/tests/test_landsat_read.py
@@ -0,0 +1,358 @@
+import gc
+import sys
+import threading
+import time
+import weakref
+
+import numpy as np
+import pandas
+import pytest
+import rasterio
+from rasterio.transform import from_origin
+
+from fmask import constant as C
+from fmask import satellite, utils
+from fmask.satellite import Landsat
+
+SHAPE = (48, 64)
+BAND_IDS = {
+    "COASTAL": "1", "BLUE": "2", "GREEN": "3", "RED": "4", "NIR": "5",
+    "SWIR1": "6", "SWIR2": "7", "CIRRUS": "9", "TIRS1": "10", "TIRS2": "11",
+}
+SPECTRAL = ["coastal", "blue", "green", "red", "nir", "swir1", "swir2", "tirs1", "tirs2", "cirrus"]
+
+
+def _write(path, data, dtype, nodata):
+    with rasterio.open(
+        path, "w", driver="GTiff", height=SHAPE[0], width=SHAPE[1], count=1, dtype=dtype,
+        crs="EPSG:32633", transform=from_origin(500000, 4000000, 30, 30), nodata=nodata,
+    ) as dst:
+        dst.write(data, 1)
+
+
+@pytest.fixture
+def scene(tmp_path):
+    rng = np.random.default_rng(0)
+    record, radiometric, thermal = {}, {}, {}
+    for offset, (name, band_id) in enumerate(BAND_IDS.items()):
+        data = rng.integers(5000, 50000, SHAPE, dtype=np.uint16)
+        data[offset, :] = 0
+        data[:, offset] = 0
+        data[-1, -1 - offset] = np.iinfo(np.uint16).max
+        filename = f"band{band_id}.tif"
+        _write(tmp_path / filename, data, "uint16", 0)
+        record[f"FILE_NAME_BAND_{band_id}"] = filename
+        radiometric[f"REFLECTANCE_MULT_BAND_{band_id}"] = str(2.0e-5 + offset * 1e-7)
+        radiometric[f"REFLECTANCE_ADD_BAND_{band_id}"] = str(-0.1 - offset * 1e-3)
+        radiometric[f"RADIANCE_MULT_BAND_{band_id}"] = str(3.342e-4)
+        radiometric[f"RADIANCE_ADD_BAND_{band_id}"] = str(0.1 + offset * 1e-2)
+        thermal[f"K1_CONSTANT_BAND_{band_id}"] = str(774.8853 + offset)
+        thermal[f"K2_CONSTANT_BAND_{band_id}"] = str(1321.0789 + offset)
+    angle = rng.integers(2000, 7000, SHAPE, dtype=np.int16)
+    _write(tmp_path / "sza.tif", angle, "int16", None)
+    record["FILE_NAME_ANGLE_SOLAR_ZENITH_BAND_4"] = "sza.tif"
+    return {
+        "folder": str(tmp_path),
+        "metadata": {
+            "LEVEL1_PROCESSING_RECORD": record,
+            "LEVEL1_RADIOMETRIC_RESCALING": radiometric,
+            "LEVEL1_THERMAL_CONSTANTS": thermal,
+        },
+        "bands": pandas.DataFrame(
+            {"NAME": list(BAND_IDS), "ID": list(BAND_IDS.values())}
+        ),
+    }
+
+
+def _make(cls, scene):
+    image = object.__new__(cls)
+    image.folder = scene["folder"]
+    image.metadata = scene["metadata"]
+    image.bands = scene["bands"]
+    image.obsmask = None
+    image.profile = None
+    image.profilefull = None
+    return image
+
+
+def _assert_same(actual, expected):
+    if expected is None:
+        assert actual is None
+        return
+    assert actual.dtype == expected.dtype
+    assert actual.shape == expected.shape
+    assert actual.tobytes() == expected.tobytes()
+
+
+def test_values_do_not_touch_self(scene):
+    image = _make(Landsat, scene)
+    sentinel_mask = np.ones(SHAPE, dtype=bool)
+    image.obsmask = sentinel_mask
+    image.profile = "profile"
+    image.profilefull = "profilefull"
+    _, _, valid = image._read_band_values("blue", "TOA")
+    assert image.obsmask is sentinel_mask
+    assert image.profile == "profile"
+    assert image.profilefull == "profilefull"
+    assert valid.dtype == np.bool_
+    assert not valid[1, :].any()
+
+
+def test_read_bands_order(scene, monkeypatch):
+    image = _make(Landsat, scene)
+    release = threading.Event()
+
+    def fake_values(self, band, **kwargs):
+        if band == "first":
+            assert release.wait(timeout=10)
+        else:
+            release.set()
+        return band
+
+    monkeypatch.setattr(Landsat, "_read_band_values", fake_values)
+    requests = [{"band": "first", "level": "TOA"}, {"band": "second", "level": "TOA"}]
+    assert list(image._read_bands(requests)) == ["first", "second"]
+
+
+def test_read_bands_limits_threads(scene, monkeypatch):
+    image = _make(Landsat, scene)
+    lock = threading.Lock()
+    active = [0, 0]
+    barrier = threading.Barrier(satellite.READ_THREADS)
+
+    def fake_values(self, band, **kwargs):
+        with lock:
+            active[0] += 1
+            active[1] = max(active)
+        barrier.wait(timeout=10)
+        with lock:
+            active[0] -= 1
+        return band
+
+    monkeypatch.setattr(Landsat, "_read_band_values", fake_values)
+    count = satellite.READ_THREADS * 3
+    list(image._read_bands([{"band": str(i), "level": "TOA"} for i in range(count)]))
+    assert active[1] == satellite.READ_THREADS
+
+
+def test_read_bands_propagates_errors(scene, monkeypatch):
+    image = _make(Landsat, scene)
+
+    def fake_values(self, band, **kwargs):
+        raise OSError("bad band")
+
+    monkeypatch.setattr(Landsat, "_read_band_values", fake_values)
+    with pytest.raises(OSError, match="bad band"):
+        list(image._read_bands([{"band": "blue", "level": "TOA"}]))
+
+
+def test_worker_threads_use_caller_error_state(scene, monkeypatch):
+    image = _make(Landsat, scene)
+    monkeypatch.setattr(Landsat, "_read_band_values", lambda self, band, **kw: np.geterr()["invalid"])
+    with np.errstate(invalid="ignore"):
+        assert list(image._read_bands([{"band": "blue", "level": "TOA"}])) == ["ignore"]
+    with np.errstate(invalid="raise"):
+        assert list(image._read_bands([{"band": "blue", "level": "TOA"}])) == ["raise"]
+
+
+def test_read_bands_releases_each_result_before_the_next(scene, monkeypatch):
+    image = _make(Landsat, scene)
+
+    class Value:
+        pass
+
+    refs = {}
+
+    def fake_values(self, band, **kwargs):
+        value = Value()
+        refs[band] = weakref.ref(value)
+        return value
+
+    monkeypatch.setattr(Landsat, "_read_band_values", fake_values)
+    results = image._read_bands([{"band": str(i), "level": "TOA"} for i in range(3)])
+    next(results)
+    gc.collect()
+    value = next(results)
+    gc.collect()
+    assert refs["0"]() is None
+    del value
+    list(results)
+
+
+def test_read_bands_limits_results_waiting_behind_a_slow_band(scene, monkeypatch):
+    image = _make(Landsat, scene)
+    started = []
+    started_during_slow_band = []
+
+    def fake_values(self, band, **kwargs):
+        started.append(band)
+        if band == "0":
+            time.sleep(0.3)
+            started_during_slow_band.append(len(started))
+        return band
+
+    monkeypatch.setattr(Landsat, "_read_band_values", fake_values)
+    count = satellite.READ_THREADS * 3
+    requests = [{"band": str(i), "level": "TOA"} for i in range(count)]
+    assert list(image._read_bands(requests)) == [str(i) for i in range(count)]
+    assert started_during_slow_band == [satellite.READ_THREADS]
+
+
+def test_read_datacube_releases_read_inputs_before_indices(scene, monkeypatch):
+    original_hot = satellite.hot
+    seen = {}
+
+    def spy(*args):
+        seen.update(sys._getframe(1).f_locals)
+        return original_hot(*args)
+
+    monkeypatch.setattr(satellite, "hot", spy)
+    _make(Landsat, scene).read_datacube(SPECTRAL + ["hot"])
+    assert seen
+    assert not {"request", "requests", "cos_sza", "saturated", "image", "satu", "valid"} & set(seen)
+
+
+def _raw(scene, band):
+    filename = scene["metadata"]["LEVEL1_PROCESSING_RECORD"][f"FILE_NAME_BAND_{BAND_IDS[band.upper()]}"]
+    with rasterio.open(f"{scene['folder']}/{filename}") as src:
+        return src.read(1)
+
+
+def _cos_sza(scene):
+    with rasterio.open(f"{scene['folder']}/sza.tif") as src:
+        return np.cos(np.deg2rad(src.read(1).astype(np.float32) / 100))
+
+
+def _expected_toa(scene, band):
+    band_id = BAND_IDS[band.upper()]
+    rescaling = scene["metadata"]["LEVEL1_RADIOMETRIC_RESCALING"]
+    image = _raw(scene, band).astype(np.float32)
+    image = image * float(rescaling[f"REFLECTANCE_MULT_BAND_{band_id}"])
+    image = image + float(rescaling[f"REFLECTANCE_ADD_BAND_{band_id}"])
+    return image / _cos_sza(scene)
+
+
+def _expected_bt(scene, band):
+    band_id = BAND_IDS[band.upper()]
+    rescaling = scene["metadata"]["LEVEL1_RADIOMETRIC_RESCALING"]
+    thermal = scene["metadata"]["LEVEL1_THERMAL_CONSTANTS"]
+    image = _raw(scene, band).astype(np.float32)
+    image = image * float(rescaling[f"RADIANCE_MULT_BAND_{band_id}"])
+    image = image + float(rescaling[f"RADIANCE_ADD_BAND_{band_id}"])
+    image = image + C.EPS
+    image = float(thermal[f"K2_CONSTANT_BAND_{band_id}"]) / np.log(
+        float(thermal[f"K1_CONSTANT_BAND_{band_id}"]) / image + 1
+    )
+    return image - 273.15
+
+
+def _expected_band(scene, band):
+    return _expected_bt(scene, band) if band.startswith("tirs") else _expected_toa(scene, band)
+
+
+def _expected_obsmask(scene, bands):
+    return np.logical_and.reduce([_raw(scene, band) != 0 for band in bands])
+
+
+@pytest.mark.parametrize("pass_sza", [True, False], ids=["sza_passed", "sza_read"])
+def test_read_band_converts_to_toa(scene, pass_sza):
+    image = _make(Landsat, scene)
+    sza = image.read_angle(angle="SOLAR_ZENITH", unit="RADIAN") if pass_sza else None
+    result, satu = image.read_band("red", level="TOA", sza=sza)
+    _assert_same(result, _expected_toa(scene, "red"))
+    assert satu is None
+
+
+def test_read_band_converts_to_brightness_temperature(scene):
+    result, satu = _make(Landsat, scene).read_band("tirs1", level="BT")
+    _assert_same(result, _expected_bt(scene, "tirs1"))
+    assert satu is None
+
+
+def test_read_band_returns_saturated_pixels(scene):
+    _, satu = _make(Landsat, scene).read_band("green", saturation=True)
+    _assert_same(satu, _raw(scene, "green") == np.iinfo(np.uint16).max)
+    assert satu.any()
+
+
+def test_read_band_overlaps_the_observation_mask(scene):
+    image = _make(Landsat, scene)
+    image.read_band("coastal")
+    _assert_same(image.obsmask, _raw(scene, "coastal") != 0)
+    image.read_band("blue")
+    _assert_same(image.obsmask, _expected_obsmask(scene, ["coastal", "blue"]))
+
+
+def test_read_band_can_leave_the_observation_mask(scene):
+    image = _make(Landsat, scene)
+    image.read_band("blue", maskupdate=False)
+    assert image.obsmask is None
+
+
+def test_read_band_profile(scene):
+    image = _make(Landsat, scene)
+    image.read_band("blue", profile=True)
+    with rasterio.open(f"{scene['folder']}/band2.tif") as src:
+        assert image.profile == src.meta
+    assert image.profilefull.name.endswith("band2.tif")
+    assert image.profilefull.closed
+
+
+def test_obsmask_all_bands(scene):
+    image = _make(Landsat, scene)
+    image.read_datacube(list(SPECTRAL))
+    _assert_same(image.obsmask, _expected_obsmask(scene, SPECTRAL))
+    assert not image.obsmask.all()
+    assert image.obsmask.any()
+
+
+@pytest.mark.parametrize("maskupdate", [True, False])
+@pytest.mark.parametrize("predictors", [SPECTRAL, list(reversed(SPECTRAL))], ids=["forward", "reversed"])
+def test_read_datacube_matches_each_band(scene, predictors, maskupdate):
+    image = _make(Landsat, scene)
+    data, _ = image.read_datacube(list(predictors), maskupdate=maskupdate)
+    expected = np.stack([_expected_band(scene, band) for band in predictors])
+    if maskupdate:
+        expected[:, ~_expected_obsmask(scene, SPECTRAL)] = 0
+    else:
+        assert image.obsmask is None
+    _assert_same(data.data, expected)
+    assert data.bands == list(predictors)
+
+
+def test_read_datacube_sets_saturation_bits_in_blue_green_red_order(scene):
+    _, statu = _make(Landsat, scene).read_datacube(list(SPECTRAL))
+    assert statu.index == 2
+    for bit, band in enumerate(["blue", "green", "red"]):
+        _assert_same(statu.get(bit), _raw(scene, band) == np.iinfo(np.uint16).max)
+    assert not np.array_equal(statu.get(0), statu.get(1))
+
+
+def test_read_datacube_computes_indices_from_the_bands_read(scene):
+    data, _ = _make(Landsat, scene).read_datacube(SPECTRAL + ["hot", "ndvi"])
+    mask = _expected_obsmask(scene, SPECTRAL)
+    blue, red, nir = (_expected_toa(scene, band) for band in ("blue", "red", "nir"))
+    for name, expected in (("hot", utils.hot(blue, red)), ("ndvi", utils.ndvi(red, nir))):
+        expected = expected.astype(np.float32).copy()
+        expected[~mask] = 0
+        _assert_same(data.data[data.bands.index(name)], expected)
+
+
+def test_read_datacube_replaces_nan_values(scene, monkeypatch):
+    original = Landsat._read_band_values
+
+    def with_nan(self, band, **kwargs):
+        image, satu, valid = original(self, band, **kwargs)
+        image[20, 20] = np.nan
+        return image, satu, valid
+
+    monkeypatch.setattr(Landsat, "_read_band_values", with_nan)
+    data, _ = _make(Landsat, scene).read_datacube(list(SPECTRAL))
+    assert (data.data[:, 20, 20] == np.float32(C.EPS)).all()
+
+
+@pytest.mark.parametrize("missing", ["blue", "green", "red"])
+def test_read_datacube_without_visible_band_raises(scene, missing):
+    predictors = [band for band in SPECTRAL if band != missing]
+    with pytest.raises(ValueError, match=missing):
+        _make(Landsat, scene).read_datacube(predictors)

From 2f5267ccf7ba94e4ebfba886f752cd6c69383828 Mon Sep 17 00:00:00 2001
From: Chris Holden <ceholden@gmail.com>
Date: Wed, 4 Nov 2026 12:53:07 -0500
Subject: [PATCH 05/11] perf: reuse UNet batch buffer and use inference mode

Byte-for-byte comparison with v5.0.1-2 on synthetic data passed on
perf/equivalence-tests; final masks and timing are to be evaluated on AWS Batch.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/fmask/unetlib.py        | 31 +++++++++++++++++++++----
 tests/test_unet_classify.py | 46 ++++++++++++++++++++++++++++++++++++-
 2 files changed, 71 insertions(+), 6 deletions(-)

diff --git a/src/fmask/unetlib.py b/src/fmask/unetlib.py
index 61c1b0c..2738e60 100644
--- a/src/fmask/unetlib.py
+++ b/src/fmask/unetlib.py
@@ -159,6 +159,23 @@ def __len__(self):
         return len(self.ids)
 
 
+def fill_batch(buffer, datacube, anchors, patch_size):
+    """Copy the datacube patches at the anchors into a buffer and return the filled part.
+
+    Args:
+        buffer (tensor): Reused tensor of (batch, channel, height, width), with room for all anchors.
+        datacube (tensor): The datacube of (channel, height, width).
+        anchors (list): The (row, column, _, _) of each patch.
+        patch_size (int): The height and width of a patch.
+
+    Returns:
+        tensor: The part of the buffer holding this batch.
+    """
+    for index, (row, col, _, _) in enumerate(anchors):
+        buffer[index] = datacube[:, row:row + patch_size, col:col + patch_size]
+    return buffer[: len(anchors)]
+
+
 class UNet(object):
     """_summary_"""
 
@@ -580,15 +597,19 @@ def classify(self, probability="none", normalize=True, shift=True) -> tuple:
             image_prob = np.zeros((h, w), dtype=np.float32)
             image_prob_full = np.zeros((h, w), dtype=np.float32)
         self.model.eval()
-        
+
+        batch_buffer = torch.empty(
+            (min(self.batch_size, len(offanchors)), datacube.shape[0], self.patch_size, self.patch_size),
+            dtype=datacube.dtype,
+            device=datacube.device,
+        )
         for ibatch in range(0, len(offanchors), self.batch_size):
             # obtain the batch anchors
             batch_anchors = offanchors[ibatch : np.min((ibatch + self.batch_size, len(offanchors)))]  # maximum batch size
-            # obtain the batch tensor of the datacube
-            # batch_tensor = torch.stack([datacube[:, r:r+self.patch_size, c:c+self.patch_size] for r, c, _, _ in batch_anchors])
             # make the classification
-            with torch.no_grad(): # one line to finish the classification, in order to save the memory
-                batch_probs = torch.nn.functional.softmax(self.model(torch.stack([datacube[:, r:r+self.patch_size, c:c+self.patch_size] for r, c, _, _ in batch_anchors])), dim=1)
+            batch_tensor = fill_batch(batch_buffer, datacube, batch_anchors, self.patch_size)
+            with torch.inference_mode(): # one line to finish the classification, in order to save the memory
+                batch_probs = torch.nn.functional.softmax(self.model(batch_tensor), dim=1)
             # batch processing the class and prob layer
             if prob_index == -2: # -2 means the top class will be extracted by this model, but not the probability
                 _, batch_classes = batch_probs.topk(1, dim=1)  # top label with top prob.
diff --git a/tests/test_unet_classify.py b/tests/test_unet_classify.py
index 97889bd..17ecc37 100644
--- a/tests/test_unet_classify.py
+++ b/tests/test_unet_classify.py
@@ -6,7 +6,7 @@
 
 from fmask import utils
 from fmask.satellite import Data
-from fmask.unetlib import UNet
+from fmask.unetlib import UNet, fill_batch
 
 BANDS = ["blue", "green", "red", "nir", "swir1"]
 PREDICTORS = ["nir", "blue", "swir1"]
@@ -64,3 +64,47 @@ def test_classify_matches_per_pixel_scores(probability, normalize, batch_size):
         expected = probs.max(axis=0) if probability == "default" else probs[1]
         assert np.allclose(image_prob[mask], expected[mask], rtol=1e-5, atol=0)
     assert stub.image.data.data.tobytes() == cube_before.tobytes()
+
+
+def _datacube():
+    return torch.from_numpy(
+        np.random.default_rng(1).random((3, 70, 83)).astype(np.float32)
+    ).to("cpu", dtype=torch.float)
+
+
+def _stack(datacube, anchors, size):
+    return torch.stack([datacube[:, r:r + size, c:c + size] for r, c, _, _ in anchors])
+
+
+def _anchors(count):
+    return [(i * 2, i * 3, 0, 0) for i in range(count)]
+
+
+def test_fill_batch_matches_stack():
+    datacube = _datacube()
+    anchors = _anchors(4)
+    buffer = torch.empty((4, 3, 32, 32), dtype=datacube.dtype)
+    actual = fill_batch(buffer, datacube, anchors, 32)
+    expected = _stack(datacube, anchors, 32)
+    assert actual.dtype == expected.dtype
+    assert actual.shape == expected.shape
+    assert actual.numpy().tobytes() == expected.numpy().tobytes()
+
+
+def test_short_last_batch():
+    datacube = _datacube()
+    buffer = torch.empty((4, 3, 32, 32), dtype=datacube.dtype)
+    fill_batch(buffer, datacube, _anchors(4), 32)
+    anchors = _anchors(2)
+    actual = fill_batch(buffer, datacube, anchors, 32)
+    expected = _stack(datacube, anchors, 32)
+    assert actual.shape == (2, 3, 32, 32)
+    assert actual.numpy().tobytes() == expected.numpy().tobytes()
+
+
+def test_fill_batch_does_not_change_datacube():
+    datacube = _datacube()
+    before = datacube.clone()
+    buffer = torch.empty((4, 3, 32, 32), dtype=datacube.dtype)
+    fill_batch(buffer, datacube, _anchors(4), 32)
+    assert torch.equal(datacube, before)

From 075c94429466ecad63a0a048246077885c62c87a Mon Sep 17 00:00:00 2001
From: Chris Holden <ceholden@gmail.com>
Date: Wed, 4 Nov 2026 16:35:15 -0500
Subject: [PATCH 06/11] perf: read Sentinel-2 bands concurrently

Spectral bands are read and converted in up to four threads, and saturation is
read from the DN in the same step. Observation mask, band and saturation updates
stay in the main thread in the original order.

Before rebasing on the saturation fix, a byte-for-byte comparison with v5.0.1-2
on synthetic data passed on perf/equivalence-tests. After it, read_datacube output
matched the fix branch byte for byte on a synthetic scene with saturated pixels.
Final masks and timing are to be evaluated on AWS Batch.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/fmask/satellite.py       | 153 +++++++++++-----------------
 tests/test_sentinel2_read.py | 186 +++++++++++++++++++++++++++++++++++
 2 files changed, 246 insertions(+), 93 deletions(-)
 create mode 100644 tests/test_sentinel2_read.py

diff --git a/src/fmask/satellite.py b/src/fmask/satellite.py
index da3c33f..6c67c27 100644
--- a/src/fmask/satellite.py
+++ b/src/fmask/satellite.py
@@ -57,6 +57,22 @@
 
 READ_THREADS = 4
 
+SENTINEL2_SPECTRAL_READS = (
+    ("coastal", False),
+    ("blue", True),
+    ("green", True),
+    ("red", True),
+    ("vre1", False),
+    ("vre2", False),
+    ("vre3", False),
+    ("wnir", False),
+    ("nir", False),
+    ("wv", False),
+    ("swir1", False),
+    ("swir2", False),
+    ("cirrus", False),
+)
+
 LANDSAT_SPECTRAL_READS = (
     ("coastal", "TOA", False),
     ("blue", "TOA", True),
@@ -957,15 +973,40 @@ def read_band(self, band, level="TOA", sza=None, maskupdate=True, saturation = F
 
         if C.MSG_FULL:
             print(f">>> loading {band.lower()} in {level.lower()}")
+        image, satu, valid = self._read_band_values(band, level, saturation)
+        if maskupdate:
+            self._update_obsmask(valid)
+        # append the profile of the geotiff, which will be used as the reference to auxiliary data warping and mask saving.
+        if profile:
+            with rasterio.open(self._band_path(band.upper())) as src:
+                self.profile = src.meta.copy()
+                self.profilefull = src  # save the full profile for the further use
+        return image, satu
+
+    def _band_path(self, band):
+        """Return the path of the band named in uppercase."""
+        band_file = self.bands[self.bands["NAME"] == band]["FILE"].values[0]
+        return os.path.join(self.folder, band_file) + ".jp2"
+
+    def _read_band_values(self, band, level="TOA", saturation=False):
+        """Read and convert a band without changing this object, so it can run in a thread.
+
+        Args:
+            band (str): The name of the band to read.
+            level (str, optional): The level of the image data. Defaults to "TOA".
+            saturation (bool, optional): Return the saturated layer of this band. Defaults to False.
+
+        Returns:
+            tuple: image of the band, saturated pixels or None, and valid pixels before conversion
+        """
         # The No Data value is represented as 0 in the L1C product. see https://sentinels.copernicus.eu/web/sentinel/technical-guides/sentinel-2-msi/level-1c/product-formatting
         src_nodata = self.nodata_index
         # specify the band information
         band = band.upper()  # uppercase
-        band_file = self.bands[self.bands["NAME"] == band]["FILE"].values[0]
         res_trgt = self.profilefull.res[
             0
         ]  # resolution of processing images, such as 20 m
-        with rasterio.open(os.path.join(self.folder, band_file) + ".jp2") as src:
+        with rasterio.open(self._band_path(band)) as src:
             res_band = src.res[0]
             if res_trgt == res_band:  # band's resolution == our target
                 image = src.read(1)  # only one band, and orginal band
@@ -977,28 +1018,13 @@ def read_band(self, band, level="TOA", sza=None, maskupdate=True, saturation = F
                 image = src.read(
                     1, out_shape=self.profilefull.shape, resampling=Resampling.nearest
                 )  # only one band, and nearest is used to resample to smaller resolution
+        valid = image != src_nodata  # should be 0 when we test the Landsat Collection 2.
         satu = image == self.saturation_index if saturation else None
-        if maskupdate:
-            # self.obsmask = np.logical_and(self.obsmask, image != src.nodata) # should be 0 when we test the Landsat Collection 2.
-            if self.obsmask is None:  # initlize the mask of observation extent
-                self.obsmask = (
-                    image != src_nodata
-                )  # should be 0 when we test the Landsat Collection 2.
-            else:  # update by overlapping it with new band's extent
-                self.obsmask = np.logical_and(
-                    self.obsmask, image != src_nodata
-                )  # should be 0 when we test the Landsat Collection 2.
-
-        # append the profile of the geotiff, which will be used as the reference to auxiliary data warping and mask saving.
-        if profile:
-            self.profile = src.meta.copy()
-            self.profilefull = src  # save the full profile for the further use
 
         # return final result
         if level == "TOA":
-            return self.convert_dn_to_toa(band, image), satu
-        else:
-            return image, satu
+            return self.convert_dn_to_toa(band, image), satu, valid
+        return image, satu, valid
 
     def convert_dn_to_toa(self, band, data):
         """
@@ -1141,72 +1167,20 @@ def read_datacube(self, predictors, maskupdate=True):
             dtype=np.float32,
         )
         # Spectral bands
-        predictor = "coastal"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], _ = self.read_band(
-                band=predictor, level="TOA", maskupdate=maskupdate
-            )
-        predictor = "blue"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], satu_blue = self.read_band(
-                band=predictor, level="TOA", maskupdate=maskupdate, saturation=True
-            )
-        predictor = "green"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], satu_green = self.read_band(
-                band=predictor, level="TOA", maskupdate=maskupdate, saturation=True
-            )
-        predictor = "red"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], satu_red = self.read_band(
-                band=predictor, level="TOA", maskupdate=maskupdate, saturation=True
-            )
-        predictor = "vre1"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], _ = self.read_band(
-                band=predictor, level="TOA", maskupdate=maskupdate
-            )
-        predictor = "vre2"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], _ = self.read_band(
-                band=predictor, level="TOA", maskupdate=maskupdate
-            )
-        predictor = "vre3"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], _ = self.read_band(
-                band=predictor, level="TOA", maskupdate=maskupdate
-            )
-
-        predictor = "wnir"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], _ = self.read_band(
-                band=predictor, level="TOA", maskupdate=maskupdate
-            )
-        predictor = "nir"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], _ = self.read_band(
-                band=predictor, level="TOA", maskupdate=maskupdate
-            )
-        predictor = "wv"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], _ = self.read_band(
-                band=predictor, level="TOA", maskupdate=maskupdate
-            )
-        predictor = "swir1"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], _ = self.read_band(
-                band=predictor, level="TOA", maskupdate=maskupdate
-            )
-        predictor = "swir2"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], _ = self.read_band(
-                band=predictor, level="TOA", maskupdate=maskupdate
-            )
-        predictor = "cirrus"
-        if predictor in predictors:
-            bands[predictors.index(predictor), :, :], _ = self.read_band(
-                band=predictor, level="TOA", maskupdate=maskupdate
-            )
+        requests = [
+            {"band": band, "level": "TOA", "saturation": saturation}
+            for band, saturation in SENTINEL2_SPECTRAL_READS
+            if band in predictors
+        ]
+        # saturation for visible bands after the calculation of whiteness
+        statu = BitLayer(bands.shape[1:])
+        for request, (image, satu, valid) in zip(requests, self._read_bands(requests)):
+            if maskupdate:
+                self._update_obsmask(valid)
+            bands[predictors.index(request["band"]), :, :] = image
+            if satu is not None:
+                statu.append(satu) # set bits 0, 1 and 2 for blue, green and red, in request order
+        del image, satu, valid
 
         # spectral index
         predictor = "hot"
@@ -1218,13 +1192,6 @@ def read_datacube(self, predictors, maskupdate=True):
                 bands[predictors.index("red"), :, :],
             )
 
-        # saturation for visible bands after the calculation of whiteness
-        statu = BitLayer(bands.shape[1:])
-        statu.append(satu_blue)
-        statu.append(satu_green)
-        statu.append(satu_red)
-        del satu_blue, satu_green, satu_red
-
         predictor = "whiteness"
         if predictor in predictors:
             if C.MSG_FULL:
diff --git a/tests/test_sentinel2_read.py b/tests/test_sentinel2_read.py
new file mode 100644
index 0000000..71775ba
--- /dev/null
+++ b/tests/test_sentinel2_read.py
@@ -0,0 +1,186 @@
+import numpy as np
+import pandas
+import pytest
+import rasterio
+from rasterio.enums import Resampling
+from rasterio.transform import from_origin
+
+from fmask import constant as C
+from fmask.satellite import Sentinel2
+
+TARGET = 96
+RES = 20
+# (name, file, resolution): 20 m is equal to the target, 10 m is averaged, 60 m is nearest.
+BANDS = [
+    ("COASTAL", "B01", 60), ("BLUE", "B02", 10), ("GREEN", "B03", 10), ("RED", "B04", 10),
+    ("VRE1", "B05", 20), ("VRE2", "B06", 20), ("VRE3", "B07", 20), ("WNIR", "B08", 10),
+    ("NIR", "B8A", 20), ("WV", "B09", 60), ("CIRRUS", "B10", 60), ("SWIR1", "B11", 20),
+    ("SWIR2", "B12", 20),
+]
+SPECTRAL = [
+    "coastal", "blue", "green", "red", "vre1", "vre2", "vre3",
+    "wnir", "nir", "wv", "swir1", "swir2", "cirrus",
+]
+OFFSETS = {name: -1000.0 - index for index, (name, _, _) in enumerate(BANDS)}
+QUANTIFICATION = 10000.0
+
+
+def _write(path, data, res, driver):
+    # One resolution level means no overviews, so reads really use the requested resampling.
+    options = {"QUALITY": 100, "REVERSIBLE": "YES", "RESOLUTIONS": 1} if driver == "JP2OpenJPEG" else {}
+    with rasterio.open(
+        path, "w", driver=driver, height=data.shape[0], width=data.shape[1], count=1,
+        dtype="uint16", crs="EPSG:32633", transform=from_origin(500000, 4000000, res, res),
+        **options,
+    ) as dst:
+        dst.write(data, 1)
+
+
+@pytest.fixture
+def scene(tmp_path):
+    rng = np.random.default_rng(0)
+    for offset, (_, filename, res) in enumerate(BANDS):
+        size = TARGET * RES // res
+        data = rng.integers(1000, 12000, (size, size), dtype=np.uint16)
+        data[offset, :] = 0
+        data[:, offset] = 0
+        _write(tmp_path / f"{filename}.jp2", data, res, "JP2OpenJPEG")
+    reference = tmp_path / "reference.tif"
+    _write(reference, np.ones((TARGET, TARGET), dtype=np.uint16), RES, "GTiff")
+    with rasterio.open(reference) as src:
+        profile = src.meta.copy()
+        profilefull = src
+    return {
+        "folder": str(tmp_path),
+        "bands": pandas.DataFrame(
+            {
+                "NAME": [name for name, _, _ in BANDS],
+                "FILE": [filename for _, filename, _ in BANDS],
+                "RADIO_ADD_OFFSET": [OFFSETS[name] for name, _, _ in BANDS],
+                "QUANTIFICATION_VALUE": [QUANTIFICATION] * len(BANDS),
+            }
+        ),
+        "profile": profile,
+        "profilefull": profilefull,
+    }
+
+
+def _make(scene):
+    image = object.__new__(Sentinel2)
+    image.folder = scene["folder"]
+    image.bands = scene["bands"]
+    image.profile = scene["profile"]
+    image.profilefull = scene["profilefull"]
+    image.obsmask = None
+    image.nodata_index = 0.0
+    image.saturation_index = 65535.0
+    image.lat_north = 0.0
+    image.lat_south = 0.0
+    return image
+
+
+def _assert_same(actual, expected):
+    assert actual.dtype == expected.dtype
+    assert actual.shape == expected.shape
+    assert actual.tobytes() == expected.tobytes()
+
+
+def _raw(scene, band):
+    name, filename, res = next(entry for entry in BANDS if entry[0] == band.upper())
+    with rasterio.open(f"{scene['folder']}/{filename}.jp2") as src:
+        if res == RES:
+            return src.read(1)
+        resampling = Resampling.average if res < RES else Resampling.nearest
+        return src.read(1, out_shape=(TARGET, TARGET), resampling=resampling)
+
+
+def _expected_toa(scene, band):
+    return (_raw(scene, band) + OFFSETS[band.upper()]) / QUANTIFICATION
+
+
+def _expected_obsmask(scene, bands):
+    return np.logical_and.reduce([_raw(scene, band) != 0 for band in bands])
+
+
+@pytest.mark.parametrize("band", ["vre1", "blue", "coastal"], ids=["equal", "average", "nearest"])
+def test_read_band_converts_to_toa(scene, band):
+    result, satu = _make(scene).read_band(band)
+    _assert_same(result, _expected_toa(scene, band))
+    assert satu is None
+
+
+@pytest.mark.parametrize("band", ["vre1", "blue", "coastal"], ids=["equal", "average", "nearest"])
+def test_read_band_returns_dn_for_other_levels(scene, band):
+    result, _ = _make(scene).read_band(band, level="DN")
+    _assert_same(result, _raw(scene, band))
+    assert result.shape == (TARGET, TARGET)
+
+
+def test_read_band_overlaps_the_observation_mask(scene):
+    image = _make(scene)
+    image.read_band("vre1")
+    _assert_same(image.obsmask, _raw(scene, "vre1") != 0)
+    image.read_band("blue")
+    _assert_same(image.obsmask, _expected_obsmask(scene, ["vre1", "blue"]))
+
+
+def test_read_band_can_leave_the_observation_mask(scene):
+    image = _make(scene)
+    image.read_band("blue", maskupdate=False)
+    assert image.obsmask is None
+
+
+def test_read_band_profile(scene):
+    image = _make(scene)
+    image.read_band("vre1", profile=True)
+    with rasterio.open(f"{scene['folder']}/B05.jp2") as src:
+        assert image.profile == src.meta
+    assert image.profilefull.name.endswith("B05.jp2")
+    assert image.profilefull.closed
+
+
+def test_values_do_not_touch_self(scene):
+    image = _make(scene)
+    mask = np.ones((TARGET, TARGET), dtype=bool)
+    profile, profilefull = image.profile, image.profilefull
+    image.obsmask = mask
+    _, _, valid = image._read_band_values("vre1", "TOA")
+    assert image.obsmask is mask
+    assert image.profile is profile
+    assert image.profilefull is profilefull
+    assert valid.dtype == np.bool_
+
+
+def test_obsmask_all_bands(scene):
+    image = _make(scene)
+    image.read_datacube(list(SPECTRAL))
+    _assert_same(image.obsmask, _expected_obsmask(scene, SPECTRAL))
+    assert not image.obsmask.all()
+    assert image.obsmask.any()
+
+
+@pytest.mark.parametrize("maskupdate", [True, False])
+@pytest.mark.parametrize("predictors", [SPECTRAL, list(reversed(SPECTRAL))], ids=["forward", "reversed"])
+def test_read_datacube_matches_each_band(scene, predictors, maskupdate):
+    image = _make(scene)
+    data, _ = image.read_datacube(list(predictors), maskupdate=maskupdate)
+    expected = np.stack([_expected_toa(scene, band) for band in predictors]).astype(np.float32)
+    if maskupdate:
+        expected[:, ~_expected_obsmask(scene, SPECTRAL)] = 0
+    else:
+        assert image.obsmask is None
+    _assert_same(data.data, expected)
+    assert data.bands == list(predictors)
+
+
+def test_read_datacube_replaces_nan_values(scene, monkeypatch):
+    original = Sentinel2._read_band_values
+
+    def with_nan(self, band, **kwargs):
+        image, satu, valid = original(self, band, **kwargs)
+        image[20, 20] = np.nan
+        return image, satu, valid
+
+    monkeypatch.setattr(Sentinel2, "_read_band_values", with_nan)
+    data, _ = _make(scene).read_datacube(list(SPECTRAL))
+    assert (data.data[:, 20, 20] == np.float32(C.EPS)).all()

From 38449eb3970c34836c375fc5c7ca843e86c240b8 Mon Sep 17 00:00:00 2001
From: Chris Holden <ceholden@gmail.com>
Date: Wed, 4 Nov 2026 17:26:58 -0500
Subject: [PATCH 07/11] test: use reflectance-scale inputs in the UNet classify
 test

With raw values up to 3000 and normalize=False, logits reach about 1800, where one
float32 step can move a probability by more than the 1e-5 tolerance. Batches of
different shapes can round differently on x86, so the test failed in CI on Linux
while passing on macOS.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 tests/test_unet_classify.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_unet_classify.py b/tests/test_unet_classify.py
index 17ecc37..2a7dea7 100644
--- a/tests/test_unet_classify.py
+++ b/tests/test_unet_classify.py
@@ -16,7 +16,7 @@
 
 def _stub(batch_size=3):
     torch.manual_seed(0)
-    cube = (np.random.default_rng(0).random((len(BANDS), *SHAPE)) * 3000).astype(np.float32)
+    cube = (np.random.default_rng(0).random((len(BANDS), *SHAPE)) * 3).astype(np.float32)
     obsmask = np.ones(SHAPE, dtype=bool)
     obsmask[:6] = False
     obsmask[:, -5:] = False

From 1c09f948d9213f113c94bf18958ccf7b6c982048 Mon Sep 17 00:00:00 2001
From: Chris Holden <ceholden@gmail.com>
Date: Wed, 4 Nov 2026 17:29:06 -0500
Subject: [PATCH 08/11] chore: remove unused UNet.classify_backup

Nothing in the repository calls it, and it had drifted from classify.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/fmask/unetlib.py | 162 -------------------------------------------
 1 file changed, 162 deletions(-)

diff --git a/src/fmask/unetlib.py b/src/fmask/unetlib.py
index 2738e60..c891ae1 100644
--- a/src/fmask/unetlib.py
+++ b/src/fmask/unetlib.py
@@ -712,168 +712,6 @@ def classify(self, probability="none", normalize=True, shift=True) -> tuple:
             image_prob[missing_mask] = image_prob_full[missing_mask]
         return image_class, image_prob
     
-    def classify_backup(self, probability="none", normalize=True, shift=True) -> tuple:
-        """Classify the image by the model and return the class and probability
-
-        Args:
-            probability (srt, optional): One item of the classes, or "none" (means not to extract the prob layer), or "default" (highest score for the classified results). Defaults to "none".
-            normalize (bool, optional): Whether to normalize the datacube. Defaults to True.
-            shift (bool, optional): Whether to shift the patch off anchors. Defaults to True.
-        Returns:
-            image_class, image_prob: class result (0: noncloud, 1: cloud) and probability result
-        """
-        # check the probability index
-        prob_index = (
-            -1 if probability == "default" else
-            -2 if probability == "none" else
-            self.classes.index(probability)
-        )
-
-        # rescale the date cube
-        datacube = self.image.data.get(self.predictors)
-        if normalize:
-            datacube = normalize_datacube(datacube, obsmask=self.image.obsmask)
-    
-        if C.MSG_FULL:
-            print(">>> classifying image by unet")
-
-        # Convert the datacube to tensor if it is not
-        if not isinstance(datacube, torch.Tensor):
-            datacube = torch.from_numpy(datacube).to(
-                self.device, dtype=torch.float
-            )  # (batch, chanel, height, width)
-        datacube = torch.unsqueeze(
-            datacube, dim=0
-        )  # add one more dimension for the batch
-        offanchors = init_patch_offanchors(self.image.obsmask, size=self.patch_size, stride=self.patch_stride_classify, shift=shift)
-        # Change the r c off
-        wincut = int((self.patch_size - self.patch_stride_classify) / 2)
-   
-        # initlized image_class and image_prob
-
-        h, w = datacube.shape[2], datacube.shape[3]
-        image_class = np.full((h, w), 255, dtype=np.uint8) # 255 means no data, will be filled later using image_class_full
-        image_class_full = np.zeros((h, w), dtype=np.uint8)
-        if prob_index == -2:
-            image_prob = None  # no need to initialize
-        else:
-            image_prob = np.zeros((h, w), dtype=np.float32)
-            image_prob_full = np.zeros((h, w), dtype=np.float32)
-        self.model.eval()
-        with torch.no_grad():  # disabled gradient calculation
-            for offanchor in offanchors:
-                r_off, c_off = offanchor[0], offanchor[1]
-                # classify the pacth
-                patchprob = torch.nn.functional.softmax(
-                    self.model(
-                        datacube[
-                            :,
-                            :,
-                            r_off : r_off + self.patch_size,
-                            c_off : c_off + self.patch_size,
-                        ]
-                    ),
-                    dim=1,
-                )
-
-                # extract the prob and class
-                if (
-                    prob_index >= 0
-                ):  # prob_class_index is the index of the class that we want to extract the probability
-                    _, patchclass = patchprob.topk(
-                        1, dim=1
-                    )  # top label with top prob.
-                    # select the prob of the defined class and index based on the classes, but keep same shape
-                    patchprob = patchprob[:, prob_index, :, :].to("cpu").numpy()
-                    # patchprob = patchprob.to("cpu").numpy()
-                    image_class[
-                        r_off + wincut: r_off + self.patch_size - wincut, c_off + wincut: c_off + self.patch_size - wincut
-                    ] = patchclass[
-                        0,
-                        0,
-                        wincut: self.patch_size - wincut,
-                        wincut: self.patch_size - wincut,
-                    ]
-                    image_class_full[
-                        r_off : r_off + self.patch_size, c_off : c_off + self.patch_size
-                    ] = patchclass[
-                        0,
-                        0,
-                        :,
-                        :
-                    ]
-                    image_prob[
-                        r_off + wincut: r_off + self.patch_size - wincut, c_off + wincut: c_off + self.patch_size - wincut
-                    ] = patchprob[
-                        0,
-                        0,
-                        wincut: self.patch_size - wincut,
-                        wincut: self.patch_size - wincut,
-                    ]
-                    image_prob_full[
-                        r_off: r_off + self.patch_size, c_off: c_off + self.patch_size
-                    ] = patchprob[
-                        0,
-                        0,
-                        :,
-                        :
-                    ]
-                elif (
-                    prob_index == -1
-                ):  # -1 means only the predicted class will be extracted by this model
-                    patchprob, patchclass = patchprob.topk(
-                        1, dim=1
-                    ).to("cpu").numpy()  # top label with top prob.
-                    # patchprob = patchprob.to("cpu").numpy()
-                    image_prob[
-                        r_off + wincut: r_off + self.patch_size - wincut, c_off + wincut: c_off + self.patch_size - wincut
-                    ] = patchprob[
-                        0,
-                        0,
-                        wincut: self.patch_size - wincut,
-                        wincut: self.patch_size - wincut,
-                    ]
-                    image_prob_full[
-                        r_off: r_off + self.patch_size, c_off: c_off + self.patch_size
-                    ] = patchprob[
-                        0,
-                        0,
-                        :,
-                        :
-                    ]
-                elif (
-                    prob_index == -2
-                ):  # -2 means the top class will be extracted by this model, but not the probability
-                    _, patchclass = patchprob.topk(
-                        1, dim=1
-                    )  # top label with top prob.
-
-                    # which will be extracted all the time
-                    patchclass = patchclass.to("cpu").numpy()
-                    image_class[
-                        r_off + wincut: r_off + self.patch_size - wincut, c_off + wincut: c_off + self.patch_size - wincut
-                    ] = patchclass[
-                        0,
-                        0,
-                        wincut: self.patch_size - wincut,
-                        wincut: self.patch_size - wincut,
-                    ]
-                    image_class_full[
-                        r_off : r_off + self.patch_size, c_off : c_off + self.patch_size
-                    ] = patchclass[
-                        0,
-                        0,
-                        :,
-                        :
-                    ]
-        # fill the 255
-        # Fill missing values from full map
-        missing_mask = (image_class == 255)
-        image_class[missing_mask] = image_class_full[missing_mask]
-        if image_prob is not None:
-            image_prob[missing_mask] = image_prob_full[missing_mask]
-        return image_class, image_prob
-
     def check_device(self) -> None:
         """examine device, CPU or GPU
 

From 1028d52810f1c0d1692935dff785dc9d835a3a95 Mon Sep 17 00:00:00 2001
From: Chris Holden <ceholden@gmail.com>
Date: Thu, 5 Nov 2026 09:34:35 -0500
Subject: [PATCH 09/11] test: organize test files by source module

Each test file now follows the module it covers, so coverage is easier to reason
about: test_satellite.py for Data and the shared band reader, with
test_satellite_landsat.py and test_satellite_sentinel.py for each sensor, then
test_utils.py, test_reconstruct.py, test_lightgbmlib.py and test_unetlib.py.

The shared band reader tests now run against the Satellite base class instead of
a Landsat scene. The same 178 tests are collected as before.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 ...est_data_sample.py => test_lightgbmlib.py} |  27 ---
 tests/test_reconstruct.py                     |  45 +----
 tests/test_satellite.py                       | 154 ++++++++++++++++++
 ...dsat_read.py => test_satellite_landsat.py} | 104 ------------
 ...el2_read.py => test_satellite_sentinel.py} |  99 ++++++++++-
 tests/test_sentinel2_saturation.py            | 143 ----------------
 ...{test_unet_classify.py => test_unetlib.py} |   0
 tests/{test_normalize.py => test_utils.py}    |  48 ++++++
 8 files changed, 294 insertions(+), 326 deletions(-)
 rename tests/{test_data_sample.py => test_lightgbmlib.py} (70%)
 create mode 100644 tests/test_satellite.py
 rename tests/{test_landsat_read.py => test_satellite_landsat.py} (73%)
 rename tests/{test_sentinel2_read.py => test_satellite_sentinel.py} (66%)
 delete mode 100644 tests/test_sentinel2_saturation.py
 rename tests/{test_unet_classify.py => test_unetlib.py} (100%)
 rename tests/{test_normalize.py => test_utils.py} (77%)

diff --git a/tests/test_data_sample.py b/tests/test_lightgbmlib.py
similarity index 70%
rename from tests/test_data_sample.py
rename to tests/test_lightgbmlib.py
index 81166ac..1636e24 100644
--- a/tests/test_data_sample.py
+++ b/tests/test_lightgbmlib.py
@@ -14,33 +14,6 @@ def _data():
     return Data(cube, BANDS)
 
 
-def _coords():
-    rng = np.random.default_rng(1)
-    rows, cols = rng.integers(0, 30, 200), rng.integers(0, 40, 200)
-    rows[:20], cols[:20] = rows[20:40], cols[20:40]
-    return rows, cols
-
-
-def test_sample_matches_get():
-    data = _data()
-    bands = ["nir", "blue", "red"]
-    rows, cols = _coords()
-    expected = data.get(bands)[:, rows, cols]
-    actual = data.sample(bands, rows, cols)
-    assert actual.dtype == expected.dtype
-    assert actual.shape == (3, 200)
-    assert np.array_equal(actual, expected)
-    assert actual.T.tobytes() == expected.T.tobytes()
-
-
-def test_sample_is_copy():
-    data = _data()
-    before = data.data.copy()
-    rows, cols = _coords()
-    data.sample(["nir", "blue"], rows, cols)[:] = -1
-    assert np.array_equal(data.data, before)
-
-
 class _RecordingModel:
     classes_ = np.array([0, 1])
 
diff --git a/tests/test_reconstruct.py b/tests/test_reconstruct.py
index fa3d62d..2762780 100644
--- a/tests/test_reconstruct.py
+++ b/tests/test_reconstruct.py
@@ -2,7 +2,7 @@
 import pytest
 from skimage.morphology import reconstruction
 
-from fmask import _reconstruct, utils
+from fmask import _reconstruct
 from fmask._reconstruct import reconstruct_erosion
 
 
@@ -169,46 +169,3 @@ def spy(seed, mask, method):
     seed = _hole_seed(mask)
     _assert_same(reconstruct_erosion(seed, mask), reconstruction(seed, mask, method="erosion"))
     assert calls == ["erosion"]
-
-
-@pytest.mark.parametrize("fill_value", [None, -1.0])
-@pytest.mark.parametrize("obs_shape", ["full", "partial", "single"])
-def test_imfill_matches_skimage(obs_shape, fill_value):
-    data = _random(np.float32, 11, (31, 37)) * 10
-    obsmask = np.ones(data.shape, dtype=bool)
-    if obs_shape == "partial":
-        obsmask[:5] = False
-        obsmask[:, -4:] = False
-    elif obs_shape == "single":
-        obsmask[:] = False
-        obsmask[15, 15] = True
-    max_val = data[obsmask].max()
-    filled_input = data.copy()
-    filled_input[~obsmask] = max_val if fill_value is None else fill_value
-    seed = filled_input.copy()
-    seed[1:-1, 1:-1] = max_val
-
-    result = utils.imfill(data, obsmask, fill_value=fill_value)
-
-    _assert_same(result, reconstruction(seed, filled_input, method="erosion"))
-    _assert_same(data, filled_input)
-
-
-def test_imfill_fills_a_closed_hole():
-    data = np.full((5, 5), 5.0, dtype=np.float32)
-    data[2, 2] = 1.0
-    result = utils.imfill(data, np.ones((5, 5), dtype=bool))
-    assert result[2, 2] == 5.0
-
-
-def test_imfill_keeps_a_depression_open_to_the_border():
-    data = np.full((5, 5), 5.0, dtype=np.float32)
-    data[2, :3] = 1.0
-    result = utils.imfill(data, np.ones((5, 5), dtype=bool))
-    assert (result[2, :3] == 1.0).all()
-
-
-def test_imfill_empty_observation_mask_raises():
-    data = _random(np.float32, 12, (8, 8))
-    with pytest.raises(ValueError):
-        utils.imfill(data, np.zeros(data.shape, dtype=bool))
diff --git a/tests/test_satellite.py b/tests/test_satellite.py
new file mode 100644
index 0000000..160cd24
--- /dev/null
+++ b/tests/test_satellite.py
@@ -0,0 +1,154 @@
+import gc
+import inspect
+import threading
+import time
+import weakref
+
+import numpy as np
+import pytest
+
+from fmask import satellite
+from fmask.satellite import Data, Landsat, Satellite, Sentinel2
+
+BANDS = ["blue", "green", "red", "nir"]
+
+
+class _Reader(Satellite):
+    def _read_band_values(self, band, **kwargs):
+        raise NotImplementedError
+
+
+def _data():
+    cube = np.random.default_rng(0).random((4, 30, 40)).astype(np.float32)
+    return Data(cube, BANDS)
+
+
+def _coords():
+    rng = np.random.default_rng(1)
+    rows, cols = rng.integers(0, 30, 200), rng.integers(0, 40, 200)
+    rows[:20], cols[:20] = rows[20:40], cols[20:40]
+    return rows, cols
+
+
+def test_sample_matches_get():
+    data = _data()
+    bands = ["nir", "blue", "red"]
+    rows, cols = _coords()
+    expected = data.get(bands)[:, rows, cols]
+    actual = data.sample(bands, rows, cols)
+    assert actual.dtype == expected.dtype
+    assert actual.shape == (3, 200)
+    assert np.array_equal(actual, expected)
+    assert actual.T.tobytes() == expected.T.tobytes()
+
+
+def test_sample_is_copy():
+    data = _data()
+    before = data.data.copy()
+    rows, cols = _coords()
+    data.sample(["nir", "blue"], rows, cols)[:] = -1
+    assert np.array_equal(data.data, before)
+
+
+def test_read_bands_order(monkeypatch):
+    image = object.__new__(_Reader)
+    release = threading.Event()
+
+    def fake_values(self, band, **kwargs):
+        if band == "first":
+            assert release.wait(timeout=10)
+        else:
+            release.set()
+        return band
+
+    monkeypatch.setattr(_Reader, "_read_band_values", fake_values)
+    requests = [{"band": "first", "level": "TOA"}, {"band": "second", "level": "TOA"}]
+    assert list(image._read_bands(requests)) == ["first", "second"]
+
+
+def test_read_bands_limits_threads(monkeypatch):
+    image = object.__new__(_Reader)
+    lock = threading.Lock()
+    active = [0, 0]
+    barrier = threading.Barrier(satellite.READ_THREADS)
+
+    def fake_values(self, band, **kwargs):
+        with lock:
+            active[0] += 1
+            active[1] = max(active)
+        barrier.wait(timeout=10)
+        with lock:
+            active[0] -= 1
+        return band
+
+    monkeypatch.setattr(_Reader, "_read_band_values", fake_values)
+    count = satellite.READ_THREADS * 3
+    list(image._read_bands([{"band": str(i), "level": "TOA"} for i in range(count)]))
+    assert active[1] == satellite.READ_THREADS
+
+
+def test_read_bands_propagates_errors(monkeypatch):
+    image = object.__new__(_Reader)
+
+    def fake_values(self, band, **kwargs):
+        raise OSError("bad band")
+
+    monkeypatch.setattr(_Reader, "_read_band_values", fake_values)
+    with pytest.raises(OSError, match="bad band"):
+        list(image._read_bands([{"band": "blue", "level": "TOA"}]))
+
+
+def test_worker_threads_use_caller_error_state(monkeypatch):
+    image = object.__new__(_Reader)
+    monkeypatch.setattr(_Reader, "_read_band_values", lambda self, band, **kw: np.geterr()["invalid"])
+    with np.errstate(invalid="ignore"):
+        assert list(image._read_bands([{"band": "blue", "level": "TOA"}])) == ["ignore"]
+    with np.errstate(invalid="raise"):
+        assert list(image._read_bands([{"band": "blue", "level": "TOA"}])) == ["raise"]
+
+
+def test_read_bands_releases_each_result_before_the_next(monkeypatch):
+    image = object.__new__(_Reader)
+
+    class Value:
+        pass
+
+    refs = {}
+
+    def fake_values(self, band, **kwargs):
+        value = Value()
+        refs[band] = weakref.ref(value)
+        return value
+
+    monkeypatch.setattr(_Reader, "_read_band_values", fake_values)
+    results = image._read_bands([{"band": str(i), "level": "TOA"} for i in range(3)])
+    next(results)
+    gc.collect()
+    value = next(results)
+    gc.collect()
+    assert refs["0"]() is None
+    del value
+    list(results)
+
+
+def test_read_bands_limits_results_waiting_behind_a_slow_band(monkeypatch):
+    image = object.__new__(_Reader)
+    started = []
+    started_during_slow_band = []
+
+    def fake_values(self, band, **kwargs):
+        started.append(band)
+        if band == "0":
+            time.sleep(0.3)
+            started_during_slow_band.append(len(started))
+        return band
+
+    monkeypatch.setattr(_Reader, "_read_band_values", fake_values)
+    count = satellite.READ_THREADS * 3
+    requests = [{"band": str(i), "level": "TOA"} for i in range(count)]
+    assert list(image._read_bands(requests)) == [str(i) for i in range(count)]
+    assert started_during_slow_band == [satellite.READ_THREADS]
+
+
+def test_read_band_arguments_match_landsat():
+    assert inspect.signature(Sentinel2.read_band) == inspect.signature(Landsat.read_band)
diff --git a/tests/test_landsat_read.py b/tests/test_satellite_landsat.py
similarity index 73%
rename from tests/test_landsat_read.py
rename to tests/test_satellite_landsat.py
index ea6392b..9cde80b 100644
--- a/tests/test_landsat_read.py
+++ b/tests/test_satellite_landsat.py
@@ -1,8 +1,4 @@
-import gc
 import sys
-import threading
-import time
-import weakref
 
 import numpy as np
 import pandas
@@ -98,106 +94,6 @@ def test_values_do_not_touch_self(scene):
     assert not valid[1, :].any()
 
 
-def test_read_bands_order(scene, monkeypatch):
-    image = _make(Landsat, scene)
-    release = threading.Event()
-
-    def fake_values(self, band, **kwargs):
-        if band == "first":
-            assert release.wait(timeout=10)
-        else:
-            release.set()
-        return band
-
-    monkeypatch.setattr(Landsat, "_read_band_values", fake_values)
-    requests = [{"band": "first", "level": "TOA"}, {"band": "second", "level": "TOA"}]
-    assert list(image._read_bands(requests)) == ["first", "second"]
-
-
-def test_read_bands_limits_threads(scene, monkeypatch):
-    image = _make(Landsat, scene)
-    lock = threading.Lock()
-    active = [0, 0]
-    barrier = threading.Barrier(satellite.READ_THREADS)
-
-    def fake_values(self, band, **kwargs):
-        with lock:
-            active[0] += 1
-            active[1] = max(active)
-        barrier.wait(timeout=10)
-        with lock:
-            active[0] -= 1
-        return band
-
-    monkeypatch.setattr(Landsat, "_read_band_values", fake_values)
-    count = satellite.READ_THREADS * 3
-    list(image._read_bands([{"band": str(i), "level": "TOA"} for i in range(count)]))
-    assert active[1] == satellite.READ_THREADS
-
-
-def test_read_bands_propagates_errors(scene, monkeypatch):
-    image = _make(Landsat, scene)
-
-    def fake_values(self, band, **kwargs):
-        raise OSError("bad band")
-
-    monkeypatch.setattr(Landsat, "_read_band_values", fake_values)
-    with pytest.raises(OSError, match="bad band"):
-        list(image._read_bands([{"band": "blue", "level": "TOA"}]))
-
-
-def test_worker_threads_use_caller_error_state(scene, monkeypatch):
-    image = _make(Landsat, scene)
-    monkeypatch.setattr(Landsat, "_read_band_values", lambda self, band, **kw: np.geterr()["invalid"])
-    with np.errstate(invalid="ignore"):
-        assert list(image._read_bands([{"band": "blue", "level": "TOA"}])) == ["ignore"]
-    with np.errstate(invalid="raise"):
-        assert list(image._read_bands([{"band": "blue", "level": "TOA"}])) == ["raise"]
-
-
-def test_read_bands_releases_each_result_before_the_next(scene, monkeypatch):
-    image = _make(Landsat, scene)
-
-    class Value:
-        pass
-
-    refs = {}
-
-    def fake_values(self, band, **kwargs):
-        value = Value()
-        refs[band] = weakref.ref(value)
-        return value
-
-    monkeypatch.setattr(Landsat, "_read_band_values", fake_values)
-    results = image._read_bands([{"band": str(i), "level": "TOA"} for i in range(3)])
-    next(results)
-    gc.collect()
-    value = next(results)
-    gc.collect()
-    assert refs["0"]() is None
-    del value
-    list(results)
-
-
-def test_read_bands_limits_results_waiting_behind_a_slow_band(scene, monkeypatch):
-    image = _make(Landsat, scene)
-    started = []
-    started_during_slow_band = []
-
-    def fake_values(self, band, **kwargs):
-        started.append(band)
-        if band == "0":
-            time.sleep(0.3)
-            started_during_slow_band.append(len(started))
-        return band
-
-    monkeypatch.setattr(Landsat, "_read_band_values", fake_values)
-    count = satellite.READ_THREADS * 3
-    requests = [{"band": str(i), "level": "TOA"} for i in range(count)]
-    assert list(image._read_bands(requests)) == [str(i) for i in range(count)]
-    assert started_during_slow_band == [satellite.READ_THREADS]
-
-
 def test_read_datacube_releases_read_inputs_before_indices(scene, monkeypatch):
     original_hot = satellite.hot
     seen = {}
diff --git a/tests/test_sentinel2_read.py b/tests/test_satellite_sentinel.py
similarity index 66%
rename from tests/test_sentinel2_read.py
rename to tests/test_satellite_sentinel.py
index 71775ba..3b9e131 100644
--- a/tests/test_sentinel2_read.py
+++ b/tests/test_satellite_sentinel.py
@@ -22,6 +22,7 @@
     "wnir", "nir", "wv", "swir1", "swir2", "cirrus",
 ]
 OFFSETS = {name: -1000.0 - index for index, (name, _, _) in enumerate(BANDS)}
+PIXELS = {"blue": (50, 50), "green": (40, 40), "red": (30, 30)}
 QUANTIFICATION = 10000.0
 
 
@@ -36,27 +37,27 @@ def _write(path, data, res, driver):
         dst.write(data, 1)
 
 
-@pytest.fixture
-def scene(tmp_path):
+def _write_scene(folder, offsets, nodata_lines):
     rng = np.random.default_rng(0)
     for offset, (_, filename, res) in enumerate(BANDS):
         size = TARGET * RES // res
         data = rng.integers(1000, 12000, (size, size), dtype=np.uint16)
-        data[offset, :] = 0
-        data[:, offset] = 0
-        _write(tmp_path / f"{filename}.jp2", data, res, "JP2OpenJPEG")
-    reference = tmp_path / "reference.tif"
+        if nodata_lines:
+            data[offset, :] = 0
+            data[:, offset] = 0
+        _write(folder / f"{filename}.jp2", data, res, "JP2OpenJPEG")
+    reference = folder / "reference.tif"
     _write(reference, np.ones((TARGET, TARGET), dtype=np.uint16), RES, "GTiff")
     with rasterio.open(reference) as src:
         profile = src.meta.copy()
         profilefull = src
     return {
-        "folder": str(tmp_path),
+        "folder": str(folder),
         "bands": pandas.DataFrame(
             {
                 "NAME": [name for name, _, _ in BANDS],
                 "FILE": [filename for _, filename, _ in BANDS],
-                "RADIO_ADD_OFFSET": [OFFSETS[name] for name, _, _ in BANDS],
+                "RADIO_ADD_OFFSET": [offsets[name] for name, _, _ in BANDS],
                 "QUANTIFICATION_VALUE": [QUANTIFICATION] * len(BANDS),
             }
         ),
@@ -79,6 +80,17 @@ def _make(scene):
     return image
 
 
+@pytest.fixture
+def scene(tmp_path):
+    return _write_scene(tmp_path, OFFSETS, nodata_lines=True)
+
+
+@pytest.fixture(params=[-1000.0, -1001.0, 0.0], ids=["offset_1000", "offset_1001", "no_offset"])
+def clear_image(tmp_path, request):
+    offsets = {name: request.param for name, _, _ in BANDS}
+    return _make(_write_scene(tmp_path, offsets, nodata_lines=False))
+
+
 def _assert_same(actual, expected):
     assert actual.dtype == expected.dtype
     assert actual.shape == expected.shape
@@ -184,3 +196,74 @@ def with_nan(self, band, **kwargs):
     monkeypatch.setattr(Sentinel2, "_read_band_values", with_nan)
     data, _ = _make(scene).read_datacube(list(SPECTRAL))
     assert (data.data[:, 20, 20] == np.float32(C.EPS)).all()
+
+
+def _set_dn(image, band, row, col, value):
+    """Set one target pixel in the band file, covering every native pixel inside it."""
+    _, filename, res = next(entry for entry in BANDS if entry[0] == band.upper())
+    path = f"{image.folder}/{filename}.jp2"
+    with rasterio.open(path) as src:
+        data = src.read(1)
+    scale = RES // res
+    data[row * scale:(row + 1) * scale, col * scale:(col + 1) * scale] = value
+    _write(path, data, res, "JP2OpenJPEG")
+
+
+@pytest.mark.parametrize("level", ["TOA", "DN"])
+def test_read_band_returns_saturated_pixels(clear_image, level):
+    _set_dn(clear_image, "blue", 50, 50, 65535)
+    _set_dn(clear_image, "blue", 40, 40, 65534)
+
+    _, satu = clear_image.read_band("blue", level=level, saturation=True)
+
+    expected = np.zeros((TARGET, TARGET), dtype=bool)
+    expected[50, 50] = True
+    assert np.array_equal(satu, expected)
+
+
+def test_read_band_leaves_saturation_out_by_default(clear_image):
+    _set_dn(clear_image, "blue", 50, 50, 65535)
+
+    _, satu = clear_image.read_band("blue")
+
+    assert satu is None
+
+
+def test_read_datacube_sets_saturation_bits_in_blue_green_red_order(clear_image):
+    for band, (row, col) in PIXELS.items():
+        _set_dn(clear_image, band, row, col, 65535)
+
+    _, statu = clear_image.read_datacube(list(SPECTRAL))
+
+    assert statu.index == 2
+    for bit, band in enumerate(PIXELS):
+        expected = np.zeros((TARGET, TARGET), dtype=bool)
+        expected[PIXELS[band]] = True
+        assert np.array_equal(statu.get(bit), expected)
+
+
+def test_read_datacube_does_not_flag_the_dn_below_saturation(clear_image):
+    for band in PIXELS:
+        _set_dn(clear_image, band, 50, 50, 65534)
+
+    _, statu = clear_image.read_datacube(list(SPECTRAL))
+
+    assert not (statu.get(0) | statu.get(1) | statu.get(2)).any()
+
+
+@pytest.mark.parametrize("missing", list(PIXELS))
+def test_read_datacube_without_visible_band_raises(clear_image, missing):
+    predictors = [band for band in SPECTRAL if band != missing]
+    with pytest.raises(ValueError, match=missing):
+        clear_image.read_datacube(predictors)
+
+
+@pytest.mark.parametrize("band", list(PIXELS))
+def test_read_datacube_zeroes_whiteness_of_saturated_pixels(clear_image, band):
+    _set_dn(clear_image, band, 50, 50, 65535)
+
+    data, _ = clear_image.read_datacube(SPECTRAL + ["whiteness"])
+
+    expected_nonzero = np.ones((TARGET, TARGET), dtype=bool)
+    expected_nonzero[50, 50] = False
+    assert np.array_equal(data.data[data.bands.index("whiteness")] != 0, expected_nonzero)
diff --git a/tests/test_sentinel2_saturation.py b/tests/test_sentinel2_saturation.py
deleted file mode 100644
index 9648eae..0000000
--- a/tests/test_sentinel2_saturation.py
+++ /dev/null
@@ -1,143 +0,0 @@
-import inspect
-
-import numpy as np
-import pandas
-import pytest
-import rasterio
-from rasterio.transform import from_origin
-
-from fmask.satellite import Landsat, Sentinel2
-
-TARGET = 64
-RES = 20
-# (name, file, resolution): 20 m is equal to the target, 10 m is averaged, 60 m is nearest.
-BANDS = [
-    ("COASTAL", "B01", 60), ("BLUE", "B02", 10), ("GREEN", "B03", 10), ("RED", "B04", 10),
-    ("VRE1", "B05", 20), ("VRE2", "B06", 20), ("VRE3", "B07", 20), ("WNIR", "B08", 10),
-    ("NIR", "B8A", 20), ("WV", "B09", 60), ("CIRRUS", "B10", 60), ("SWIR1", "B11", 20),
-    ("SWIR2", "B12", 20),
-]
-SPECTRAL = [
-    "coastal", "blue", "green", "red", "vre1", "vre2", "vre3",
-    "wnir", "nir", "wv", "swir1", "swir2", "cirrus",
-]
-PIXELS = {"blue": (50, 50), "green": (40, 40), "red": (30, 30)}
-QUANTIFICATION = 10000.0
-
-
-def _write(path, data, res, driver):
-    # One resolution level means no overviews, so reads really use the requested resampling.
-    options = {"QUALITY": 100, "REVERSIBLE": "YES", "RESOLUTIONS": 1} if driver == "JP2OpenJPEG" else {}
-    with rasterio.open(
-        path, "w", driver=driver, height=data.shape[0], width=data.shape[1], count=1,
-        dtype="uint16", crs="EPSG:32633", transform=from_origin(500000, 4000000, res, res),
-        **options,
-    ) as dst:
-        dst.write(data, 1)
-
-
-@pytest.fixture(params=[-1000.0, -1001.0, 0.0], ids=["offset_1000", "offset_1001", "no_offset"])
-def scene(tmp_path, request):
-    rng = np.random.default_rng(0)
-    for _, filename, res in BANDS:
-        size = TARGET * RES // res
-        data = rng.integers(1000, 12000, (size, size), dtype=np.uint16)
-        _write(tmp_path / f"{filename}.jp2", data, res, "JP2OpenJPEG")
-    reference = tmp_path / "reference.tif"
-    _write(reference, np.ones((TARGET, TARGET), dtype=np.uint16), RES, "GTiff")
-    with rasterio.open(reference) as src:
-        profile = src.meta.copy()
-        profilefull = src
-    image = object.__new__(Sentinel2)
-    image.folder = str(tmp_path)
-    image.bands = pandas.DataFrame(
-        {
-            "NAME": [name for name, _, _ in BANDS],
-            "FILE": [filename for _, filename, _ in BANDS],
-            "RADIO_ADD_OFFSET": [request.param] * len(BANDS),
-            "QUANTIFICATION_VALUE": [QUANTIFICATION] * len(BANDS),
-        }
-    )
-    image.profile = profile
-    image.profilefull = profilefull
-    image.obsmask = None
-    image.nodata_index = 0.0
-    image.saturation_index = 65535.0
-    image.lat_north = 0.0
-    image.lat_south = 0.0
-    return image
-
-
-def _set_dn(image, band, row, col, value):
-    """Set one target pixel in the band file, covering every native pixel inside it."""
-    _, filename, res = next(entry for entry in BANDS if entry[0] == band.upper())
-    path = f"{image.folder}/{filename}.jp2"
-    with rasterio.open(path) as src:
-        data = src.read(1)
-    scale = RES // res
-    data[row * scale:(row + 1) * scale, col * scale:(col + 1) * scale] = value
-    _write(path, data, res, "JP2OpenJPEG")
-
-
-def test_read_band_arguments_match_landsat():
-    assert inspect.signature(Sentinel2.read_band) == inspect.signature(Landsat.read_band)
-
-
-@pytest.mark.parametrize("level", ["TOA", "DN"])
-def test_read_band_returns_saturated_pixels(scene, level):
-    _set_dn(scene, "blue", 50, 50, 65535)
-    _set_dn(scene, "blue", 40, 40, 65534)
-
-    _, satu = scene.read_band("blue", level=level, saturation=True)
-
-    expected = np.zeros((TARGET, TARGET), dtype=bool)
-    expected[50, 50] = True
-    assert np.array_equal(satu, expected)
-
-
-def test_read_band_leaves_saturation_out_by_default(scene):
-    _set_dn(scene, "blue", 50, 50, 65535)
-
-    _, satu = scene.read_band("blue")
-
-    assert satu is None
-
-
-def test_read_datacube_sets_saturation_bits_in_blue_green_red_order(scene):
-    for band, (row, col) in PIXELS.items():
-        _set_dn(scene, band, row, col, 65535)
-
-    _, statu = scene.read_datacube(list(SPECTRAL))
-
-    assert statu.index == 2
-    for bit, band in enumerate(PIXELS):
-        expected = np.zeros((TARGET, TARGET), dtype=bool)
-        expected[PIXELS[band]] = True
-        assert np.array_equal(statu.get(bit), expected)
-
-
-def test_read_datacube_does_not_flag_the_dn_below_saturation(scene):
-    for band in PIXELS:
-        _set_dn(scene, band, 50, 50, 65534)
-
-    _, statu = scene.read_datacube(list(SPECTRAL))
-
-    assert not (statu.get(0) | statu.get(1) | statu.get(2)).any()
-
-
-@pytest.mark.parametrize("missing", list(PIXELS))
-def test_read_datacube_without_visible_band_raises(scene, missing):
-    predictors = [band for band in SPECTRAL if band != missing]
-    with pytest.raises(ValueError, match=missing):
-        scene.read_datacube(predictors)
-
-
-@pytest.mark.parametrize("band", list(PIXELS))
-def test_read_datacube_zeroes_whiteness_of_saturated_pixels(scene, band):
-    _set_dn(scene, band, 50, 50, 65535)
-
-    data, _ = scene.read_datacube(SPECTRAL + ["whiteness"])
-
-    expected_nonzero = np.ones((TARGET, TARGET), dtype=bool)
-    expected_nonzero[50, 50] = False
-    assert np.array_equal(data.data[data.bands.index("whiteness")] != 0, expected_nonzero)
diff --git a/tests/test_unet_classify.py b/tests/test_unetlib.py
similarity index 100%
rename from tests/test_unet_classify.py
rename to tests/test_unetlib.py
diff --git a/tests/test_normalize.py b/tests/test_utils.py
similarity index 77%
rename from tests/test_normalize.py
rename to tests/test_utils.py
index 7c3de12..f0f7f0c 100644
--- a/tests/test_normalize.py
+++ b/tests/test_utils.py
@@ -1,5 +1,6 @@
 import numpy as np
 import pytest
+from skimage.morphology import reconstruction
 
 from fmask import constant as C
 from fmask import utils
@@ -162,3 +163,50 @@ def test_predictor_order_matches_data_get():
         utils.normalize_bands(data.data, [data.bands.index(p) for p in predictors], obsmask=mask),
         _bands(data.get(predictors), obsmask=mask),
     )
+
+
+def _random(dtype, seed=0, shape=(40, 57)):
+    return np.random.default_rng(seed).random(shape).astype(dtype)
+
+
+@pytest.mark.parametrize("fill_value", [None, -1.0])
+@pytest.mark.parametrize("obs_shape", ["full", "partial", "single"])
+def test_imfill_matches_skimage(obs_shape, fill_value):
+    data = _random(np.float32, 11, (31, 37)) * 10
+    obsmask = np.ones(data.shape, dtype=bool)
+    if obs_shape == "partial":
+        obsmask[:5] = False
+        obsmask[:, -4:] = False
+    elif obs_shape == "single":
+        obsmask[:] = False
+        obsmask[15, 15] = True
+    max_val = data[obsmask].max()
+    filled_input = data.copy()
+    filled_input[~obsmask] = max_val if fill_value is None else fill_value
+    seed = filled_input.copy()
+    seed[1:-1, 1:-1] = max_val
+
+    result = utils.imfill(data, obsmask, fill_value=fill_value)
+
+    _assert_same(result, reconstruction(seed, filled_input, method="erosion"))
+    _assert_same(data, filled_input)
+
+
+def test_imfill_fills_a_closed_hole():
+    data = np.full((5, 5), 5.0, dtype=np.float32)
+    data[2, 2] = 1.0
+    result = utils.imfill(data, np.ones((5, 5), dtype=bool))
+    assert result[2, 2] == 5.0
+
+
+def test_imfill_keeps_a_depression_open_to_the_border():
+    data = np.full((5, 5), 5.0, dtype=np.float32)
+    data[2, :3] = 1.0
+    result = utils.imfill(data, np.ones((5, 5), dtype=bool))
+    assert (result[2, :3] == 1.0).all()
+
+
+def test_imfill_empty_observation_mask_raises():
+    data = _random(np.float32, 12, (8, 8))
+    with pytest.raises(ValueError):
+        utils.imfill(data, np.zeros(data.shape, dtype=bool))

From b9c0b8f1a3026f5452fbcf9ce2a56fea856596a6 Mon Sep 17 00:00:00 2001
From: Chris Holden <ceholden@gmail.com>
Date: Thu, 5 Nov 2026 10:13:37 -0500
Subject: [PATCH 10/11] test: make the signed zero case reach the skimage
 fallback

With values from 0 to 2, zero is the smallest value, so removing the signbit guard
in reconstruct_erosion left the results unchanged. Adding negative values puts
signed zeros mid-range, where the order of equal values matters, so the case now
fails without the guard.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 tests/test_reconstruct.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_reconstruct.py b/tests/test_reconstruct.py
index 2762780..428c083 100644
--- a/tests/test_reconstruct.py
+++ b/tests/test_reconstruct.py
@@ -30,7 +30,7 @@ def _with_nan(dtype):
 
 def _signed_zeros(dtype):
     rng = np.random.default_rng(3)
-    mask = rng.integers(0, 3, (40, 57)).astype(dtype)
+    mask = rng.integers(-1, 3, (40, 57)).astype(dtype)
     mask[mask == 0] = rng.choice([-0.0, 0.0], size=int((mask == 0).sum()))
     return mask
 

From 45ca94649d8cbe0343b4450ce08da7d096241159 Mon Sep 17 00:00:00 2001
From: Chris Holden <ceholden@gmail.com>
Date: Thu, 5 Nov 2026 10:14:01 -0500
Subject: [PATCH 11/11] test: remove duplicate and low-value tests

Removes tests that repeat another test or pin incidental details, and trims
parametrization that adds cases without adding behavior, from 178 cases to 106.
Test_read_bands_limits_results_waiting_behind_a_slow_band is replaced by a test
that counts submitted reads after each result, so it needs no sleep.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 tests/test_reconstruct.py        | 54 +++++---------------------
 tests/test_satellite.py          | 64 +++++++++----------------------
 tests/test_satellite_landsat.py  | 33 +++-------------
 tests/test_satellite_sentinel.py | 62 +++++++++++-------------------
 tests/test_unetlib.py            | 33 ++++++----------
 tests/test_utils.py              | 65 ++++++++------------------------
 6 files changed, 83 insertions(+), 228 deletions(-)

diff --git a/tests/test_reconstruct.py b/tests/test_reconstruct.py
index 428c083..b404daa 100644
--- a/tests/test_reconstruct.py
+++ b/tests/test_reconstruct.py
@@ -39,12 +39,6 @@ def _only_negative_zero(dtype):
     return np.full((9, 11), -0.0, dtype=dtype)
 
 
-def _positive_zero_and_values(dtype):
-    mask = _ties(dtype)
-    mask[0, 0] = 0.0
-    return mask
-
-
 def _infinities(dtype):
     mask = _random(dtype, 4)
     mask[3, 3] = -np.inf
@@ -63,7 +57,6 @@ def _negative_values(dtype):
     "nan": _with_nan,
     "signed_zeros": _signed_zeros,
     "only_negative_zero": _only_negative_zero,
-    "positive_zero": _positive_zero_and_values,
     "infinities": _infinities,
     "negative_values": _negative_values,
     "one_by_one": lambda dtype: np.array([[1.5]], dtype=dtype),
@@ -89,8 +82,15 @@ def test_matches_skimage(case, dtype):
     )
 
 
-@pytest.mark.parametrize("seed_dtype", [np.float32, np.float64])
-@pytest.mark.parametrize("mask_dtype", [np.uint8, np.int16, np.bool_, np.float16])
+@pytest.mark.parametrize(
+    "mask_dtype,seed_dtype",
+    [
+        (np.uint8, np.float32),
+        (np.int16, np.float64),
+        (np.bool_, np.float32),
+        (np.float16, np.float64),
+    ],
+)
 def test_matches_skimage_integer_and_small_types(mask_dtype, seed_dtype):
     mask = np.random.default_rng(8).integers(0, 2, (20, 23)).astype(mask_dtype)
     seed = _hole_seed(mask.astype(seed_dtype))
@@ -100,26 +100,6 @@ def test_matches_skimage_integer_and_small_types(mask_dtype, seed_dtype):
     )
 
 
-def test_matches_skimage_on_imfill_input():
-    mask = _random(np.float32, 9, (61, 47)) * 100
-    mask[10:30, 10:30] -= 50
-    seed = mask.copy()
-    seed[1:-1, 1:-1] = mask.max()
-    _assert_same(
-        reconstruct_erosion(seed, mask),
-        reconstruction(seed, mask, method="erosion"),
-    )
-
-
-def test_leaves_inputs_unchanged():
-    mask = _random(np.float32, 10)
-    seed = _hole_seed(mask)
-    mask_before, seed_before = mask.copy(), seed.copy()
-    reconstruct_erosion(seed, mask)
-    _assert_same(mask, mask_before)
-    _assert_same(seed, seed_before)
-
-
 def test_seed_below_mask_raises_same_error():
     mask = np.ones((5, 5), dtype=np.float32)
     seed = np.zeros((5, 5), dtype=np.float32)
@@ -130,22 +110,8 @@ def test_seed_below_mask_raises_same_error():
     assert str(actual.value) == str(expected.value)
 
 
-def test_shape_mismatch_raises_same_error():
-    with pytest.raises(AssertionError):
-        reconstruction(np.ones((5, 5)), np.ones((4, 5)), method="erosion")
-    with pytest.raises(AssertionError):
-        reconstruct_erosion(np.ones((5, 5)), np.ones((4, 5)))
-
-
-def test_private_api_available():
+def test_mirror_is_available():
     assert _reconstruct.MIRROR_AVAILABLE
-    from skimage.morphology._grayreconstruct import reconstruction_loop
-
-    ranks = np.array([1, 1, 1, 1, 1, 1, 1, 1, 0], dtype=np.uint32)
-    prev = np.full(9, -1, np.int32)
-    next_ = np.full(9, -1, np.int32)
-    strides = np.array([1, -1], np.int32)
-    reconstruction_loop(ranks, prev, next_, strides, 8, 3)
 
 
 @pytest.mark.parametrize(
diff --git a/tests/test_satellite.py b/tests/test_satellite.py
index 160cd24..0b2c781 100644
--- a/tests/test_satellite.py
+++ b/tests/test_satellite.py
@@ -1,7 +1,6 @@
 import gc
 import inspect
 import threading
-import time
 import weakref
 
 import numpy as np
@@ -42,14 +41,6 @@ def test_sample_matches_get():
     assert actual.T.tobytes() == expected.T.tobytes()
 
 
-def test_sample_is_copy():
-    data = _data()
-    before = data.data.copy()
-    rows, cols = _coords()
-    data.sample(["nir", "blue"], rows, cols)[:] = -1
-    assert np.array_equal(data.data, before)
-
-
 def test_read_bands_order(monkeypatch):
     image = object.__new__(_Reader)
     release = threading.Event()
@@ -66,27 +57,6 @@ def fake_values(self, band, **kwargs):
     assert list(image._read_bands(requests)) == ["first", "second"]
 
 
-def test_read_bands_limits_threads(monkeypatch):
-    image = object.__new__(_Reader)
-    lock = threading.Lock()
-    active = [0, 0]
-    barrier = threading.Barrier(satellite.READ_THREADS)
-
-    def fake_values(self, band, **kwargs):
-        with lock:
-            active[0] += 1
-            active[1] = max(active)
-        barrier.wait(timeout=10)
-        with lock:
-            active[0] -= 1
-        return band
-
-    monkeypatch.setattr(_Reader, "_read_band_values", fake_values)
-    count = satellite.READ_THREADS * 3
-    list(image._read_bands([{"band": str(i), "level": "TOA"} for i in range(count)]))
-    assert active[1] == satellite.READ_THREADS
-
-
 def test_read_bands_propagates_errors(monkeypatch):
     image = object.__new__(_Reader)
 
@@ -131,24 +101,26 @@ def fake_values(self, band, **kwargs):
     list(results)
 
 
-def test_read_bands_limits_results_waiting_behind_a_slow_band(monkeypatch):
+def test_read_band_arguments_match_landsat():
+    assert inspect.signature(Sentinel2.read_band) == inspect.signature(Landsat.read_band)
+
+
+def test_read_bands_starts_a_read_only_as_results_are_consumed(monkeypatch):
     image = object.__new__(_Reader)
-    started = []
-    started_during_slow_band = []
+    submitted = []
 
-    def fake_values(self, band, **kwargs):
-        started.append(band)
-        if band == "0":
-            time.sleep(0.3)
-            started_during_slow_band.append(len(started))
-        return band
+    class CountingPool(satellite.ThreadPoolExecutor):
+        def submit(self, *args, **kwargs):
+            submitted.append(args)
+            return super().submit(*args, **kwargs)
 
-    monkeypatch.setattr(_Reader, "_read_band_values", fake_values)
+    monkeypatch.setattr(satellite, "ThreadPoolExecutor", CountingPool)
+    monkeypatch.setattr(_Reader, "_read_band_values", lambda self, band, **kwargs: band)
     count = satellite.READ_THREADS * 3
-    requests = [{"band": str(i), "level": "TOA"} for i in range(count)]
-    assert list(image._read_bands(requests)) == [str(i) for i in range(count)]
-    assert started_during_slow_band == [satellite.READ_THREADS]
-
+    results = image._read_bands([{"band": str(i), "level": "TOA"} for i in range(count)])
 
-def test_read_band_arguments_match_landsat():
-    assert inspect.signature(Sentinel2.read_band) == inspect.signature(Landsat.read_band)
+    assert next(results) == "0"
+    assert len(submitted) == satellite.READ_THREADS
+    assert next(results) == "1"
+    assert len(submitted) == satellite.READ_THREADS + 1
+    assert list(results) == [str(i) for i in range(2, count)]
diff --git a/tests/test_satellite_landsat.py b/tests/test_satellite_landsat.py
index 9cde80b..70294d2 100644
--- a/tests/test_satellite_landsat.py
+++ b/tests/test_satellite_landsat.py
@@ -1,5 +1,3 @@
-import sys
-
 import numpy as np
 import pandas
 import pytest
@@ -7,7 +5,7 @@
 from rasterio.transform import from_origin
 
 from fmask import constant as C
-from fmask import satellite, utils
+from fmask import utils
 from fmask.satellite import Landsat
 
 SHAPE = (48, 64)
@@ -94,20 +92,6 @@ def test_values_do_not_touch_self(scene):
     assert not valid[1, :].any()
 
 
-def test_read_datacube_releases_read_inputs_before_indices(scene, monkeypatch):
-    original_hot = satellite.hot
-    seen = {}
-
-    def spy(*args):
-        seen.update(sys._getframe(1).f_locals)
-        return original_hot(*args)
-
-    monkeypatch.setattr(satellite, "hot", spy)
-    _make(Landsat, scene).read_datacube(SPECTRAL + ["hot"])
-    assert seen
-    assert not {"request", "requests", "cos_sza", "saturated", "image", "satu", "valid"} & set(seen)
-
-
 def _raw(scene, band):
     filename = scene["metadata"]["LEVEL1_PROCESSING_RECORD"][f"FILE_NAME_BAND_{BAND_IDS[band.upper()]}"]
     with rasterio.open(f"{scene['folder']}/{filename}") as src:
@@ -194,16 +178,11 @@ def test_read_band_profile(scene):
     assert image.profilefull.closed
 
 
-def test_obsmask_all_bands(scene):
-    image = _make(Landsat, scene)
-    image.read_datacube(list(SPECTRAL))
-    _assert_same(image.obsmask, _expected_obsmask(scene, SPECTRAL))
-    assert not image.obsmask.all()
-    assert image.obsmask.any()
-
-
-@pytest.mark.parametrize("maskupdate", [True, False])
-@pytest.mark.parametrize("predictors", [SPECTRAL, list(reversed(SPECTRAL))], ids=["forward", "reversed"])
+@pytest.mark.parametrize(
+    "predictors,maskupdate",
+    [(SPECTRAL, True), (list(reversed(SPECTRAL)), False)],
+    ids=["forward_with_mask", "reversed_without_mask"],
+)
 def test_read_datacube_matches_each_band(scene, predictors, maskupdate):
     image = _make(Landsat, scene)
     data, _ = image.read_datacube(list(predictors), maskupdate=maskupdate)
diff --git a/tests/test_satellite_sentinel.py b/tests/test_satellite_sentinel.py
index 3b9e131..25363a7 100644
--- a/tests/test_satellite_sentinel.py
+++ b/tests/test_satellite_sentinel.py
@@ -85,10 +85,14 @@ def scene(tmp_path):
     return _write_scene(tmp_path, OFFSETS, nodata_lines=True)
 
 
-@pytest.fixture(params=[-1000.0, -1001.0, 0.0], ids=["offset_1000", "offset_1001", "no_offset"])
-def clear_image(tmp_path, request):
-    offsets = {name: request.param for name, _, _ in BANDS}
-    return _make(_write_scene(tmp_path, offsets, nodata_lines=False))
+def _clear_image(folder, offset):
+    offsets = {name: offset for name, _, _ in BANDS}
+    return _make(_write_scene(folder, offsets, nodata_lines=False))
+
+
+@pytest.fixture
+def clear_image(tmp_path):
+    return _clear_image(tmp_path, -1000.0)
 
 
 def _assert_same(actual, expected):
@@ -121,10 +125,9 @@ def test_read_band_converts_to_toa(scene, band):
     assert satu is None
 
 
-@pytest.mark.parametrize("band", ["vre1", "blue", "coastal"], ids=["equal", "average", "nearest"])
-def test_read_band_returns_dn_for_other_levels(scene, band):
-    result, _ = _make(scene).read_band(band, level="DN")
-    _assert_same(result, _raw(scene, band))
+def test_read_band_returns_dn_for_other_levels(scene):
+    result, _ = _make(scene).read_band("blue", level="DN")
+    _assert_same(result, _raw(scene, "blue"))
     assert result.shape == (TARGET, TARGET)
 
 
@@ -163,16 +166,11 @@ def test_values_do_not_touch_self(scene):
     assert valid.dtype == np.bool_
 
 
-def test_obsmask_all_bands(scene):
-    image = _make(scene)
-    image.read_datacube(list(SPECTRAL))
-    _assert_same(image.obsmask, _expected_obsmask(scene, SPECTRAL))
-    assert not image.obsmask.all()
-    assert image.obsmask.any()
-
-
-@pytest.mark.parametrize("maskupdate", [True, False])
-@pytest.mark.parametrize("predictors", [SPECTRAL, list(reversed(SPECTRAL))], ids=["forward", "reversed"])
+@pytest.mark.parametrize(
+    "predictors,maskupdate",
+    [(SPECTRAL, True), (list(reversed(SPECTRAL)), False)],
+    ids=["forward_with_mask", "reversed_without_mask"],
+)
 def test_read_datacube_matches_each_band(scene, predictors, maskupdate):
     image = _make(scene)
     data, _ = image.read_datacube(list(predictors), maskupdate=maskupdate)
@@ -209,26 +207,19 @@ def _set_dn(image, band, row, col, value):
     _write(path, data, res, "JP2OpenJPEG")
 
 
-@pytest.mark.parametrize("level", ["TOA", "DN"])
-def test_read_band_returns_saturated_pixels(clear_image, level):
-    _set_dn(clear_image, "blue", 50, 50, 65535)
-    _set_dn(clear_image, "blue", 40, 40, 65534)
+@pytest.mark.parametrize("offset", [-1000.0, -1001.0, 0.0])
+def test_read_band_returns_saturated_pixels(tmp_path, offset):
+    image = _clear_image(tmp_path, offset)
+    _set_dn(image, "blue", 50, 50, 65535)
+    _set_dn(image, "blue", 40, 40, 65534)
 
-    _, satu = clear_image.read_band("blue", level=level, saturation=True)
+    _, satu = image.read_band("blue", saturation=True)
 
     expected = np.zeros((TARGET, TARGET), dtype=bool)
     expected[50, 50] = True
     assert np.array_equal(satu, expected)
 
 
-def test_read_band_leaves_saturation_out_by_default(clear_image):
-    _set_dn(clear_image, "blue", 50, 50, 65535)
-
-    _, satu = clear_image.read_band("blue")
-
-    assert satu is None
-
-
 def test_read_datacube_sets_saturation_bits_in_blue_green_red_order(clear_image):
     for band, (row, col) in PIXELS.items():
         _set_dn(clear_image, band, row, col, 65535)
@@ -242,15 +233,6 @@ def test_read_datacube_sets_saturation_bits_in_blue_green_red_order(clear_image)
         assert np.array_equal(statu.get(bit), expected)
 
 
-def test_read_datacube_does_not_flag_the_dn_below_saturation(clear_image):
-    for band in PIXELS:
-        _set_dn(clear_image, band, 50, 50, 65534)
-
-    _, statu = clear_image.read_datacube(list(SPECTRAL))
-
-    assert not (statu.get(0) | statu.get(1) | statu.get(2)).any()
-
-
 @pytest.mark.parametrize("missing", list(PIXELS))
 def test_read_datacube_without_visible_band_raises(clear_image, missing):
     predictors = [band for band in SPECTRAL if band != missing]
diff --git a/tests/test_unetlib.py b/tests/test_unetlib.py
index 2a7dea7..2fc344b 100644
--- a/tests/test_unetlib.py
+++ b/tests/test_unetlib.py
@@ -44,9 +44,17 @@ def _whole_image_probs(stub, normalize):
         return torch.softmax(stub.model(torch.from_numpy(cube)[None]), dim=1)[0].numpy()
 
 
-@pytest.mark.parametrize("batch_size", [1, 3, 500], ids=["single", "short_last_batch", "one_batch"])
-@pytest.mark.parametrize("normalize", [True, False])
-@pytest.mark.parametrize("probability", ["none", "default", "cloud"])
+@pytest.mark.parametrize(
+    "probability,normalize,batch_size",
+    [
+        ("none", True, 3),
+        ("default", True, 3),
+        ("cloud", True, 3),
+        ("cloud", True, 1),
+        ("cloud", True, 500),
+        ("cloud", False, 3),
+    ],
+)
 def test_classify_matches_per_pixel_scores(probability, normalize, batch_size):
     stub = _stub(batch_size)
     cube_before = stub.image.data.data.copy()
@@ -80,17 +88,6 @@ def _anchors(count):
     return [(i * 2, i * 3, 0, 0) for i in range(count)]
 
 
-def test_fill_batch_matches_stack():
-    datacube = _datacube()
-    anchors = _anchors(4)
-    buffer = torch.empty((4, 3, 32, 32), dtype=datacube.dtype)
-    actual = fill_batch(buffer, datacube, anchors, 32)
-    expected = _stack(datacube, anchors, 32)
-    assert actual.dtype == expected.dtype
-    assert actual.shape == expected.shape
-    assert actual.numpy().tobytes() == expected.numpy().tobytes()
-
-
 def test_short_last_batch():
     datacube = _datacube()
     buffer = torch.empty((4, 3, 32, 32), dtype=datacube.dtype)
@@ -100,11 +97,3 @@ def test_short_last_batch():
     expected = _stack(datacube, anchors, 32)
     assert actual.shape == (2, 3, 32, 32)
     assert actual.numpy().tobytes() == expected.numpy().tobytes()
-
-
-def test_fill_batch_does_not_change_datacube():
-    datacube = _datacube()
-    before = datacube.clone()
-    buffer = torch.empty((4, 3, 32, 32), dtype=datacube.dtype)
-    fill_batch(buffer, datacube, _anchors(4), 32)
-    assert torch.equal(datacube, before)
diff --git a/tests/test_utils.py b/tests/test_utils.py
index f0f7f0c..42283cf 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -4,7 +4,6 @@
 
 from fmask import constant as C
 from fmask import utils
-from fmask.satellite import Data
 
 SHAPE = (4, 23, 31)
 
@@ -63,20 +62,18 @@ def test_percentiles_ignore_pixels_outside_the_mask():
     _assert_same(_bands(with_outliers, obsmask=mask), _bands(cube, obsmask=mask))
 
 
-@pytest.mark.parametrize("srange", [(-1, 1), (0, 1), (-3, 7)])
-def test_pixels_outside_the_mask_get_the_range_minimum(srange):
-    result = _bands(_cube(3), obsmask=_partial_mask(), srange=srange)
-    assert (result[:, ~_partial_mask()] == min(srange)).all()
+def test_pixels_outside_the_mask_get_the_range_minimum():
+    result = _bands(_cube(3), obsmask=_partial_mask(), srange=(-3, 7))
+    assert (result[:, ~_partial_mask()] == -3).all()
 
 
-@pytest.mark.parametrize("percentiles,srange", [((1, 99), (-1, 1)), ((5, 50), (-3, 7))])
-def test_values_inside_the_mask_stay_in_range(percentiles, srange):
+def test_values_inside_the_mask_stay_in_range():
     mask = _partial_mask()
-    result = _bands(_cube(4), obsmask=mask, percentiles=percentiles, srange=srange)
-    assert result.min() >= srange[0]
-    assert result.max() <= srange[1]
-    assert result[:, mask].min() == srange[0]
-    assert result[:, mask].max() == srange[1]
+    result = _bands(_cube(4), obsmask=mask, percentiles=(5, 50), srange=(-3, 7))
+    assert result.min() >= -3
+    assert result.max() <= 7
+    assert result[:, mask].min() == -3
+    assert result[:, mask].max() == 7
 
 
 def test_constant_band_gets_the_range_minimum():
@@ -105,19 +102,17 @@ def test_nan_outside_the_mask_is_replaced():
 
 
 @pytest.mark.parametrize("use_mask", [True, False])
-@pytest.mark.parametrize("block_rows", [1, 7, 1000])
-def test_block_size_does_not_change_results(block_rows, use_mask):
+def test_block_size_does_not_change_results(use_mask):
     cube = _cube(8)
     obsmask = _partial_mask() if use_mask else None
     _assert_same(
-        _bands(cube, obsmask=obsmask, block_rows=block_rows),
+        _bands(cube, obsmask=obsmask, block_rows=7),
         _bands(cube, obsmask=obsmask, block_rows=SHAPE[1]),
     )
 
 
-@pytest.mark.parametrize("dtype", [np.float32, np.float64])
-def test_result_keeps_the_input_dtype(dtype):
-    assert _bands(_cube(9, dtype)).dtype == dtype
+def test_result_keeps_the_input_dtype():
+    assert _bands(_cube(9, np.float64)).dtype == np.float64
 
 
 def test_normalize_datacube_matches_normalize_bands():
@@ -126,11 +121,6 @@ def test_normalize_datacube_matches_normalize_bands():
     _assert_same(utils.normalize_datacube(cube), _bands(cube))
 
 
-def test_empty_obsmask_raises():
-    with pytest.raises(IndexError):
-        _bands(_cube(11), obsmask=np.zeros(SHAPE[1:], dtype=bool))
-
-
 @pytest.mark.filterwarnings("ignore::RuntimeWarning")
 def test_input_not_mutated():
     mask = _partial_mask()
@@ -142,11 +132,6 @@ def test_input_not_mutated():
     _assert_same(mask, mask_before)
 
 
-def test_result_does_not_share_memory_with_input():
-    cube = _cube(13)
-    assert not np.shares_memory(_bands(cube), cube)
-
-
 def test_bands_subset_follows_the_requested_order():
     cube, mask = _cube(14), _partial_mask()
     _assert_same(
@@ -155,31 +140,19 @@ def test_bands_subset_follows_the_requested_order():
     )
 
 
-def test_predictor_order_matches_data_get():
-    cube, mask = _cube(15), _partial_mask()
-    data = Data(cube, ["blue", "green", "red", "nir"])
-    predictors = ["nir", "blue", "red"]
-    _assert_same(
-        utils.normalize_bands(data.data, [data.bands.index(p) for p in predictors], obsmask=mask),
-        _bands(data.get(predictors), obsmask=mask),
-    )
-
-
 def _random(dtype, seed=0, shape=(40, 57)):
     return np.random.default_rng(seed).random(shape).astype(dtype)
 
 
-@pytest.mark.parametrize("fill_value", [None, -1.0])
-@pytest.mark.parametrize("obs_shape", ["full", "partial", "single"])
+@pytest.mark.parametrize(
+    "obs_shape,fill_value", [("full", None), ("partial", None), ("partial", -1.0)]
+)
 def test_imfill_matches_skimage(obs_shape, fill_value):
     data = _random(np.float32, 11, (31, 37)) * 10
     obsmask = np.ones(data.shape, dtype=bool)
     if obs_shape == "partial":
         obsmask[:5] = False
         obsmask[:, -4:] = False
-    elif obs_shape == "single":
-        obsmask[:] = False
-        obsmask[15, 15] = True
     max_val = data[obsmask].max()
     filled_input = data.copy()
     filled_input[~obsmask] = max_val if fill_value is None else fill_value
@@ -204,9 +177,3 @@ def test_imfill_keeps_a_depression_open_to_the_border():
     data[2, :3] = 1.0
     result = utils.imfill(data, np.ones((5, 5), dtype=bool))
     assert (result[2, :3] == 1.0).all()
-
-
-def test_imfill_empty_observation_mask_raises():
-    data = _random(np.float32, 12, (8, 8))
-    with pytest.raises(ValueError):
-        utils.imfill(data, np.zeros(data.shape, dtype=bool))