Make quality_mask safe for the EMIT L2A Mask V002 product

[thc1006]

Summary

quality_mask in python/modules/emit_tools.py silently produces a wrong mask for the EMIT L2A Mask V002 product. It hard-codes bands 5 and 6 (the V001 AOD550 / H2O data bands) as the only non-flag layers:

if any(x in quality_bands for x in [5, 6]):
    raise AttributeError("Selected flags include a data band (5 or 6) not just flag bands")
else:
    qmask = np.sum(mask_ds["mask"][:, :, quality_bands].values, axis=-1)
    qmask[qmask > 1] = 1

The V002 mask adds continuous layers and reorders the flags, so any V002 continuous band other than 5/6 falls through to the else branch and is treated as a binary flag: its values are summed and clipped >1 -> 1. The result is not a {0, 1} mask, and downstream qmask == 1 masking then behaves incorrectly, with no error raised.

Verified V002 band layout

Confirmed by reading sensor_band_parameters/mask_bands from a real granule (EMIT_L2A_MASK_002_20220810T034103_2222203_001.nc, EMITL2AMASK v002). V002 has 11 bands: 7 binary flags and 4 continuous data layers.

index	band name	type
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 flags were reordered relative to V001 (e.g. Dilated Cloud Flag moved to index 4), so band indices are not interchangeable between versions and index-based logic is fragile.

Reproduction (silent failure)

import sys; sys.path.append("python/modules")
from emit_tools import quality_mask

# band 8 = "SpecTf-Cloud Probability", a continuous probability layer
mask = quality_mask("EMIT_L2A_MASK_002_20220810T034103_2222203_001.nc", [8])

On the real granule the current code returns a float32 array with 31 distinct values in the range [0.2637, 1.0] (30 of them fractional) instead of a {0, 1} mask. It never raises. SpecTf-Buffer Distance (index 10) is similarly accepted; it happens to be 0/1 in this scene but is a pixel-distance layer that will hold values >1 in cloudier scenes and be silently clipped.

Fix

There is no explicit "is-a-flag" metadata field in the mask file (the sensor_band_parameters group has no such attribute, and the single mask variable shares one units), so the band name is the only per-layer discriminator. quality_mask is now metadata-driven and hardened:

• reads sensor_band_parameters/mask_bands and classifies each requested layer by name;
• treats a layer as a binary flag only when its name ends with Flag (case-insensitive, matching both V001 and V002);
• rejects any continuous data layer (AOD550, H2O (g cm-2), SpecTf-Cloud Probability, SpecTf-Buffer Distance) with a ValueError that names the band and lists the file's actual flag bands;
• adds an optional threshold argument so a probability layer can be turned into a mask explicitly (probability >= threshold);
• as defense in depth, refuses a Flag-named layer whose finite values are not binary rather than clipping it.

For binary flags the combination is a logical OR, which reproduces the previous sum-then-clip result exactly.

Non-finite (fill / no-data) handling — fail-closed

The mask variable declares _FillValue = -9999, which xarray's default mask_and_scale decodes to NaN, so a layer can contain non-finite (no-data) pixels. The real granule above contains no non-finite pixels, but fill/no-data is a documented part of the product schema, so the code must handle it. A no-data pixel is not a clean observation, so the mask is now built fail-closed: any non-finite value is excluded (mask = 1), never returned as 0 (clear). This applies to both the flag path and the threshold (probability) path.

Concretely, for a flag layer with pixel values [0, 1, NaN, +Inf, -Inf], quality_mask now returns [0, 1, 1, 1, 1], whereas the naive layer > 0 returned [0, 1, 0, 1, 0] (NaN and -Inf silently treated as clear). The threshold path behaves the same way.

Resource management, and additional validation

• The two datasets that quality_mask opens are now opened in a with block, so both file handles are always closed, including when the function raises (the previous code, and upstream, leaked them). The returned array is a detached numpy array, so closing is safe.
• Boolean indices are rejected explicitly (bool is a subclass of int, and np.bool_ would otherwise be used as a boolean index deep inside xarray).
• Out-of-range band indices are rejected with a message listing the file's bands.
• A mismatch between the mask's band axis length and the number of mask_bands names raises a clear ValueError instead of a later IndexError or a mis-slice.

Verification on the real granule

Running the previous function and the new one side by side on EMIT_L2A_MASK_002_20220810T034103_2222203_001.nc:

• Binary-flag path is value-identical old vs new (e.g. bands [0, 1, 4] -> same {0, 1} mask, 1588851 pixels masked in both).
• Band 8 (SpecTf-Cloud Probability): old returns the non-binary float array described above; new raises ValueError: Band 8 ('SpecTf-Cloud Probability') is a continuous data layer, not a binary quality flag ....
• Band 10 (SpecTf-Buffer Distance): new raises naming the band.
• quality_mask(fp, [8], threshold=0.5) returns a uint8 {0, 1} mask equal to a directly-computed numpy.where(isfinite, raw >= 0.5, True) reference.
• Bands 5/6 remain rejected (now ValueError naming AOD550 / H2O (g cm-2)).

Backward compatibility

• Signature is quality_mask(filepath, quality_bands, threshold=None); existing two-argument calls are unchanged, and legacy V001 flag combinations return the same mask as before.
• The rejection of a data band now raises ValueError (was AttributeError) with a clearer, band-naming message.
• The returned array dtype is uint8 {0, 1} (was float32 {0.0, 1.0}); finite values are unchanged and qmask == 1 masking is unaffected.

Tests

Adds python/modules/tests/ (network-free pytest) with small synthetic V001 (8-band) and V002 (11-band) mask files whose mask_bands metadata mirrors the real products. Coverage: legacy V001 flags still combine correctly, V002 binary flags combine correctly, continuous V002 bands now raise instead of silently returning a wrong mask, the threshold path, non-finite handling (NaN, +Inf, -Inf, and a declared _FillValue) is fail-closed, bool/out-of-range/axis-mismatch indices raise, and the opened datasets are closed on both the success and error paths. An optional test runs against a real granule when EMIT_L2A_MASK_V002 points to one and checks the threshold result against a directly-computed reference.

The fixtures use only xarray + h5netcdf (already required by emit_tools), so the tests add no netCDF backend dependency; python/modules/tests/requirements-test.txt lists the test dependencies (pytest is the only addition). The tests stub emit_tools' optional heavy geospatial imports (GDAL, geopandas, rasterio, ...) when they are absent so the suite runs in a minimal environment, and stub only the specific missing module so a real import/ABI failure is not masked.

A short note describing the V002 mask layout was added to python/how-tos/How_to_use_EMIT_Quality_data.ipynb.

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.