From 2268e319ce330d44bcc85259d53f8822f49df1c5 Mon Sep 17 00:00:00 2001 From: thc1006 <84045975+thc1006@users.noreply.github.com> Date: Mon, 20 Jul 2026 01:21:17 +0800 Subject: [PATCH 1/2] Make quality_mask safe for the EMIT L2A Mask V002 product quality_mask hard-coded bands 5 and 6 (the V001 AOD/H2O data bands) as the only non-flag layers. The V002 mask adds three continuous layers -- SpecTf-Cloud Probability and SpecTf-Buffer Distance -- and reorders the flags, so requesting one of the new data bands was silently treated as a binary flag: its values were summed and clipped to {0, 1}, producing a semantically meaningless mask. Classify each requested layer from its sensor_band_parameters/mask_bands name instead. There is no separate flag-vs-data metadata field in either product, so the band name is the only per-layer discriminator: - binary flags ("... Flag") are OR-combined, exactly as before; - continuous data layers (AOD550, H2O, SpecTf-Cloud Probability, SpecTf-Buffer Distance) are rejected with an error naming the band; - a probability layer can be turned into a mask via the new optional threshold argument (probability >= threshold); - band indices are validated and the returned mask is asserted to be {0, 1}. Verified against a real EMITL2AMASK V002 granule: the binary-flag path is value-identical to the previous implementation, while the new SpecTf data bands are now rejected instead of silently accepted. Add pytest regression tests covering the V001 (8-band) and V002 (11-band) layouts, including the previously-silent failure. Co-authored-by: junnncct1106 --- CHANGE_LOG.md | 10 + python/modules/emit_tools.py | 118 ++++++++- python/modules/tests/conftest.py | 42 ++++ python/modules/tests/test_quality_mask.py | 281 ++++++++++++++++++++++ 4 files changed, 439 insertions(+), 12 deletions(-) create mode 100644 python/modules/tests/conftest.py create mode 100644 python/modules/tests/test_quality_mask.py diff --git a/CHANGE_LOG.md b/CHANGE_LOG.md index d67990e..62aa9f1 100644 --- a/CHANGE_LOG.md +++ b/CHANGE_LOG.md @@ -5,6 +5,16 @@ The format is based on [Keep a Changelog](http://keepachangelog.com/) and this project adheres to [Semantic Versioning](http://semver.org/). _________________________________________________________________________ +## 2026-07-19 + +> ### Changed +> +> - Made `quality_mask` in `python/modules/emit_tools.py` metadata-driven so it is safe for the EMIT L2A Mask V002 product. It now classifies each requested layer from its `sensor_band_parameters/mask_bands` name and refuses to build a mask from a continuous data layer (`AOD550`, `H2O (g cm-2)`, `SpecTf-Cloud Probability`, `SpecTf-Buffer Distance`). Previously only bands 5 and 6 were hard-coded as data bands, so the new V002 continuous bands were silently treated as flags. A continuous probability layer can be turned into a mask with the new optional `threshold` argument. + +> ### Added +> +> - `python/modules/tests/` with pytest regression tests for `quality_mask` covering the V001 (8-band) and V002 (11-band) mask layouts. + ## 2024-06-17 > ### Added diff --git a/python/modules/emit_tools.py b/python/modules/emit_tools.py index e4acf45..df7d639 100644 --- a/python/modules/emit_tools.py +++ b/python/modules/emit_tools.py @@ -264,33 +264,127 @@ def ortho_xr(ds, GLT_NODATA_VALUE=0, fill_value=-9999): return out_xr -def quality_mask(filepath, quality_bands): +def quality_mask(filepath, quality_bands, threshold=None): """ - This function builds a single layer mask to apply based on the bands selected from an EMIT L2A Mask file. + This function builds a single binary mask to apply, based on the flag bands selected from an EMIT L2A Mask file. + + The layers in an EMIT L2A Mask file fall into two categories: + + - binary quality flags (e.g. `Cloud Flag`, `Cirrus Flag`, `Water Flag`, + `Spacecraft Flag`, `Dilated Cloud Flag`, `Aggregate Flag`, and, in the + V002 product, `SpecTf-Cloud Flag`), whose values are 0 or 1 and which can + be OR-combined into a mask, and + - continuous data layers (e.g. `AOD550`, `H2O (g cm-2)`, and, in the V002 + product, `SpecTf-Cloud Probability` and `SpecTf-Buffer Distance`), which + are physical quantities, not masks. + + The mask file distinguishes these two categories only through the band + *names* stored in `sensor_band_parameters/mask_bands`; there is no separate + "is-a-flag" metadata field in either the V001 or the V002 product. Summing a + continuous layer and clipping values > 1 to 1, as if it were a flag, silently + produces a semantically meaningless mask. This function therefore classifies + each requested band from its name and refuses to treat a continuous data + layer as a flag. The prior implementation hard-coded bands 5 and 6 (the V001 + AOD/H2O data bands) as the only non-flag layers, which silently accepted the + new V002 continuous bands (`SpecTf-Cloud Probability`, `SpecTf-Buffer + Distance`). Parameters: filepath: an EMIT L2A Mask netCDF file. - quality_bands: a list of bands (quality flags only) from the mask file that should be used in creation of mask. + quality_bands: an index or list of band indices (quality flags only) from the mask file that should be used in creation of the mask. Indices refer to the order of `sensor_band_parameters/mask_bands`. + threshold: optional float in [0, 1]. When provided, a continuous *probability* layer (e.g. `SpecTf-Cloud Probability`) included in `quality_bands` is converted to a binary mask using `probability >= threshold`. Non-probability data layers (AOD, water vapor, buffer distance) are never accepted. Returns: - qmask: a numpy array that can be used with the emit_xarray function to apply a quality mask. + qmask: a uint8 numpy array of {0, 1} that can be used with the emit_xarray function to apply a quality mask. """ # Open Dataset mask_ds = xr.open_dataset(filepath, engine="h5netcdf") - # Open Sensor band Group + # Open Sensor band Group, which names every layer in the mask file mask_parameters_ds = xr.open_dataset( filepath, engine="h5netcdf", group="sensor_band_parameters" ) + band_names = [ + str(b) for b in np.asarray(mask_parameters_ds["mask_bands"].data).ravel() + ] + n_bands = len(band_names) + + # Accept a single index or an iterable of indices + if isinstance(quality_bands, (int, np.integer)): + quality_bands = [quality_bands] + quality_bands = list(quality_bands) + + # Validate the requested indices against this file's actual bands + for b in quality_bands: + if not isinstance(b, (int, np.integer)) or b < 0 or b >= n_bands: + raise ValueError( + f"quality_bands index {b!r} is out of range for this mask file, " + f"which has {n_bands} bands: " + f"{[f'{i}: {name}' for i, name in enumerate(band_names)]}" + ) + + if threshold is not None and not (0.0 <= float(threshold) <= 1.0): + raise ValueError(f"threshold must be within [0, 1]; got {threshold!r}") + # Print Flags used - flags_used = mask_parameters_ds["mask_bands"].data[quality_bands] + flags_used = [band_names[b] for b in quality_bands] print(f"Flags used: {flags_used}") - # Check for data bands and build mask - if any(x in quality_bands for x in [5, 6]): - err_str = f"Selected flags include a data band (5 or 6) not just flag bands" - raise AttributeError(err_str) + + def _is_flag(name): + # Binary quality-flag layers are named "... Flag" in both V001 and V002. + return name.strip().lower().endswith("flag") + + def _is_probability(name): + return "probability" in name.lower() + + flag_band_hint = [ + f"{i}: {name}" for i, name in enumerate(band_names) if _is_flag(name) + ] + + # Build one binary layer per requested band + layers = [] + for b in quality_bands: + name = band_names[b] + layer = mask_ds["mask"][:, :, b].values + if _is_flag(name): + # Defense in depth: a flag layer must be binary. Never silently clip. + finite = layer[np.isfinite(layer)] + uniq = np.unique(finite) + if uniq.size and not np.all(np.isin(uniq, (0.0, 1.0))): + raise ValueError( + f"Band {b} ('{name}') is named as a flag but contains non-binary " + f"values (e.g. {uniq[:5]}); refusing to build a mask from it." + ) + layers.append((layer > 0).astype(np.uint8)) + elif _is_probability(name) and threshold is not None: + layers.append((layer >= float(threshold)).astype(np.uint8)) + else: + # Continuous / data band: refuse rather than silently clip it to {0, 1}. + hint = "" + if _is_probability(name): + hint = ( + " Pass threshold= to convert this probability " + "layer into a binary mask (probability >= threshold)." + ) + raise ValueError( + f"Band {b} ('{name}') is a continuous data layer, not a binary " + f"quality flag, so it cannot be combined into a mask.{hint} " + f"Flag bands available in this file: {flag_band_hint}." + ) + + # Combine the binary layers. Logical OR reproduces the previous + # sum-then-clip behavior exactly for binary flags. + if layers: + qmask = np.zeros_like(layers[0]) + for layer in layers: + qmask |= layer else: - qmask = np.sum(mask_ds["mask"][:, :, quality_bands].values, axis=-1) - qmask[qmask > 1] = 1 + qmask = np.zeros(mask_ds["mask"].shape[:2], dtype=np.uint8) + qmask[qmask > 1] = 1 + + # The returned mask must be binary; fail loudly if that invariant is broken. + assert set(int(v) for v in np.unique(qmask)).issubset( + {0, 1} + ), "quality_mask produced a non-binary mask" return qmask diff --git a/python/modules/tests/conftest.py b/python/modules/tests/conftest.py new file mode 100644 index 0000000..30092b0 --- /dev/null +++ b/python/modules/tests/conftest.py @@ -0,0 +1,42 @@ +"""Test configuration for the ``emit_tools`` module tests. + +``emit_tools`` imports the full geospatial stack (``osgeo.gdal``, ``geopandas``, +``rasterio``, ``rioxarray``, ``scikit-image``, ``spectral`` ...) at module load +time. None of that is needed to exercise ``quality_mask``, which only uses +``numpy`` and ``xarray``. To let these unit tests run in a minimal environment +(e.g. CI without GDAL installed), any of those optional dependencies that are +not importable are replaced with lightweight stand-ins before ``emit_tools`` is +imported. When the full environment is present, the real modules are used and +nothing is stubbed. +""" + +import importlib +import sys +from pathlib import Path +from unittest import mock + +# Make ``import emit_tools`` resolve, matching the notebooks' convention of +# adding ``python/modules`` to ``sys.path``. +MODULES_DIR = Path(__file__).resolve().parents[1] +if str(MODULES_DIR) not in sys.path: + sys.path.insert(0, str(MODULES_DIR)) + +# Optional heavy dependencies that ``emit_tools`` imports but ``quality_mask`` +# does not need. Stub only the ones that are actually missing. +_OPTIONAL_DEPENDENCIES = ( + "osgeo", + "spectral", + "spectral.io", + "skimage", + "geopandas", + "rasterio", + "rioxarray", + "rioxarray.merge", + "s3fs", +) + +for _name in _OPTIONAL_DEPENDENCIES: + try: + importlib.import_module(_name) + except Exception: # ImportError or any import-time failure + sys.modules[_name] = mock.MagicMock() diff --git a/python/modules/tests/test_quality_mask.py b/python/modules/tests/test_quality_mask.py new file mode 100644 index 0000000..624ea36 --- /dev/null +++ b/python/modules/tests/test_quality_mask.py @@ -0,0 +1,281 @@ +"""Regression tests for :func:`emit_tools.quality_mask`. + +These tests are network-free: they build small synthetic EMIT L2A Mask netCDF +files whose ``sensor_band_parameters/mask_bands`` metadata mirrors the real +V001 (8-band) and V002 (11-band) products. The V002 band names and ordering +were verified against a real granule +(``EMIT_L2A_MASK_002_20220810T034103_2222203_001.nc``): + + 0 Cloud Flag (binary flag) + 1 Cirrus Flag (binary flag) + 2 Water Flag (binary flag) + 3 Spacecraft Flag (binary flag) + 4 Dilated Cloud Flag (binary flag) + 5 AOD550 (continuous data) + 6 H2O (g cm-2) (continuous data) + 7 Aggregate Flag (binary flag) + 8 SpecTf-Cloud Probability (continuous data, probability) + 9 SpecTf-Cloud Flag (binary flag) + 10 SpecTf-Buffer Distance (continuous data, distance) + +The bug being guarded against: the previous implementation hard-coded bands 5 +and 6 as the only non-flag layers, so it silently accepted the new V002 +continuous bands (8 and 10) and returned a semantically meaningless "mask". + +An optional test also runs against a real granule when the environment variable +``EMIT_L2A_MASK_V002`` points to one. +""" + +import os + +import numpy as np +import netCDF4 as nc +import pytest + +from emit_tools import quality_mask + + +# --- Verified real band layouts ------------------------------------------------- + +V002_MASK_BANDS = [ + "Cloud Flag", + "Cirrus Flag", + "Water Flag", + "Spacecraft Flag", + "Dilated Cloud Flag", + "AOD550", + "H2O (g cm-2)", + "Aggregate Flag", + "SpecTf-Cloud Probability", + "SpecTf-Cloud Flag", + "SpecTf-Buffer Distance", +] + +# Representative EMIT L2A Mask V001 layout (8 bands). V001 mixed the +# capitalization of "flag"/"Flag", which the classifier must handle +# case-insensitively. AOD550/H2O sit at indices 5/6, matching the constants the +# previous implementation hard-coded. +V001_MASK_BANDS = [ + "Cloud flag", + "Dilated Cloud Flag", + "Cirrus flag", + "Water flag", + "Spacecraft Flag", + "AOD550", + "H2O (g cm-2)", + "Aggregate Flag", +] + + +def _write_mask_nc(path, band_names, mask_array): + """Write a minimal EMIT-like L2A Mask netCDF file. + + Structure mirrors the real product: a root ``mask`` variable with dims + ``(downtrack, crosstrack, bands)`` and a ``sensor_band_parameters`` group + holding the ``mask_bands`` string variable. + """ + downtrack, crosstrack, n_bands = mask_array.shape + assert n_bands == len(band_names) + with nc.Dataset(path, "w", format="NETCDF4") as ds: + ds.createDimension("downtrack", downtrack) + ds.createDimension("crosstrack", crosstrack) + ds.createDimension("bands", n_bands) + mask_var = ds.createVariable( + "mask", "f4", ("downtrack", "crosstrack", "bands") + ) + mask_var[:] = mask_array.astype("f4") + mask_var.long_name = "Masks" + mask_var.units = "unitless" + grp = ds.createGroup("sensor_band_parameters") + band_var = grp.createVariable("mask_bands", str, ("bands",)) + for i, name in enumerate(band_names): + band_var[i] = name + band_var.long_name = "Mask Band Names" + return str(path) + + +def _binary(pattern): + return np.array(pattern, dtype="f4") + + +# Reusable binary flag patterns on a 3x4 grid. +_A = _binary([[1, 0, 0, 0], [0, 0, 1, 0], [0, 0, 0, 0]]) +_B = _binary([[0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]) +_C = _binary([[1, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 0]]) +_D = _binary([[1, 0, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]) +_E = _binary([[0, 1, 0, 1], [0, 1, 1, 0], [1, 0, 0, 0]]) +_ZERO = np.zeros((3, 4), dtype="f4") + +# Continuous layers (deliberately include values > 1 and fractional values). +_AOD = _binary([[-0.04, 0.20, 0.50, 1.10], [0.30, 0.80, 0.95, 1.05], [0.00, 0.10, 0.90, 0.99]]) +_H2O = _binary([[0.05, 1.20, 2.00, 3.60], [0.50, 1.10, 0.90, 2.20], [0.30, 0.40, 1.50, 0.80]]) +_PROB = _binary([[0.30, 0.95, 0.50, 1.00], [0.26, 0.70, 0.99, 0.40], [0.80, 0.35, 0.60, 0.45]]) +# Buffer distance is conceptually continuous (pixel distance); include 2 and 3. +_DIST = _binary([[0, 1, 2, 3], [0, 0, 1, 2], [3, 2, 1, 0]]) + + +@pytest.fixture +def v002_mask(tmp_path): + layers = [_A, _B, _ZERO, _ZERO, _C, _AOD, _H2O, _D, _PROB, _E, _DIST] + arr = np.stack(layers, axis=-1) + path = _write_mask_nc(tmp_path / "v002_mask.nc", V002_MASK_BANDS, arr) + return {"path": path, "names": V002_MASK_BANDS, "array": arr} + + +@pytest.fixture +def v001_mask(tmp_path): + layers = [_A, _B, _C, _ZERO, _ZERO, _AOD, _H2O, _D] + arr = np.stack(layers, axis=-1) + path = _write_mask_nc(tmp_path / "v001_mask.nc", V001_MASK_BANDS, arr) + return {"path": path, "names": V001_MASK_BANDS, "array": arr} + + +def _expected_or(arr, bands): + """Reference mask: logical OR of the given (binary) band layers.""" + out = np.zeros(arr.shape[:2], dtype=np.uint8) + for b in bands: + out |= (arr[:, :, b] > 0).astype(np.uint8) + return out + + +def _assert_binary_mask(qmask): + assert isinstance(qmask, np.ndarray) + assert qmask.dtype == np.uint8 + assert set(int(v) for v in np.unique(qmask)).issubset({0, 1}) + + +# --- V001: legacy binary flags still work -------------------------------------- + +def test_v001_legacy_flags_unchanged(v001_mask): + # Mixed-case flag names (Cloud flag / Dilated Cloud Flag / Cirrus flag) + # must all be recognized as flags. + bands = [0, 1, 2] + qmask = quality_mask(v001_mask["path"], bands) + _assert_binary_mask(qmask) + np.testing.assert_array_equal(qmask, _expected_or(v001_mask["array"], bands)) + + +def test_v001_single_aggregate_flag(v001_mask): + qmask = quality_mask(v001_mask["path"], [7]) + _assert_binary_mask(qmask) + np.testing.assert_array_equal(qmask, _expected_or(v001_mask["array"], [7])) + + +@pytest.mark.parametrize("band, name", [(5, "AOD550"), (6, "H2O (g cm-2)")]) +def test_v001_data_bands_rejected(v001_mask, band, name): + with pytest.raises(ValueError) as exc: + quality_mask(v001_mask["path"], [band]) + assert name in str(exc.value) + + +def test_v001_mixed_flag_and_data_rejected(v001_mask): + with pytest.raises(ValueError) as exc: + quality_mask(v001_mask["path"], [0, 5]) + assert "AOD550" in str(exc.value) + + +# --- V002: binary flags work --------------------------------------------------- + +def test_v002_binary_flags_combined(v002_mask): + # Includes the new V002 binary SpecTf-Cloud Flag at index 9. + bands = [0, 1, 4, 7, 9] + qmask = quality_mask(v002_mask["path"], bands) + _assert_binary_mask(qmask) + np.testing.assert_array_equal(qmask, _expected_or(v002_mask["array"], bands)) + + +def test_v002_all_zero_flag(v002_mask): + # Water Flag (idx 2) is all zeros in this fixture. + qmask = quality_mask(v002_mask["path"], [2]) + _assert_binary_mask(qmask) + assert qmask.sum() == 0 + + +# --- V002: the silent-failure case now raises ---------------------------------- + +@pytest.mark.parametrize( + "band, name", + [(8, "SpecTf-Cloud Probability"), (10, "SpecTf-Buffer Distance")], +) +def test_v002_continuous_band_rejected(v002_mask, band, name): + with pytest.raises(ValueError) as exc: + quality_mask(v002_mask["path"], [band]) + assert name in str(exc.value) + + +def test_v002_old_behavior_was_silently_wrong(v002_mask): + """Document the bug: the previous sum-then-clip logic silently accepted a + continuous band and returned a non-binary array, while the fixed function + refuses it.""" + arr = v002_mask["array"] + # Reproduce the previous implementation on the SpecTf-Cloud Probability band. + legacy = np.sum(arr[:, :, [8]], axis=-1) + legacy[legacy > 1] = 1 + legacy_values = set(float(v) for v in np.unique(legacy)) + # The "mask" the old code returned contained fractional probabilities, i.e. + # it was never a valid {0, 1} mask. + assert not legacy_values.issubset({0.0, 1.0}) + # The fixed function refuses the same band instead. + with pytest.raises(ValueError): + quality_mask(v002_mask["path"], [8]) + + +# --- V002: probability accepted only via explicit threshold -------------------- + +def test_v002_probability_threshold(v002_mask): + qmask = quality_mask(v002_mask["path"], [8], threshold=0.5) + _assert_binary_mask(qmask) + expected = (v002_mask["array"][:, :, 8] >= 0.5).astype(np.uint8) + np.testing.assert_array_equal(qmask, expected) + + +def test_v002_probability_threshold_combined_with_flag(v002_mask): + qmask = quality_mask(v002_mask["path"], [0, 8], threshold=0.5) + _assert_binary_mask(qmask) + expected = (v002_mask["array"][:, :, 0] > 0).astype(np.uint8) + expected |= (v002_mask["array"][:, :, 8] >= 0.5).astype(np.uint8) + np.testing.assert_array_equal(qmask, expected) + + +def test_threshold_only_applies_to_probability(v002_mask): + # Buffer Distance is continuous but not a probability: threshold must not + # silently accept it. + with pytest.raises(ValueError): + quality_mask(v002_mask["path"], [10], threshold=0.5) + + +@pytest.mark.parametrize("bad", [-0.1, 1.5]) +def test_threshold_out_of_range(v002_mask, bad): + with pytest.raises(ValueError): + quality_mask(v002_mask["path"], [8], threshold=bad) + + +# --- Argument validation ------------------------------------------------------- + +def test_out_of_range_index_rejected(v002_mask): + with pytest.raises(ValueError): + quality_mask(v002_mask["path"], [11]) # only 0..10 exist + + +def test_single_int_accepted(v002_mask): + qmask = quality_mask(v002_mask["path"], 0) + _assert_binary_mask(qmask) + np.testing.assert_array_equal(qmask, _expected_or(v002_mask["array"], [0])) + + +# --- Optional: verify against a real downloaded granule ------------------------ + +@pytest.mark.skipif( + not os.environ.get("EMIT_L2A_MASK_V002"), + reason="set EMIT_L2A_MASK_V002 to a real EMIT L2A MASK V002 granule to run", +) +def test_real_v002_granule(): + fp = os.environ["EMIT_L2A_MASK_V002"] + # Binary flags combine into a valid mask. + qmask = quality_mask(fp, [0, 1, 4]) + _assert_binary_mask(qmask) + # SpecTf-Cloud Probability (index 8) is rejected unless a threshold is given. + with pytest.raises(ValueError): + quality_mask(fp, [8]) + qmask_thr = quality_mask(fp, [8], threshold=0.5) + _assert_binary_mask(qmask_thr) From a0e7839d959bf292611d01e328bae620bb50b0e9 Mon Sep 17 00:00:00 2001 From: thc1006 <84045975+thc1006@users.noreply.github.com> Date: Mon, 20 Jul 2026 02:20:45 +0800 Subject: [PATCH 2/2] Harden quality_mask: fail-closed non-finite handling, close datasets, stricter validation Follow-up to the V002 metadata-driven quality_mask fix, addressing review feedback: - Non-finite (fill / no-data) pixels are now fail-closed: any NaN/+Inf/-Inf is excluded (mask = 1), never returned as 0 (clear), in both the flag path and the threshold path. The mask variable declares _FillValue = -9999, which xarray decodes to NaN, so this is a real product concern. - Open the mask dataset and sensor-band group in a `with` block so both file handles are always closed, including on error (the original code leaked them). - Reject boolean indices explicitly (bool is a subclass of int) and validate that the mask band axis length matches the number of mask_bands names. - Use named indexing (isel) and incremental OR; the result is uint8 {0, 1} by construction, so the redundant final assert was removed. - Docstring: add a V002 note (11 bands, reordered indices, print mask_bands rather than hard-coding); add the same note to the quality how-to notebook. Tests: switch the synthetic fixtures to the xarray/h5netcdf backend (no netCDF4 test dependency), add a requirements-test.txt, and add cases for non-finite (NaN, +/-Inf, and a declared _FillValue) fail-closed behavior, boolean and axis-mismatch rejection, and dataset closure on success and error. The conftest now stubs only the specific missing optional dependency instead of masking any import error. Co-authored-by: junnncct1106 --- CHANGE_LOG.md | 4 +- .../How_to_use_EMIT_Quality_data.ipynb | 11 + python/modules/emit_tools.py | 187 ++++++++------- python/modules/tests/conftest.py | 12 +- python/modules/tests/requirements-test.txt | 9 + python/modules/tests/test_quality_mask.py | 214 ++++++++++++++++-- 6 files changed, 332 insertions(+), 105 deletions(-) create mode 100644 python/modules/tests/requirements-test.txt diff --git a/CHANGE_LOG.md b/CHANGE_LOG.md index 62aa9f1..28740ad 100644 --- a/CHANGE_LOG.md +++ b/CHANGE_LOG.md @@ -10,10 +10,12 @@ _________________________________________________________________________ > ### Changed > > - Made `quality_mask` in `python/modules/emit_tools.py` metadata-driven so it is safe for the EMIT L2A Mask V002 product. It now classifies each requested layer from its `sensor_band_parameters/mask_bands` name and refuses to build a mask from a continuous data layer (`AOD550`, `H2O (g cm-2)`, `SpecTf-Cloud Probability`, `SpecTf-Buffer Distance`). Previously only bands 5 and 6 were hard-coded as data bands, so the new V002 continuous bands were silently treated as flags. A continuous probability layer can be turned into a mask with the new optional `threshold` argument. +> - Hardened `quality_mask`: non-finite (fill / no-data) pixels are now excluded from the mask (fail-closed) rather than silently returned as clear; the opened datasets are always closed; boolean indices, out-of-range indices, and a band-axis / `mask_bands` length mismatch are rejected with clear errors. > ### Added > -> - `python/modules/tests/` with pytest regression tests for `quality_mask` covering the V001 (8-band) and V002 (11-band) mask layouts. +> - `python/modules/tests/` with network-free pytest regression tests for `quality_mask` covering the V001 (8-band) and V002 (11-band) mask layouts, non-finite handling, and resource cleanup, plus a `requirements-test.txt`. +> - A note in `python/how-tos/How_to_use_EMIT_Quality_data.ipynb` describing the V002 mask layout. ## 2024-06-17 diff --git a/python/how-tos/How_to_use_EMIT_Quality_data.ipynb b/python/how-tos/How_to_use_EMIT_Quality_data.ipynb index 1c0b7d1..49839fb 100644 --- a/python/how-tos/How_to_use_EMIT_Quality_data.ipynb +++ b/python/how-tos/How_to_use_EMIT_Quality_data.ipynb @@ -196,6 +196,17 @@ "Some of these bands are direct masks (Cloud, Dilated, Currus, Water, Spacecraft), and some (AOD550 and H2O (g cm-2)) are information calculated during the L2A reflectance retrieval that may be used as additional screening, depending on the application. The final mask that the EMIT mission will use for its minerological applications is shown as the Aggreged Flag - but not all users might want this particular mask." ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Note on the L2A Mask V002 product\n", + "\n", + "The reprocessed **V002** mask has **11 bands** (7 binary flags + 4 continuous data layers), and the flag ordering differs from V001, so band indices are **not** interchangeable between versions. Always read `mask_bands` for the file you are working with (as shown above) instead of hard-coding indices.\n", + "\n", + "The four continuous data layers — `AOD550`, `H2O (g cm-2)`, `SpecTf-Cloud Probability`, and `SpecTf-Buffer Distance` — are physical quantities, not flags. `quality_mask` rejects them so they cannot be silently mistaken for a binary mask. To turn the cloud-probability layer into a mask, pass a probability cutoff, e.g. `quality_mask(fp_mask, [], threshold=0.5)`." + ] + }, { "cell_type": "markdown", "metadata": {}, diff --git a/python/modules/emit_tools.py b/python/modules/emit_tools.py index df7d639..8ece8eb 100644 --- a/python/modules/emit_tools.py +++ b/python/modules/emit_tools.py @@ -289,102 +289,131 @@ def quality_mask(filepath, quality_bands, threshold=None): new V002 continuous bands (`SpecTf-Cloud Probability`, `SpecTf-Buffer Distance`). + V002 note: the L2A Mask V002 product has 11 bands (7 flags + 4 data layers) + and the flag ordering differs from V001, so band indices are not portable + between versions. Inspect `sensor_band_parameters/mask_bands` for the file at + hand (e.g. `print(xr.open_dataset(fp, group="sensor_band_parameters").mask_bands.values)`) + rather than hard-coding indices. `SpecTf-Cloud Probability` requires the + `threshold` argument; `SpecTf-Buffer Distance`, `AOD550` and `H2O` are data + layers and are rejected. + + Non-finite handling: the `mask` variable declares a `_FillValue` (-9999) that + xarray decodes to NaN, so a layer can contain non-finite (no-data) pixels. + Because a no-data pixel is not a clean observation, the mask is built + fail-closed: any non-finite value is excluded (mask = 1), never returned as 0 + (clear). + Parameters: filepath: an EMIT L2A Mask netCDF file. - quality_bands: an index or list of band indices (quality flags only) from the mask file that should be used in creation of the mask. Indices refer to the order of `sensor_band_parameters/mask_bands`. + quality_bands: an integer index or list of band indices (quality flags only) from the mask file that should be used in creation of the mask. Indices refer to the order of `sensor_band_parameters/mask_bands`. threshold: optional float in [0, 1]. When provided, a continuous *probability* layer (e.g. `SpecTf-Cloud Probability`) included in `quality_bands` is converted to a binary mask using `probability >= threshold`. Non-probability data layers (AOD, water vapor, buffer distance) are never accepted. Returns: qmask: a uint8 numpy array of {0, 1} that can be used with the emit_xarray function to apply a quality mask. """ - # Open Dataset - mask_ds = xr.open_dataset(filepath, engine="h5netcdf") - # Open Sensor band Group, which names every layer in the mask file - mask_parameters_ds = xr.open_dataset( + # Open the mask dataset and the sensor band group (which names every layer) + # in a context manager so both file handles are always closed, including on + # error. ``qmask`` is a detached numpy array, so closing the files is safe. + with xr.open_dataset(filepath, engine="h5netcdf") as mask_ds, xr.open_dataset( filepath, engine="h5netcdf", group="sensor_band_parameters" - ) - band_names = [ - str(b) for b in np.asarray(mask_parameters_ds["mask_bands"].data).ravel() - ] - n_bands = len(band_names) - - # Accept a single index or an iterable of indices - if isinstance(quality_bands, (int, np.integer)): - quality_bands = [quality_bands] - quality_bands = list(quality_bands) - - # Validate the requested indices against this file's actual bands - for b in quality_bands: - if not isinstance(b, (int, np.integer)) or b < 0 or b >= n_bands: + ) as mask_parameters_ds: + band_names = [ + str(b) for b in np.asarray(mask_parameters_ds["mask_bands"].data).ravel() + ] + n_bands = len(band_names) + + # Guard against a layout where the band axis and the band-name list + # disagree, instead of silently mis-slicing or raising a later IndexError. + mask_var = mask_ds["mask"] + if mask_var.ndim != 3 or mask_var.shape[-1] != n_bands: raise ValueError( - f"quality_bands index {b!r} is out of range for this mask file, " - f"which has {n_bands} bands: " - f"{[f'{i}: {name}' for i, name in enumerate(band_names)]}" + f"Unexpected mask layout: 'mask' has shape {tuple(mask_var.shape)} " + f"but 'sensor_band_parameters/mask_bands' lists {n_bands} bands; " + f"expected the last axis to index the mask bands " + f"(dims (downtrack, crosstrack, bands))." ) + band_dim = mask_var.dims[-1] - if threshold is not None and not (0.0 <= float(threshold) <= 1.0): - raise ValueError(f"threshold must be within [0, 1]; got {threshold!r}") - - # Print Flags used - flags_used = [band_names[b] for b in quality_bands] - print(f"Flags used: {flags_used}") - - def _is_flag(name): - # Binary quality-flag layers are named "... Flag" in both V001 and V002. - return name.strip().lower().endswith("flag") - - def _is_probability(name): - return "probability" in name.lower() + # Accept a single index or an iterable of indices + if isinstance(quality_bands, (int, np.integer, bool, np.bool_)): + quality_bands = [quality_bands] + quality_bands = list(quality_bands) - flag_band_hint = [ - f"{i}: {name}" for i, name in enumerate(band_names) if _is_flag(name) - ] - - # Build one binary layer per requested band - layers = [] - for b in quality_bands: - name = band_names[b] - layer = mask_ds["mask"][:, :, b].values - if _is_flag(name): - # Defense in depth: a flag layer must be binary. Never silently clip. - finite = layer[np.isfinite(layer)] - uniq = np.unique(finite) - if uniq.size and not np.all(np.isin(uniq, (0.0, 1.0))): + # Validate the requested indices. Reject bool explicitly: bool is a + # subclass of int and would otherwise be used as a boolean index. + for b in quality_bands: + if isinstance(b, (bool, np.bool_)) or not isinstance(b, (int, np.integer)): raise ValueError( - f"Band {b} ('{name}') is named as a flag but contains non-binary " - f"values (e.g. {uniq[:5]}); refusing to build a mask from it." + f"quality_bands must contain integer band indices; got {b!r} " + f"of type {type(b).__name__}." ) - layers.append((layer > 0).astype(np.uint8)) - elif _is_probability(name) and threshold is not None: - layers.append((layer >= float(threshold)).astype(np.uint8)) - else: - # Continuous / data band: refuse rather than silently clip it to {0, 1}. - hint = "" - if _is_probability(name): - hint = ( - " Pass threshold= to convert this probability " - "layer into a binary mask (probability >= threshold)." + if b < 0 or b >= n_bands: + raise ValueError( + f"quality_bands index {int(b)} is out of range for this mask " + f"file, which has {n_bands} bands: " + f"{[f'{i}: {name}' for i, name in enumerate(band_names)]}" ) - raise ValueError( - f"Band {b} ('{name}') is a continuous data layer, not a binary " - f"quality flag, so it cannot be combined into a mask.{hint} " - f"Flag bands available in this file: {flag_band_hint}." - ) - # Combine the binary layers. Logical OR reproduces the previous - # sum-then-clip behavior exactly for binary flags. - if layers: - qmask = np.zeros_like(layers[0]) - for layer in layers: - qmask |= layer - else: - qmask = np.zeros(mask_ds["mask"].shape[:2], dtype=np.uint8) - qmask[qmask > 1] = 1 + if threshold is not None and not (0.0 <= float(threshold) <= 1.0): + raise ValueError(f"threshold must be within [0, 1]; got {threshold!r}") + + # Print Flags used + flags_used = [band_names[b] for b in quality_bands] + print(f"Flags used: {flags_used}") + + def _is_flag(name): + # Binary quality-flag layers are named "... Flag" in both V001 and V002. + return name.strip().lower().endswith("flag") + + def _is_probability(name): + return "probability" in name.lower() + + flag_band_hint = [ + f"{i}: {name}" for i, name in enumerate(band_names) if _is_flag(name) + ] + + # Combine one binary layer per requested band. Non-finite (fill / no-data) + # values are always excluded (mask = 1), never returned as 0 (clear), in + # both the flag and threshold paths. + qmask = None + for b in quality_bands: + name = band_names[b] + layer = mask_var.isel({band_dim: b}).values + finite = np.isfinite(layer) + if _is_flag(name): + # Defense in depth: a flag layer's finite values must be binary. + finite_vals = np.unique(layer[finite]) + if finite_vals.size and not np.all(np.isin(finite_vals, (0.0, 1.0))): + raise ValueError( + f"Band {b} ('{name}') is named as a flag but contains " + f"non-binary values (e.g. {finite_vals[:5]}); refusing to " + f"build a mask from it." + ) + binary = np.where(finite, layer > 0, True).astype(np.uint8) + elif _is_probability(name) and threshold is not None: + binary = np.where( + finite, layer >= float(threshold), True + ).astype(np.uint8) + else: + # Continuous / data band: refuse rather than silently clip to {0, 1}. + hint = "" + if _is_probability(name): + hint = ( + " Pass threshold= to convert this probability " + "layer into a binary mask (probability >= threshold)." + ) + raise ValueError( + f"Band {b} ('{name}') is a continuous data layer, not a binary " + f"quality flag, so it cannot be combined into a mask.{hint} " + f"Flag bands available in this file: {flag_band_hint}." + ) + # Incremental OR keeps the result strictly {0, 1} by construction. + qmask = binary if qmask is None else (qmask | binary) + + if qmask is None: + # No bands requested: nothing is masked. + qmask = np.zeros(mask_var.shape[:2], dtype=np.uint8) - # The returned mask must be binary; fail loudly if that invariant is broken. - assert set(int(v) for v in np.unique(qmask)).issubset( - {0, 1} - ), "quality_mask produced a non-binary mask" return qmask diff --git a/python/modules/tests/conftest.py b/python/modules/tests/conftest.py index 30092b0..bf95a8d 100644 --- a/python/modules/tests/conftest.py +++ b/python/modules/tests/conftest.py @@ -38,5 +38,13 @@ for _name in _OPTIONAL_DEPENDENCIES: try: importlib.import_module(_name) - except Exception: # ImportError or any import-time failure - sys.modules[_name] = mock.MagicMock() + except ModuleNotFoundError as _exc: + # Only stub the specific optional dependency that is genuinely absent. + # A ModuleNotFoundError naming an unrelated module (e.g. a broken + # transitive dependency or an ABI/version mismatch) is re-raised so it is + # not silently masked as a passing test run. + _missing = _exc.name or "" + if _missing == _name or _name.startswith(_missing + "."): + sys.modules[_name] = mock.MagicMock() + else: + raise diff --git a/python/modules/tests/requirements-test.txt b/python/modules/tests/requirements-test.txt new file mode 100644 index 0000000..878bed1 --- /dev/null +++ b/python/modules/tests/requirements-test.txt @@ -0,0 +1,9 @@ +# Dependencies for running the emit_tools unit tests (python/modules/tests). +# numpy, xarray and h5netcdf are already required by emit_tools itself; pytest is +# the only additional dependency. The tests are network-free (they synthesize +# small mask files) and do not require the geospatial stack (GDAL, geopandas, +# rasterio, ...) that emit_tools imports for its other functions. +pytest +numpy +xarray +h5netcdf diff --git a/python/modules/tests/test_quality_mask.py b/python/modules/tests/test_quality_mask.py index 624ea36..e430566 100644 --- a/python/modules/tests/test_quality_mask.py +++ b/python/modules/tests/test_quality_mask.py @@ -29,8 +29,8 @@ import os import numpy as np -import netCDF4 as nc import pytest +import xarray as xr from emit_tools import quality_mask @@ -67,31 +67,57 @@ ] -def _write_mask_nc(path, band_names, mask_array): - """Write a minimal EMIT-like L2A Mask netCDF file. +def _write_mask_nc(path, band_names, mask_array, fill_value=None): + """Write a minimal EMIT-like L2A Mask netCDF file with the h5netcdf backend. Structure mirrors the real product: a root ``mask`` variable with dims ``(downtrack, crosstrack, bands)`` and a ``sensor_band_parameters`` group - holding the ``mask_bands`` string variable. + holding the ``mask_bands`` string variable. Uses only xarray + h5netcdf (both + already required by emit_tools), so the tests add no netCDF backend + dependency. When ``fill_value`` is given, the ``mask`` variable declares that + ``_FillValue`` (the real V002 product declares ``_FillValue = -9999``, which + xarray decodes to NaN on read). """ - downtrack, crosstrack, n_bands = mask_array.shape - assert n_bands == len(band_names) - with nc.Dataset(path, "w", format="NETCDF4") as ds: - ds.createDimension("downtrack", downtrack) - ds.createDimension("crosstrack", crosstrack) - ds.createDimension("bands", n_bands) - mask_var = ds.createVariable( - "mask", "f4", ("downtrack", "crosstrack", "bands") - ) - mask_var[:] = mask_array.astype("f4") - mask_var.long_name = "Masks" - mask_var.units = "unitless" - grp = ds.createGroup("sensor_band_parameters") - band_var = grp.createVariable("mask_bands", str, ("bands",)) - for i, name in enumerate(band_names): - band_var[i] = name - band_var.long_name = "Mask Band Names" - return str(path) + path = str(path) + root = xr.Dataset( + {"mask": (("downtrack", "crosstrack", "bands"), mask_array.astype("float32"))} + ) + root["mask"].attrs = {"long_name": "Masks", "units": "unitless"} + encoding = None if fill_value is None else {"mask": {"_FillValue": fill_value}} + root.to_netcdf(path, engine="h5netcdf", encoding=encoding) + sbp = xr.Dataset({"mask_bands": (("bands",), np.array(band_names))}) + sbp["mask_bands"].attrs = {"long_name": "Mask Band Names"} + sbp.to_netcdf(path, engine="h5netcdf", group="sensor_band_parameters", mode="a") + return path + + +class _DatasetSpy: + """Transparent proxy around an xarray Dataset that records ``close()``. + + Used to prove ``quality_mask`` closes the file handles it opens (including + on error) without depending on private xarray internals. + """ + + def __init__(self, ds): + self._ds = ds + self.closed = False + + def __getattr__(self, item): + return getattr(self._ds, item) + + def __getitem__(self, item): + return self._ds[item] + + def __enter__(self): + return self + + def __exit__(self, *exc): + self.close() + return False + + def close(self): + self.closed = True + return self._ds.close() def _binary(pattern): @@ -263,6 +289,130 @@ def test_single_int_accepted(v002_mask): np.testing.assert_array_equal(qmask, _expected_or(v002_mask["array"], [0])) +@pytest.mark.parametrize("bad", [True, False, np.bool_(True), np.bool_(False)]) +def test_bool_index_rejected(v002_mask, bad): + # bool is a subclass of int; it must not be accepted as a band index. + with pytest.raises(ValueError): + quality_mask(v002_mask["path"], [bad]) + + +def test_bool_scalar_index_rejected(v002_mask): + with pytest.raises(ValueError): + quality_mask(v002_mask["path"], True) + + +def test_axis_metadata_mismatch_raises(tmp_path): + # 'mask' has 11 bands but 'mask_bands' lists only 8 -> inconsistent layout. + path = str(tmp_path / "mismatch.nc") + root = xr.Dataset( + {"mask": (("downtrack", "crosstrack", "bands"), np.zeros((3, 4, 11), "float32"))} + ) + root.to_netcdf(path, engine="h5netcdf") + sbp = xr.Dataset({"mask_bands": (("mask_band",), np.array(V001_MASK_BANDS))}) + sbp.to_netcdf(path, engine="h5netcdf", group="sensor_band_parameters", mode="a") + with pytest.raises(ValueError): + quality_mask(path, [0]) + + +# --- Non-finite (fill / no-data) handling -------------------------------------- + +def test_flag_nan_pixel_is_masked_not_clear(tmp_path): + # A no-data (NaN) pixel in a flag layer must be EXCLUDED (1), never returned + # as a clean observation (0). _A[2, 3] is 0 (clear) before we blank it. + cloud = _A.copy() + cloud[2, 3] = np.nan + layers = [cloud, _B, _ZERO, _ZERO, _C, _AOD, _H2O, _D, _PROB, _E, _DIST] + arr = np.stack(layers, axis=-1) + path = _write_mask_nc(tmp_path / "v002_nan.nc", V002_MASK_BANDS, arr) + + qmask = quality_mask(path, [0]) + _assert_binary_mask(qmask) + assert qmask[2, 3] == 1 + # The previous `(layer > 0)` logic would have returned 0 (clear) here. + assert bool((np.nan_to_num(arr[:, :, 0]) > 0)[2, 3]) is False + + +def test_flag_fillvalue_pixel_is_masked(tmp_path): + # Same, but via a declared _FillValue (-9999) as in the real product, which + # xarray decodes to NaN on read. _A[1, 0] is 0 (clear) before we blank it. + cloud = _A.copy() + cloud[1, 0] = -9999.0 + layers = [cloud, _B, _ZERO, _ZERO, _C, _AOD, _H2O, _D, _PROB, _E, _DIST] + arr = np.stack(layers, axis=-1) + path = _write_mask_nc( + tmp_path / "v002_fill.nc", V002_MASK_BANDS, arr, fill_value=-9999.0 + ) + + qmask = quality_mask(path, [0]) + _assert_binary_mask(qmask) + assert qmask[1, 0] == 1 + + +def test_flag_infinite_pixels_are_masked(tmp_path): + # Both +Inf and -Inf are non-finite and must be excluded (1). Naive + # `layer > 0` would return -Inf as 0 (clear) -- the fail-closed bug. + # _A[0, 1] and _A[2, 0] are both 0 (clear) before we blank them. + cloud = _A.copy() + cloud[0, 1] = np.inf + cloud[2, 0] = -np.inf + layers = [cloud, _B, _ZERO, _ZERO, _C, _AOD, _H2O, _D, _PROB, _E, _DIST] + arr = np.stack(layers, axis=-1) + path = _write_mask_nc(tmp_path / "v002_inf.nc", V002_MASK_BANDS, arr) + + qmask = quality_mask(path, [0]) + _assert_binary_mask(qmask) + assert qmask[0, 1] == 1 # +inf excluded + assert qmask[2, 0] == 1 # -inf excluded + + +def test_threshold_nan_pixel_is_masked(tmp_path): + # A NaN in the probability layer must also be excluded, not clear. + # _PROB[0, 0] is 0.30, which is < 0.5 (clear) before we blank it. + prob = _PROB.copy() + prob[0, 0] = np.nan + layers = [_A, _B, _ZERO, _ZERO, _C, _AOD, _H2O, _D, prob, _E, _DIST] + arr = np.stack(layers, axis=-1) + path = _write_mask_nc(tmp_path / "v002_prob_nan.nc", V002_MASK_BANDS, arr) + + qmask = quality_mask(path, [8], threshold=0.5) + _assert_binary_mask(qmask) + assert qmask[0, 0] == 1 + # `(layer >= 0.5)` alone would have returned 0 (clear) here. + assert bool((np.nan_to_num(arr[:, :, 8]) >= 0.5)[0, 0]) is False + + +# --- Resource management: opened datasets are closed --------------------------- + +def _install_open_spy(monkeypatch): + import emit_tools + + real_open = emit_tools.xr.open_dataset + spies = [] + + def spy_open(*args, **kwargs): + spy = _DatasetSpy(real_open(*args, **kwargs)) + spies.append(spy) + return spy + + monkeypatch.setattr(emit_tools.xr, "open_dataset", spy_open) + return spies + + +def test_datasets_closed_on_success(v002_mask, monkeypatch): + spies = _install_open_spy(monkeypatch) + quality_mask(v002_mask["path"], [0, 1]) + assert len(spies) == 2 + assert all(s.closed for s in spies) + + +def test_datasets_closed_on_error(v002_mask, monkeypatch): + spies = _install_open_spy(monkeypatch) + with pytest.raises(ValueError): + quality_mask(v002_mask["path"], [8]) # continuous band -> raises inside `with` + assert len(spies) == 2 + assert all(s.closed for s in spies) + + # --- Optional: verify against a real downloaded granule ------------------------ @pytest.mark.skipif( @@ -271,11 +421,29 @@ def test_single_int_accepted(v002_mask): ) def test_real_v002_granule(): fp = os.environ["EMIT_L2A_MASK_V002"] + + # The verified V002 SpecTf layers sit at indices 8, 9, 10. + with xr.open_dataset(fp, engine="h5netcdf", group="sensor_band_parameters") as sbp: + names = [str(x) for x in np.asarray(sbp["mask_bands"].data).ravel()] + assert len(names) == 11 + assert names[8] == "SpecTf-Cloud Probability" + assert names[9] == "SpecTf-Cloud Flag" + assert names[10] == "SpecTf-Buffer Distance" + # Binary flags combine into a valid mask. qmask = quality_mask(fp, [0, 1, 4]) _assert_binary_mask(qmask) + # SpecTf-Cloud Probability (index 8) is rejected unless a threshold is given. with pytest.raises(ValueError): quality_mask(fp, [8]) - qmask_thr = quality_mask(fp, [8], threshold=0.5) + + # With a threshold, the result equals a directly-computed fail-closed + # reference over the decoded probability layer. + t = 0.5 + with xr.open_dataset(fp, engine="h5netcdf") as ds: + raw8 = np.asarray(ds["mask"].isel(bands=8).values) + reference = np.where(np.isfinite(raw8), raw8 >= t, True).astype(np.uint8) + qmask_thr = quality_mask(fp, [8], threshold=t) _assert_binary_mask(qmask_thr) + np.testing.assert_array_equal(qmask_thr, reference)