feat: add script to apply ESA quality mask

[ceholden]

What I'm changing

This PR is intended to address this ticket, NASA-IMPACT/hls-sentinel#155

Specifically this is intended to replace a utility written in C to apply a binary quality mask to Sentinel-2 L1C data that masks lost or degraded MSI pixel values. The original PR was, NASA-IMPACT/hls-sentinel#152

This utility script is intended to fit into the rest of the HLS Sentinel-2 related processing. To apply the script we would need to do a few extra steps,

1. Merge this PR and create a new release in Github to point to (e.g., hls-utilities@v1.10)
2. Update the installation of hls-utilities in the hls-sentinel container to include this update (https://github.com/NASA-IMPACT/hls-sentinel/blob/20180e0d0cd286ffefe61bd6aeeaeff056874102/Dockerfile#L123)
3. Update the hls-sentinel processing steps to include running this script as part of sentinel_granule.sh (e.g., around here https://github.com/NASA-IMPACT/hls-sentinel/blob/20180e0d0cd286ffefe61bd6aeeaeff056874102/scripts/sentinel_granule.sh#L44)
4. Test
5. Update the hls-sentinel image used in our production pipelines by tagging our new container appropriately (i.e., "latest")

How I did it

The ESA quality mask is an 8 band single bit image that encodes 0/1 for 8 quality related attributes. The bands we're interested in are bands 3 & 4 (counting from 1 in GDAL notation) which encode if the pixel retrieval was either lost or degraded.

I followed the original C code implementation with a few modifications...

1. Finds the mask + image pairs
   • The original C code will exit 1 if not all of the images have been found, or would exit 0 if the granule doesn't have the expected image/mask pair (e.g., prior to baseline 04.00).
   • I wanted to ensure we don't hard fail if we can't apply the mask, and favored doing a "no-op" if the granule file structure isn't as we expect
2. Reads the 3rd and 4th bands of the mask
   • In the original implementation the JPEG2000 is first converted to an ENVI file, a binary in "band sequential" (BSQ) format. The C code reads all bands for the entire image and fills the imagery band if EITHER mask band 3 or 4 have a value of 1 for a pixel row/column.
   • In our implementation we can use GDAL to only read the 3rd and 4th bands
3. Updates the imagery in-place by assigning HLS_REFL_FILLVAL (-9999) to any pixel where the mask was TRUE for either band.
   • The -9999 value is not suitable for us because the original L1C images in JPEG2000 format use 0 as the nodata value (since -9999 is out of data type range)
   • The 0 value is used as a no data value for L1C imagery and this is noted in the metadata
   • I see in the hls-sentinel that some utilities modify in place while others create new data. If it'd be useful for a "debug mode" it would be easy to rewrite the images to new files

How you can test it

Install the test dependencies and run unit tests for this script,

pip install -e .[test]
pytest tests/test_apply_s2_quality_mask.py

The unit tests are what helped me discover that JPEG2000 driver will apply lossy compression when we update the file (r+ access mode) even though the original file had lossless compression in creation options 🙃

I ran this for one of the granules mentioned in the original C code PR and recorded a tiny demo,

Note the "inspect tool" that shows the imagery values on the right hand side. The S2 L1C imagery don't have a No Data Value set, so QGIS shows the 0 as valid pixels. Still we can see the individual pixel values are correct (0 instead of ~1,000) and the image is stretched differently because of the different minimum values.

[junchangju]

Checking the quality mask for every mask could be very slow, but there is no need to do that. The quality mask is to be applied only if the metadata indicates data loss. For example, the following xml file
S2A_MSIL1C_20240830T092031_N0511_R093_T35UMQ_20240830T112048.SAFE/GRANULE/L1C_T35UMQ_A047995_20240830T092120/QI_DATA/GENERAL_QUALITY.xml
indicates data loss by having a line:
There is data loss in this tile

So we can grep "There is data loss in this tile" on GENERAL_QUALITY.xml.
This xml is easy to find because the GRANULE directory contains only one subdirectory, which is can be referred to as "*" in bash. This odd directory structure of GRANULE is a vestige of the original SAFE format which did include dozens of granules in one SAFE.

PS: I don't know why the xml markup tags were removed once I posted the comment. But the message remains.

[ceholden]

Thanks for the great suggestion and background information @junchangju! It sounds like it will save a lot of time so I'll add that in as a pre-check.

I found this metadata in the examples I've been using and saw that it describes the impacted bands as "Affected_Bands". I suspect this will probably be all bands, but it's easy enough to get the list,

...
<Earth_Explorer_File>
    <Data_Block>
        <report>
            <checkList>
                ...
                <check>
                    <inspection creation="2022-11-21T07:37:42.488Z" duration="11940" execution="2022-11-21T07:37:42.495Z" id="Data_Loss" item="S2A_OPER_MSI_L1C_TL_ATOS_20221121T053331_A038726_T45TXF_N04.00" itemURL="/mnt/nfs-l1-02/l1-processing/processing-id-17769-2022-11-21-06-03-12/TaskTable_20221121T060312/FORMAT_METADATA_TILE_L1C_ADAPT/output/PDI_DS_TILE_LIST/S2A_OPER_MSI_L1C_TL_ATOS_20221121T053331_A038726_T45TXF_N04.00/" name="Check TECQUA for data loss " priority="5" processingStatus="done" status="FAILED"/>
                    <message contentType="text/plain">There is data loss in this tile</message>
                    <extraValues>
                        <value name="Affected_Bands">B01 B02 B03 B04 B05 B06 B07 B08 B09 B10 B11 B12 B8A</value>
                    </extraValues>
                </check>
            </checkList>
        </report>
    </Data_Block>
</Earth_Explorer_File>

[ceholden]

Ah actually it's not always all of the bands. I had this in the 3rd example I checked,

                <check>
                    <inspection creation="2023-03-05T08:36:45.870Z" duration="23651" execution="2023-03-05T08:36:45.877Z" id="Data_Loss" item="S2B_OPER_MSI_L1C_TL_2BPS_20230305T072429_A031305_T43QDG_N05.09" itemURL="/mnt/local/l1ctile/work/3be743a1-28c3-4010-84c9-8ccca1586310/ipf_output_L1CTile_20230305082123/L1CTile/TaskTable_20230305T082127/FORMAT_METADATA_TILE_L1C/output/PDI_DS_TILE_LIST/S2B_OPER_MSI_L1C_TL_2BPS_20230305T072429_A031305_T43QDG_N05.09/" name="Check TECQUA for data loss " priority="5" processingStatus="done" status="FAILED"/>
                    <message contentType="text/plain">There is data loss in this tile</message>
                    <extraValues>
                        <value name="Affected_Bands">B01 B02 B03 B09 B10 B11 B12 B8A</value>
                    </extraValues>
                </check>

[junchangju]

Ah actually it's not always all of the bands. I had this in the 3rd example I checked,

                <check>
                    <inspection creation="2023-03-05T08:36:45.870Z" duration="23651" execution="2023-03-05T08:36:45.877Z" id="Data_Loss" item="S2B_OPER_MSI_L1C_TL_2BPS_20230305T072429_A031305_T43QDG_N05.09" itemURL="/mnt/local/l1ctile/work/3be743a1-28c3-4010-84c9-8ccca1586310/ipf_output_L1CTile_20230305082123/L1CTile/TaskTable_20230305T082127/FORMAT_METADATA_TILE_L1C/output/PDI_DS_TILE_LIST/S2B_OPER_MSI_L1C_TL_2BPS_20230305T072429_A031305_T43QDG_N05.09/" name="Check TECQUA for data loss " priority="5" processingStatus="done" status="FAILED"/>
                    <message contentType="text/plain">There is data loss in this tile</message>
                    <extraValues>
                        <value name="Affected_Bands">B01 B02 B03 B09 B10 B11 B12 B8A</value>
                    </extraValues>
                </check>

Then loop over all bands should be fine.

[ceholden]

Closing this PR because it isn't running tests in CI. I now have access to this repo and can submit a PR from a branch. See #15 for followup