Merge pull request #65 from lukasValentin/dev

Draft: Fix different grid sizes in Mapper class (#64)

.gitignore

@@ -239,5 +239,6 @@ Temporary Items
 
 eodal.code-workspace
 .vscode/extensions.json
+.vscode/
 
 *.xml

CHANGELOG.rst

@@ -10,6 +10,19 @@ The format is based on `Keep a Changelog`_, and this project adheres to `Semanti
 
 Categories for changes are: Added, Changed, Deprecated, Removed, Fixed, Security.
 
+Version `0.2.1 < https://github.com/EOA-team/eodal/releases/tag/v0.2.1>`__
+--------------------------------------------------------------------------------
+
+Release date: 2023-06-05
+
+- Fixed: Small bugs when loading the EOdal SceneCollection from a pickled binary object (thanks to @atoparseks) (#55).
+- Fixed: The color of no-data values of RasterCollection objects is now set to black (instead of white) when plotting the data (#56).
+- Added: The user can now specify a custom color map when plotting RasterCollection objects (#56).
+- Changed: The user can now specify a custom directory for writing log files to (#52).
+- Added: The MapperConfigs record the data source from which satellite data was read (#62).
+- Fixed: A work-around using HTTP retries has been implemented to surpass HTTP 500 errors when connecting to MS Planetary Computer (#58).
+- Fixed: All scenes in a SceneCollection returned from the mapper now share the same spatial extent and are aligned on a common reference grid (important for re-projections) (#64).
+
 Version `0.2.0 < https://github.com/EOA-team/eodal/releases/tag/v0.2.0>`__
 --------------------------------------------------------------------------------
 

eodal/__meta__.py

@@ -15,4 +15,4 @@
 description = "Earth Observation Data Analysis Library"  # One-liner
 url = "https://github.com/EOA-team/eodal"  # your project home-page
 license = "GNU General Public License version 3"  # See https://choosealicense.com
-version = "0.2.0"
+version = "0.2.1"

eodal/core/band.py

@@ -176,7 +176,7 @@ class GeoInfo(object):
 
     def __init__(
         self,
-        epsg: int,
+        epsg: int | CRS,
         ulx: Union[int, float],
         uly: Union[int, float],
         pixres_x: Union[int, float],
@@ -205,10 +205,13 @@ def __init__(
             system.
         """
         # make sure the EPSG code is valid
-        try:
-            CRS.from_epsg(epsg)
-        except Exception as e:
-            raise ValueError(e)
+        if isinstance(epsg, CRS):
+            epsg = epsg.to_epsg()
+        else:
+            try:
+                CRS.from_epsg(epsg)
+            except Exception as e:
+                raise ValueError(e)
 
         object.__setattr__(self, "epsg", epsg)
         object.__setattr__(self, "ulx", ulx)
@@ -1918,9 +1921,7 @@ def resample(
     def reproject(
         self,
         target_crs: Union[int, CRS],
-        dst_transform: Optional[Affine] = None,
         interpolation_method: Optional[int] = Resampling.nearest,
-        num_threads: Optional[int] = 1,
         inplace: Optional[bool] = False,
         **kwargs,
     ):
@@ -1936,13 +1937,12 @@ def reproject(
         :param interpolation_method:
             interpolation method to use for interpolating grid cells after
             reprojection. Default is neares neighbor interpolation.
-        :param num_threads:
-            number of threads to use for the operation. Uses a single thread by
-            default.
         :param inplace:
             if False (default) returns a copy of the ``Band`` instance
             with the changes applied. If True overwrites the values
             in the current instance.
+        :param kwargs:
+            optional keyword arguments to pass to ``rasterio.warp.reproject``.
         :returns:
             ``Band`` instance if `inplace` is False, None instead.
         """
@@ -1953,9 +1953,7 @@ def reproject(
             "src_transform": self.transform,
             "dst_crs": target_crs,
             "src_nodata": self.nodata,
-            "resampling": interpolation_method,
-            "num_threads": num_threads,
-            "dst_transform": dst_transform,
+            "resampling": interpolation_method
         }
         reprojection_options.update(kwargs)
 
@@ -1971,13 +1969,11 @@ def reproject(
 
         try:
             # set destination array in case dst_transfrom is provided
-            if (
-                "dst_transform" in reprojection_options.keys()
-                and reprojection_options.get("dst_transfrom") is not None
-            ):
+            if reprojection_options.get("dst_transfrom") is not None:
                 if "destination" not in reprojection_options.keys():
-                    dst = np.zeros_like(band_data)
-                    reprojection_options.update({"destination": dst})
+                    raise ValueError(
+                        '"destination" must be provided ' +
+                        'alongside "dst_transform"')
 
             out_data, out_transform = reproject_raster_dataset(
                 raster=band_data, **reprojection_options
@@ -1986,14 +1982,20 @@ def reproject(
             raise ReprojectionError(f"Could not re-project band {self.band_name}: {e}")
 
         # cast array back to original dtype
-        out_data = out_data[0, :, :].astype(self.values.dtype)
+        if len(out_data.shape) == 2:
+            out_data = out_data.astype(self.values.dtype)
+        elif len(out_data.shape) == 3:
+            out_data = out_data[0, :, :].astype(self.values.dtype)
 
         # reproject the mask separately
         if self.is_masked_array:
             out_mask, _ = reproject_raster_dataset(
                 raster=band_mask, **reprojection_options
             )
-            out_mask = out_mask[0, :, :].astype(bool)
+            if len(out_mask.shape) == 2:
+                out_mask = out_mask.astype(bool)
+            elif len(out_mask.shape) == 3:
+                out_mask = out_mask[0, :, :].astype(bool)
             # mask also those pixels which were set to nodata after reprojection
             # due to the raster alignment
             nodata = reprojection_options.get("src_nodata", 0)

eodal/mapper/mapper.py

@@ -36,6 +36,7 @@
 from datetime import datetime
 from pathlib import Path
 from sqlalchemy.exc import DatabaseError
+from rasterio import Affine
 from shapely.geometry import box
 from typing import Any, Callable, Dict, List, Optional
 
@@ -221,6 +222,7 @@ def from_yaml(cls, fpath: str | Path):
             # overwrite data_source attribute with content from
             # file in case it is available
             configs.data_source = yaml_content.get('data_source', 'UNKNOWN')
+            return configs
         except KeyError as e:
             raise ValueError(f"Missing keys in yaml file: {e}")
 
@@ -438,7 +440,7 @@ def _process_scene(
         scene.reproject(
             target_crs=item.target_epsg,
             interpolation_method=reprojection_method,
-            inplace=True,
+            inplace=True
         )
 
         return scene
@@ -462,6 +464,24 @@ def _load_scenes_collection(
         Mosaicing operations are handled on the fly so calling programs will
         always receive analysis-ready data.
 
+        IMPORTANT:
+
+            When mosaicing is necessary, it **must** be ensured that all bands
+            have the spatial resolution. Thus, either read only bands with the
+            same spatial resolution or use the `scene_modifier` function to
+            resample the bands into a common spatial resolution. If this
+            requirement is not fulfilled the mosaicing will fail!
+
+        ..versionadd:: 0.2.1
+
+            To make sure all scenes in the collection share the same extent and
+            grid, a post-processing step is carried out re-projecting all scenes
+            into a common reference grid using nearest neighbor interpolation.
+            This ensures not only all scenes having the same spatial extent which
+            is convenient for stacking scenes in time, but also makes sure there
+            are no pixel shifts due to re-projections from different spatial
+            reference systems such as different UTM zones.
+
         :param scene_constructor:
             Callable used to read the scenes found into `RasterCollection` fulfilling
             the `is_scene` criterion (i.e., a time stamp is available). The callable is
@@ -616,7 +636,59 @@ def _load_scenes_collection(
 
         # sort scenes by their timestamps and save as data attribute
         # to mapper instance
-        self.data = scoll.sort()
+        scoll = scoll.sort()
+
+        # ..versionadd:: 0.2.1
+        # we have to make sure that the scenes in self.data all share the same
+        # extent and are aligned onto the same grid.
+        # We use the boundary of all scenes to update the upper left corner
+        bounds_list = []
+        ncols = []
+        nrows = []
+        for _, scene in scoll:
+            bounds = scene[scene.band_names[0]].bounds
+            ncols.append(scene[scene.band_names[0]].ncols)
+            nrows.append(scene[scene.band_names[0]].nrows)
+            crs = scene[scene.band_names[0]].crs
+            bounds_list.append(gpd.GeoDataFrame(
+                geometry=[bounds],
+                crs=crs
+            ))
+        bounds_gdf = pd.concat(bounds_list)
+        total_bounds = bounds_gdf.total_bounds
+
+        # 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]]
+                geo_info = reference_band.geo_info
+                # get the transform of the destination extent
+                minx, maxy = total_bounds[0], total_bounds[-1]
+                dst_transform = Affine(
+                    a=geo_info.pixres_x,
+                    b=0,
+                    c=minx,
+                    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
+                    )
+                    band.reproject(
+                        inplace=True,
+                        target_crs=reference_band.crs,
+                        dst_transform=dst_transform,
+                        destination=destination)
+
+        self.data = scoll
+
         logger.info(f"Finished extraction of {self.sensor} scenes")
 
     def _load_pixels(

scripts/sentinel2_mapper.py

@@ -36,16 +36,17 @@
 import geopandas as gpd
 
 from datetime import datetime
+from pathlib import Path
+from shapely.geometry import box
+from typing import List
+
 from eodal.config import get_settings
 from eodal.core.scene import SceneCollection
 from eodal.core.sensors.sentinel2 import Sentinel2
 from eodal.mapper.feature import Feature
 from eodal.mapper.filter import Filter
 from eodal.mapper.mapper import Mapper, MapperConfigs
 
-from pathlib import Path
-from typing import List
-
 
 Settings = get_settings()
 # set to False to use a local data archive
@@ -79,16 +80,26 @@ def preprocess_sentinel2_scenes(
 
 if __name__ == '__main__':
 
+    import os
+    cwd = Path(__file__).absolute().parent.parent
+    os.chdir(cwd)
+
     # user-inputs
     # -------------------------- Collection -------------------------------
     collection: str = 'sentinel2-msi'
 
     # ------------------------- Time Range ---------------------------------
-    time_start: datetime = datetime(2022, 6, 1)  		# year, month, day (incl.)
-    time_end: datetime = datetime(2022, 6, 30)   		# year, month, day (incl.)
+    time_start: datetime = datetime(2022, 6, 10)  		# year, month, day (incl.)
+    time_end: datetime = datetime(2022, 6, 15)   		# year, month, day (incl.)
 
     # ---------------------- Spatial Feature  ------------------------------
-    geom: Path = Path('data/sample_polygons/lake_lucerne.gpkg')
+    bbox = box(*[6.5738, 46.4565, 7.2628, 47.2190])
+    feature = Feature(
+        name='lake_neuchatel',
+        geometry=bbox,
+        epsg=4326,
+        attributes={}
+    )
 
     # ------------------------- Metadata Filters ---------------------------
     metadata_filters: List[Filter] = [
@@ -97,7 +108,6 @@ def preprocess_sentinel2_scenes(
     ]
 
     # query the scenes available (no I/O of scenes, this only fetches metadata)
-    feature = Feature.from_geoseries(gpd.read_file(geom).geometry)
     mapper_configs = MapperConfigs(
         collection=collection,
         time_start=time_start,
@@ -115,16 +125,11 @@ def preprocess_sentinel2_scenes(
     # the metadata is loaded into a GeoPandas GeoDataFrame
     mapper.metadata
 
-    # get the least cloudy scene
-    mapper.metadata = mapper.metadata[
-          mapper.metadata.cloudy_pixel_percentage ==
-          mapper.metadata.cloudy_pixel_percentage.min()].copy()
-
     # load the least cloudy scene available from STAC
     scene_kwargs = {
         'scene_constructor': Sentinel2.from_safe,
-        'scene_constructor_kwargs': {'band_selection':
-                                     ["B02", "B03", "B04", "B05", "B8A"]}}
+        'scene_constructor_kwargs': {'band_selection': [
+            "B02", "B03", "B04"], 'read_scl': False}}
 
     mapper.load_scenes(scene_kwargs=scene_kwargs)
     # the data loaded into `mapper.data` as a EOdal SceneCollection

setup.py

@@ -140,7 +140,7 @@ def run(self):
     name='eodal',
     # setup_requires=['setuptools_scm'],
     # use_scm_version={'version_scheme': 'python-simplified-semver'},
-    version='0.2.0',
+    version='0.2.1',
     description='The Earth Observation Data Analysis Library EOdal',
     long_description=long_description,
     long_description_content_type='text/markdown',
@@ -158,9 +158,12 @@ def run(self):
      "Natural Language :: English",
      "Programming Language :: Python :: 3",
      "Operating System :: OS Independent",
+     "Programming Language :: Python :: 3.6",
+     "Programming Language :: Python :: 3.7",
      "Programming Language :: Python :: 3.8",
      "Programming Language :: Python :: 3.9",
-     "Programming Language :: Python :: 3.10"
+     "Programming Language :: Python :: 3.10",
+     "Programming Language :: Python :: 3.11",
     ],
     # Could also include keywords, download_url, project_urls, etc.
     # Custom commands

tests/core/test_scene_collection.py

@@ -1,7 +1,5 @@
 '''
-Created on Nov 24, 2022
-
-@author: graflu
+Tests for the SceneCollection class.
 '''
 
 import pytest
@@ -48,7 +46,7 @@ def _generate_random(number: int, polygon: Polygon) -> List[Point]:
     return _generate_random
 
 
-def test_raster_is_scene(apply_settings, get_bandstack):
+def test_raster_is_scene(get_bandstack):
     """test the is_scene attribute of RasterCollections"""
 
     fpath_raster = get_bandstack()