Merge pull request #85 from lukasValentin/fix_landsat_mapper

Fix Landsat no-data issues

eodal/core/band.py

@@ -789,6 +789,18 @@ def from_rasterio(
             else:
                 area_or_point = kwargs["area_or_point"]
                 kwargs.pop("area_or_point")
+            # check no data value
+            if "nodata" not in kwargs.keys():
+                if src.nodata is not None:
+                    nodata = np.array([src.nodata]).astype(
+                        src.dtypes[0])[0]
+                else:
+                    raise TypeError(
+                        'Nodata must not be None. Please pass ' +
+                        '`nodata` argument')
+            else:
+                nodata = kwargs["nodata"]
+                kwargs.pop("nodata")
 
             # overwrite band_idx if band_name_src is provided
             band_names = list(src.descriptions)
@@ -854,13 +866,6 @@ def from_rasterio(
             if units is not None:
                 unit = units[band_idx - 1]
 
-        if "nodata" in kwargs.keys():
-            nodata = kwargs["nodata"]
-            kwargs.pop("nodata")
-        else:
-            nodata, nodata_vals = None, attrs.get("nodatavals", None)
-            if nodata_vals is not None:
-                nodata = nodata_vals[band_idx - 1]
         # make sure the nodata type matches the datatype of the
         # band values
         if band_data.dtype.kind not in 'fc' and np.isnan(nodata):

eodal/core/sensors/landsat.py

@@ -450,22 +450,26 @@ def from_usgs(
 
             # get alias name of band
             band_alias = None
+            nodata = None
             if band_name in landsat_band_mapping[sensor].values():
                 band_alias = [
                     k for k, v in landsat_band_mapping[sensor].items()
                     if v == band_name][0]
-            elif band_name in landsat_band_mapping['quality_flags'].keys() \
+            elif band_name in landsat_band_mapping['quality_flags'] \
                     or band_name in landsat_band_mapping['quality_flags'].values():
                 band_alias = band_name.lower()
-            elif band_name in landsat_band_mapping['atmospheric_correction'].keys() \
+                nodata = 0
+            elif band_name in landsat_band_mapping['atmospheric_correction'] \
                     or band_name in landsat_band_mapping['atmospheric_correction'
                                                          ].values():
                 band_alias = band_name.lower()
+                nodata = 0
 
             # read band
             try:
-                if band_name in landsat_band_mapping[sensor].values():
-                    kwargs.update({'scale': gain, 'offset': offset})
+                kwargs.update({'scale': gain, 'offset': offset})
+                if nodata is not None:
+                    kwargs.update({'nodata': nodata})
                 landsat.add_band(
                     Band.from_rasterio,
                     fpath_raster=band_fpath,

tests/core/test_band.py

@@ -110,6 +110,13 @@ def test_band_from_rasterio(get_test_band, get_bandstack):
         Band.from_rasterio(
             fpath_raster=fpath_raster,
             band_idx=22,
+            nodata=0
+        )
+    # nodata value no specified
+    with pytest.raises(TypeError):
+        Band.from_rasterio(
+            fpath_raster=fpath_raster,
+            band_idx=22
         )
 
 

tests/mapper/test_landsat_mapper.py

@@ -0,0 +1,145 @@
+'''
+Tests for the Mapper class for Sentinel-2.
+
+.. versionadded:: 0.2.0
+
+Copyright (C) 2023 Lukas Valentin Graf
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+'''
+
+import geopandas as gpd
+import numpy as np
+import pytest
+
+from datetime import datetime
+from shapely.geometry import box, Polygon
+
+from eodal.config import get_settings
+from eodal.core.raster import RasterCollection
+from eodal.core.scene import SceneCollection
+from eodal.core.sensors import Landsat
+from eodal.mapper.feature import Feature
+from eodal.mapper.filter import Filter
+from eodal.mapper.mapper import Mapper, MapperConfigs
+
+Settings = get_settings()
+
+
+def preprocess_landsat_scene(
+        ds: Landsat
+) -> Landsat:
+    """
+    Mask clouds and cloud shadows in a Landsat scene based
+    on the 'qa_pixel' band.
+
+    :param ds:
+        Landsat scene before cloud mask applied.
+    :return:
+        Landsat scene with clouds and cloud shadows masked.
+    """
+    ds.mask_clouds_and_shadows(inplace=True)
+    return ds
+
+
+@pytest.mark.parametrize(
+    'collection,time_start,time_end,geom,metadata_filters,apply_scaling',
+    [(
+        'landsat-c2-l2',
+        datetime(2022,7,1),
+        datetime(2022,7,15),
+        Polygon(
+            [[7.04229, 47.01202],
+            [7.08525, 47.01202],
+            [7.08525, 46.96316],
+            [7.04229, 46.96316],
+            [7.04229, 47.01202]]
+        ),
+        [Filter('eo:cloud_cover','<', 80)],
+        False
+    ),(
+        'landsat-c2-l2',
+        datetime(2022,7,1),
+        datetime(2022,7,15),
+        Polygon(
+            [[7.04229, 47.01202],
+            [7.08525, 47.01202],
+            [7.08525, 46.96316],
+            [7.04229, 46.96316],
+            [7.04229, 47.01202]]
+        ),
+        [Filter('eo:cloud_cover','<', 80)],
+        True
+    ),
+    (
+        'landsat-c2-l1',
+        datetime(1972, 9, 1),
+        datetime(1972, 10, 31),
+        box(*[8.4183, 47.2544, 8.7639, 47.5176]),
+        [Filter('eo:cloud_cover','<', 80)],
+        False
+    )
+    ]
+)
+def test_landsat_mapper(
+        collection, time_start, time_end, geom, metadata_filters,
+        apply_scaling):
+    """
+    Test the mapper class for handling Landsat collection 2 data
+    (level 1 and 2).
+    """
+    Settings.USE_STAC = True
+    feature = Feature(
+        name='Test Area',
+        geometry=geom,
+        epsg=4326,
+        attributes={'id': 1}
+    )
+    mapper_configs = MapperConfigs(
+        collection=collection,
+        time_start=time_start,
+        time_end=time_end,
+        feature=feature,
+        metadata_filters=metadata_filters
+    )
+    mapper = Mapper(mapper_configs)
+    mapper.query_scenes()
+    assert isinstance(mapper.metadata, gpd.GeoDataFrame), 'expected a GeoDataFrame'
+    assert not mapper.metadata.empty, 'metadata must not be empty'
+
+    scene_kwargs = {
+        'scene_constructor': Landsat.from_usgs,
+        'scene_constructor_kwargs': {
+            'band_selection': ['red', 'nir08'],
+            'read_qa': True,
+            'apply_scaling': apply_scaling
+        },
+        'scene_modifier': preprocess_landsat_scene
+    }
+    mapper.load_scenes(scene_kwargs=scene_kwargs)
+
+    assert isinstance(mapper.data, SceneCollection), 'expected a SceneCollection'
+    assert mapper.metadata.target_epsg.nunique() == 1, 'there must not be more than a single target EPSG'
+    scenes_crs = [x['red'].crs.to_epsg() for _, x in mapper.data]
+    assert (scenes_crs == mapper.metadata.target_epsg.unique()[0]).all(), \
+        'all scenes must be projected into the target EPSG'
+    assert mapper.metadata.shape[0] == len(mapper.data), \
+        'mis-match between length of metadata and data'
+
+    for _, scene in mapper.data:
+        dtype = scene['red'].values.dtype
+        if apply_scaling:
+            assert dtype in [np.float32, np.float64]
+        else:
+            assert dtype in [np.uint8, np.uint16]

tests/mapper/test_sentinel2_mapper.py

@@ -1,5 +1,5 @@
 '''
-Tests for the filter class.
+Tests for the Mapper class for Sentinel-2.
 
 .. versionadded:: 0.2.0
 
@@ -20,6 +20,7 @@
 '''
 
 import geopandas as gpd
+import numpy as np
 import pytest
 
 from datetime import datetime
@@ -36,7 +37,7 @@
 Settings = get_settings()
 
 @pytest.mark.parametrize(
-    'collection,time_start,time_end,geom,metadata_filters',
+    'collection,time_start,time_end,geom,metadata_filters,apply_scaling',
     [(
         'sentinel2-msi',
         datetime(2022,7,1),
@@ -49,9 +50,25 @@
             [7.04229, 47.01202]]
         ),
         [Filter('cloudy_pixel_percentage','<', 80), Filter('processing_level', '==', 'Level-2A')],
-    ),]
+        False
+    ),
+        (
+        'sentinel2-msi',
+        datetime(2022,7,1),
+        datetime(2022,7,15),
+        Polygon(
+            [[7.04229, 47.01202],
+            [7.08525, 47.01202],
+            [7.08525, 46.96316],
+            [7.04229, 46.96316],
+            [7.04229, 47.01202]]
+        ),
+        [Filter('cloudy_pixel_percentage','<', 80), Filter('processing_level', '==', 'Level-2A')],
+        True
+    )]
 )
-def test_sentinel2_mapper(collection, time_start, time_end, geom, metadata_filters):
+def test_sentinel2_mapper(collection, time_start, time_end, geom, metadata_filters,
+                          apply_scaling):
     """
     Test the mapper class for handling Sentinel-2 data including mosaicing data
     from two different Sentinel-2 tiles
@@ -80,7 +97,10 @@ def resample(ds: RasterCollection, **kwargs):
 
     scene_kwargs = {
         'scene_constructor': Sentinel2.from_safe,
-        'scene_constructor_kwargs': {'band_selection': ['B04', 'B08']},
+        'scene_constructor_kwargs': {
+            'band_selection': ['B04', 'B08'],
+            'apply_scaling': apply_scaling
+        },
         'scene_modifier': resample,
         'scene_modifier_kwargs': {'target_resolution': 10}
     }
@@ -93,3 +113,10 @@ def resample(ds: RasterCollection, **kwargs):
         'all scenes must be projected into the target EPSG'
     assert mapper.metadata.shape[0] == len(mapper.data), \
         'mis-match between length of metadata and data'
+
+    for _, scene in mapper.data:
+        dtype = scene['red'].values.dtype
+        if apply_scaling:
+            assert dtype in [np.float32, np.float64]
+        else:
+            assert dtype in [np.uint8, np.uint16]