Merge pull request #88 from lukasValentin/dev

Generation of binary masks from Landsat and Sentinel-2 scenes

.coveragerc

@@ -1,6 +1,6 @@
 [run]
 branch = True
-omit = */tests/*
+omit = */tests/*, eodal/downloader/, eodal/metadata/database/
 
 [report]
 exclude_lines =

eodal/core/sensors/landsat.py

@@ -328,12 +328,15 @@ def _preprocess_band_selection(
                     k for k, v in sensor_bands.items() if v == 'red'][0]
                 band_fpath = in_dir.glob(f"*_{sensor_specific_band_name}.TIF")[0]
             band_res = None
-            if band_name in sensor_bands.values():
-                band_res = band_resolution[sensor][band_name]
-            elif read_qa or band_name in qa_bands:
-                band_res = band_resolution['quality_flags'][band_name]
-            elif read_atcor or band_name in atcor_bands:
-                band_res = band_resolution['atmospheric_correction'][band_name]
+            try:
+                if band_name in sensor_bands.values():
+                    band_res = band_resolution[sensor][band_name]
+                elif read_qa or band_name in qa_bands:
+                    band_res = band_resolution['quality_flags'][band_name]
+                elif read_atcor or band_name in atcor_bands:
+                    band_res = band_resolution['atmospheric_correction'][band_name]
+            except KeyError:
+                continue
 
             item = {
                 'band_name': band_name,
@@ -553,14 +556,77 @@ def calc_si(
             self, si_name=si_name, inplace=inplace,
             band_mapping=band_mapping)
 
+    def get_cloud_and_shadow_mask(
+            self,
+            mask_band: Optional[str] = 'qa_pixel',
+            inplace: Optional[bool] = False,
+            name_cloud_mask: Optional[str] = 'cloud_mask',
+            cloud_classes: list[int] = [1, 2, 3, 5]
+    ) -> Band | None:
+        """
+        Get a cloud and shadow mask from the pixel qualuty band (using bits
+        [1, 2, 3, 5] by default)
+
+        NOTE:
+            Since the `mask_band` can be set to *any* `Band` it is also
+            possible to use a different cloud mask, e.g., from
+            a custom classifier as long as the bit map used for the classes
+            is the same.
+            See also
+            https://www.usgs.gov/media/images/landsat-collection-2-pixel-quality-assessment-bit-index
+            for more details.
+
+        :param mask_band:
+            The name of the band to use for masking. Default is 'qa_pixel'.
+        :param inplace:
+            Whether to return a new `Band` instance or add a `Band` to the
+            current `Landsat` instance. Default is False.
+        :param name_cloud_mask:
+            The name of the cloud mask band. Default is 'cloud_mask'.
+        :param cloud_classes:
+            Classes that define "clouds". Default is [1, 2, 3, 5].
+        :return:
+            A `Band` instance containing the cloud mask or `None` if
+            `inplace` is True.
+        """
+        mask_list = []
+        for cloud_class in cloud_classes:
+            # construct the bit range
+            bit_range = (cloud_class, cloud_class)
+            # get the binary mask as boolean array
+            mask = self.mask_from_qa_bits(
+                bit_range=bit_range,
+                band_name=mask_band).astype('bool')
+            # add the mask to the list
+            mask_list.append(mask)
+
+        # combine the masks into a single cloud mask
+        cloud_mask = np.any(mask_list, axis=0)
+        # update cloud mask with the mask of the area of interest,
+        if self[mask_band].is_masked_array:
+            cloud_mask = np.logical_and(
+                cloud_mask, ~self[mask_band].values.mask)
+        cloud_mask_band = Band(
+            band_name=name_cloud_mask,
+            values=cloud_mask,
+            geo_info=self[mask_band].geo_info,
+            nodata=0,
+        )
+        if inplace:
+            self.add_band(cloud_mask_band)
+            return
+        else:
+            return cloud_mask_band
+
     def get_water_mask(
             self,
             mask_band: Optional[str] = 'qa_pixel',
             inplace: Optional[bool] = False,
-            name_water_mask: Optional[str] = 'water_mask'
+            name_water_mask: Optional[str] = 'water_mask',
+            water_class: int = 7
     ) -> Band | None:
         """
-        Get a water mask from the pixel quality band (bit 7).
+        Get a water mask from the pixel quality band (bit 7 by default).
 
         NOTE:
             Since the `mask_band` can be set to *any* `Band` it is also
@@ -578,14 +644,20 @@ def get_water_mask(
             current `Landsat` instance. Default is False.
         :param name_water_mask:
             The name of the water mask band. Default is 'water_mask'.
+        :param water_class:
+            Bit number of the water class. 7 by default.
         :return:
             A `Band` instance containing the water mask or `None` if
             `inplace` is True.
         """
         water_mask = self.mask_from_qa_bits(
-            bit_range=(7, 7),
+            bit_range=(water_class, water_class),
             band_name=mask_band
-        )
+        ).astype('bool')
+        # update water mask with the mask of the area of interest,
+        if self[mask_band].is_masked_array:
+            water_mask = np.logical_and(
+                water_mask, ~self[mask_band].values.mask)
         water_mask_band = Band(
                 band_name=name_water_mask,
                 values=water_mask,
@@ -594,7 +666,7 @@ def get_water_mask(
             )
         if inplace:
             self.add_band(water_mask_band)
-            return None
+            return
         else:
             return water_mask_band
 
@@ -678,31 +750,21 @@ def mask_clouds_and_shadows(
         :param kwargs:
             optional kwargs to pass to `~eodal.core.raster.RasterCollection.mask`
         :returns:
-            depending on `inplace` (passed in the kwargs) a new `Sentinel2` instance
-            or None
+            depending on `inplace` (passed in the kwargs) a new `Landsat`
+            instance or None
         """
         if bands_to_mask is None:
             bands_to_mask = self.band_names
         # the mask band should never be masked
         if mask_band in bands_to_mask:
             bands_to_mask.remove(mask_band)
         try:
-            mask_list = []
-            for cloud_class in cloud_classes:
-                # construct the bit range
-                bit_range = (cloud_class, cloud_class)
-                # get the binary mask as boolean array
-                mask = self.mask_from_qa_bits(
-                    bit_range=bit_range,
-                    band_name=mask_band).astype('bool')
-                # add the mask to the list
-                mask_list.append(mask)
-
-            # combine the masks
-            cloud_mask = np.any(mask_list, axis=0)
+            cloud_mask = self.get_cloud_and_shadow_mask(
+                mask_band=mask_band,
+                cloud_classes=cloud_classes)
             # mask the bands
             return self.mask(
-                mask=cloud_mask,
+                mask=cloud_mask.values,
                 bands_to_mask=bands_to_mask,
                 **kwargs
             )

eodal/core/sensors/sentinel2.py

@@ -650,6 +650,54 @@ def mask_clouds_and_shadows(
         except Exception as e:
             raise Exception(f"Could not mask clouds and shadows: {e}")
 
+    def get_cloud_and_shadow_mask(
+        self,
+        mask_band: Optional[str] = 'scl',
+        inplace: Optional[bool] = False,
+        name_cloud_mask: Optional[str] = 'cloud_mask',
+        cloud_classes: Optional[list[int]] = [1, 2, 3, 7, 8, 9, 10, 11]
+    ) -> Sentinel2 | None:
+        """
+        Generate a binary cloud mask from the Scene Classification
+        Layer (SCL) that is part of the standard ESA product.
+
+        NOTE:
+            You can specify *any* band instance instead of SCL as
+            long as the band has an integer data type.
+
+        :param mask_band:
+            The name of the band to use for masking. Default is 'scl'.
+        :param inplace:
+            Whether to return a new `Band` instance or add a `Band` to the
+            current `Sentinel2` instance. Default is False.
+        :param name_cloud_mask:
+            The name of the water mask band. Default is 'water_mask'.
+        :param cloud_classes:
+            list of SCL values to be considered as clouds/shadows and snow.
+            By default, all three cloud classes and cloud shadows are considered
+            plus snow.
+        :return:
+            A `Band` instance containing the cloud mask or `None` if
+            `inplace` is True.
+        """
+        cloud_mask = np.isin(self[mask_band].values, cloud_classes)
+        # update cloud mask with the mask of the area of interest
+        if self[mask_band].is_masked_array:
+            cloud_mask = np.logical_and(
+                cloud_mask, ~self[mask_band].values.mask)
+
+        cloud_mask_band = Band(
+            band_name=name_cloud_mask,
+            values=cloud_mask,
+            geo_info=self[mask_band].geo_info,
+            nodata=0,
+        )
+        if inplace:
+            self.add_band(cloud_mask_band)
+            return
+        else:
+            return cloud_mask_band
+
     def get_scl_stats(self) -> pd.DataFrame:
         """
         Returns a ``DataFrame`` with the number of pixel for each

tests/core/test_landsat.py

@@ -3,6 +3,7 @@
 '''
 
 import geopandas as gpd
+import numpy as np
 import pytest
 
 from datetime import datetime
@@ -52,3 +53,48 @@ def test_landsat_from_usgs():
         for band_name in landsat.band_names[:-1]:
             assert 0 < landsat[band_name].values.min() <= 1, 'wrong value'
             assert 0 < landsat[band_name].values.max() <= 1, 'wrong value'
+
+        # get the cloud mask -> this should fail because we didn't
+        # read the qa bands
+        with pytest.raises(KeyError):
+            cloud_mask = landsat.get_cloud_and_shadow_mask()
+
+        # repeat the reading of the scene WITH the qa bands
+        band_selection = ['blue', 'green', 'red', 'nir08', 'swir12']
+        landsat = Landsat.from_usgs(
+            in_dir=landsat_scene_item['assets'],
+            vector_features=gdf,
+            read_qa=True,
+            read_atcor=False,
+            band_selection=band_selection
+        )
+        assert 'qa_pixel' in landsat.band_names, 'missing quality band'
+
+        # now the generation of a binary cloud mask should work
+        cloud_mask = landsat.get_cloud_and_shadow_mask()
+        assert cloud_mask.values.dtype == bool, 'expected boolean'
+
+        water_mask = landsat.get_water_mask()
+        assert water_mask.values.dtype == 'bool', 'expected boolean'
+
+        # mask clouds and shadows -> check if the mask has an effect
+        landsat.mask_clouds_and_shadows(inplace=True)
+        assert landsat['blue'].is_masked_array, 'expected as masked array'
+
+        # calculate the NDVI
+        landsat.calc_si(si_name='NDVI', inplace=True)
+        assert 'NDVI' in landsat.band_names
+        assert -1 <= landsat['NDVI'].values.min() <= 1
+        assert -1 <= landsat['NDVI'].values.max() <= 1
+        assert np.isnan(landsat['NDVI'].nodata)
+
+        # repeat the reading of the scene WITH the atcor bands
+        band_selection = ['blue', 'green', 'red', 'nir08', 'swir12']
+        landsat = Landsat.from_usgs(
+            in_dir=landsat_scene_item['assets'],
+            vector_features=gdf,
+            read_qa=True,
+            read_atcor=True,
+            band_selection=band_selection
+        )
+        assert 'lwir' in landsat.band_names

tests/core/test_sentinel2.py

@@ -179,6 +179,10 @@ def test_read_from_safe_with_mask_l2a(datadir, get_s2_safe_l2a, get_polygons, ge
     # make sure meta information was saved correctly
     assert handler['scl'].meta['dtype'] == 'uint8', 'wrong data type for SCL in meta'
 
+    # get a binary cloud mask
+    cloud_mask = handler.get_cloud_and_shadow_mask()
+    assert cloud_mask.values.dtype == bool, 'expected boolean'
+
     # to_xarray should fail because of different spatial resolutions
     with pytest.raises(ValueError):
         handler.to_xarray()
@@ -413,7 +417,7 @@ def test_read_from_safe_l2a(datadir, get_s2_safe_l2a):
         inplace=True
     )
 
-    assert reader.is_bandstack(), 'data should fulfill bandstack criteria but doesnot'
+    assert reader.is_bandstack(), 'data should fulfill bandstack criteria but does not'
     assert reader['blue'].crs == 32633, 'projection was not updated'
 
     # try writing bands to output file