[ceholden]

Description

Closes https://github.com/NASA-IMPACT/hls_development/issues/441

The Fmask code from GERSL has no library versions included in the installation details, but anecodtally it seems like it was developed on very old versions, including numpy<2. NumPy v2 was released in June, 2024, so anyone following their conda-based install details in recent years would get a version of NumPy v2.

This major version difference is problematic because of a change in NumPy v2 regarding type promotion rules following NEP-50. The code does a comparison of band data in float32 against a (Python-dtype) float, so in NumPy v2 a type promotion occurs and saturation comparison checks fail.

To give a brief example:

cube_value = np.float32(6.4535)  # stored in the float32 cube, rounded
scalar = np.float64(6.4535)      # the comparison value, not rounded

cube_value == scalar
# NumPy 1: scalar is cast down to float32, so True
# NumPy 2: cube_value is widened to float64, keeps its rounding error, so False

This doesn't happen in the Landsat code path because the saturation check was done in units of DNs (integers, before scaling) rather than in (float) units of TOA reflectance.

What I did

I fixed this by mimicking the Landsat code path:

1. Add saturation: bool = False to the Sentinel2.read_band method and when True return a boolean array of the saturation check done in units of DNs
2. Run the saturation check for the 3 bands used for Sentinel-2 checks during read_band calls

I added unit tests for this use case and set up the repo to run the unit tests in CI for PRs and merges to main.

---

Once this is merged I'll create an issue upstream and perhaps include a PR if the merge isn't too gnarly. Alas the GERSL repository has switched off of the main branch and created an entirely new Git history in their develop branch, making synchronization more difficult than it needs to be.