Merge pull request #91 from lukasValentin/fix-reference-scene-heuristic

FIX: reference scene heuristic and uniform pixel size in `SceneCollection` returned from `mapper`

eodal/config/settings.py

@@ -7,7 +7,7 @@
 The ``Settings`` class uses ``pydantic``. This means all attributes of the class can
 be **overwritten** using environmental variables or a `.env` file.
 
-Copyright (C) 2022 Lukas Valentin Graf
+Copyright (C) 2022, 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

eodal/mapper/mapper.py

@@ -453,11 +453,18 @@ def _process_scene(
             scene = scene_modifier.__call__(scene, **scene_modifier_kwargs)
 
         # reproject scene if necessary
-        scene.reproject(
-            target_crs=item.target_epsg,
-            interpolation_method=reprojection_method,
-            inplace=True
-        )
+        # ..versionadd 0.2.3 check for the EPSG to make sure no unnecessary
+        # operations are undertaken to save runtime and avoid floating
+        # point inaccuracies
+        epsg_scene = [b.geo_info.epsg for _, b in scene]
+        intersection = set(epsg_scene).intersection(set([item.target_epsg]))
+        # we need to reproject only if the intersection returns an empty set.
+        if len(intersection) == 0: 
+            scene.reproject(
+                target_crs=item.target_epsg,
+                interpolation_method=reprojection_method,
+                inplace=True
+            )
 
         return scene
 
@@ -705,12 +712,34 @@ def _load_scenes_collection(
             ))
         bounds_gdf = pd.concat(bounds_list)
         total_bounds = bounds_gdf.total_bounds
+
+        # ..versionadd 0.2.3 (https://github.com/EOA-team/eodal/issues/90)
+        # we need to ensure that all scenes have not only the same extent but
+        # also the same pixel size and, hence, number of rows and columns
+        # The heuristic is:
+        # 1) is there a since that was not reprojected (_duplicated is False)
+        # 2) If no: use the first scene as reference scene
+        # 3) If yes: use the first scene that was not reprojected
+
+        # is there at a least a single scene that was not reprojected
+        if not self.metadata._duplicated.all():
+            reference_time = self.metadata[
+                ~self.metadata._duplicated][self.time_column].values[0]
+            timestamps = scoll.timestamps
+            reference_timestamp_idx = np.argmin(
+                [pd.to_datetime(x) - reference_time for x in timestamps])
+            reference_scene = scoll[
+                timestamps[reference_timestamp_idx]]
+        else:
+            reference_scene = scoll[scoll.timestamps[0]]
 
         # loop over scenes and project them onto the total bounds
         for _, scene in scoll:
-            if scene.is_bandstack():
-                # get a band of the SceneCollection
-                reference_band = scene[scene.band_names[0]]
+            # loop over the bands and reproject them
+            # onto the target grid
+            for band_name, band in scene:
+                # get the reference band
+                reference_band = reference_scene[band_name]
                 geo_info = reference_band.geo_info
                 # get the transform of the destination extent
                 minx, maxy = total_bounds[0], total_bounds[-1]
@@ -721,28 +750,25 @@ def _load_scenes_collection(
                     d=0,
                     e=geo_info.pixres_y,
                     f=maxy)
-
-                # loop over the bands and reproject them
-                # onto the target grid
-                for _, band in scene:
-                    dst_shape = (max(nrows), max(ncols))
-                    destination = np.zeros(
-                        dst_shape,
-                        dtype=band.values.dtype
-                    )
-                    # determine nodata
-                    if not np.isnan(band.nodata):
-                        dst_nodata = band.nodata
-                    else:
-                        # if band nodata is NaN we set
-                        # no-data to None (rasterio default)
-                        dst_nodata = None
-                    band.reproject(
-                        inplace=True,
-                        target_crs=reference_band.crs,
-                        dst_transform=dst_transform,
-                        destination=destination,
-                        dst_nodata=dst_nodata)
+
+                dst_shape = (max(nrows), max(ncols))
+                destination = np.zeros(
+                    dst_shape,
+                    dtype=band.values.dtype
+                )
+                # determine nodata
+                if not np.isnan(band.nodata):
+                    dst_nodata = band.nodata
+                else:
+                    # if band nodata is NaN we set
+                    # no-data to None (rasterio default)
+                    dst_nodata = None
+                band.reproject(
+                    inplace=True,
+                    target_crs=reference_band.crs,
+                    dst_transform=dst_transform,
+                    destination=destination,
+                    dst_nodata=dst_nodata)
 
         self.data = scoll
 

tests/conftest.py

@@ -254,6 +254,21 @@ def _get_polygons():
     return _get_polygons
 
 
+@pytest.fixture()
+def get_polygons_utm_zones(get_project_root_path):
+    """
+    Returns path to an area of interest at the boundary of two UTM zones
+    """
+    def _get_polygons():
+
+        testdata_dir = get_project_root_path.joinpath('data')
+        testdata_polys = testdata_dir.joinpath(
+            Path('sample_polygons').joinpath('utm_zones_western-ch.gpkg')
+        )
+        return testdata_polys
+    return _get_polygons
+
+
 @pytest.fixture()
 def get_scene_collection(get_bandstack):
     """fixture returing a SceneCollection with three scenes"""

tests/mapper/test_grid_alignment.py

@@ -0,0 +1,126 @@
+'''
+Ensure that all scenes in a SceneCollection have the same spatial
+extent, pixel size and number of rows and columns. We test this
+behavior for Sentinel-2 scenes at the boundary of two UTM zones.
+'''
+
+import geopandas as gpd
+import pytest
+
+from datetime import datetime
+from eodal.config import get_settings
+from eodal.core.raster import RasterCollection
+from eodal.core.sensors import Sentinel2
+from eodal.mapper.feature import Feature
+from eodal.mapper.filter import Filter
+from eodal.mapper.mapper import Mapper, MapperConfigs
+
+settings = get_settings()
+
+
+def test_grid_alignment(get_polygons_utm_zones):
+    """test grid alignment for Sentinel-2"""
+
+    settings.USE_STAC = True
+
+    time_start = datetime(2018, 3, 2)
+    time_end = datetime(2018, 4, 6)
+    metadata_filters = [
+        Filter('cloudy_pixel_percentage', '<', 70),
+        Filter('processing_level', '==', 'Level-2A')
+    ]
+    feature = Feature.from_geoseries(
+        gds=gpd.read_file(get_polygons_utm_zones()).geometry)
+    mapper_configs = MapperConfigs(
+        collection='sentinel2-msi',
+        time_start=time_start,
+        time_end=time_end,
+        feature=feature,
+        metadata_filters=metadata_filters
+    )
+
+    mapper = Mapper(mapper_configs)
+    mapper.query_scenes()
+
+    def resample(ds: RasterCollection, **kwargs):
+        return ds.resample(inplace=False, **kwargs)
+
+    scene_kwargs = {
+        'scene_constructor': Sentinel2.from_safe,
+        'scene_constructor_kwargs': {
+            'band_selection': ['red', 'red_edge_1'],
+            'apply_scaling': False
+        },
+        'scene_modifier': resample,
+        'scene_modifier_kwargs': {'target_resolution': 10}
+    }
+    mapper.load_scenes(scene_kwargs=scene_kwargs)
+
+    # all scenes should have the same extent and, since, we resampled
+    # all bands to a common spatial extent, the same pixel size and number
+    # of rows and columns
+    scoll = mapper.data
+
+    pixres_x, pixres_y, ulx, uly, nrows, ncols = [], [], [], [], [], []
+    for _, scene in scoll:
+        for _, band in scene:
+            geo_info = band.geo_info
+            pixres_x.append(geo_info.pixres_x)
+            pixres_y.append(geo_info.pixres_y)
+            ulx.append(geo_info.ulx)
+            uly.append(geo_info.uly)
+            nrows.append(band.nrows)
+            ncols.append(band.ncols)
+
+    assert len(set(pixres_x)) == 1
+    assert pixres_x[0] == 10
+    assert len(set(pixres_y)) == 1
+    assert pixres_y[0] == -10
+    assert len(set(ulx)) == 1
+    assert len(set(uly)) == 1
+    assert len(set(nrows)) == 1
+    assert len(set(ncols)) == 1
+
+    # this should also work when the bands are not bandstacked
+    # i.e., have a different pixel size
+    scene_kwargs = {
+        'scene_constructor': Sentinel2.from_safe,
+        'scene_constructor_kwargs': {
+            'band_selection': ['red', 'red_edge_1'],
+            'apply_scaling': False
+        }
+    }
+    mapper.load_scenes(scene_kwargs=scene_kwargs)
+
+    scoll = mapper.data
+
+    band_names = scoll[scoll.timestamps[0]].band_names
+    # the upper left corner must be always the same
+    ulx, uly = [], []
+    for band_name in band_names:
+        pixres_x, pixres_y, nrows, ncols = [], [], [], []
+        for _, scene in scoll:
+            band = scene[band_name]
+            geo_info = band.geo_info
+            pixres_x.append(geo_info.pixres_x)
+            pixres_y.append(geo_info.pixres_y)
+            ulx.append(geo_info.ulx)
+            uly.append(geo_info.uly)
+            nrows.append(band.nrows)
+            ncols.append(band.ncols)
+
+        assert len(set(pixres_x)) == 1
+        assert len(set(pixres_y)) == 1
+        assert len(set(nrows)) == 1
+        assert len(set(ncols)) == 1
+        # only the red band should have a spatial resolution of 10 m
+        # as we did not resample the 20 m bands
+        if band_name == 'B04':
+            assert pixres_x[0] == 10.
+            assert pixres_y[0] == -10.
+        else:
+            assert pixres_x[0] == 20.
+            assert pixres_y[0] == -20.
+
+    assert len(set(ulx)) == 1
+    assert len(set(uly)) == 1