diff --git a/eodal/core/band.py b/eodal/core/band.py index cf8556c..83ce314 100644 --- a/eodal/core/band.py +++ b/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): diff --git a/eodal/core/sensors/landsat.py b/eodal/core/sensors/landsat.py index 4d73d95..3ff321e 100644 --- a/eodal/core/sensors/landsat.py +++ b/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, diff --git a/tests/core/test_band.py b/tests/core/test_band.py index 1ab7580..84fbfec 100644 --- a/tests/core/test_band.py +++ b/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 ) diff --git a/tests/mapper/test_landsat_mapper.py b/tests/mapper/test_landsat_mapper.py new file mode 100644 index 0000000..19cf027 --- /dev/null +++ b/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 . +''' + +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] diff --git a/tests/mapper/test_sentinel2_mapper.py b/tests/mapper/test_sentinel2_mapper.py index dc91c6a..5eb0cbe 100644 --- a/tests/mapper/test_sentinel2_mapper.py +++ b/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]