Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix Landsat no-data issues #85

Merged
merged 6 commits into from
Oct 21, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 12 additions & 7 deletions eodal/core/band.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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):
Expand Down
12 changes: 8 additions & 4 deletions eodal/core/sensors/landsat.py
Original file line number Diff line number Diff line change
Expand Up @@ -450,22 +450,26 @@

# 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

Check warning on line 466 in eodal/core/sensors/landsat.py

View check run for this annotation

Codecov / codecov/patch

eodal/core/sensors/landsat.py#L466

Added line #L466 was not covered by tests

# 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,
Expand Down
7 changes: 7 additions & 0 deletions tests/core/test_band.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
)


Expand Down
145 changes: 145 additions & 0 deletions tests/mapper/test_landsat_mapper.py
Original file line number Diff line number Diff line change
@@ -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]
37 changes: 32 additions & 5 deletions tests/mapper/test_sentinel2_mapper.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
'''
Tests for the filter class.
Tests for the Mapper class for Sentinel-2.

.. versionadded:: 0.2.0

Expand All @@ -20,6 +20,7 @@
'''

import geopandas as gpd
import numpy as np
import pytest

from datetime import datetime
Expand All @@ -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),
Expand All @@ -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
Expand Down Expand Up @@ -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}
}
Expand All @@ -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]