Repository: EOA-team/eodal (language: python)

Issue:
Sentinel2Mapper unique_id_attribute error
If unique_id_attribute = 'id' points to field of type 'string', everything works as expected.

If unique_id_attribute = 'id' points to field of type 'int', following example happens:

Example:
# get a new mapper instance
mapper = Sentinel2Mapper(
	date_start=date_start,
	date_end=date_end,
	processing_level=processing_level,
	cloud_cover_threshold=scene_cloud_cover_threshold,
	mapper_configs=mapper_configs,
	feature_collection=aoi,
    unique_id_attribute = 'id'
)


File "/home/miguel/code/eodal_fork/eodal/eodal/operational/mapping/sentinel2.py", line 428, in get_complete_timeseries
    if isinstance(res, gpd.GeoDataFrame):

UnboundLocalError: local variable 'res' referenced before assignment
Mapper routine for polygon-features: Different dimensions [xsize,ysize] returned for dates with two (and maybe more) acquisitions
sorry, no code-example at hand. Maybe different polygon extraction approaches are implemented (if a merging does somehow happen).

Creates issues in terms of data inconsistency and stacking of arrays from returned list of raster-collections.


Discussion / hints:
I suggest to make this issue part of the "next generation of the Mapper class" (see #10 )
added the `wontfix` label since the `mapper` class will undergo substantial reworking (targeted for EOdal version 0.**3**.0)


Fix patch:
diff --git a/README.md b/README.md
new file mode 100644
index 00000000..ea549238
--- /dev/null
+++ b/README.md
@@ -0,0 +1,42 @@
+|GHA tests| |Codecov report| |pre-commit| |black|
+
+
+#E:earth_africa:dal Earth Observation Data Analysis Library
+**A truely open-source package for unified analysis of Earth Observation (EO) data**
+
+## About E:earth_africa:dal
+
+E:earth_africa:dal is a Python library enabling the acquisition, organization, and analysis of Earth observation data in a completely open-source manner within a unified framework.
+
+E:earth_africa:dal enables open-source, reproducible geo-spatial data science. At the same time, E:earth_africa:dal lowers the burden of data handling and provides access to **global satellite data archives** through **downloaders** and the fantastic **SpatioTemporalAssetsCatalogs** (STAC).
+
+E:earth_africa:dal supports working in **cloud-environments** using [STAC catalogs](https://stacspec.org/) ("online" mode) and on **local premises** using a spatial PostgreSQL/PostGIS database to organize metadata ("offline" mode).
+
+Read more about E:earth_africa:dal in [this peer reviewed article](https://doi.org/10.1016/j.compag.2022.107487).
+
+## Citing E:earth_africa:dal
+
+We put a lot of effort in developing E:earth_africa:dal. To give us proper credit please respect our [license agreement](LICENSE). When you use E:earth_africa:dal for your **research** please **cite our paper** in addition to give us proper scientific credit.
+
+.. code::latex
+
+	@article{GRAF2022107487,
+	title = {EOdal: An open-source Python package for large-scale agroecological research using Earth Observation and gridded environmental data},
+	journal = {Computers and Electronics in Agriculture},
+	volume = {203},
+	pages = {107487},
+	year = {2022},
+	issn = {0168-1699},
+	doi = {https://doi.org/10.1016/j.compag.2022.107487},
+	url = {https://www.sciencedirect.com/science/article/pii/S0168169922007955},
+	author = {Lukas Valentin Graf and Gregor Perich and Helge Aasen},
+	keywords = {Satellite data, Python, Open-source, Earth Observation, Ecophysiology},
+	abstract = {Earth Observation by means of remote sensing imagery and gridded environmental data opens tremendous opportunities for systematic capture, quantification and interpretation of plant–environment interactions through space and time. The acquisition, maintenance and processing of these data sources, however, requires a unified software framework for efficient and scalable integrated spatio-temporal analysis taking away the burden of data and file handling from the user. Existing software products either cover only parts of these requirements, exhibit a high degree of complexity, or are closed-source, which limits reproducibility of research. With the open-source Python library EOdal (Earth Observation Data Analysis Library) we propose a novel software that enables the development of fully reproducible spatial data science chains through the strict use of open-source developments. Thanks to its modular design, EOdal enables advanced data warehousing especially for remote sensing data, sophisticated spatio-temporal analysis and intersection of different data sources, as well as nearly unlimited expandability through application programming interfaces (APIs).}
+	}
+
+## Examples
+We have compiled a set of [Jupyter notebooks](https://github.com/EOA-team/eodal_notebooks) showing you the capabilities of E:earth_africa:dal and how to unlock them.
+
+## Contributing
+
+Contributions to E:earth_africa:dal are welcome. Please make sure to read the [contribution guidelines](Contributing.rst) first.
diff --git a/README.rst b/README.rst
deleted file mode 100644
index 7164e7dd..00000000
--- a/README.rst
+++ /dev/null
@@ -1,131 +0,0 @@
-|GHA tests| |Codecov report| |pre-commit| |black|
-
-""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
-E:earth_africa:dal Earth Observation Data Analysis Library
-""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
-
-................................................................................
-A truely open-source package for unified analysis of Earth Observation (EO) data
-................................................................................
-
-.. contents:: Overview
-   :depth: 3
-
-========================
-About E:earth_africa:dal
-========================
-
-E:earth_africa:dal is a Python library enabling the acquisition, organization, and analysis of Earth observation data in a completely open-source manner.
-
-E:earth_africa:dal Python allows to
-
-	* load
-	* modify
-	* analyze
-	* write
-	* and interface
-
-EO data within an unified framework. E:earth_africa:dal thus enables open-source, reproducible geo-spatial data science while lowering the burden of data handling on the user-side.
-E:earth_africa:dal supports working in **cloud-environments** using [STAC catalogs](https://stacspec.org/) ("online" mode) and
-on **local premises** using a spatial PostgreSQL/PostGIS database to organize metadata ("offline" mode).
-
-Read more about E:earth_africa:dal in [this peer reviewed article](https://doi.org/10.1016/j.compag.2022.107487).
-
-========================
-Citing E:earth_africa:dal
-========================
-
-When using EOdal not only refer to our [license agreement](LICENSE) but also cite us properly:
-
-.. code::latex
-
-	@article{GRAF2022107487,
-	title = {EOdal: An open-source Python package for large-scale agroecological research using Earth Observation and gridded environmental data},
-	journal = {Computers and Electronics in Agriculture},
-	volume = {203},
-	pages = {107487},
-	year = {2022},
-	issn = {0168-1699},
-	doi = {https://doi.org/10.1016/j.compag.2022.107487},
-	url = {https://www.sciencedirect.com/science/article/pii/S0168169922007955},
-	author = {Lukas Valentin Graf and Gregor Perich and Helge Aasen},
-	keywords = {Satellite data, Python, Open-source, Earth Observation, Ecophysiology},
-	abstract = {Earth Observation by means of remote sensing imagery and gridded environmental data opens tremendous opportunities for systematic capture, quantification and interpretation of plant–environment interactions through space and time. The acquisition, maintenance and processing of these data sources, however, requires a unified software framework for efficient and scalable integrated spatio-temporal analysis taking away the burden of data and file handling from the user. Existing software products either cover only parts of these requirements, exhibit a high degree of complexity, or are closed-source, which limits reproducibility of research. With the open-source Python library EOdal (Earth Observation Data Analysis Library) we propose a novel software that enables the development of fully reproducible spatial data science chains through the strict use of open-source developments. Thanks to its modular design, EOdal enables advanced data warehousing especially for remote sensing data, sophisticated spatio-temporal analysis and intersection of different data sources, as well as nearly unlimited expandability through application programming interfaces (APIs).}
-	}
-
-==============================
-Examples of E:earth_africa:dal
-==============================
-
-The following code snippet reads spectral bands from a Sentinel-2 scene
-organized in .SAFE folder structure acquired over Southern Germany in
-Level2A (bottom-of-atmosphere reflectance). The Sentinel-2 scene can be
-downloaded `here <https://data.mendeley.com/datasets/ckcxh6jskz/1>`__ (
-S2A_MSIL2A_20190524T101031_N0212_R022_T32UPU_20190524T130304.zip):
-
-.. code:: python
-
-   import geopandas as gpd
-   from pathlib import Path
-   from shapely.geometry import Polygon
-   from eodal.core.sensors import Sentinel2
-
-   # file-path to the .SAFE dataset
-   dot_safe_dir = Path('../data/S2A_MSIL2A_20190524T101031_N0212_R022_T32UPU_20190524T130304.SAFE')
-
-   # construct a bounding box for reading a spatial subset of the scene (geographic coordinates)
-   ymin, ymax = 47.949, 48.027
-   xmin, xmax = 11.295, 11.385
-   bbox = Polygon([(xmin, ymin), (xmax, ymin), (xmax, ymax), (xmin, ymax)])
-
-   # eodal expects a vector file or a GeoDataFrame for spatial sub-setting
-   bbox_gdf = gpd.GeoDataFrame(geometry=[bbox], crs=4326)
-
-   # read data from .SAFE (all 10 and 20m bands + scene classification layer)
-   s2_ds = Sentinel2().from_safe(
-       in_dir=dot_safe_dir,
-       vector_features=bbox_gdf
-   )
-
-   # eodal support band aliasing. Thus, you can access the bands by their name ...
-   s2_ds.band_names
-
-Output
-
-.. code:: shell
-
-   >>> ['B02', 'B03', 'B04', 'B05', 'B06', 'B07', 'B08', 'B8A', 'B11', 'B12', 'SCL']
-
-.. code:: python
-
-   # ... or by their alias, i.e., their color name
-   s2_ds.band_aliases
-
-.. code:: shell
-
-   >>> ['blue', 'green', 'red', 'red_edge_1', 'red_edge_2', 'red_edge_3', 'nir_1', 'nir_2', 'swir_1', 'swir_2', 'scl']
-
-.. code:: python
-
-   # plot false-color infrared preview
-   s2_ds.plot_multiple_bands(band_selection=['nir_1','red','green'])
-
-.. image:: img/eodal_Sentinel-2_NIR.png
-  :width: 400
-  :alt: Sentinel-2 False-Color NIR
-
-.. code:: python
-
-   # plot scene classification layer
-   s2_ds.plot_scl()
-
-.. image:: img/eodal_Sentinel-2_SCL.png
-  :width: 400
-  :alt: Sentinel-2 Scene classification layer
-
-
-============
-Contributing
-============
-
-Contributions to E:earth_africa:dal are welcome. Please make sure to read the [contribution guidelines](Contributing.rst) first.
diff --git a/data/.gitignore b/data/.gitignore
index 14efea52..a119923f 100644
--- a/data/.gitignore
+++ b/data/.gitignore
@@ -5,3 +5,5 @@
 /MCARI.tif
 20220414_101133_47_227b/
 
+/sample_mapper_call.yaml
+/sample_mapper_data.pkl
diff --git a/eodal/__meta__.py b/eodal/__meta__.py
index 5db25e70..4b3ad75e 100644
--- a/eodal/__meta__.py
+++ b/eodal/__meta__.py
@@ -14,4 +14,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.1.1"
+version = "0.2.0"
diff --git a/eodal/config/__init__.py b/eodal/config/__init__.py
index 348e9b9e..17751ea1 100644
--- a/eodal/config/__init__.py
+++ b/eodal/config/__init__.py
@@ -1,3 +1,2 @@
 from .settings import get_settings
-from .sentinel2 import Sentinel2
 from .stac_providers import STAC_Providers
diff --git a/eodal/config/sentinel2.py b/eodal/config/sentinel2.py
deleted file mode 100644
index 45cbb4d5..00000000
--- a/eodal/config/sentinel2.py
+++ /dev/null
@@ -1,61 +0,0 @@
-"""
-Global settings for Sentinel-2 that can be customized if required.
-
-The ``Sentinel2`` 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
-
-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/>.
-
-"""
-
-from __future__ import annotations
-
-from pydantic.main import BaseModel
-from typing import List
-
-
-class Sentinel2(BaseModel):
-    """
-    base class defining Sentinel-2 product, archive and sensor details
-    """
-
-    PROCESSING_LEVELS: List[str] = ["L1C", "L2A"]
-
-    SPATIAL_RESOLUTIONS: dict = {
-        60.0: ["B01", "B09", "B10"],
-        10.0: ["B02", "B03", "B04", "B08"],
-        20.0: ["B05", "B06", "B07", "B8A", "B11", "B12", "SCL"],
-    }
-    BAND_INDICES: dict = {
-        "B01": 0,
-        "B02": 1,
-        "B03": 2,
-        "B04": 3,
-        "B05": 4,
-        "B06": 5,
-        "B07": 6,
-        "B08": 7,
-        "B8A": 8,
-        "B09": 9,
-        "B10": 10,
-        "B11": 11,
-        "B12": 12,
-    }
-
-    # define nodata values for ESA Sentinel-2 spectral bands (reflectance) and
-    # scene classification layer (SCL)
-    NODATA_REFLECTANCE: int = 64537
-    NODATA_SCL: int = 254
diff --git a/eodal/config/settings.py b/eodal/config/settings.py
index 17db6203..aa8f0454 100644
--- a/eodal/config/settings.py
+++ b/eodal/config/settings.py
@@ -63,6 +63,8 @@ class Settings(BaseSettings):
     # define CREODIAS username and password
     CREODIAS_USER: str = ""
     CREODIAS_PASSWORD: str = ""
+    # maximum number of records per request: 2000 (CREODIAS currently does not allow more)
+    CREODIAS_MAX_RECORDS: int = 2000
 
     # define Planet-API token
     PLANET_API_KEY: str = ""
diff --git a/eodal/config/stac_providers.py b/eodal/config/stac_providers.py
index 663075fd..c2907911 100644
--- a/eodal/config/stac_providers.py
+++ b/eodal/config/stac_providers.py
@@ -27,8 +27,8 @@ class AWS:
         """Amazon Web Services"""
 
         URL: str = "https://earth-search.aws.element84.com/v0"
-        S2L1C: str = "sentinel-s2-l1c"
-        S2L2A: str = "sentinel-s2-l2a"
+        S2Level_1C: str = "sentinel-s2-l1c"
+        S2Level_2A: str = "sentinel-s2-l2a"
 
         class Sentinel2:
             product_uri: str = "sentinel:product_id"
@@ -48,8 +48,8 @@ class MSPC:
         """Microsoft Planetary Computer"""
 
         URL: str = "https://planetarycomputer.microsoft.com/api/stac/v1"
-        S2L1C: str = "sentinel-2-l1c"
-        S2L2A: str = "sentinel-2-l2a"
+        S2Level_1C: str = "sentinel-2-l1c"
+        S2Level_2A: str = "sentinel-2-l2a"
         S1RTC: str = "sentinel-1-rtc" # radiometric and terrain corrected using PlanetDEM (IW mode)
         S1GRD: str = "sentinel-1-grd" # corrected to ground range using ellipsoid model WGS84
         
diff --git a/eodal/core/algorithms.py b/eodal/core/algorithms.py
index 7a4efe7c..09e4622f 100644
--- a/eodal/core/algorithms.py
+++ b/eodal/core/algorithms.py
@@ -19,6 +19,7 @@
 '''
 from __future__ import annotations
 
+import eodal
 import os
 import geopandas as gpd
 import uuid
@@ -80,7 +81,7 @@ def merge_datasets(
         to merge into a single raster
     :param out_file:
         name of the resulting raster dataset (optional). If None (default)
-        returns a new ``SatDataHandler`` instance otherwise writes the data
+        returns a new ``RasterCollection`` instance otherwise writes the data
         to disk as new raster dataset.
     :param target_crs:
         optional target spatial coordinate reference system in which the output
@@ -135,7 +136,7 @@ def merge_datasets(
     # when out_file was provided the merged data is written to file directly
     if out_file is not None:
         return
-    # otherwise, create new SatDataHandler instance from merged datasets
+    # otherwise, create new RasterCollection instance from merged datasets
     # add scene properties if available
     if sensor is None:
         raster = RasterCollection(scene_properties=scene_properties)
@@ -164,46 +165,43 @@ def merge_datasets(
         offset = band_attrs.get("offsets")
         if isinstance(offset, tuple):
             offset = offset[0]
-        description = band_attrs.get("descriptions")
-        if isinstance(description, tuple):
-            description = description[0]
         unit = band_attrs.get("units")
         if isinstance(unit, tuple):
             unit = unit[0]
+
+        # get band name and alias if provided
+        band_name = band_options.get('band_names', f'B{idx+1}')
+        if isinstance(band_name, list):
+            if len(band_name) == n_bands:
+                band_name = band_name[idx]
+            else:
+                band_name = f'B{idx+1}'
+        band_alias = band_options.get('band_aliases', f'B{idx+1}')
+        if isinstance(band_alias, list):
+            if len(band_alias) == n_bands:
+                band_alias = band_alias[idx]
+            else:
+                band_alias = f'B{idx+1}'
+
         raster.add_band(
             band_constructor=Band,
-            band_name=f"B{idx+1}",
+            band_name=band_name,
             values=out_ds[idx, :, :],
             geo_info=geo_info,
             is_tiled=is_tiled,
             scale=scale,
             offset=offset,
-            band_alias=description,
+            band_alias=band_alias,
             unit=unit,
         )
 
-    # clip raster collection if required to vector_features
+    # clip raster collection if required to vector_features to keep consistency
+    # of masks - this ensures that is_masked_array is set to True
     if vector_features is not None:
-        tmp_dir = Settings.TEMP_WORKING_DIR
-        fname_tmp = tmp_dir.joinpath(f"{uuid.uuid4()}.tif")
-        raster.to_rasterio(fname_tmp)
-        if sensor is None:
-            expr = "RasterCollection"
-        else:
-            expr = f"eodal.core.sensors.{sensor.lower()}.{sensor[0].upper() + sensor[1::]}()"
-        if band_options is None:
-            band_options = {}
-        raster = eval(
-            f"""{expr}.from_multi_band_raster(
-            fpath_raster=fname_tmp,
-            vector_features=vector_features,
-            full_bounding_box_only=False,
-            **band_options
-        )"""
-        )
-        # set scene properties if available
-        if scene_properties is not None:
-            raster.scene_properties = scene_properties
-        os.remove(fname_tmp)
+        raster.clip_bands(inplace=True, clipping_bounds=vector_features)
+
+    # set scene properties
+    if scene_properties is not None:
+        raster.scene_properties = scene_properties
 
     return raster
diff --git a/eodal/core/band.py b/eodal/core/band.py
index b6f8fb70..feda7962 100644
--- a/eodal/core/band.py
+++ b/eodal/core/band.py
@@ -55,6 +55,7 @@
 from shapely.geometry import box, MultiPolygon, Point, Polygon
 from typing import Any, Dict, List, Optional, Tuple, Union
 
+from eodal.config import get_settings
 from eodal.core.operators import Operator
 from eodal.core.utils.geometry import check_geometry_types, convert_3D_2D
 from eodal.core.utils.raster import get_raster_attributes, bounds_window
@@ -67,6 +68,8 @@
 )
 from eodal.utils.reprojection import reproject_raster_dataset, check_aoi_geoms
 
+Settings = get_settings()
+
 class BandOperator(Operator):
     """
     Band operator supporting basic algebraic operations on `Band` objects
@@ -636,7 +639,7 @@ def _get_pixel_geometries(
     @classmethod
     def from_rasterio(
         cls,
-        fpath_raster: Path,
+        fpath_raster: Path | Dict,
         band_idx: Optional[int] = 1,
         band_name_src: Optional[str] = "",
         band_name_dst: Optional[str] = "B1",
@@ -655,7 +658,11 @@ def from_rasterio(
             For Point features refer to the `read_pixels` method.
 
         :param fpath_raster:
-            file-path to the raster file from which to read a band
+            file-path to the raster file from which to read a band or
+
+            .. versionadd:: 0.2.0
+                can be also an `assets` dictionary returned from a STAC query
+
         :param band_idx:
             band index of the raster band to read (starting with 1). If not
             provided the first band will be always read. Ignored if
@@ -689,14 +696,21 @@ def from_rasterio(
         :returns:
             new ``Band`` instance from a ``rasterio`` dataset.
         """
+        _fpath_raster = deepcopy(fpath_raster)
+        # check if fpath_raster is STAC item or file system path
+        if isinstance(fpath_raster, dict):
+            if band_name_src is not None:
+                _fpath_raster = _fpath_raster[band_name_src]['href']
+            else:
+                _fpath_raster = _fpath_raster[list(_fpath_raster.keys())[band_idx]]['href']
+
         # check vector features if provided
         masking = False
         if vector_features is not None:
-
             masking = True
             gdf_aoi = check_aoi_geoms(
                 in_dataset=vector_features,
-                fname_raster=fpath_raster,
+                fname_raster=_fpath_raster,
                 full_bounding_box_only=full_bounding_box_only,
             )
             # check for third dimension (has_z) and flatten it to 2d
@@ -709,7 +723,7 @@ def from_rasterio(
             )
 
         # read data using rasterio
-        with rio.open(fpath_raster, "r") as src:
+        with rio.open(_fpath_raster, "r") as src:
 
             # parse image attributes
             attrs = get_raster_attributes(riods=src)
@@ -1118,7 +1132,7 @@ def copy(self):
 
     def clip(
         self,
-        clipping_bounds: Path | gpd.GeoDataFrame | Tuple[float,float,float,float] | Polygon  | MultiPolygon,
+        clipping_bounds: Path | gpd.GeoDataFrame | gpd.GeoSeries | Tuple[float,float,float,float] | Polygon  | MultiPolygon,
         full_bounding_box_only: Optional[bool] = False,
         inplace: Optional[bool] = False
     ):
@@ -1134,7 +1148,7 @@ def clip(
 
         :param clipping_bounds:
             spatial bounds to clip the Band to. Can be either a vector file, a shapely
-            `Polygon` or `MultiPolygon`, a `GeoDataFrame` or a coordinate tuple with
+            `Polygon` or `MultiPolygon`, a `GeoDataFrame`, `GeoSeries` or a coordinate tuple with
             (xmin, ymin, xmax, ymax).
             Vector files and `GeoDataFrame` are reprojected into the bands' coordinate
             system if required, while the coordinate tuple and shapely geometry **MUST**
@@ -1143,6 +1157,7 @@ def clip(
             if False (default), clips to the bounding box of the geometry and masks values
             outside the actual geometry boundaries. To obtain all values within the
             bounding box set to True.
+            .. versionadded:: 0.1.1
         :param inplace:
             if False (default) returns a copy of the ``Band`` instance
             with the changes applied. If True overwrites the values
@@ -1150,23 +1165,28 @@ def clip(
         :returns:
             clipped band instance.
         """
+        # prepare geometries
         if isinstance(clipping_bounds, Path):
             clipping_bounds = gpd.read_file(clipping_bounds)
-
-        # check inputs
+        if isinstance(clipping_bounds, gpd.GeoSeries):
+            clipping_bounds = gpd.GeoDataFrame(geometry=clipping_bounds)
         if isinstance(clipping_bounds, tuple):
             if len(clipping_bounds) != 4:
                 raise ValueError('Expected four coordinates (xmin, ymin, xmax, ymax)')
             xmin, ymin, xmax, ymax = clipping_bounds
             actual_geom = box(*clipping_bounds)
-        elif isinstance(clipping_bounds, gpd.GeoDataFrame):
-            # get the bounding box of the FIRST feature
+
+        # get bounding box
+        if isinstance(clipping_bounds, gpd.GeoDataFrame):
+            # reproject GeoDataFrame if necessary
             _clipping_bounds = clipping_bounds.copy()
-            _clipping_bounds = _clipping_bounds.bounds
-            xmin, ymin, xmax, ymax = _clipping_bounds.values[0]
+            _clipping_bounds.to_crs(epsg=self.geo_info.epsg, inplace=True)
+            # get the bounding box of the FIRST feature
+            _clipping_bounds_boundaries = _clipping_bounds.bounds
+            xmin, ymin, xmax, ymax = _clipping_bounds_boundaries.values[0]
             # the actual geometries are dissolved in case there is more than one record
             # and converted to a shapely object
-            actual_geom = clipping_bounds.dissolve().geometry.values[0]
+            actual_geom = _clipping_bounds.dissolve().geometry.values[0]
         elif (isinstance(clipping_bounds, Polygon) or isinstance(clipping_bounds, MultiPolygon)):
             xmin, ymin, xmax, ymax = clipping_bounds.bounds
             actual_geom = clipping_bounds
diff --git a/eodal/core/raster.py b/eodal/core/raster.py
index 46088d46..6e7239c5 100644
--- a/eodal/core/raster.py
+++ b/eodal/core/raster.py
@@ -95,6 +95,7 @@
 from typing import Optional
 from typing import Union
 
+from eodal.config import get_settings
 from eodal.core.band import Band
 from eodal.core.operators import Operator
 from eodal.core.spectral_indices import SpectralIndices
@@ -102,6 +103,8 @@
 from eodal.utils.decorators import check_band_names
 from eodal.utils.exceptions import BandNotFoundError
 
+Settings = get_settings()
+
 class SceneProperties(object):
     """
     A class for storing scene-relevant properties
@@ -285,6 +288,7 @@ def calc(
         cls.check_operator(operator=operator)
         # make a copy of a to avoid overwriting the original values
         _a = deepcopy(a)
+
         # if `other` is a Band object get its values
         if isinstance(other, Band):
             _other = deepcopy(other)
@@ -596,18 +600,23 @@ def get_band_alias(self, band_name: str) -> Union[Dict[str, str], None]:
 
     @staticmethod
     def _bands_from_selection(
-        fpath_raster: Path,
+        fpath_raster: Path | Dict,
         band_idxs: Optional[List[int]] = None,
         band_names_src: Optional[List[str]] = None,
         band_names_dst: Optional[List[str]] = None,
     ) -> Dict[str, Union[str, int]]:
         """
         Selects bands in a multi-band raster dataset based on a custom
-        selection of band indices or band names
+        selection of band indices or band names.
+
+        .. versionadd:: 0.2.0
+            works also with a dictionary of hrefs returned from a
+            STAC query
 
         :param fpath_raster:
             file-path to the raster file (technically spoken, this
-            can also have just a single band)
+            can also have just a single band) **or** when `USE_STAC` is True
+            the `assets` dictionary returned from a STAC call.
         :param band_idxs:
             optional list of band indices in the raster dataset
             to read. If not provided (default) all bands are loaded.
@@ -628,12 +637,16 @@ def _bands_from_selection(
         # check band selection
         band_names, band_count = None, None
         if band_idxs is None:
-            try:
-                with rio.open(fpath_raster, "r") as src:
-                    band_names = list(src.descriptions)
-                    band_count = src.count
-            except Exception as e:
-                raise IOError(f"Could not read {fpath_raster}: {e}")
+            if isinstance(fpath_raster, dict):
+                band_names = list(fpath_raster.keys())
+                band_count = len(band_names)
+            else:
+                try:
+                    with rio.open(fpath_raster, "r") as src:
+                        band_names = list(src.descriptions)
+                        band_count = src.count
+                except Exception as e:
+                    raise IOError(f"Could not read {fpath_raster}: {e}")
             # use default band names if not provided in data set
             if len(band_names) == 0:
                 band_names_src = [f"B{idx+1}" for idx in range(band_count)]
diff --git a/eodal/core/scene.py b/eodal/core/scene.py
index dd176cdc..0bfaa659 100644
--- a/eodal/core/scene.py
+++ b/eodal/core/scene.py
@@ -71,8 +71,7 @@ class methods (e.g., `RasterCollection.from_multi_band_raster`)
         self.collection = dict()
         self._frozen = True
         self._is_sorted = True
-
-        object.__setattr__(self, 'indexed_by_timestamps', indexed_by_timestamps)
+        self._indexed_by_timestamps = indexed_by_timestamps
 
         self._identifiers = []
         if scene_constructor is not None:
@@ -234,6 +233,10 @@ def __repr__(self) -> str:
                 f'# Scenes:    {len(self)}\nTimestamps:    {timestamps}\n' +  \
                 f'Scene Identifiers:    {", ".join(self.identifiers)}'
 
+    @property
+    def indexed_by_timestamps(self) -> bool:
+        return self._indexed_by_timestamps
+
     @staticmethod
     def _sort_keys(
         sort_direction: str,
diff --git a/eodal/core/sensors/planet_scope.py b/eodal/core/sensors/planet_scope.py
index ff6cb737..c56e3a82 100644
--- a/eodal/core/sensors/planet_scope.py
+++ b/eodal/core/sensors/planet_scope.py
@@ -124,8 +124,8 @@ def mask_non_clear_pixels(self, confidence_threshold: Optional[int] = 100):
     @classmethod
     def read_pixels(
         cls,
-        vector_features: Union[Path, gpd.GeoDataFrame],
         in_dir: Path,
+        vector_features: Union[Path, gpd.GeoDataFrame],
         band_selection: Optional[List[str]] = None,
         read_ql: Optional[bool] = True,
         apply_scaling: Optional[bool] = True,    
@@ -135,12 +135,12 @@ def read_pixels(
         vector geometry features read from a vector file (e.g., ESRI shapefile) or
         ``GeoDataFrame``.
 
-        :param point_features:
-            vector file (e.g., ESRI shapefile or geojson) or ``GeoDataFrame``
-            defining point locations for which to extract pixel values
         :param in_dir:
             Planet Scope SuperDove scene from which to extract
             pixel values at the provided point locations
+        :param point_features:
+            vector file (e.g., ESRI shapefile or geojson) or ``GeoDataFrame``
+            defining point locations for which to extract pixel values
         :param band_selection:
             list of bands to read. Per default all raster bands available are read.
         :param read_ql:
diff --git a/eodal/core/sensors/sentinel1.py b/eodal/core/sensors/sentinel1.py
index 8de12165..f8c22aeb 100644
--- a/eodal/core/sensors/sentinel1.py
+++ b/eodal/core/sensors/sentinel1.py
@@ -30,9 +30,8 @@
 import pandas as pd
 import planetary_computer
 
-from datetime import date
 from pathlib import Path
-from typing import Dict, List, Optional
+from typing import Any, Dict, List, Optional
 
 from eodal.config import get_settings
 from eodal.core.band import Band
@@ -156,8 +155,8 @@ def from_safe(
     @prepare_point_features
     def read_pixels_from_safe(
             cls,
+            in_dir: Dict[str, Any] | Path,
             vector_features: Path | gpd.GeoDataFrame,
-            in_dir: Path | Dict[str, str],
             polarizations: Optional[List[str]] = ['VV', 'VH']
         ) -> gpd.GeoDataFrame:
         """
@@ -169,12 +168,13 @@ def read_pixels_from_safe(
             A point is dimension-less, therefore, the raster grid cell (pixel) closest
             to the point is returned if the point lies within the raster.
 
+        :param in_dir:
+            Sentinel-1 scene in .SAFE structure from which to extract pixel values at
+            the provided point locations (GRD or RTC).Can be either a dictionary
+            item returned from STAC or a physical file-path
         :param vector_features:
             vector file (e.g., ESRI shapefile or geojson) or ``GeoDataFrame``
             defining point locations for which to extract pixel values
-        :param in_dir:
-            Sentinel-1 scene in .SAFE structure from which to extract
-            pixel values at the provided point locations (GRD or RTC)
         :param polarizations:
             selection of polarization to read. 'VV' and 'VH' by default.
         :returns:
diff --git a/eodal/core/sensors/sentinel2.py b/eodal/core/sensors/sentinel2.py
index 1318c9d1..23ad91dd 100644
--- a/eodal/core/sensors/sentinel2.py
+++ b/eodal/core/sensors/sentinel2.py
@@ -37,9 +37,9 @@
 from pathlib import Path
 from rasterio.mask import raster_geometry_mask
 from shapely.geometry import box
-from typing import Dict, Optional, List, Tuple, Union
+from typing import Any, Dict, Optional, List, Tuple, Union
 
-from eodal.core.band import Band, WavelengthInfo, GeoInfo
+from eodal.core.band import Band, WavelengthInfo
 from eodal.core.raster import RasterCollection, SceneProperties
 from eodal.utils.constants.sentinel2 import (
     band_resolution,
@@ -91,7 +91,7 @@ def is_blackfilled(self) -> bool:
         else:
             band_name = self.band_names[0]
             if self[band_name].is_masked_array:
-                return (self[band_name].values.date == 0).all()
+                return (self[band_name].values.data == 0).all()
             elif self[band_name].is_ndarray:
                 return (self[band_name].values == 0).all()
             elif self[band_name].is_zarr:
@@ -306,6 +306,12 @@ def from_safe(
                         vector_features_df = gpd.read_file(vector_features)
                     elif isinstance(vector_features, gpd.GeoDataFrame):
                         vector_features_df = vector_features.copy()
+                    elif isinstance(vector_features, gpd.GeoSeries):
+                        vector_features_df = gpd.GeoDataFrame(vector_features.copy())
+                    else:
+                        raise TypeError(
+                            'Geometry must be vector file, GeoSeries or GeoDataFrame'
+                        )
 
                     # drop Nones in geometry column
                     none_idx = vector_features_df[
@@ -351,6 +357,8 @@ def from_safe(
                     # remember to mask the feature after clipping the data
                     if not kwargs.get("full_bounding_box_only", False):
                         masking_after_read_required = True
+                    # update the vector_features entry
+                    kwargs.update({"vector_features": bounds_df})
 
         # determine platform (S2A or S2B)
         try:
@@ -392,7 +400,6 @@ def from_safe(
 
         # loop over bands and add them to the collection of bands
         sentinel2 = cls(scene_properties=scene_properties)
-        kwargs_orig = deepcopy(kwargs)
         for band_name in list(band_df_safe.band_name):
 
             # get entry from dataframe with file-path of band
@@ -422,17 +429,6 @@ def from_safe(
 
             # read band
             try:
-                if align_shapes:
-                    # if the current band already has the lowest resolution
-                    # reading the mask directly is possible. Otherwise we use
-                    # the bounding of the coarsest resolution and apply the masking
-                    # later
-                    if band_safe.band_resolution.values != lowest_resolution:
-                        kwargs.update({"vector_features": bounds_df})
-                    else:
-                        kwargs.update(
-                            {"vector_features": kwargs_orig.get("vector_features")}
-                        )
                 sentinel2.add_band(
                     Band.from_rasterio,
                     fpath_raster=band_fpath,
@@ -476,8 +472,8 @@ def from_safe(
     @prepare_point_features
     def read_pixels_from_safe(
         cls,
+        in_dir: Dict[str, Any] | Path,
         vector_features: Union[Path, gpd.GeoDataFrame],
-        in_dir: Path,
         band_selection: Optional[List[str]] = None,
         read_scl: Optional[bool] = True,
         apply_scaling: Optional[bool] = True,
@@ -505,12 +501,13 @@ def read_pixels_from_safe(
             to do spatial resampling! The underlying ``rasterio.sample`` function always
             snaps to the closest pixel in the current spectral band.
 
+        :param in_dir:
+            Sentinel-2 scene in .SAFE structure from which to extract
+            pixel values at the provided point locations. Can be either a dictionary
+            item returned from STAC or a physical file-path
         :param vector_features:
             vector file (e.g., ESRI shapefile or geojson) or ``GeoDataFrame``
             defining point locations for which to extract pixel values
-        :param in_dir:
-            Sentinel-2 scene in .SAFE structure from which to extract
-            pixel values at the provided point locations
         :param band_selection:
             list of bands to read. Per default all raster bands available are read.
         :param read_scl:
@@ -641,23 +638,28 @@ def plot_scl(self, colormap: Optional[str] = "") -> Figure:
 
     def mask_clouds_and_shadows(
         self,
-        bands_to_mask: List[str],
+        bands_to_mask: Optional[List[str]] = None,
         cloud_classes: Optional[List[int]] = [1, 2, 3, 7, 8, 9, 10, 11],
+        mask_band: Optional[str] = 'SCL',
         **kwargs,
-    ):
+    ) -> Sentinel2:
         """
         A Wrapper around the inherited ``mask`` method to mask clouds,
-        shadows, water and snow based on the SCL band. Works therefore on L2A data,
-        only.
+        shadows, water and snow based on (by default) the SCL band.
+        Works therefore on L2A data, only.
 
-        Masks the SCL classes 1, 2, 3, 7, 8, 9, 10, 11.
+        NOTE:
+            Since the `mask_band` can be set to *any* `Band` it is also
+            possible to use a different cloud/shadow etc. mask, e.g., from
+            a custom classifier.
 
-        If another class selection is desired consider using the mask function
-        from `eodal.utils.io` directly or change the default
-        ``cloud_classes``.
+        NOTE:
+            You might also use the mask function from `eodal.core.raster.RasterCollection`
+            directly.
 
         :param bands_to_mask:
-            list of bands on which to apply the SCL mask
+            list of bands on which to apply the SCL mask. If not specified all bands
+            are masked.
         :param cloud_classes:
             list of SCL values to be considered as clouds/shadows and snow.
             By default, all three cloud classes and cloud shadows are considered
@@ -668,7 +670,12 @@ def mask_clouds_and_shadows(
             depending on `inplace` (passed in the kwargs) a new `Sentinel2` instance
             or None
         """
-        mask_band = "SCL"
+        if bands_to_mask is None:
+            bands_to_mask = self.band_names
+        # the mask band should never be masked as otherwise the SCL functions
+        # might not work as expected
+        if mask_band in bands_to_mask:
+            bands_to_mask.remove('SCL')
         try:
             return self.mask(
                 mask=mask_band,
@@ -692,12 +699,15 @@ class of the scene classification layer. Works for data in
             class occurences.
         """
         # check if SCL is available
-        if not "SCL" in self.band_names:
+        if not 'scl' in self.band_names and not 'SCL' in self.band_names:
             raise BandNotFoundError(
                 "Could not find scene classification layer. Is scene L2A?"
             )
 
-        scl = self.get_band("SCL")
+        try:
+            scl = self.get_band('SCL')
+        except BandNotFoundError:
+            scl = self.get_band('scl')
         # if the scl array is a masked array consider only those pixels
         # not masked out
         if scl.is_masked_array:
diff --git a/eodal/core/utils/geometry.py b/eodal/core/utils/geometry.py
index a47a77c9..788831de 100644
--- a/eodal/core/utils/geometry.py
+++ b/eodal/core/utils/geometry.py
@@ -34,11 +34,15 @@ def read_geometries(in_dataset: Union[Path, gpd.GeoDataFrame]) -> gpd.GeoDataFra
     Returns a geodataframe containing vector features
 
     :param in_dataset:
-        path-like object or ``GeoDataFrame``
+        path-like object ``GeoSeries`` or ``GeoDataFrame``
+    :returns:
+        ``GeoDataFrame`` representation of vector features
     """
 
     if isinstance(in_dataset, gpd.GeoDataFrame):
-        return in_dataset
+        return in_dataset.copy()
+    elif isinstance(in_dataset, gpd.GeoSeries):
+        return gpd.GeoDataFrame(geometry=in_dataset.copy())
     elif isinstance(in_dataset, Path):
         try:
             return gpd.read_file(in_dataset)
diff --git a/eodal/core/utils/raster.py b/eodal/core/utils/raster.py
index 108fc0a2..7db70388 100644
--- a/eodal/core/utils/raster.py
+++ b/eodal/core/utils/raster.py
@@ -71,6 +71,8 @@ def spatial_to_image_coordinates(x: float, y: float, affine: Affine, op=math.flo
     """
     Convert spatial x and y coordinat to image coordinates (row + column)
 
+    .. versionadded:: 0.1.1
+
     taken from:
         `rasterstats` package under BSD 3-Clause "New" or "Revised" License
         Copyright (c) 2013 Matthew Perry
@@ -85,6 +87,8 @@ def bounds_window(bounds: Tuple[float,float,float,float], affine: Affine
     """
     Create a full cover rasterio-style window
 
+    .. versionadded:: 0.1.1
+
     taken from:
         `rasterstats` package under BSD 3-Clause "New" or "Revised" License
         Copyright (c) 2013 Matthew Perry
diff --git a/eodal/mapper/feature.py b/eodal/mapper/feature.py
index 3371d788..b654e6ad 100644
--- a/eodal/mapper/feature.py
+++ b/eodal/mapper/feature.py
@@ -1,6 +1,8 @@
 '''
 Module defining geographic features for mapping.
 
+.. versionadded:: 0.1.1
+
 Copyright (C) 2022 Lukas Valentin Graf
 
 This program is free software: you can redistribute it and/or modify
@@ -65,7 +67,7 @@ def __init__(
         if name == '':
             raise ValueError(f'Empty feature names are not allowed')
         if type(geometry) not in allowed_geom_types:
-            raise ValueError(f'geometry must of type {",".join(allowed_geom_types)}')
+            raise ValueError(f'geometry must of type {allowed_geom_types}')
         if type(epsg) != int or epsg <= 0:
             raise ValueError('EPSG code must be a positive integer value')
         if not isinstance(attributes, pd.Series) and not isinstance(attributes, dict):
@@ -144,7 +146,6 @@ def from_dict(cls, dictionary: Dict[str, Any]):
             )
         except WKTReadingError as e:
             raise ValueError(f'Invalid Geometry: {e}')
-        
 
     def to_epsg(self, epsg: int):
         """
@@ -171,8 +172,7 @@ def to_geoseries(self) -> gpd.GeoSeries:
         gds = gpd.GeoSeries([self.geometry], crs=f'EPSG:{self.epsg}')
         # add attributes from Feature
         gds.attrs = self.attributes
-        # set name of Feature to GeoSeries
-        gds.name = self.name
+        gds.name = 'geometry'
         return gds
 
     def to_dict(self) -> Dict[str, Any]:
@@ -188,42 +188,3 @@ def to_dict(self) -> Dict[str, Any]:
         feature_dict['geometry'] = self.geometry.wkt
         feature_dict['attributes'] = self.attributes
         return feature_dict
-
-if __name__ == '__main__':
-
-    # working constructor calls
-    geom = Point([49,11])
-    epsg = 4326
-    name = 'Test Point'
-    feature = Feature(name, geom, epsg)
-
-    assert feature.geometry == geom, 'geometry differs'
-    assert feature.epsg == epsg, 'EPSG code differs'
-    assert feature.name == name, 'name differs'
-    assert feature.attributes == {}, 'attributes must be empty'
-
-    attributes = {'key': 'value'}
-    feature = Feature(name, geom, epsg, attributes)
-    assert feature.attributes == attributes, 'attributes differ'
-
-    attributes = pd.Series({'key1': 'value1', 'key2': 'value2'})
-    feature = Feature(name, geom, epsg, attributes)
-    assert feature.attributes == attributes.to_dict(), 'attributes differ'
-
-    gds = feature.to_geoseries()
-    assert gds.name == feature.name, 'name differs'
-    assert gds.crs.to_epsg() == feature.epsg, 'EPSG differs'
-    assert gds.attrs == feature.attributes, 'attributes differ'
-
-    # from_geoseries class method
-    gds.attrs = {}
-    feature = Feature.from_geoseries(gds)
-    assert gds.name == feature.name, 'name differs'
-    assert gds.crs.to_epsg() == feature.epsg, 'EPSG differs'
-    assert gds.attrs == feature.attributes, 'attributes differ'
-
-    # project into another spatial reference system
-    feature_utm = feature.to_epsg(epsg=32632)
-    assert feature_utm.epsg == 32632, 'projection had no effect'
-    assert feature_utm.name == feature.name, 'name got lost'
-    assert feature_utm.attributes == feature.attributes, 'attributes got lost'
diff --git a/eodal/mapper/filter.py b/eodal/mapper/filter.py
index 77c2ce17..c0e23b21 100644
--- a/eodal/mapper/filter.py
+++ b/eodal/mapper/filter.py
@@ -19,42 +19,77 @@
 
 from __future__ import annotations
 
+from typing import Any
+
+operators = ['<', '<=', '==', '!=', '>', '>=']
+
 class Filter:
     """
-    The generic filter class.
+    The generic filter class. A filter is used to query data catalogs.
+
+    Each filter follows the following structure:
+    
+    <entity> <operator> <value>
+
+    For instance, `cloudy_pixel_percentage < 10`, where
+
+        * `cloudy_pixel_percentage` is the entity to filter
+        * `<` is the operator
+        * `10` is the value
 
     :attrib entity:
         metadata entity to use for filtering
-    :attrib condition:
-        condition that must be met to keep a metadata item in the
-        selection
+    :attrib operator:
+            comparison operator to use, e.g., "gt" for "greater than" (>)
+    :attrib value:
+        value on the right-hand side of the filter expression
     """
-    def __init__(self, entity: str, condition: str):
+    def __init__(self, entity: str, operator: str, value: Any):
         """
         Constructor method
 
         :param entity:
             metadata entity to use for filtering
-        :param condition:
-            condition that must be met to keep a metadata item in the
-            selection
+        :param operator:
+            comparison operator to use, e.g., ">" for "greater than" (value)
+        :param value:
+            value on the right-hand side of the filter expression
         """
+        # check inputs
+        if not isinstance(entity, str):
+            raise TypeError('Entity argument must be a string')
+        if entity == '':
+            raise ValueError('Entity argument must not be an empty string')
+        if not isinstance(operator, str):
+            raise TypeError('Operator argument must be a string')
+        if operator not in operators:
+            raise ValueError('Operator must be one of: ' + ','.join(operators))
+        if value is None:
+            raise ValueError('Value cannot be None')
+
         self._entity = entity
-        self._condition = condition
+        self._operator = operator
+        self._value = value
 
     def __repr__(self) -> str:
-        return f'Filter by {self.entity} {self.condition}'
+        return self.expression
 
     @property
     def entity(self) -> str:
+        """metadata entity used to filter"""
         return self._entity
 
     @property
-    def condition(self) -> str:
-        return self._condition
+    def operator(self) -> str:
+        """filter operator"""
+        return self._operator
 
-if __name__ == '__main__':
-    
-    cc_filter = Filter(entity='cloudy_pixel_percentage', condition='<30')
-    cc_filter
-    
\ No newline at end of file
+    @property
+    def value(self) -> Any:
+        """right-side value of the filter"""
+        return self._value
+
+    @property
+    def expression(self) -> str:
+        """returns the filter expression as string"""
+        return f'{self.entity} {self.operator} {self.value}'
diff --git a/eodal/mapper/mapper.py b/eodal/mapper/mapper.py
index 1d716db3..fcc2565d 100644
--- a/eodal/mapper/mapper.py
+++ b/eodal/mapper/mapper.py
@@ -1,19 +1,56 @@
 '''
-Created on Dec 12, 2022
+The EOdal `Mapper` class allows to extract and handle EO data in space and time
+and bring the data into Analysis-Ready-Format (ARD).
 
-@author: graflu
+.. versionadded:: 0.2.0
+
+Copyright (C) 2022 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/>.
 '''
 
 from __future__ import annotations
 
+import cv2
+import eodal
+import geopandas as gpd
+import getpass
+import numpy as np
+import pandas as pd
+import warnings
+import uuid
 import yaml
 
+from copy import deepcopy
 from datetime import datetime
 from pathlib import Path
-from typing import List, Optional
+from sqlalchemy.exc import DatabaseError
+from shapely.geometry import box
+from typing import Any, Callable, Dict, List, Optional
 
+from eodal.config import get_settings
+from eodal.core.algorithms import merge_datasets
+from eodal.core.raster import RasterCollection
+from eodal.core.scene import SceneCollection
 from eodal.mapper.feature import Feature
 from eodal.mapper.filter import Filter
+from eodal.metadata.database.querying import find_raw_data_by_bbox
+from eodal.metadata.utils import reconstruct_path
+from eodal.utils.exceptions import STACError
+
+settings = get_settings()
+logger = settings.logger
 
 class MapperConfigs:
     """
@@ -47,8 +84,10 @@ def __init__(
         default class constructor
 
         :param collection:
-            name of the collection (<platform>-<sensor>-<processing level>) to use.
-            E.g., "Sentinel2-MSI-L2A
+            name of the collection (<platform>-<sensor>) to use.
+            E.g., "sentinel2-msi". <sensor> is optional and can be omitted
+            if a platform does not carry more than a single sensor. I.e.,
+            one could also pass "sentinel2" instead.
         :param feature:
             geographic feature(s) for which to extract data from collection
         :param time_start:
@@ -63,18 +102,19 @@ def __init__(
         # check inputs
         if not isinstance(collection, str):
             raise TypeError('Collection must be a string')
-        if len(collection) < 6:
-            raise ValueError('Collections must have at least six characters')
-        if collection.count('-') != 2:
+
+        if len(collection) < 3:
+            raise ValueError('Collections must have at least 3 characters')
+        if collection.count('-') > 2:
             raise ValueError(
-                f'Collections must obey the format <platform>-<sensor>-<prcoessing_level>'
+                f'Collections must obey the format <platform>-<sensor> where <sesor> is optional'
             )
         if not isinstance(feature, Feature):
             raise TypeError('Expected a Feature object')
         if not isinstance(time_start, datetime) and not isinstance(time_end, datetime):
             raise TypeError('Expected datetime objects')
         if metadata_filters is not None:
-            if not [isinstance(x, Filter) for x in metadata_filters].all():
+            if not np.array([isinstance(x, Filter) for x in metadata_filters]).all():
                 raise TypeError('All filters must be instances of the Filter class')
 
         self._collection = collection
@@ -83,6 +123,11 @@ def __init__(
         self._time_end = time_end
         self._metadata_filters = metadata_filters
 
+    def __repr__(self) -> str:
+        return f'EOdal MapperConfig\n------------------\nCollection: {self.collection}' + \
+            f'\nTime Range: {self.time_start} - {self.time_end}\nFeature:\n' + \
+            f'{self.feature.__repr__()}\nMetadata Filters: {str(self.metadata_filters)}'
+
     @property
     def collection(self) -> str:
         return self._collection
@@ -93,11 +138,10 @@ def platform(self) -> str:
 
     @property
     def sensor(self) -> str:
-        return self._collection.split('-')[1]
-
-    @property
-    def processing_level(self) -> str:
-        return self._collection.split('-')[-1]
+        try:
+            return self._collection.split('-')[1]
+        except IndexError:
+            return ''
 
     @property
     def feature(self) -> Feature:
@@ -116,9 +160,14 @@ def metadata_filters(self) -> List[Filter] | None:
         return self._metadata_filters
 
     @classmethod
-    def from_yaml(cls, fpath: str | Path):
+    def from_yaml(cls, fpath: str | Path) -> cls:
         """
         Load mapping configurations from YAML file
+
+        :param fpath:
+            file-path to yaml with Mapper configurations
+        :returns:
+            new `MapperConfigs` instance
         """
         with open(fpath, 'r') as f:
             try:
@@ -129,6 +178,12 @@ def from_yaml(cls, fpath: str | Path):
         if 'feature' not in yaml_content.keys():
             raise ValueError('"feature" attribute is required"')
         feature_yaml = yaml_content['feature']
+        # reconstruct the filter objects
+        filter_list = []
+        if 'metadata_filters' in yaml_content.keys():
+            filters_yaml = yaml_content['metadata_filters']
+            for filter_yaml in filters_yaml:
+                filter_list.append(Filter(*filter_yaml.split()))
         try:
             feature = Feature.from_dict(feature_yaml)
             return cls(
@@ -136,14 +191,14 @@ def from_yaml(cls, fpath: str | Path):
                 feature=feature,
                 time_start=yaml_content['time_start'],
                 time_end=yaml_content['time_end'],
-                metadata_filters=yaml_content['metadata_filters']
+                metadata_filters=filter_list
             )
         except KeyError as e:
-            raise ValueError(f'IMissing keys in yaml file: {e}')
+            raise ValueError(f'Missing keys in yaml file: {e}')
 
     def to_yaml(self, fpath: str | Path) -> None:
         """
-        save MapperConfig to YAML file (*.yml)
+        save MapperConfig and some meta-information to YAML file (*.yml)
 
         :param fpath:
             file-path where saving the Feature instance to
@@ -153,7 +208,12 @@ def to_yaml(self, fpath: str | Path) -> None:
         mapper_configs_dict['feature'] = self.feature.to_dict()
         mapper_configs_dict['time_start'] = self.time_start
         mapper_configs_dict['time_end'] = self.time_end
-        mapper_configs_dict['metadata_filters'] = self.metadata_filters
+        mapper_configs_dict['metadata_filters'] = [x.__repr__() for x in self.metadata_filters]
+        # add some meta-information about creation time and current user
+        mapper_configs_dict['created_at'] = datetime.now()
+        mapper_configs_dict['created_by'] = getpass.getuser()
+        mapper_configs_dict['eodal_version'] = eodal.__version__
+
         with open(fpath, 'w+') as f:
             yaml.dump(mapper_configs_dict, f, allow_unicode=True)
 
@@ -164,35 +224,453 @@ class Mapper:
 
     The mapper class takes over searching for EO scenes, merging them and
     filling eventually occurring black-fill (no-data regions).
+
+    :attrib mapper_configs:
+        `MapperConfig` instance defining search criteria
+    :attrib data:
+        loaded scene data as EOdal `SceneCollection`
+    :attrib metadata:
+        corresponding scene metadata as `GeoDataFrame`.
     """
 
-    def __init__(self, mapper_configs: MapperConfigs):
-        pass
-    
-if __name__ == '__main__':
-
-    from shapely.geometry import Point
-
-    collection = 'sentinel2-msi-l2a'
-    geom = Point([49,11])
-    epsg = 4326
-    name = 'Test Point'
-    attributes = {'this': 'is a test', 'a': 123}
-    feature = Feature(name, geom, epsg, attributes)
-    time_start = datetime(2022,12,1)
-    time_end = datetime.now()
-
-    mapper_configs = MapperConfigs(collection, feature, time_start, time_end)
-    fpath = '/mnt/ides/Lukas/test.yml'
-    mapper_configs.to_yaml(fpath)
-
-    mapper_configs_rl = MapperConfigs.from_yaml(fpath)
-    assert mapper_configs.feature.name == mapper_configs_rl.feature.name, 'wrong feature name'
-    assert mapper_configs.feature.epsg == mapper_configs_rl.feature.epsg, 'wrong feature EPSG'
-    assert mapper_configs.feature.geometry == mapper_configs_rl.feature.geometry, 'wrong feature geometry'
-    assert set(mapper_configs.feature.attributes) == set(mapper_configs_rl.feature.attributes), \
-        'wrong feature attributes'
-    assert mapper_configs.time_start == mapper_configs_rl.time_start, 'wrong start time'
-    assert mapper_configs.time_end == mapper_configs_rl.time_end, 'wrong end time'
-    assert mapper_configs.collection == mapper_configs_rl.collection, 'wrong collection'
-        
\ No newline at end of file
+    def __init__(
+        self,
+        mapper_configs: MapperConfigs,
+        time_column: Optional[str] = 'sensing_time'
+    ):
+        """
+        Class constructor
+
+        :param mapper_configs:
+            `MapperConfig` instance defining search criteria
+        :param time_column:
+            name of the metadata column denoting the time stamps of the scenes.
+            `sensing_time` by default.
+        """
+        if not isinstance(mapper_configs, MapperConfigs):
+            raise TypeError(f'Expected a MapperConfigs instance')
+        self._mapper_configs = mapper_configs
+        self._time_column = time_column
+        self._metadata = None
+        self._sensor = self.mapper_configs.collection.split('-')[0]
+        self._data = None
+        self._geoms_are_points = False
+
+    def __repr__(self) -> str:
+        return f'EOdal Mapper\n============\n{self.mapper_configs.__repr__()}'
+
+    @property
+    def data(self) -> None | SceneCollection | gpd.GeoDataFrame:
+        """
+        SceneCollection with scenes found or GeoDataFrame in case
+        single pixel values are extracted
+        """
+        return self._data
+
+    @data.setter
+    def data(self, scoll: Optional[SceneCollection] = None):
+        if scoll is not None:
+            if not isinstance(scoll, SceneCollection) and \
+            not isinstance(scoll, gpd.GeoDataFrame):
+                raise TypeError('Expected a EOdal SceneCollection or GeoDataFrame')
+        self._data = scoll
+
+    @property
+    def mapper_configs(self) -> MapperConfigs:
+        """mapper configurations"""
+        return self._mapper_configs
+
+    @property
+    def metadata(self) -> None | gpd.GeoDataFrame:
+        """scene metadata found"""
+        return self._metadata
+
+    @metadata.setter
+    def metadata(self, values: Optional[gpd.GeoDataFrame] = None):
+        """set scene metadata"""
+        if values is not None:
+            if not isinstance(values, gpd.GeoDataFrame):
+                raise TypeError('Expected a GeoDataFrame')
+        self._metadata = deepcopy(values)
+
+    @property
+    def sensor(self) -> str | None:
+        return self._sensor
+
+    @property
+    def time_column(self) -> str:
+        return self._time_column
+
+    def query_scenes(self) -> None:
+        """
+        Query available scenes for the current `MapperConfigs` and loads
+        into the `observations` attribute.
+
+        Depending on the settings of EOdal, this method either makes a
+        query to a STAC resource or to a PostgreSQL/PostGIS database.
+
+        NOTE:
+            This method only queries a metadata catalog without reading data.
+        """
+        # check for point geometries
+        self._geoms_are_points = self.mapper_configs.feature.geometry.geom_type in \
+            ['Point', 'MultiPoint']
+
+        # determine bounding box of the feature using
+        # its representation in geographic coordinates (WGS84, EPSG: 4326)
+        feature_wgs84 = self.mapper_configs.feature.to_epsg(4326)
+        bbox = box(*feature_wgs84.geometry.bounds)
+
+        # determine platform
+        platform = self.mapper_configs.platform
+
+        # put kwargs together
+        kwargs = {
+            'platform': platform,
+            'time_start': self.mapper_configs.time_start,
+            'time_end': self.mapper_configs.time_end,
+            'bounding_box': bbox,
+            'metadata_filters': self.mapper_configs.metadata_filters
+        }
+
+        # query the metadata catalog (STAC or database depending on settings)
+        if settings.USE_STAC:
+            try:
+                exec(f'from eodal.metadata.stac import {platform}')
+                scenes_df = eval(f'{platform}(**kwargs)')
+            except Exception as e:
+                raise STACError(f"Querying STAC catalog failed: {e}")
+        else:
+            try:
+                scenes_df = find_raw_data_by_bbox(**kwargs)
+            except Exception as e:
+                raise DatabaseError(f"Querying metadata DB failed: {e}")
+
+        # populate the metadata attribute
+        self.metadata = scenes_df
+
+    def _process_scene(
+        self,
+        item: pd.Series,
+        scene_constructor: Callable[...,RasterCollection],
+        scene_constructor_kwargs: Dict[str, Any],
+        scene_modifier: Callable[...,RasterCollection],
+        scene_modifier_kwargs: Dict[str, Any],
+        reprojection_method: int
+    ) -> RasterCollection:
+        """
+        Pre-process a scene so it can be added to a SceneCollection.
+        This includes applying user-defined reading and modification functions.
+        In addition, the `Scene` is projected into the target spatial reference
+        system by using the reference system the majority of the scenes in a
+        collection has in common.
+
+        IMPORTANT:
+            The reprojection step into the target spatial reference system (if
+            scene is not projected in it already) is **always** done!
+
+        :param item:
+            metadata item of the scene including its file-path or URL
+        :param scene_constructor:
+            Callable used to read the scenes found into `RasterCollection` fulfilling
+            the `is_scene` criterion (i.e., a time stamp is available).
+        :param scene_constructor_kwargs:
+            keyword-arguments to pass to `scene_constructor`. `fpath_raster`
+            and `vector_features` are filled in by the `Mapper` instance automatically,
+        :param scene_modifier:
+            optional Callable modifying a `RasterCollection` or returning a new
+            `RasterCollection`.
+        :param scene_modifier_kwargs:
+            keyword arguments for `scene_modifier`
+        :returns:
+            `Scene` with all pre-processing steps applied.
+        """
+        scene_constructor_kwargs.update({
+            'vector_features': self.mapper_configs.feature.to_geoseries()
+        })
+        try:
+            # call scene constructor. The file-path (or URL) goes first
+            scene = scene_constructor.__call__(item.real_path, **scene_constructor_kwargs)
+            scene.scene_properties.sensing_time = item[self.time_column]
+        except Exception as e:
+            raise ValueError(f'Could not load scene:  {e}')
+
+        # apply scene modifier callable if available
+        if scene_modifier is not None:
+            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
+        )
+
+        return scene
+
+    def _load_scenes_collection(
+        self,
+        reprojection_method: Optional[int] = cv2.INTER_NEAREST_EXACT,
+        scene_constructor: Optional[Callable[...,RasterCollection]] = RasterCollection.from_multi_band_raster,
+        scene_constructor_kwargs: Optional[Dict[str, Any]] = {},
+        scene_modifier: Optional[Callable[...,RasterCollection]] = None,
+        scene_modifier_kwargs: Optional[Dict[str, Any]] = {}
+    ) -> None:
+        """
+        Auxiliary method to handle EOdal scenes and store them into a SceneCollection.
+
+        This method is called when the geometries used for calling the `Mapper`
+        instance are of type `Polygon` or `MultiPolygon`.
+
+        Mosaicing operations are handled on the fly so calling programs will
+        always receive analysis-ready data.
+
+        :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
+            applied to all scenes found in the metadata query call.
+            By default the standard class-method call `~RasterCollection.from_multi_band_raster`
+            is used. It can be replaced, however, with a custom-written callable that
+            can be of any design except that it **MUST** accept a keyword argument
+            `fpath_raster` used for reading the Scene data and `vector_features` for
+            cropping the data to the spatial extent of the Mapper instance.
+        :param scene_constructor_kwargs:
+            optional keyword-arguments to pass to `scene_constructor`. `fpath_raster`
+            and `vector_features` are filled in by the `Mapper` instance automatically,
+            i.e., any custom values passed will be overwritten.
+        :param scene_modifier:
+            optional Callable modifying a `RasterCollection` or returning a new
+            `RasterCollection`. The Callable is applied to all scenes in the
+            `SceneCollection` when loaded by the `Mapper`. Can be used, e.g.,
+            to calculate spectral indices on the fly or for applying masks.
+        :param scene_modifier_kwargs:
+            optional keyword arguments for `scene_modifier` (if any).
+        """
+        # open a SceneCollection for storing the data
+        scoll = SceneCollection()
+        logger.info(f'Starting extraction of {self.sensor} scenes')
+        # filter out datasets where mosaicing is necessary (time stamp is the same)
+        self.metadata['_duplicated'] = self.metadata[self.time_column].duplicated(keep=False)
+        # datasets where the 'duplicated' entry is False are truely unqiue
+        _metadata_unique = self.metadata[~self.metadata._duplicated].copy()
+        _metadata_nonunique = self.metadata[self.metadata._duplicated].copy()
+
+        # mosaic the non-unique datasets first
+        if not _metadata_nonunique.empty:
+            _metadata_nonunique.sort_values(by=self.time_column, inplace=True)
+            unique_time_stamps = _metadata_nonunique[self.time_column].unique()
+            # loop over unique time stamps. In the end there should be a single
+            # scene per time stamp
+            update_scene_properties_list = []
+            for unique_time_stamp in unique_time_stamps:
+                scenes = _metadata_nonunique[
+                    _metadata_nonunique[self.time_column].dt.strftime('%Y-%m-%d %H:%M') == \
+                    pd.to_datetime(unique_time_stamp).strftime('%Y-%m-%d %H:%M')[0:16]
+                ]
+                # read the datasets one by one, save them into a temporary directory
+                # and merge them using rasterio
+                dataset_list = []
+                scene_properties_list = []
+                for _, item in scenes.iterrows():
+                    _scene = self._process_scene(
+                        item=item,
+                        scene_constructor=scene_constructor,
+                        scene_constructor_kwargs=scene_constructor_kwargs,
+                        reprojection_method=reprojection_method,
+                        scene_modifier=scene_modifier,
+                        scene_modifier_kwargs=scene_modifier_kwargs
+                    )
+                    fname_scene = settings.TEMP_WORKING_DIR.joinpath(f"{uuid.uuid4()}.tif")
+                    _scene.to_rasterio(fname_scene)
+                    dataset_list.append(fname_scene)
+                    scene_properties_list.append(_scene.scene_properties)
+                # merge datasets using rasterio and read results back into a scene
+                band_options={'band_names': _scene.band_names, 'band_aliases': _scene.band_aliases}
+                scene = merge_datasets(
+                    datasets=dataset_list,
+                    target_crs=self.metadata.target_epsg.unique()[0],
+                    vector_features=self.mapper_configs.feature.to_geoseries(),
+                    sensor=self.sensor,
+                    band_options=band_options
+                )
+                # handle scene properties. They need to be merged as well
+                merged_scene_properties = scene_properties_list[0]
+                for other_scene_properties in scene_properties_list[1::]:
+                    scene_props_keys = list(merged_scene_properties.__dict__.keys())
+                    for scene_prop in scene_props_keys:
+                        first_val = eval(f'merged_scene_properties.{scene_prop}')
+                        this_val =  eval(f'other_scene_properties.{scene_prop}')
+                        if first_val != this_val:
+                            # only string values can be merged (connected by '&&')
+                            if isinstance(first_val, str):
+                                new_val = first_val + '&&' + this_val
+                                exec(f'merged_scene_properties.{scene_prop} = new_val')
+
+                scene.scene_properties = merged_scene_properties
+                # because ESA has some mess with the naming of their file names and
+                # metadata, there might be duplicated seems not detected by the mapper
+                # since the time stamps slightly differ (few seconds in some cases) but
+                # their IDs are the same
+                try:
+                    scoll.add_scene(scene)
+                except KeyError:
+                    logger.warn(
+                        f'Scene with ID {merged_scene_properties.product_uri} ' + \
+                        f'already added to SceneCollection - continue'
+                    )
+                    continue
+                update_scene_properties_list.append(merged_scene_properties)
+
+            # update the metadata entries to avoid mis-matches in the number of
+            # scenes and their URIs
+            self.metadata.drop_duplicates(subset=[self.time_column], keep='first', inplace=True)
+            for updated_scene_properties in update_scene_properties_list:
+                # use the time stamp for finding the correct metadata records. There might be
+                # some disagreement in the milliseconds because of different precision levels
+                # therefore, an offset of less than 1 second is tolerated
+                idx = self.metadata[
+                    abs(self.metadata[self.time_column] - pd.to_datetime(updated_scene_properties.acquisition_time)
+                    ) < pd.Timedelta(60, unit='minutes')
+                ].index
+                for scene_property in updated_scene_properties.__dict__:
+                    self.metadata.loc[idx, scene_property] = \
+                        eval(f'updated_scene_properties.{scene_property}')
+
+        # then add those scenes that are unique, i.e., no mosaicing is required
+        if not _metadata_unique.empty:
+            for _, item in _metadata_unique.iterrows():
+                scene = self._process_scene(
+                    item=item,
+                    scene_constructor=scene_constructor,
+                    scene_constructor_kwargs=scene_constructor_kwargs,
+                    scene_modifier=scene_modifier,
+                    scene_modifier_kwargs=scene_modifier_kwargs,
+                    reprojection_method=reprojection_method
+                )
+                # because ESA has some mess with the naming of their file names and
+                # metadata, there might be duplicated seems not detected by the mapper
+                # since the time stamps slightly differ (few seconds in some cases) but
+                # their IDs are the same
+                try:
+                    scoll.add_scene(scene)
+                except KeyError:
+                    logger.warn(
+                        f'Scene with ID {merged_scene_properties.product_uri} ' + \
+                        f'already added to SceneCollection - continue'
+                    )
+                    continue
+
+        # sort scenes by their timestamps and save as data attribute
+        # to mapper instance
+        self.data = scoll.sort()
+        logger.info(f'Finished extraction of {self.sensor} scenes')
+
+    def _load_pixels(
+        self,
+        pixel_reader: Optional[Callable[...,gpd.GeoDataFrame]] = RasterCollection.read_pixels,
+        pixel_reader_kwargs: Optional[Dict[str, Any]] = {}
+    ) -> None:
+        """
+        Load pixel values from 0:N scenes into a `GeoDataFrame`
+
+        :param pixel_reader:
+            Callable to read the pixels from scenes. `RasterCollection.read_pixels` by
+            default. Any other callable can be provided as long as it returns a
+            GeoDataFrame and accepts `vector_features` and `fpath_raster` as input
+            keyword arguments.
+        :param pixel_reader_kwargs:
+            optional keyword arguments to pass on to `pixel_reader`.
+        """
+        # loop over scenes and read the pixel values. Carry out reprojection where
+        # necessary
+        pixel_scene_list = []
+        for _, item in self.metadata.iterrows():
+            pixel_reader_kwargs.update({
+                'vector_features': self.mapper_configs.feature.to_geoseries(),
+            })
+            try:
+                pixels = pixel_reader.__call__(item.real_path, **pixel_reader_kwargs)
+                pixels[self.time_column] = item[self.time_column]
+            except Exception as e:
+                raise ValueError(f'Could not read pixel data: {e}')
+
+            # reproject pixels if necessary
+            pixels.to_crs(epsg=item.target_epsg, inplace=True)
+            pixel_scene_list.append(pixels)
+
+        self.data = pd.concat(pixel_scene_list)
+
+    def load_scenes(
+        self,
+        scene_kwargs: Optional[Dict[str, Any]] = None,
+        pixel_kwargs: Optional[Dict[str, Any]] = None
+    ) -> None:
+        """
+        Load scenes from `~Mapper.query_scenes` result into a `SceneCollection`
+        (Polygon and MultiPolygon geometries) or into a `GeoDataFrame` (Point
+        and MultiPoint geometries).
+
+        Mosaicing operations are handled on the fly so calling programs will
+        always receive analysis-ready data (e.g., when working across image
+        tiles).
+
+        :param scene_kwargs:
+            key-word arguments to pass to `~Mapper._load_scenes_collection`
+            for handling EOdal scenes. These arguments *MUST* be provided when
+            using `Polygon` or `MulitPolygon` geometries in the `Mapper` call.
+        :param pixel_kwargs:
+            key-word arguments to pass to `~Mapper._load_pixels` for handling
+            single pixel values. These arguments *MUST* be provided when
+            using `Point` or `MulitPoint` geometries in the `Mapper` call.
+        """
+        # check if the correct keyword arguments have been passed
+        if not self._geoms_are_points:
+            if scene_kwargs is None:
+                raise ValueError(
+                    'Since Polygon/MultiPolygon geometries are provided ' + \
+                    '`pixel_kwargs` must not be empty'
+                )
+
+        # check if scenes have been queried and found
+        if self.metadata is None:
+            warnings.warn(
+                'No scenes are available - have you already executed Mapper.query_scenes()?'
+            )
+            return
+        if self.metadata.empty:
+            warnings.warn(
+                'No scenes were found - consider modifying your search criteria'
+            )
+            return
+
+        # check the spatial reference system of the scenes found. The mapper class
+        # will ensure that all of them will be available in the same CRS. The CRS
+        # is selected which most of the scenes already have to keep the reprojection
+        # costs as small as possible (and also to avoid resampling induced errors)
+        try:
+            self.metadata['target_epsg'] = self.metadata.epsg.mode().values[0]
+        except KeyError as e:
+            raise ValueError(f'Could not determine CRS of scenes: {e}')
+
+        # check if mosaicing scenes is required. This is done by checking the sensing_time
+        # time stamps. If there are multiple scenes with the same time stamp they must be
+        # mosaiced into a single scene
+        self.metadata['mosaicing'] = False
+        duplicated_idx = self.metadata[self.metadata.duplicated([self.time_column])].index
+        self.metadata.loc[duplicated_idx, 'mosaicing'] = True
+
+        # provide paths to raster data. Depending on th settings, this is a path on the
+        # file system or a URL
+        self.metadata['real_path'] = ''
+        if settings.USE_STAC:
+            self.metadata['real_path'] = self.metadata['assets']
+        else:
+            self.metadata['real_path'] = self.metadata.apply(
+                lambda x: reconstruct_path(record=x), axis=1
+            )
+
+        # load the data depending on the geometry type of the feature(s)
+        if self._geoms_are_points:
+            self._load_pixels(**pixel_kwargs)
+        else:
+            self._load_scenes_collection(**scene_kwargs)
diff --git a/eodal/metadata/database/querying.py b/eodal/metadata/database/querying.py
index 7acd8f05..d18dc8ed 100644
--- a/eodal/metadata/database/querying.py
+++ b/eodal/metadata/database/querying.py
@@ -21,10 +21,16 @@
 
 import geopandas as gpd
 
-from sqlalchemy import create_engine
+from datetime import datetime
+from geoalchemy2.functions import ST_Intersects
+from geoalchemy2.functions import ST_GeomFromText
+from shapely.geometry import Polygon
+from sqlalchemy import asc, and_, create_engine, inspect
 from sqlalchemy.orm.session import sessionmaker
+from typing import List
 
 from eodal.config import get_settings
+from eodal.mapper.filter import Filter
 from eodal.metadata.database import Regions, S2_Raw_Metadata
 from eodal.utils.exceptions import DataNotFoundError, RegionNotFoundError
 
@@ -35,6 +41,12 @@
 engine = create_engine(DB_URL, echo=Settings.ECHO_DB)
 session = sessionmaker(bind=engine)()
 
+# map platforms to metadata models
+platform_mapping = {
+    'sentinel1': 'S1_Raw_Metadata',
+    'sentinel2': 'S2_Raw_Metadata',
+    'ps_superdove': 'PS_SuperDove_Metadata'
+}
 
 def get_region(region: str) -> gpd.GeoDataFrame:
     """
@@ -57,6 +69,62 @@ def get_region(region: str) -> gpd.GeoDataFrame:
     except Exception as e:
         raise RegionNotFoundError(f"{region} not found: {e}")
 
+def find_raw_data_by_bbox(
+    platform: str,
+    time_start: datetime,
+    time_end: datetime,
+    bounding_box: Polygon,
+    metadata_filters: List[Filter]
+) -> gpd.GeoDataFrame:
+    """
+    Query the metadata DB by a geographic bounding box, time period
+    and custom filters
+    """
+    # get the DB table for a given platform
+    inspector = inspect(engine)
+    db_tables = inspector.get_table_names(schema=Settings.DEFAULT_SCHEMA)
+    # the table names follow the scheme <platform>_raw_metadata
+    platform_table = f'{platform}_raw_metadata'
+    if platform_table not in db_tables:
+        raise ValueError(f'{platform} was not found in database')
+
+    # loop over metadata filters and construct the sqlalchemy expressions
+    filter_str = ''
+    for metadata_filter in metadata_filters:
+        value = metadata_filter.value
+        if type(metadata_filter.value) == str:
+            value = f'"{metadata_filter.value}"'
+        filter_str += f'.filter(db_model.{metadata_filter.entity} ' + \
+            f'{metadata_filter.operator} {value})'
+
+    # build the query using "db_model" as a placeholder for the actual db-model
+    # of the platform
+    bounding_box_wkt = f'SRID=4326;{bounding_box.wkt}'
+    query_statement_exc = \
+        f"""
+        (session.query(db_model)
+        .filter(ST_Intersects(db_model.geom, ST_GeomFromText(bounding_box_wkt)))
+        .filter(
+            and_(
+                db_model.sensing_time <= time_end,
+                db_model.sensing_time >= time_start,
+            )
+        )
+        """ + filter_str + \
+        """
+        .order_by(db_model.sensing_date.asc())
+        ).statement
+        """
+    exec(f'from eodal.metadata.database import {platform_mapping[platform]} as db_model')
+    query_statement = eval(query_statement_exc.replace('\n','').replace(' ',''))
+
+    # read returned records as a GeoDataFrame and return
+    try:
+        return gpd.read_postgis(query_statement, session.bind)
+    except Exception as e:
+        raise DataNotFoundError(f"Could not find {platform} data by bounding box: {e}")
+    
+
 def get_s2_tile_footprint(tile_name: str) -> gpd.GeoDataFrame:
     """
     Queries the geographic extent of a Sentinel-2 tile
diff --git a/eodal/metadata/sentinel2/parsing.py b/eodal/metadata/sentinel2/parsing.py
index bb783275..6bb57df7 100644
--- a/eodal/metadata/sentinel2/parsing.py
+++ b/eodal/metadata/sentinel2/parsing.py
@@ -35,13 +35,12 @@
 from datetime import date
 
 from eodal.config import get_settings
-from eodal.config.sentinel2 import Sentinel2
+from eodal.utils.constants.sentinel2 import s2_band_mapping
 from eodal.utils.exceptions import UnknownProcessingLevel
 from eodal.utils.exceptions import InputError
 from eodal.utils.warnings import NothingToDo
 
 logger = get_settings().logger
-S2 = Sentinel2()
 
 
 def parse_MTD_DS(in_file: Path) -> Dict[str, Any]:
@@ -64,7 +63,7 @@ def parse_MTD_DS(in_file: Path) -> Dict[str, Any]:
 
     # now, the values of some relevant tags can be extracted:
     metadata = dict()
-    band_names = list(S2.BAND_INDICES.keys())
+    band_names = list(s2_band_mapping.keys())
 
     datatakeIdentifier_xml = xmldoc.getElementsByTagName("Datatake_Info")
     element = datatakeIdentifier_xml[0]
diff --git a/eodal/metadata/stac/client.py b/eodal/metadata/stac/client.py
index 4d42c829..4ade8bf8 100644
--- a/eodal/metadata/stac/client.py
+++ b/eodal/metadata/stac/client.py
@@ -4,24 +4,26 @@
 
 from __future__ import annotations
 
+import geopandas as gpd
 import pandas as pd
+import warnings
 
-from datetime import date, datetime
+from datetime import datetime
 from pystac_client import Client
 from pystac_client.stac_api_io import StacApiIO
 from shapely.geometry import box, Polygon
-from typing import Any, Dict, List, Optional
+from typing import Any, Dict, List
 
 from eodal.config import get_settings, STAC_Providers
+from eodal.mapper.filter import Filter
 from eodal.utils.decorators import prepare_bbox
 from eodal.utils.reprojection import infer_utm_zone
-from eodal.utils.sentinel2 import ProcessingLevels
 
 Settings = get_settings()
 
 def query_stac(
-    date_start: date,
-    date_end: date,
+    time_start: datetime,
+    time_end: datetime,
     collection: str,
     bounding_box: Polygon,
 ) -> List[Dict[str, Any]]:
@@ -31,10 +33,10 @@ def query_stac(
 
     This is a sensor-agnostic function called by sensor-specific ones.
 
-    :param date_start:
-        start date of the time period
-    :param date_end:
-        end date of the time period
+    :param time_start:
+        start of the time period
+    :param time_end:
+        end of the time period
     :param collection:
         name of the collection
         (e.g., sentinel-2-l2a for Sentinel-2 L2A data)
@@ -49,7 +51,7 @@ def query_stac(
     stac_api_io.session.verify = Settings.STAC_API_IO_CA_BUNDLE
     cat = Client.from_file(Settings.STAC_BACKEND.URL, stac_io=stac_api_io)
     # transform dates into the required format (%Y-%m-%d)
-    datestr = f'{date_start.strftime("%Y-%m-%d")}/{date_end.strftime("%Y-%m-%d")}'
+    datestr = f'{time_start.strftime("%Y-%m-%d")}/{time_end.strftime("%Y-%m-%d")}'
     # search datasets on catalog
     search = cat.search(
         collections=collection,
@@ -65,33 +67,68 @@ def query_stac(
     scenes = item_json["features"]
     return scenes
 
+def _filter_criteria_fulfilled(
+    metadata_dict: Dict[str, Any],
+    metadata_filters: List[Filter]
+) -> bool:
+    """
+    Check if a scene fulfills the metadata filter criteria
+
+    :param metadata_dict:
+        scene metadata returned from STAC item
+    :param metadata_filters:
+        scene metadata filters to apply on metadata_dict
+    :returns:
+        `True` if all criteria passed, `False` if a single
+        criterion was not met.
+    """
+    criteria_fulfilled = True
+    for filter in metadata_filters:
+        if filter.entity in ['processing_level', 'product_type']:
+            continue
+        if filter.entity not in metadata_dict.keys():
+            warnings.warn(
+                f'{filter.entity} could not be retrieved from STAC -> skipping filter'
+            )
+        # check if the filter condition is met
+        condition_met = eval(
+            f'metadata_dict["{filter.entity}"] {filter.operator} {filter.value}'
+        )
+        if not condition_met:
+            criteria_fulfilled = False
+            break
+    return criteria_fulfilled
+
 @prepare_bbox
 def sentinel2(
-    cloud_cover_threshold: float, processing_level: ProcessingLevels, **kwargs
-) -> pd.DataFrame:
+    metadata_filters: List[Filter], **kwargs
+) -> gpd.GeoDataFrame:
     """
     Sentinel-2 specific STAC query allows filtering by scene-wide cloudy pixel
     percentage.
 
-    :param cloud_cover_threshold:
-        optional cloud cover threshold to filter datasets by scene cloud coverage.
-        Must be provided as number between 0 and 100%.
-    :param processing_level:
-        Sentinel-2 processing level
+    :param metadata_filters:
+        custom filters for filtering scenes in catalog by some metadata attributes
     :param kwargs:
         keyword arguments to pass to `query_stac` function
     :returns:
-        dataframe with references to found Sentinel-2 mapper
+        dataframe with references to found Sentinel-2 scenes
     """
     # check for processing level of the data and set the collection accordingly
-    processing_level_stac = eval(f"Settings.STAC_BACKEND.S2{processing_level.name}")
+    filter_entities = [x.entity for x in metadata_filters]
+    if 'processing_level' in filter_entities:
+        processing_level = [x.value for x in metadata_filters if x.entity == 'processing_level'][0]
+        processing_level = processing_level.replace('-','_')
+    processing_level_stac = eval(f"Settings.STAC_BACKEND.S2{processing_level}")
     kwargs.update({"collection": processing_level_stac})
 
     # query STAC catalog
-    scenes = query_stac(**kwargs)
+    stac_kwargs = kwargs.copy()
+    del stac_kwargs['platform']
+    scenes = query_stac(**stac_kwargs)
     # get STAC provider specific naming conventions
     s2 = Settings.STAC_BACKEND.Sentinel2
-    # loop over found mapper and check their cloud cover
+    # loop over scenes found and apply the Filters provided
     metadata_list = []
     for scene in scenes:
         # extract scene metadata required for Sentinel-2
@@ -115,6 +152,9 @@ def sentinel2(
             scene_id = props[s2.scene_id]
         except KeyError:
             scene_id = scene[s2.scene_id]
+        # TODO: think about a more generic way to do this. The problem is:
+        # we need to map the different STAC provider settings into EOdals
+        # metadata model to avoid having the user to think about it
         meta_dict = {
             "product_uri": product_uri,
             "scene_id": scene_id,
@@ -129,52 +169,77 @@ def sentinel2(
                 props[s2.sensing_time], s2.sensing_time_fmt
             ),
             "sun_azimuth_angle": props[s2.sun_azimuth_angle],
-            "sun_zenith_angle": props[s2.sun_zenith_angle]
+            "sun_zenith_angle": props[s2.sun_zenith_angle],
+            "geom": Polygon(scene['geometry']['coordinates'][0])
         }
         # get links to actual Sentinel-2 bands
         meta_dict["assets"] = scene["assets"]
-        # only keep scene if the cloudy pixel percentage is not above
-        # the user-defined threshold (in theory, this could also be directly
-        # passed to the STAC API but then the function is less generic
-        if meta_dict["cloudy_pixel_percentage"] <= cloud_cover_threshold:
+
+        # apply filters
+        append_scene = _filter_criteria_fulfilled(meta_dict, metadata_filters)
+        if append_scene:
             metadata_list.append(meta_dict)
 
-    # create pandas DataFrame out of scene metadata records
-    return pd.DataFrame(metadata_list)
+    # create geppandas GeoDataFrame out of scene metadata records
+    df = pd.DataFrame(metadata_list)
+    df.sort_values(by='sensing_time', inplace=True)
+    return gpd.GeoDataFrame(df, geometry='geom', crs=4326)
 
 @prepare_bbox
-def sentinel1(collection: Optional[str] = 'sentinel-1-rtc', **kwargs) -> pd.DataFrame:
+def sentinel1(
+    metadata_filters: List[Filter], **kwargs
+) -> gpd.GeoDataFrame:
     """
-    Sentinel-1 specific STAC query function to retrieve mapper from MSPC
+    Sentinel-1 specific STAC query function to retrieve mapper from MSPC.
 
-    :param collection:
-        Sentinel-1 collection to use. Must be one of 'sentinel-1-grd' (ground
-        range detected), 'sentinel-1-rtc' (radiometrically terrain corrected)
+    IMPORTANT:
+        Returns the RTC product by default if not stated otherwise (using
+        `EOdal.mapper.filter` on `product_type`
+
+    :param metadata_filters:
+        custom filters for filtering scenes in catalog by some metadata attributes
     :param kwargs:
         :param kwargs:
         keyword arguments to pass to `query_stac` function
     :returns:
-        dataframe with references to found Sentinel-1 mapper
+        dataframe with references to found Sentinel-1 scenes
     """
 
     if Settings.STAC_BACKEND != STAC_Providers.MSPC:
-        raise ValueError('This method requires Microsoft Planetary Computer')
+        raise ValueError('This method currently requires Microsoft Planetary Computer')
 
-    # set collection to sentinel1-rtc
-    kwargs.update({'collection': collection})
+    # construct collection string (defaults to S1RTC if not stated otherwise in the
+    # metadata filters)
+    collection = 'S1'
+    if len([x.entity for x in metadata_filters if x.entity == 'product_type']) == 1:
+        collection += [x.value for x in metadata_filters if x.entity == 'product_type'][0]
+    else:
+        collection += 'RTC'
 
     # query the catalog
-    scenes = query_stac(**kwargs)
+    stac_kwargs = kwargs.copy()
+    del stac_kwargs['platform']
+    stac_kwargs.update({'collection': eval(f"Settings.STAC_BACKEND.{collection}")})
+    scenes = query_stac(**stac_kwargs)
     metadata_list = []
     for scene in scenes:
         metadata_dict = scene['properties']
         metadata_dict['assets'] = scene['assets']
+        metadata_dict['sensing_time'] = metadata_dict['datetime']
         metadata_dict['sensing_date'] = datetime.strptime(
-                metadata_dict['datetime'].split("T")[0], "%Y-%m-%d"
-            ).date()
-        metadata_list.append(metadata_dict)
+            metadata_dict['sensing_time'].split("T")[0], "%Y-%m-%d"
+        ).date()
+        del metadata_dict['datetime']
         # infer EPSG code of the scene in UTM coordinates from its bounding box
         bbox = box(*scene['bbox'])
         metadata_dict['epsg'] = infer_utm_zone(bbox)
-
-    return pd.DataFrame(metadata_list)
+        metadata_dict['geom'] = bbox
+        # apply filters
+        append_scene = _filter_criteria_fulfilled(metadata_dict, metadata_filters)
+        if append_scene:
+            metadata_list.append(metadata_dict)
+
+    # create geppandas GeoDataFrame out of scene metadata records
+    df = pd.DataFrame(metadata_list)
+    df.sort_values(by='sensing_time', inplace=True)
+    return gpd.GeoDataFrame(df, geometry='geom', crs=4326)
diff --git a/eodal/operational/__init__.py b/eodal/operational/__init__.py
deleted file mode 100644
index e69de29b..00000000
diff --git a/eodal/operational/archive/__init__.py b/eodal/operational/archive/__init__.py
deleted file mode 100644
index e69de29b..00000000
diff --git a/eodal/operational/archive/sentinel1.py b/eodal/operational/archive/sentinel1.py
deleted file mode 100644
index 4584696a..00000000
--- a/eodal/operational/archive/sentinel1.py
+++ /dev/null
@@ -1,155 +0,0 @@
-"""
-Function to keep a local Sentinel-1 archive up-to-date by comparing data available
-in your local archive for a given geographic region, time period, product type and
-sensor mode.
-
-The function not only checks if CREODIAS has new datasets available, it also automatically
-downloads them into a user-defined location.
-
-IMPORTANT: In order to receive results the region (i.e., geographic extent) for which
-to download data must be defined in the metadata DB.
-
-IMPORTANT: CREODIAS does not allow more than 2000 records (each Sentinel-1 scene is a record)
-to be queried at once. If you might exceed this threshold (e.g., your region is large and/or
-your time period is long) split your query into smaller chunks by using, e.g., shorter time
-periods for querying.
-
-Copyright (C) 2022 Gregor Perich & 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/>.
-"""
-
-from __future__ import annotations
-
-import pandas as pd
-import numpy as np
-
-from datetime import date
-from pathlib import Path
-from typing import Optional
-
-from eodal.downloader.sentinel1.creodias import query_creodias
-from eodal.downloader.utils.creodias import download_datasets
-
-from eodal.downloader.utils import unzip_datasets
-from eodal.metadata.database.querying import get_region
-from eodal.metadata.sentinel1.database.querying import find_raw_data_by_bbox
-from eodal.config import get_settings
-
-logger = get_settings()
-
-def pull_from_creodias(
-    date_start: date,
-    date_end: date,
-    path_out: Path,
-    region: str,
-    unzip: Optional[bool] = True,
-    overwrite_existing_zips: Optional[bool] = False,
-    **kwargs
-) -> pd.DataFrame:
-    """
-    Checks if CREODIAS has Sentinel-1 datasets not yet available locally
-    and downloads these datasets from CREODIAS.
-
-    IMPORTANT:
-        Unless specified otherwise, GRD (Ground Range Detected) products
-        acquired in IW (Interferometric Wide Swath) mode are considered.
-
-    NOTE:
-        CREODIAS limits the maximum amount of datasets to download within a
-        single query to somewhat around 2000. We therefore recommend to split
-        the query into smaller ones if the query is likely to hit this limit.
-
-    :param date start:
-        Start date of the database & creodias query
-    :param date end:
-        End date of the database & creodias query
-    :param path_out:
-        Out directory where the additional data from CREODIAS should be
-        downloaded to
-    :param region:
-        Region identifier of the Sentinel-2 archive. By region we mean a
-        geographic extent (bounding box) in which the data is organized. The bounding
-        box extent is taken from the metadata DB based on the region identifier.
-    :param unzip:
-        if True (default) datasets are unzipped and zip archives are deleted
-    :param overwrite_existing_zips:
-        if False (default) overwrites eventually existing zip files. If the download
-        process was interrupted (e.g., due to a connection timeout) setting the flag
-        to True can save time because datasets already downloaded are ignored. NOTE:
-        The function does **not** check if a dataset was downloaded completely!
-    :param kwargs:
-        optional key-word arguments to pass to 
-        `~eodal.downloader.sentinel1.creodias.query_creodias` such sensor_mode and
-        product_type
-    :return:
-        dataframe with references to downloaded datasets
-    """
-
-    # query database to get the bounding box of the selected region
-    try:
-        region_gdf = get_region(region)
-    except Exception as e:
-        logger.error(f"Failed to query region: {e}")
-        return
-
-    # parse the region's geometry as extended well-known-text
-    bounding_box = region_gdf.geometry.iloc[0]
-    bounding_box_ewkt = f"SRID=4326;{bounding_box.wkt}"
-
-    # local database query to check what is already available locally
-    try:
-        product_type = kwargs.get('product_type', 'GRD')
-        sensor_mode = kwargs.get('sensor_mode', 'IW')
-        meta_db_df = find_raw_data_by_bbox(
-            date_start=date_start,
-            date_end=date_end,
-            product_type=product_type,
-            sensor_mode=sensor_mode,
-            bounding_box=bounding_box_ewkt,
-        )
-    except Exception as e:
-        logger.error(f"Failed to query local datasets: {e}")
-        return
-
-    # set max_records to 2000 (CREODIAS currently does not allow more)
-    max_records = 2000
-
-    # check for available datasets
-    datasets = query_creodias
-
-    # get .SAFE datasets from CREODIAS
-    datasets["product_uri"] = datasets.properties.apply(
-        lambda x: Path(x["productIdentifier"]).name
-    )
-
-    # compare with records from local metadata DB and keep those records
-    # not available locally
-    missing_datasets = np.setdiff1d(
-        datasets["product_uri"].values, meta_db_df["product_uri"].values
-    )
-    datasets_filtered = datasets[datasets.product_uri.isin(missing_datasets)]
-
-    # download those mapper not available in the local database from Creodias
-    download_datasets(
-        datasets=datasets_filtered,
-        download_dir=path_out,
-        overwrite_existing_zips=overwrite_existing_zips,
-    )
-
-    # unzip datasets
-    if unzip:
-        unzip_datasets(download_dir=path_out)
-
-    return datasets_filtered
diff --git a/eodal/operational/archive/sentinel2.py b/eodal/operational/archive/sentinel2.py
deleted file mode 100644
index c2cfcc43..00000000
--- a/eodal/operational/archive/sentinel2.py
+++ /dev/null
@@ -1,159 +0,0 @@
-"""
-Function to keep a local Sentinel-2 archive up-to-date by comparing data available
-in your local archive for a given geographic region, time period and processing level
-with records available from CREODIAS (user credentials required).
-
-The function not only checks if CREODIAS has new datasets available, it also automatically
-downloads them into a user-defined location.
-
-IMPORTANT: In order to receive results the region (i.e., geographic extent) for which
-to download data must be defined in the metadata DB.
-
-IMPORTANT: CREODIAS does not allow more than 2000 records (each Sentinel-2 scene is a record)
-to be queried at once. If you might exceed this threshold (e.g., your region is large and/or
-your time period is long) split your query into smaller chunks by using, e.g., shorter time
-periods for querying.
-
-Copyright (C) 2022 Gregor Perich & 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/>.
-"""
-
-from __future__ import annotations
-
-import pandas as pd
-import numpy as np
-
-from datetime import date
-from pathlib import Path
-from typing import Optional
-
-from eodal.utils.constants import ProcessingLevels
-from eodal.downloader.sentinel2.creodias import query_creodias
-from eodal.downloader.utils.creodias import download_datasets
-
-from eodal.downloader.utils import unzip_datasets
-from eodal.metadata.database.querying import get_region
-from eodal.metadata.sentinel2.database.querying import find_raw_data_by_bbox
-from eodal.config import get_settings
-
-logger = get_settings()
-
-
-def pull_from_creodias(
-    date_start: date,
-    date_end: date,
-    processing_level: ProcessingLevels,
-    path_out: Path,
-    region: str,
-    cloud_cover_threshold: Optional[int] = 100,
-    unzip: Optional[bool] = True,
-    overwrite_existing_zips: Optional[bool] = False,
-) -> pd.DataFrame:
-    """
-    Checks if CREODIAS has Sentinel-2 datasets not yet available locally
-    and downloads these datasets from CREODIAS.
-
-    NOTE:
-        CREODIAS limits the maximum amount of datasets to download within a
-        single query to somewhat around 2000. We therefore recommend to split
-        the query into smaller ones if the query is likely to hit this limit.
-
-    :param date start:
-        Start date of the database & creodias query
-    :param date end:
-        End date of the database & creodias query
-    :param processing_level:
-        Select S2 processing level L1C or L2A
-    :param path_out:
-        Out directory where the additional data from CREODIAS should be
-        downloaded to
-    :param region:
-        Region identifier of the Sentinel-2 archive. By region we mean a
-        geographic extent (bounding box) in which the data is organized. The bounding
-        box extent is taken from the metadata DB based on the region identifier.
-    :param cloud_cover_threshold:
-        cloud cover threshold (0-100%) to be used for CREODIAS query. Defaults
-        to 100% (i.e., also completely cloudy mapper are downloaded).
-    :param unzip:
-        if True (default) datasets are unzipped and zip archives are deleted
-    :param overwrite_existing_zips:
-        if False (default) overwrites eventually existing zip files. If the download
-        process was interrupted (e.g., due to a connection timeout) setting the flag
-        to True can save time because datasets already downloaded are ignored. NOTE:
-        The function does **not** check if a dataset was downloaded completely!
-    :return:
-        dataframe with references to downloaded datasets
-    """
-
-    # query database to get the bounding box of the selected region
-    try:
-        region_gdf = get_region(region)
-    except Exception as e:
-        logger.error(f"Failed to query region: {e}")
-        return
-
-    # parse the region's geometry as extended well-known-text
-    bounding_box = region_gdf.geometry.iloc[0]
-    bounding_box_ewkt = f"SRID=4326;{bounding_box.wkt}"
-
-    # local database query to check what is already available locally
-    try:
-        meta_db_df = find_raw_data_by_bbox(
-            date_start=date_start,
-            date_end=date_end,
-            processing_level=processing_level,
-            bounding_box=bounding_box_ewkt,
-        )
-    except Exception as e:
-        logger.error(f"Failed to query local datasets: {e}")
-        return
-
-    # set max_records to 2000 (CREODIAS currently does not allow more)
-    max_records = 2000
-
-    # check for available datasets
-    datasets = query_creodias(
-        start_date=date_start,
-        end_date=date_end,
-        max_records=max_records,
-        processing_level=processing_level,
-        bounding_box=bounding_box,
-        cloud_cover_threshold=cloud_cover_threshold,
-    )
-
-    # get .SAFE datasets from CREODIAS
-    datasets["product_uri"] = datasets.properties.apply(
-        lambda x: Path(x["productIdentifier"]).name
-    )
-
-    # compare with records from local metadata DB and keep those records
-    # not available locally
-    missing_datasets = np.setdiff1d(
-        datasets["product_uri"].values, meta_db_df["product_uri"].values
-    )
-    datasets_filtered = datasets[datasets.product_uri.isin(missing_datasets)]
-
-    # download those mapper not available in the local database from Creodias
-    download_datasets(
-        datasets=datasets_filtered,
-        download_dir=path_out,
-        overwrite_existing_zips=overwrite_existing_zips,
-    )
-
-    # unzip datasets
-    if unzip:
-        unzip_datasets(download_dir=path_out)
-
-    return datasets_filtered
diff --git a/eodal/operational/cli/__init__.py b/eodal/operational/cli/__init__.py
deleted file mode 100644
index d1e434eb..00000000
--- a/eodal/operational/cli/__init__.py
+++ /dev/null
@@ -1,4 +0,0 @@
-from .sentinel2 import cli_s2_pipeline_fun
-from .sentinel2 import cli_s2_creodias_update
-from .sentinel2 import cli_s2_scene_selection
-from .sentinel2 import cli_s2_sen2cor_update
diff --git a/eodal/operational/cli/planet_scope.py b/eodal/operational/cli/planet_scope.py
deleted file mode 100644
index 4b88888a..00000000
--- a/eodal/operational/cli/planet_scope.py
+++ /dev/null
@@ -1,64 +0,0 @@
-'''
-Operational Command-Line Utilities for Planet-Scope data
-'''
-
-from __future__ import annotations
-
-from pathlib import Path
-from typing import Any, Optional, Dict
-
-from eodal.config.settings import get_settings
-from eodal.metadata.planet_scope.database.querying import get_scene_metadata
-from eodal.metadata.planet_scope.parsing import parse_metadata
-from eodal.metadata.planet_scope.database.ingestion import metadata_dict_to_database
-
-logger = get_settings().logger
-
-def cli_ps_scenes_ingestion(
-    ps_raw_data_archive: Path,
-    path_options: Optional[Dict[str, Any]] = None,
-) -> None:
-    """
-    Loops over the Planet-Scope raw data archive and checks for each
-    scene if it exists already in the eodal metadata DB. If not
-    the entry is added.
-
-    :param ps_raw_data_archive:
-        Planet-Scope raw data archive to monitor.
-    :param path_options:
-        optional dictionary specifying storage_device_ip, storage_device_ip_alias
-        (if applicable) and mount point in case the data is stored on a NAS
-        and should be accessible from different operating systems or file systems
-        with different mount points of the NAS share. If not provided, the absolute
-        path of the dataset is used in the database.
-    """
-    # loop over all mapper
-    for scene_dir in ps_raw_data_archive.iterdir():
-
-        try:
-            # query database for a given scene
-            scene_metadata = parse_metadata(scene_dir)
-            scene_df = get_scene_metadata(scene_id=scene_metadata['scene_id'])
-            if scene_df.empty:
-                # some path handling if required
-                if path_options is not None:
-                    scene_metadata["storage_device_ip"] = path_options.get(
-                        "storage_device_ip", ""
-                    )
-                    scene_metadata["storage_device_ip_alias"] = path_options.get(
-                        "storage_device_ip_alias", ""
-                    )
-                    mount_point = path_options.get("mount_point", "")
-                    mount_point_replacement = path_options.get(
-                        "mount_point_replacement", ""
-                    )
-                    scene_metadata["storage_share"] = scene_metadata[
-                        "storage_share"
-                    ].replace(mount_point, mount_point_replacement)
-                # ingest it into the database
-                metadata_dict_to_database(scene_metadata)
-                logger.info(f"Ingested scene metadata for {scene_metadata['scene_id']} into DB")
-            else:
-                logger.info(f"{scene_metadata['scene_id']} already in database")
-        except Exception as e:
-            logger.error(f"{scene_metadata['scene_id']} produced an error: {e}")
diff --git a/eodal/operational/cli/sentinel2.py b/eodal/operational/cli/sentinel2.py
deleted file mode 100644
index 46f59791..00000000
--- a/eodal/operational/cli/sentinel2.py
+++ /dev/null
@@ -1,465 +0,0 @@
-"""
-Scriptable high-level function interfaces to the Sentinel-2 processing pipeline and
-related operational functionalities.
-
-NOTE:
-    All of the functions make either use of eodal's metadata DB ("offline" mode)
-    or query a STAC provider ("online" mode)
-
-Copyright (C) 2022 Gregor Perich & 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/>.
-"""
-
-from __future__ import annotations
-
-import matplotlib.pyplot as plt
-import numpy as np
-import pandas as pd
-import shutil
-
-from datetime import date
-from datetime import datetime
-from pathlib import Path
-from typing import Any
-from typing import Dict
-from typing import Optional
-from typing import Union
-
-from eodal.config import get_settings
-from eodal.metadata.sentinel2.database import meta_df_to_database
-from eodal.metadata.sentinel2.database.ingestion import metadata_dict_to_database
-from eodal.metadata.sentinel2.database.querying import find_raw_data_by_tile
-from eodal.metadata.sentinel2.database.querying import get_scene_metadata
-from eodal.metadata.sentinel2.parsing import parse_s2_scene_metadata
-from eodal.operational.archive.sentinel2 import pull_from_creodias
-from eodal.operational.resampling.sentinel2 import exec_pipeline
-from eodal.utils.constants.sentinel2 import ProcessingLevels
-
-
-logger = get_settings().logger
-
-def cli_s2_pipeline_fun(
-    processed_data_archive: Path,
-    date_start: date,
-    date_end: date,
-    tile: str,
-    processing_level: ProcessingLevels,
-    n_threads: int,
-    resampling_options: dict,
-    path_options: Optional[dict] = None,
-) -> None:
-    """
-    Function calling Sentinel-2 pre-processing pipeline for user-defined
-    time period and Sentinel-2 tile. The pre-processing pipeline brings
-    those Sentinel-2 bands with a spatial resolution of 20m into a spatial
-    resolution of 10m and generates scene quicklooks. The function accepts
-    Sentinel-2 input data in .SAFE format in L1C and L2A processing level.
-
-    :param processed_data_archive:
-        directory where to store the processed data
-    :param date_start:
-        date defining the begin of the time period to process
-    :param date_end:
-        date defining the end of the time period to process
-    :param tile:
-        Sentinel-2 tile to process (e.g., 'T32TLT')
-    :param processing_level:
-        Processing level of the Sentinel-2 data (accepts L1C and L2A)
-    :param n_threads:
-        number of threads to use for execution of the pipeline
-    :param resampling_options:
-        key-word arguments to pass to `~exex_pipeline`
-    :param path_options:
-        optional key-word arguments for handling filepaths of the resulting
-        datasets
-
-    Example
-    ------
-    # TODO!
-    """
-    # start the processing
-    metadata, failed_datasets = exec_pipeline(
-        processed_data_archive=processed_data_archive,
-        date_start=date_start,
-        date_end=date_end,
-        tile=tile,
-        processing_level=processing_level,
-        n_threads=n_threads,
-        **resampling_options,
-    )
-
-    # set storage paths
-    metadata["storage_share"] = str(processed_data_archive)
-
-    if path_options is not None:
-        metadata["storage_device_ip"] = path_options.get("storage_device_ip", "")
-        metadata["storage_device_ip_alias"] = path_options.get(
-            "storage_device_ip_alias", ""
-        )
-        metadata["path_type"] = "posix"
-        mount_point = path_options.get("mount_point", "")
-        mount_point_replacement = path_options.get("mount_point_replacement", "")
-        metadata["storage_share"] = metadata["storage_share"].apply(
-            lambda x, mount_point=mount_point, mount_point_replacement=mount_point_replacement: str(
-                x
-            ).replace(
-                mount_point, mount_point_replacement
-            )
-        )
-
-    # rename columns to match DB attributes
-    metadata = metadata.rename(columns={"rgb_preview": "preview"})
-
-    # remove storage share from file-paths (breaks otherwise the entry in the DB)
-    metadata["bandstack"] = metadata["bandstack"].apply(lambda x: Path(x).name)
-    cols = ["preview"]
-    if processing_level.name == "L2A":
-        cols.append("scl")
-    for col in cols:
-        metadata[col] = metadata[col].apply(
-            lambda x: str(Path(Path(x).parent.name).joinpath(Path(x).name))
-        )
-
-    if "scl_preview" in metadata.columns:
-        metadata.drop("scl_preview", axis=1, inplace=True)
-
-    # write to database (set raw_metadata option to False)
-    meta_df_to_database(meta_df=metadata, raw_metadata=False)
-
-    # write failed datasets to disk
-    failed_datasets.to_csv(processed_data_archive.joinpath("failed_datasets.csv"))
-
-
-def cli_s2_creodias_update(
-    s2_raw_data_archive: Path,
-    region: str,
-    processing_level: ProcessingLevels,
-    cloud_cover_threshold: Optional[int] = 100,
-    path_options: Optional[Dict[str, str]] = None,
-    overwrite_existing_zips: Optional[bool] = False,
-) -> None:
-    """
-    Loops over an existing Sentinel-2 rawdata (i.e., *.SAFE datasets) archive
-    and checks the datasets available locally with those datasets available at
-    CREODIAS for a defined Area of Interest (AOI) and time period (we store the
-    data year by year, therefore, we always check an entire year).
-
-    The missing data (if any) is downloaded into a temporary download directory
-    `temp_dl` and unzip. The data is then copied into the actual Sentinel-2
-    archive and the metadata is extracted and ingested into the metadata base.
-
-    IMPORTANT: Requires a CREODIAS user account (user name and password).
-
-    IMPORTANT: For database ingestion it is important to map the paths correctly
-    and store them in the database in a way that allows accessing the data from
-    all your system components (``file_system_options``).
-
-    * In the easiest case, this the absolute path to the datasets (or URI)
-    * If your data is stored on a NAS, you might specify the address of the NAS
-      in the variable `storage_device_ip` and provide a `mount_point`, i.e.,
-      the path (or drive on Windows) where the NAS is mounted into your local
-      file system. Also aliasing of the is supported (`storage_device_ip_alias`),
-      however, if possible try to avoid it.
-
-
-    :param s2_raw_data_archive:
-        Sentinel-2 raw data archive (containing *.SAFE datasets) to monitor.
-        Existing datasets must have been already ingested into the metadata
-        base!
-    :param region:
-        eodal's archive philosophy organizes datasets by geographic regions.
-        Each region is identified by a unique region identifier (e.g., we use
-        `CH` for Switzerland) and has a geographic extent described by a polygon
-        geometry (bounding box) in geographic coordinates (WGS84). The geometry
-        also defines the geographic dimension of the CREODIAS query. It is
-        stored as a entry in the metadata base.
-    :param processing_level:
-        Sentinel-2 processing level (L1C or L2A) to check.
-    :param cloud_cover_threshold:
-        optional cloud cover threshold to filter out to cloudy mapper as integer
-        between 0 and 100%.
-    :param path_options:
-        optional dictionary specifying storage_device_ip, storage_device_ip_alias
-        (if applicable) and mount point in case the data is stored on a NAS
-        and should be accessible from different operating systems or file systems
-        with different mount points of the NAS share. If not provided, the absolute
-        path of the dataset is used in the database.
-    :param overwrite_existing_zips:
-        if False (default) overwrites eventually existing zip files. If the download
-        process was interrupted (e.g., due to a connection timeout) setting the flag
-        to True can save time because datasets already downloaded are ignored. NOTE:
-        The function does **not** check if a dataset was downloaded completely!
-
-    Example
-    -------
-
-    .. code-block:: python
-
-    from pathlib import Path
-    from eodal.operational.cli.sentinel2 import cli_s2_creodias_update
-    from eodal.utils.constants import ProcessingLevels
-
-    # define processing level (usually L2A but also L1C works)
-    processing_level = ProcessingLevels.L2A
-
-    # specifiy region, we use Switzerland (the extent of Switzerland is defined in the database)
-    region = 'CH'
-
-    # the archive should be always mounted in the same way for each user
-    user_name = '<your_username>'
-    s2_raw_data_archive = Path(f'/home/{user_name}/public/Evaluation/Satellite_data/Sentinel-2/Rawdata')
-
-    # file-system specific handling: this allows to store the paths in the database
-    # in such way that the dataset paths can be found from Linux and Windows machines
-    file_system_options = {
-        'storage_device_ip': '<your_nas_ip>',
-        'storage_device_ip_alias': '<alternative_nas_ip_if_any',
-        'mount_point': f'/home/{user_name}/public/'
-    }
-
-    cli_s2_creodias_update(
-        s2_raw_data_archive=s2_raw_data_archive,
-        region=region,
-        processing_level=processing_level,
-        path_options=file_system_options
-    )
-
-    """
-    # since the data is stored by year (each year is a single sub-directory) we
-    # can simple loop over the sub-directories and do the check
-    for path in Path(s2_raw_data_archive).iterdir():
-
-        if path.is_dir():
-
-            # get year automatically
-            year = int(path.name)
-
-            # create temporary download directory
-            path_out = path.joinpath(f"temp_dl_{year}")
-
-            if not path_out.exists():
-                path_out.mkdir()
-
-            # download data from CREODIAS
-            downloaded_ds = pull_from_creodias(
-                date_start=date(year, 1, 1),
-                date_end=date(year, 12, 31),
-                processing_level=processing_level,
-                path_out=path_out,
-                region=region,
-                cloud_cover_threshold=cloud_cover_threshold,
-                overwrite_existing_zips=overwrite_existing_zips,
-            )
-
-            if downloaded_ds.empty:
-                logger.info(f"No new datasets found for year {year} on CREODIAS")
-                continue
-
-            # move the datasets into the actual SAT archive (on level up)
-            error_happened = False
-            errored_datasets = []
-            error_msgs = []
-            parent_dir = path_out.parent
-            for _, record in downloaded_ds.iterrows():
-                try:
-                    shutil.move(record.dataset_name, "..")
-                    logger.info(f"Moved {record.dataset_name} to {parent_dir}")
-                    # once the dataset is moved successfully parse its metadata and
-                    # ingest it into the database
-                    in_dir = path.joinpath(record.dataset_name)
-                    scene_metadata, _ = parse_s2_scene_metadata(in_dir)
-
-                    # some path handling if required
-                    if path_options != {}:
-                        scene_metadata["storage_device_ip"] = path_options.get(
-                            "storage_device_ip", ""
-                        )
-                        scene_metadata["storage_device_ip_alias"] = path_options.get(
-                            "storage_device_ip_alias", ""
-                        )
-                        mount_point = path_options.get("mount_point", "")
-                        mount_point_replacement = path_options.get(
-                            "mount_point_replacement", ""
-                        )
-                        scene_metadata["storage_share"] = scene_metadata[
-                            "storage_share"
-                        ].replace(mount_point, mount_point_replacement)
-
-                    # database insert
-                    metadata_dict_to_database(scene_metadata)
-                    logger.info(
-                        f"Ingested scene metadata for {record.dataset_name} into DB"
-                    )
-
-                except Exception as e:
-                    error_happened = True
-                    logger.error(f"{record.dataset_name} produced an error: {e}")
-                    errored_datasets.append(record.dataset_name)
-                    error_msgs.append(e)
-
-            # delete the temp_dl directory
-            shutil.rmtree(path_out)
-
-            # log the errored datasets if any error happened
-            if error_happened:
-
-                # get the timestamp (date) of the current run
-                processing_date = datetime.now().date().strftime("%Y-%m-%d")
-
-                # write to log directory
-                log_dir = path_out.parent.joinpath("logs")
-                if not log_dir.exists():
-                    log_dir.mkdir()
-                errored_logfile = log_dir.joinpath("datasets_errored.csv")
-                if not errored_logfile.exists():
-                    with open(errored_logfile, "w") as src:
-                        line = "download_date,dataset_name,error_message\n"
-                        src.writelines(line)
-                with open(errored_logfile, "a") as src:
-                    for error in list(zip(errored_datasets, error_msgs)):
-                        line = processing_date + "," + error[0] + "," + str(error[1])
-                        src.writelines(line + "\n")
-
-
-def cli_s2_sen2cor_update(
-    s2_raw_data_archive: Path,
-    path_options: Optional[Dict[str, Any]] = {},
-) -> None:
-    """
-    Loops over the Sentinel-2 raw data archive and checks for each
-    scene (*.SAFE) if it exists already in the eodal metadata DB. If not
-    the entry is added. To speed-up the procedure single years can be
-    considered only.
-
-    :param s2_raw_data_archive:
-        Sentinel-2 raw data archive (containing *.SAFE datasets) to monitor.
-    :param path_options:
-        optional dictionary specifying storage_device_ip, storage_device_ip_alias
-        (if applicable) and mount point in case the data is stored on a NAS
-        and should be accessible from different operating systems or file systems
-        with different mount points of the NAS share. If not provided, the absolute
-        path of the dataset is used in the database.
-    """
-    # loop over all mapper (S2*.SAFE)
-    for scene_dir in s2_raw_data_archive.rglob("S2*.SAFE"):
-
-        try:
-            # check if the scene exists already in the DB (using its product_uri
-            # which is the same as the scene directory name)
-            product_uri = scene_dir.name
-            # query database
-            scene_df = get_scene_metadata(product_uri)
-            if scene_df.empty:
-                scene_metadata, _ = parse_s2_scene_metadata(scene_dir)
-                # some path handling if required
-                if path_options != {}:
-                    scene_metadata["storage_device_ip"] = path_options.get(
-                        "storage_device_ip", ""
-                    )
-                    scene_metadata["storage_device_ip_alias"] = path_options.get(
-                        "storage_device_ip_alias", ""
-                    )
-                    mount_point = path_options.get("mount_point", "")
-                    mount_point_replacement = path_options.get(
-                        "mount_point_replacement", ""
-                    )
-                    scene_metadata["storage_share"] = scene_metadata[
-                        "storage_share"
-                    ].replace(mount_point, mount_point_replacement)
-                # ingest it into the database
-                metadata_dict_to_database(scene_metadata)
-                logger.info(f"Ingested scene metadata for {product_uri} into DB")
-            else:
-                logger.info(f"{product_uri} already in database")
-        except Exception as e:
-            logger.error(f"{product_uri} produced an error: {e}")
-
-
-def cli_s2_scene_selection(
-    tile: str,
-    date_start: date,
-    date_end: date,
-    processing_level: ProcessingLevels,
-    out_dir: Path,
-    cloud_cover_threshold: Optional[Union[int, float]] = 100,
-) -> pd.DataFrame:
-    """
-    Function to query the Sentinel-2 metadata using a set of search criteria, including
-    filtering by date range, cloud cover and Sentinel-2 tile.
-
-    As a result, a CSV file with those metadata entries fulfilling the criteria is
-    returned plus a plot showing the cloud cover (extracted from the scene metadata)
-    per image acquisition date.
-
-    :param tile:
-        Sentinel-2 tile for which the scene selection should be performed. For instance,
-        'T32TMT'
-    :param date_start:
-        start date for filtering the data
-    :param end_date:
-        end data for filtering the data
-    :param processing_level:
-        Sentinel-2 processing level (L1C or L2A).
-    :param out_dir:
-        directory where to store the subset metadata CSV file and the cloud cover plot.
-    :param cloud_cover_threshold:
-        optional cloud cover threshold to filter out to cloudy mapper as integer
-        between 0 and 100%.
-    :returns:
-        metadata of mapper
-    """
-
-    # query metadata from database
-    try:
-        metadata = find_raw_data_by_tile(
-            date_start=date_start,
-            date_end=date_end,
-            processing_level=processing_level,
-            tile=tile,
-            cloud_cover_threshold=cloud_cover_threshold,
-        )
-    except Exception as e:
-        logger.error(f"Metadata query for Sentinel-2 data failed: {e}")
-        return
-
-    # calculate average cloud cover for the selected scenes
-    cc_avg = metadata.cloudy_pixel_percentage.mean()
-
-    # get timestamp of query execution
-    query_time = datetime.now().strftime("%Y%m%d_%H%M%S")
-    # write out metadata of the query as CSV
-    metadata.to_csv(out_dir.joinpath(f"{query_time}_query.csv"), index=False)
-
-    # Plot available mapper for query
-    fig = plt.figure(figsize=(15, 10), dpi=300)
-    ax = fig.add_subplot(111)
-    ax.plot(
-        metadata["sensing_date"],
-        metadata["cloudy_pixel_percentage"],
-        marker="o",
-        markersize=10,
-    )
-    ax.set_xlabel("Sensing Date")
-    ax.set_ylabel("Cloud cover [%]")
-    ax.set_ylim(0.0, 100.0)
-    ax.set_title(
-        f"Tile {tile} - No. of mapper: {metadata.shape[0]}"
-        + "\n"
-        + f"Average cloud cover: {np.round(cc_avg, 2)}%"
-    )
-    plt.savefig(out_dir.joinpath(f"{query_time}_query_CCplot.png"), bbox_inches="tight")
-    plt.close(fig)
-
-    return metadata
diff --git a/eodal/operational/mapping/__init__.py b/eodal/operational/mapping/__init__.py
deleted file mode 100644
index 895f9c3b..00000000
--- a/eodal/operational/mapping/__init__.py
+++ /dev/null
@@ -1,6 +0,0 @@
-from eodal.operational.mapping.sentinel1 import Sentinel1Mapper
-from eodal.operational.mapping.sentinel2 import Sentinel2Mapper
-from eodal.operational.mapping.mapper import MapperConfigs
-from eodal.operational.mapping.mapper import Mapper
-
-from eodal.operational.mapping.utils import plot_feature
diff --git a/eodal/operational/mapping/mapper.py b/eodal/operational/mapping/mapper.py
deleted file mode 100644
index a878cdc1..00000000
--- a/eodal/operational/mapping/mapper.py
+++ /dev/null
@@ -1,587 +0,0 @@
-"""
-Generic mapper module
-
-Copyright (C) 2022 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/>.
-"""
-
-from __future__ import annotations
-
-import cv2
-import geopandas as gpd
-import pandas as pd
-import uuid
-
-from datetime import date
-from geopandas import GeoDataFrame
-from pandas import DataFrame
-from pathlib import Path
-from shapely.geometry import box, MultiPolygon, Point, Polygon
-from sqlalchemy.exc import DatabaseError
-from typing import Any, Dict, List, Optional, Tuple, Union
-
-from eodal.config import get_settings
-from eodal.core.raster import RasterCollection
-from eodal.core.sensors import Sentinel1, Sentinel2
-from eodal.metadata.stac import sentinel1, sentinel2
-from eodal.metadata.utils import reconstruct_path
-from eodal.metadata.sentinel1.database.querying import find_raw_data_by_bbox as \
-    find_raw_data_by_bbox_sentinel1
-from eodal.metadata.sentinel2.database.querying import find_raw_data_by_bbox as \
-    find_raw_data_by_bbox_sentinel2
-from eodal.operational.resampling.utils import identify_split_scenes
-from eodal.utils.exceptions import BlackFillOnlyError, DataNotFoundError, InputError, \
-    STACError
-
-settings = get_settings()
-logger = settings.logger
-
-class Feature(object):
-    """
-    Class representing a feature, e.g., an area of interest.
-
-    :attrib identifier:
-        unique identifier of the feature
-    :attrib geom:
-        geometry of the feature
-    :attrib epsg:
-        epsg code of the feature's geometry
-    :attrib properties:
-        any key-value dictionary-like mapping of feature properties
-        (e.g., its name or other attributes spoken in terms of an
-        ESRI shapefile's table of attributes)
-    """
-
-    def __init__(
-        self,
-        identifier: Any,
-        geom: Union[Point, Polygon, MultiPolygon],
-        epsg: int,
-        properties: Optional[Dict[str, Any]],
-    ):
-        """
-        Initializes a new ``Feature`` instance.
-
-        :param identifier:
-            unique identifier of the feature
-        :param geom:
-            geometry of the feature
-        :param epsg:
-            epsg code of the feature's geometry
-        :param properties:
-            any key-value dictionary-like mapping of feature properties
-            (e.g., its name or other attributes spoken in terms of an
-            ESRI shapefile's table of attributes)
-        """
-        # some checks
-        if epsg <= 0:
-            raise TypeError("EPSG codes must be >= 0")
-        if not hasattr(geom, "__geo_interface__"):
-            raise TypeError("Geometries must implement the __geo_interface__")
-        if not isinstance(properties, dict):
-            raise TypeError("Only dictionary are accepted")
-
-        object.__setattr__(self, "identifier", identifier)
-        object.__setattr__(self, "geom", geom)
-        object.__setattr__(self, "epsg", epsg)
-        object.__setattr__(self, "properties", properties)
-
-    def __setattr__(self, *args):
-        raise TypeError("Feature attributes are immutable")
-
-    def __delattr__(self, *args):
-        raise TypeError("Feature attributes are immutable")
-
-    def __repr__(self):
-        return str(self.__dict__)
-
-    def to_gdf(self) -> GeoDataFrame:
-        """
-        Returns the feature as ``GeoDataFrame``
-
-        :returns:
-            ``Feature`` instance as ``GeoDataFrame``
-        """
-        self.properties.update({"epsg": self.epsg})
-        return GeoDataFrame(
-            self.properties,
-            index=[self.identifier],
-            crs=f"epsg:{self.epsg}",
-            geometry=[self.geom],
-        )
-
-class MapperConfigs(object):
-    """
-    Class defining configurations for the ``Mapper`` class
-
-    :attrib band_names:
-        names of raster bands to process from each dataset found
-    :attrib resampling_method:
-        resampling might become necessary when the spatial resolution
-        changes. Nearest neighbor by default.
-    :attrib spatial_resolution:
-        if provided brings all raster bands into the same spatial resolution
-    :attrib reducers:
-        optional list of spatial reducers (e.g., 'mean') converting all
-        raster observations from 2d arrays to scalars.
-    :atrrib tile_selection:
-        optional selection of tile ids for sensors following a tiling scheme
-        (e.g., S2 tiles, or Landsat PathRows).
-    """
-
-    def __init__(
-        self,
-        band_names: Optional[List[str]] = None,
-        resampling_method: Optional[int] = cv2.INTER_NEAREST_EXACT,
-        spatial_resolution: Optional[Union[int, float]] = 10.0,
-        reducers: Optional[List[str]] = None,
-        tile_selection: Optional[List[str]] = None,
-    ):
-        """
-        Constructs a new ``MapperConfig`` instance.
-
-        :param band_names:
-            names of raster bands to process from each dataset found
-        :param resampling_method:
-            resampling might become necessary when the spatial resolution
-            changes. Nearest neighbor by default.
-        :param spatial_resolution:
-            if provided brings all raster bands into the same spatial resolution
-        :param reducers:
-            optional list of spatial reducers (e.g., 'mean') converting all
-            raster observations from 2d arrays to scalars.
-        """
-        object.__setattr__(self, "band_names", band_names)
-        object.__setattr__(self, "resampling_method", resampling_method)
-        object.__setattr__(self, "spatial_resolution", spatial_resolution)
-        object.__setattr__(self, "reducers", reducers)
-        object.__setattr__(self, "tile_selection", tile_selection)
-
-    def __setattr__(self, *args):
-        raise TypeError("MapperConfigs attributes are immutable")
-
-    def __delattr__(self, *args):
-        raise TypeError("MapperConfigs attributes are immutable")
-
-    def __repr__(self):
-        return str(self.__dict__)
-
-class Mapper(object):
-    """
-    Generic Mapping class to extract raster data for a selection of areas of interest
-    (AOIs) and time period.
-
-    :attrib date_start:
-        start date of the time period to consider (inclusive)
-    :attrib date_end:
-        end date of the time period to consider (inclusive)
-    :attrib feature_collection:
-        ``GeoDataFrame`` or any vector file understood by ``fiona`` with
-        geometries of type ``Point``, ``Polygon`` or ``MultiPolygon``
-        defining the Areas Of Interest (AOIs) to extract (e.g., agricultural
-        field parcels). Each feature in the collection will be processed
-        separately
-    :attrib unique_id_attribute:
-        attribute in the `polygon_features`'s attribute table making each
-        feature (AOI) uniquely identifiable. If None (default) the features
-        are labelled by a unique-identifier created on the fly.
-    :attrib mapping_configs:
-        Mapping configurations specified by `~eodal.operational.mapper.MapperConfigs`.
-        Uses default configurations if not provided.
-    :attrib observations:
-        data structure for storing DB query results per AOI.
-    """
-
-    def __init__(
-        self,
-        date_start: date,
-        date_end: date,
-        feature_collection: Union[Path, GeoDataFrame],
-        unique_id_attribute: Optional[str] = None,
-        mapper_configs: MapperConfigs = MapperConfigs(),
-        collection: Optional[str] = ''
-    ):
-        """
-        Constructs a new ``Mapper`` instance.
-
-        :param date_start:
-            start date of the time period to consider (inclusive)
-        :param date_end:
-            end date of the time period to consider (inclusive)
-        :param feature_collection:
-            ``GeoDataFrame`` or any vector file understood by ``fiona`` with
-            geometries of type ``Point``, ``Polygon`` or ``MultiPolygon``
-            defining the Areas Of Interest (AOIs) to extract (e.g., agricultural
-            field parcels). Each feature in the collection will be processed
-            separately
-        :param unique_id_attribute:
-            attribute in the `polygon_features`'s attribute table making each
-            feature (AOI) uniquely identifiable. If None (default) the features
-            are labelled by a unique-identifier created on the fly.
-        :param mapping_configs:
-            Mapping configurations specified by `~eodal.operational.mapper.MapperConfigs`.
-            Uses default configurations if not provided.
-        """
-        object.__setattr__(self, "date_start", date_start)
-        object.__setattr__(self, "date_end", date_end)
-        object.__setattr__(self, "feature_collection", feature_collection)
-        object.__setattr__(self, "unique_id_attribute", unique_id_attribute)
-        object.__setattr__(self, "mapper_configs", mapper_configs)
-        object.__setattr__(self, "collection", collection)
-
-        observations: Dict[str, DataFrame] = None
-        object.__setattr__(self, "observations", observations)
-
-        features: Dict[str, Feature] = None
-        object.__setattr__(self, "features", features)
-
-    def __setattr__(self, *args):
-        raise TypeError("Mapper attributes are immutable")
-
-    def __delattr__(self, *args):
-        raise TypeError("Mapper attributes are immutable")
-
-    def _get_scenes(self, sensor: str) -> None:
-        """
-        Method to query available mapper. Works sensor-agnostic but requires a
-        sensor to be specified to select the correct metadata queries
-
-        :param sensor:
-            name of the sensor for which to search for mapper
-        """
-        # prepare features
-        aoi_features = self._prepare_features()
-
-        scenes = {}
-        # features implement the __geo_interface__
-        features = []
-        # extract all properties from the input features to preserve them
-        cols_to_ignore = [self.unique_id_attribute, "geometry", "geometry_wgs84"]
-        property_columns = [x for x in aoi_features.columns if x not in cols_to_ignore]
-        for _, feature in aoi_features.iterrows():
-
-            feature_uuid = feature[self.unique_id_attribute]
-            properties = feature[property_columns].to_dict()
-            feature_obj = Feature(
-                identifier=feature_uuid,
-                geom=feature["geometry"],
-                epsg=aoi_features.crs.to_epsg(),
-                properties=properties,
-            )
-            features.append(feature_obj.to_gdf())
-
-            # determine bounding box of the current feature using
-            # its representation in geographic coordinates
-            bbox = box(*feature.geometry_wgs84.bounds)
-
-            # use the resulting bbox to query the bounding box
-            # there a two options: use STAC or the PostgreSQL DB
-            # prepare sensor specific kwargs
-            kwargs = {
-                'date_start': self.date_start,
-                'date_end': self.date_end,
-                'bounding_box': bbox
-            }
-            if sensor.lower() == 'sentinel2':
-                kwargs.update({
-                    'processing_level': self.processing_level,
-                    'cloud_cover_threshold': self.cloud_cover_threshold
-                })
-            elif sensor.lower() == 'sentinel1':
-                kwargs.update({
-                    'collection': self.collection
-                })
-
-            if settings.USE_STAC:
-                try:
-                    scenes_df = eval(f'{sensor}(**kwargs)')
-                except Exception as e:
-                    raise STACError(f"Querying STAC catalog failed: {e}")
-            else:
-                try:
-                    scenes_df = eval(f'find_raw_data_by_bbox_{sensor}(**kwargs)')
-                except Exception as e:
-                    raise DatabaseError(f"Querying metadata DB failed: {e}")
-
-            # filter by tile if required
-            tile_ids = self.mapper_configs.tile_selection
-            if tile_ids is not None:
-                other_tile_idx = scenes_df[~scenes_df.tile_id.isin(tile_ids)].index
-                scenes_df.drop(other_tile_idx, inplace=True)
-
-            if scenes_df.empty:
-                logger.info(
-                    f"The query for feature {feature_uuid} returned now results"
-                )
-                continue
-
-            # check if the satellite data is in different projections
-            in_single_crs = scenes_df.epsg.unique().shape[0] == 1
-
-            # check if there are several mapper available for a single sensing date
-            # in this case merging of different datasets might be necessary
-            scenes_df_split = identify_split_scenes(scenes_df)
-            scenes_df["is_split"] = False
-
-            if not scenes_df_split.empty:
-                scenes_df.loc[
-                    scenes_df.product_uri.isin(scenes_df_split.product_uri), "is_split"
-                ] = True
-
-            # in case the mapper have different projections (most likely different UTM
-            # zone numbers) figure out which will be target UTM zone. To avoid too many
-            # reprojection operations of raster data later, the target CRS is that CRS
-            # most mapper have (expressed as EPSG code)
-            scenes_df["target_crs"] = scenes_df.epsg
-            if not in_single_crs:
-                most_common_epsg = scenes_df.epsg.mode().values
-                scenes_df.loc[
-                    ~scenes_df.epsg.isin(most_common_epsg), "target_crs"
-                ] = most_common_epsg[0]
-
-            # add the scenes_df DataFrame to the dictionary that contains the data for
-            # all features
-            scenes[feature_uuid] = scenes_df
-
-        # create feature collection
-        features_gdf = pd.concat(features)
-        # append raw scene count
-        features_gdf["raw_scene_count"] = 0
-        if len(scenes) > 0:
-            features_gdf["raw_scene_count"] = features_gdf.apply(
-                lambda x, scenes=scenes: scenes[x.name].shape[0], axis=1
-            )
-        features = features_gdf.__geo_interface__
-        object.__setattr__(self, "observations", scenes)
-        object.__setattr__(self, "feature_collection", features)
-
-    def _prepare_features(self) -> pd.DataFrame:
-        """
-        Prepares the feature collection for mapper
-
-        :returns:
-            `DataFrame` with prepared features
-        """
-        # read features and loop over them to process each feature separately
-        is_file = isinstance(self.feature_collection, Path) or isinstance(
-            self.feature_collection, str
-        )
-        if is_file:
-            try:
-                aoi_features = gpd.read_file(self.feature_collection)
-            except Exception as e:
-                raise InputError(
-                    f"Could not read polygon features from file "
-                    f"{self.feature_collection}: {e}"
-                )
-        else:
-            aoi_features = self.feature_collection.copy()
-
-        # for the query, the geometries are required in geographic coordinates
-        # however, we keep the original coordinates as well to avoid to many reprojections
-        aoi_features["geometry_wgs84"] = aoi_features["geometry"].to_crs(4326)
-
-        # check if there is a unique feature column
-        # otherwise create it using uuid
-        if self.unique_id_attribute is not None:
-            if not aoi_features[self.unique_id_attribute].is_unique:
-                raise ValueError(
-                    f'"{self.unique_id_attribute}" is not unique for each feature'
-                )
-        else:
-            object.__setattr__(self, "unique_id_attribute", "_uuid4")
-            aoi_features[self.unique_id_attribute] = [
-                str(uuid.uuid4()) for _ in aoi_features.iterrows()
-            ]
-
-        return aoi_features
-
-    def get_feature_ids(self) -> List:
-        """
-        Lists feature identifiers in feature collection
-
-        :returns:
-            list of feature identifiers
-        """
-        if isinstance(self.feature_collection, Path):
-            return []
-        return [x["id"] for x in self.feature_collection["features"]]
-
-    def get_feature(self, feature_id: Any) -> Dict[str, Any]:
-        """
-        Returns a feature in its ``__geo_interface__`` representation
-        out of the feature collection
-
-        :param feature_id:
-            feature identifier to use for extraction
-        :param as_gdf:
-            return feature as dictionary (default) or ``GeoDataFrame``?
-        :returns:
-            the feature with its properties and geometry
-        """
-        if isinstance(self.feature_collection, Path):
-            return {}
-        if isinstance(self.feature_collection, GeoDataFrame):
-            gdf = self.feature_collection[
-                self.feature_collection[self.unique_id_attribute] == feature_id
-            ]
-            return {} if gdf.empty else gdf.__geo_interface__
-        else:
-            res = [
-                x for x in self.feature_collection["features"] if x["id"] == feature_id
-            ]
-            if len(res) == 0:
-                raise KeyError(f'No feature found with ID "{feature_id}"')
-            return {"type": "FeatureCollection", "features": res}
-
-    def get_feature_scenes(self, feature_identifier: Any) -> DataFrame:
-        """
-        Returns a ``DataFrame`` with all mapper found for a
-        feature in the feature collection
-
-        NOTE:
-            The scene count is termed ``raw_scene_count``. This
-            highlights that the final scene count might be
-            different due to orbit and spatial design pattern.
-
-        :param feature_identifier:
-            unique identifier of the aoi. Must be the same identifier
-            used during the database query
-        :returns:
-            ``DataFrame`` with all mapper found for a given
-            set of search parameters
-        """
-        try:
-            return self.observations[feature_identifier].copy()
-        except Exception as e:
-            raise KeyError(f"{feature_identifier} did not return any results: {e}")
-
-    def _get_observation(
-        self, feature_id: Any, sensing_date: date, sensor: str, **kwargs
-    ) -> Union[gpd.GeoDataFrame, RasterCollection, Tuple,  None]:
-        """
-        Returns the scene data (observations) for a selected feature and date.
-
-        If for the date provided no mapper are found, the data from the scene(s)
-        closest in time is returned
-
-        :param feature_id:
-            identifier of the feature for which to extract observations
-        :param sensing_date:
-            date for which to extract observations (or the closest date if
-            no observations are available for the given date)
-        :param sensor:
-            name of the sensor for which to extract the observation.
-        :param kwargs:
-            optional key-word arguments to pass on to
-            `~eodal.core.sensors` class-specific methods
-        :returns:
-            depending on the geometry type of the feature either a
-            ``GeoDataFrame`` (geometry type: ``Point``) or ``Sentinel2Handler``
-            (geometry types ``Polygon`` or ``MultiPolygon``) is returned. if
-            the observation contains nodata, only, None is returned. If multiple
-            mapper must be read to get a single observation, the status 'multiple'
-            is returned.
-        """
-        # define variable for returning results
-        res = None
-        # get available observations for the AOI feature
-        scenes_df = self.observations.get(feature_id, None)
-        if scenes_df is None:
-            raise DataNotFoundError(
-                f'Could not find any mapper for feature with ID "{feature_id}"'
-            )
-
-        # get scene(s) closest to the sensing_date provided
-        min_delta = abs((scenes_df.sensing_date - sensing_date)).min()
-        scenes_date = scenes_df[
-            abs((scenes_df.sensing_date - sensing_date)) == min_delta
-        ].copy()
-
-        # map the dataset path(s) when working locally (no STAC)
-        if not settings.USE_STAC:
-            try:
-                scenes_date["real_path"] = scenes_date.apply(
-                    lambda x: reconstruct_path(record=x), axis=1
-                )
-            except Exception as e:
-                raise DataNotFoundError(
-                    f"Cannot find the mapper on the file system: {e}"
-                )
-        # get properties and geometry of the current feature from the collection
-        feature_dict = self.get_feature(feature_id)
-        feature_gdf = gpd.GeoDataFrame.from_features(feature_dict)
-        feature_gdf.crs = feature_dict["features"][0]["properties"]["epsg"]
-        # parse feature geometry in kwargs so that only a spatial subset is read
-        # in addition parse the S2 gain factor as "scale" argument
-        kwargs.update({"vector_features": feature_gdf})
-        # multiple mapper for a single date
-        # check what to do (re-projection, merging)
-        if scenes_date.shape[0] > 1:
-            return ('multiple', scenes_date, feature_gdf)
-        else:
-            # determine scene path (local environment) or URLs (STAC)
-            if settings.USE_STAC:
-                in_dir = scenes_date["assets"].iloc[0]
-            else:
-                in_dir = scenes_date["real_path"].iloc[0]
-            # if there is only one scene all we have to do is to read
-            # read pixels in case the feature's dtype is point
-            feature_geom = feature_dict["features"][0]["geometry"]["type"]
-            if feature_geom in ["Point", "MultiPoint"]:
-                if sensor.lower() == 'sentinel1':
-                    res = Sentinel1.read_pixels_from_safe(
-                        in_dir=in_dir,
-                        polarizations=self.mapper_configs.band_names,
-                        **kwargs
-                    )
-                elif sensor.lower() == 'sentinel2':
-                    res = Sentinel2.read_pixels_from_safe(
-                        in_dir=in_dir,
-                        band_selection=self.mapper_configs.band_names,
-                        **kwargs,
-                    )
-                res["sensing_date"] = scenes_date["sensing_date"].values[0]
-                res["scene_id"] = scenes_date["scene_id"].values[0]
-                res['sensing_time'] = scenes_date['sensing_time'].values[0]
-                # make sure the result is projected into the target EPSG code, otherwise
-                # the resulting GeoDataFrame contains wrong coordinates, i.e., the
-                # coordinates were not projected into the target CRS
-                if res.crs.to_epsg() != scenes_date.target_crs.unique()[0]:
-                    res.to_crs(epsg=scenes_date.target_crs.unique()[0], inplace=True)
-            # or the feature
-            else:
-                if sensor.lower() == 'sentinel1':
-                    try:
-                        res = Sentinel1.from_safe(
-                            in_dir=in_dir,
-                            band_selection=self.mapper_configs.band_names,
-                            **kwargs
-                        )
-                    except Exception as e:
-                        raise Exception from e
-                elif sensor.lower() == 'sentinel2':
-                    try:
-                        res = Sentinel2.from_safe(
-                            in_dir=in_dir,
-                            band_selection=self.mapper_configs.band_names,
-                            **kwargs,
-                        )
-                        self._resample_s2_scene(s2_scene=res)
-                    except BlackFillOnlyError:
-                        return res
-                    except Exception as e:
-                        raise Exception from e
-        return res
diff --git a/eodal/operational/mapping/merging.py b/eodal/operational/mapping/merging.py
deleted file mode 100644
index fd1d78fe..00000000
--- a/eodal/operational/mapping/merging.py
+++ /dev/null
@@ -1,212 +0,0 @@
-"""
-Module for merging raster datasets.
-
-Copyright (C) 2022 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/>.
-"""
-from __future__ import annotations
-# TODO: delete this module once the new Mapper works
-
-import os
-import geopandas as gpd
-import uuid
-
-from pathlib import Path
-from typing import Any, Dict, List, Optional, Tuple, Union
-from rasterio.merge import merge
-from rasterio.crs import CRS
-
-import eodal
-from eodal.config import get_settings
-from eodal.core.band import Band, GeoInfo
-from eodal.core.raster import RasterCollection, SceneProperties
-
-Settings = get_settings()
-
-
-def _get_crs_and_attribs(
-    in_file: Path, **kwargs
-) -> Tuple[GeoInfo, List[Dict[str, Any]]]:
-    """
-    Returns the ``GeoInfo`` from a multi-band raster dataset
-
-    :param in_file:
-        raster datasets from which to extract the ``GeoInfo`` and
-        attributes
-    :param kwargs:
-        optional keyword-arguments to pass to
-        `~eodal.core.raster.RasterCollection.from_multi_band_raster`
-    :returns:
-        ``GeoInfo`` and metadata attributes of the raster dataset
-    """
-
-    ds = RasterCollection.from_multi_band_raster(fpath_raster=in_file, **kwargs)
-    geo_info = ds[ds.band_names[0]].geo_info
-    attrs = [ds[x].get_attributes() for x in ds.band_names]
-    return geo_info, attrs
-
-
-def merge_datasets(
-    datasets: List[Path],
-    out_file: Optional[Path] = None,
-    target_crs: Optional[Union[int, CRS]] = None,
-    vector_features: Optional[Union[Path, gpd.GeoDataFrame]] = None,
-    scene_properties: Optional[SceneProperties] = None,
-    band_options: Optional[Dict[str, Any]] = None,
-    sensor: Optional[str] = None,
-    **kwargs,
-) -> Union[None, RasterCollection]:
-    """
-    Merges a list of raster datasets using the ``rasterio.merge`` module.
-
-    The function can handle datasets in different coordinate systems by resampling
-    the data into a common spatial reference system either provided in the function
-    call or infered from the first dataset in the list.
-
-    ATTENTION:
-        All datasets must have the same number of bands and data type!
-
-    :param datasets:
-        list of datasets (as path-like objects or opened raster datasets)
-        to merge into a single raster
-    :param out_file:
-        name of the resulting raster dataset (optional). If None (default)
-        returns a new ``SatDataHandler`` instance otherwise writes the data
-        to disk as new raster dataset.
-    :param target_crs:
-        optional target spatial coordinate reference system in which the output
-        product shall be generated. Must be passed as integer EPSG code or CRS
-        instance.
-    :param vector_features:
-        optional vector features to clip the merged dataset to (full bounding box).
-    :param scene_properties:
-        optional scene properties to set to the resulting merged dataset
-    :param band_options:
-        optional sensor-specific band options to pass to the sensor's
-        ``RasterCollection`` constructor
-    :param sensor:
-        if the data is from a sensor explicitly supported by eodal such as
-        Sentinel-2 the raster data is loaded into a sensor-specific collection
-    :param kwargs:
-        kwargs to pass to ``rasterio.warp.reproject``
-    :returns:
-        depending on the kwargs passed either `None` (if output is written to file directly)
-        or a `RasterCollection` instance if the operation takes place in memory
-    """
-    # check the CRS and attributes of the datasets first
-    crs_list = []
-    attrs_list = []
-    for dataset in datasets:
-        geo_info, attrs = _get_crs_and_attribs(in_file=dataset)
-        crs_list.append(geo_info.epsg)
-        attrs_list.append(attrs)
-
-    if target_crs is None:
-        # use CRS from first dataset
-        target_crs = crs_list[0]
-    # coordinate systems are not the same -> re-projection of raster datasets
-    if len(set(crs_list)) > 1:
-        pass
-    # all datasets have one coordinate system, check if it is the desired one
-    else:
-        if target_crs is not None:
-            if crs_list[0] != target_crs:
-                # re-projection into target CRS required
-                pass
-
-    # use rasterio merge to get a new raster dataset
-    dst_kwds = {"QUALITY": "100", "REVERSIBLE": "YES"}
-    try:
-        res = merge(datasets=datasets, dst_path=out_file, dst_kwds=dst_kwds, **kwargs)
-        if res is not None:
-            out_ds, out_transform = res[0], res[1]
-    except Exception as e:
-        raise Exception(f"Could not merge datasets: {e}")
-
-    # when out_file was provided the merged data is written to file directly
-    if out_file is not None:
-        return
-    # otherwise, create new SatDataHandler instance from merged datasets
-    # add scene properties if available
-    if sensor is None:
-        raster = RasterCollection(scene_properties=scene_properties)
-    else:
-        raster = eval(f'eodal.core.sensors.{sensor.lower()}.{sensor[0].upper() + sensor[1::]}' \
-                      + '(scene_properties=scene_properties)')
-    n_bands = out_ds.shape[0]
-    # take attributes of the first dataset
-    attrs = attrs_list[0]
-    geo_info = GeoInfo(
-        epsg=target_crs,
-        ulx=out_transform.c,
-        uly=out_transform.f,
-        pixres_x=out_transform.a,
-        pixres_y=out_transform.e,
-    )
-    for idx in range(n_bands):
-        band_attrs = attrs[idx]
-        nodata = band_attrs.get("nodatavals")
-        if isinstance(nodata, tuple):
-            nodata = nodata[0]
-        is_tiled = band_attrs.get("is_tiled")
-        scale = band_attrs.get("scales")
-        if isinstance(scale, tuple):
-            scale = scale[0]
-        offset = band_attrs.get("offsets")
-        if isinstance(offset, tuple):
-            offset = offset[0]
-        description = band_attrs.get("descriptions")
-        if isinstance(description, tuple):
-            description = description[0]
-        unit = band_attrs.get("units")
-        if isinstance(unit, tuple):
-            unit = unit[0]
-        raster.add_band(
-            band_constructor=Band,
-            band_name=f"B{idx+1}",
-            values=out_ds[idx, :, :],
-            geo_info=geo_info,
-            is_tiled=is_tiled,
-            scale=scale,
-            offset=offset,
-            band_alias=description,
-            unit=unit,
-        )
-
-    # clip raster collection if required to vector_features
-    if vector_features is not None:
-        tmp_dir = Settings.TEMP_WORKING_DIR
-        fname_tmp = tmp_dir.joinpath(f"{uuid.uuid4()}.tif")
-        raster.to_rasterio(fname_tmp)
-        if sensor is None:
-            expr = "RasterCollection"
-        else:
-            expr = f"eodal.core.sensors.{sensor.lower()}.{sensor[0].upper() + sensor[1::]}()"
-        if band_options is None:
-            band_options = {}
-        raster = eval(
-            f"""{expr}.from_multi_band_raster(
-            fpath_raster=fname_tmp,
-            vector_features=vector_features,
-            full_bounding_box_only=False,
-            **band_options
-        )"""
-        )
-        # set scene properties if available
-        if scene_properties is not None:
-            raster.scene_properties = scene_properties
-        os.remove(fname_tmp)
-
-    return raster
diff --git a/eodal/operational/mapping/sentinel1.py b/eodal/operational/mapping/sentinel1.py
deleted file mode 100644
index 1bc5f707..00000000
--- a/eodal/operational/mapping/sentinel1.py
+++ /dev/null
@@ -1,139 +0,0 @@
-'''
-Created on Aug 9, 2022
-
-@author: graflu
-'''
-from __future__ import annotations
-
-import geopandas as gpd
-
-from datetime import date
-from shapely.geometry import box
-from typing import Any, Union
-
-from eodal.config import get_settings
-from eodal.core.sensors import Sentinel1
-from .mapper import Mapper
-
-settings = get_settings()
-
-class Sentinel1Mapper(Mapper):
-    """
-    Spatial mapper class for Sentinel-1 data.
-    """
-    def __init__(
-        self,
-        *args,
-        **kwargs,
-    ):
-        """
-        Initializes a new Sentinel-1 Mapper object
-
-        :param args:
-            arguments to pass to the constructor of the ``Mapper`` super class
-        :param kwargs:
-            key-word arguments to pass to the constructor of the ``Mapper``
-            super class
-        """
-        # initialize super-class
-        Mapper.__init__(self, *args, **kwargs)
-
-    def get_scenes(self) -> None:
-        """
-        Queries the Sentinel-1 metadata catalog for a selected time period and
-        feature collection.
-
-        The scene selection and processing workflow contains several steps:
-
-        1.  Query the metadata catalog for **ALL** available mapper that overlap
-            the bounding box of a given ``Polygon`` or ``MultiPolygon``
-            feature.
-        2.  Check if for a single sensing date several mapper are available
-        3.  If yes check if that's due to Sentinel-1 tiling grid design. If yes
-            flag these mapper as potential merge candidates.
-        4.  If the mapper found have different spatial coordinate systems (CRS)
-            (usually different UTM zones) flag the data accordingly. The target
-            CRS is defined as that CRS the majority of mapper shares.
-        """
-        self._get_scenes(sensor='sentinel1')
-
-    def get_observation(
-        self, feature_id: Any, sensing_date: date, **kwargs
-    ) -> Union[gpd.GeoDataFrame, Sentinel1, None]:
-        """
-        Returns the scene data (observations) for a selected feature and date.
-
-        If for the date provided no mapper are found, the data from the scene(s)
-        closest in time is returned
-
-        :param feature_id:
-            identifier of the feature for which to extract observations
-        :param sensing_date:
-            date for which to extract observations (or the closest date if
-            no observations are available for the given date)
-        :param kwargs:
-            optional key-word arguments to pass on to
-            `~eodal.core.sensors.Sentinel2.from_safe`
-        :returns:
-            depending on the geometry type of the feature either a
-            ``GeoDataFrame`` (geometry type: ``Point``) or ``Sentinel2Handler``
-            (geometry types ``Polygon`` or ``MultiPolygon``) is returned. if
-            the observation contains nodata, only, None is returned.
-        """
-        res = self._get_obervation(feature_id=feature_id, sensing_date=sensing_date,
-                                    sensor='sentinel1', **kwargs)
-        if isinstance(res, tuple):
-            _, scenes_date, feature_gdf = res
-            # TODO: merge logic according to Sentinel-2
-        return res
-
-    def get_complete_timeseries(self):
-        pass
-
-# TODO
-if __name__ == '__main__':
-
-    from eodal.metadata.stac import sentinel1 as s1_stac_client
-    from eodal.core.sensors import Sentine1 
-    import matplotlib.pyplot as plt
-    from datetime import date
-    from pathlib import Path
-
-    # define time period
-    date_start = date(2022, 5, 1)
-    date_end = date(2022, 5, 31)
-
-    # provide bounding box
-    bounding_box_fpath = Path(
-        # '/home/graflu/public/Evaluation/Projects/KP0031_lgraf_PhenomEn/02_Field-Campaigns/Strickhof/WW_2022/Bramenwies.shp'
-        '/mnt/ides/Lukas/software/eodal/data/sample_polygons/ZH_Polygon_73129_ESCH_EPSG32632.shp'
-    )
-
-    # Sentinel-1
-    res_s1 = s1_stac_client(
-        date_start=date_start,
-        date_end=date_end,
-        vector_features=bounding_box_fpath,
-        collection='sentinel-1-rtc'
-    )
-    res_s1 = res_s1.sort_values(by='datetime')
-
-
-    f, ax = plt.subplots(nrows=1, ncols=res_s1.shape[0], figsize=(20,5), sharex=True)
-
-    # loop datasets and plot CR
-    for idx, rec in res_s1.iterrows():
-        asset = rec['assets']
-    
-        s1 = Sentinel1.from_safe(in_dir=asset, vector_features=bounding_box_fpath)
-        s1.calc_si('CR', inplace=True)
-
-        s1.plot_band('CR', ax=ax[idx], colorbar_label='CR [-]', vmin=0, vmax=4)
-        if idx > 0:
-            ax[idx].set_ylabel('')
-        ax[idx].set_title(f'{s1.scene_properties.acquisition_time}\n{s1.scene_properties.platform}')
-
-    f.tight_layout()
-    f.savefig('/mnt/ides/Lukas/software/eodal/img/sentinel1_cr.png', bbox_inches='tight')
-    plt.close(f)
-            
\ No newline at end of file
diff --git a/eodal/operational/mapping/sentinel2.py b/eodal/operational/mapping/sentinel2.py
deleted file mode 100644
index 809e1267..00000000
--- a/eodal/operational/mapping/sentinel2.py
+++ /dev/null
@@ -1,443 +0,0 @@
-"""
-Mapping module for Sentinel-2 data
-
-Copyright (C) 2022 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/>.
-"""
-from __future__ import annotations
-
-import os
-import cv2
-import geopandas as gpd
-import numpy as np
-import pandas as pd
-import uuid
-
-from datetime import date
-from typing import Any, Dict, List, Optional, Union
-
-from eodal.config import get_settings
-from eodal.core.sensors import Sentinel2
-from eodal.operational.mapping.mapper import Mapper, Feature
-from eodal.operational.mapping.merging import merge_datasets
-from eodal.utils.constants.sentinel2 import ProcessingLevels
-from eodal.utils.exceptions import (
-    InputError,
-    BlackFillOnlyError,
-    DataNotFoundError
-)
-from eodal.metadata.sentinel2.utils import identify_updated_scenes
-from eodal.core.raster import SceneProperties
-
-settings = get_settings()
-logger = settings.logger
-
-class Sentinel2Mapper(Mapper):
-    """
-    Spatial mapper class for Sentinel-2 data.
-
-    :attrib processing_level:
-        Sentinel-2 data processing level (L1C or L2A)
-    :attrib cloud_cover_threshold:
-        global (scene-wide) cloud coverage threshold between 0 and 100% cloud cover.
-        Scenes with cloud coverage reported higher than the threshold are discarded.
-        To obtain *all* mapper in the archive use the default of 100%.
-    :attrib use_latest_pdgs_baseline:
-        since a scene can possible occur in different PDGS baseline numbers (the scene_id
-        and product_uri will be different, which is supported by our data model)
-        it is necessary to decide for a baseline in those cases where multiple mapper
-        from the same sensing and data take time are available (originating from the
-        same Sentinel-2 data but processed slightly differently). By default we use
-        those mapper from the latest baseline. Otherwise, it is possible to use
-        the baseline most mapper were processed with.
-    """
-
-    def __init__(
-        self,
-        processing_level: ProcessingLevels,
-        cloud_cover_threshold: Optional[Union[int, float]] = 100,
-        use_latest_pdgs_baseline: Optional[bool] = True,
-        *args,
-        **kwargs,
-    ):
-        """
-        Initializes a new Sentinel-2 Mapper object
-
-        :param processing_level:
-            Sentinel-2 processing level to query
-        :param cloud_cover_threshold:
-            cloud cover threshold in percent (0-100%). Default is 100% to
-            consider all mapper in the archive
-        :param use_latest_pdgs_baseline:
-            if True (default) forces *eodal* to use the latest processing
-            baseline in case a scene is available in different processing levels
-        :param tile_selection:
-            optional list of Sentinel-2 tiles (e.g., ['T32TMT','T32TGM']) to use for
-            filtering. Only mapper belonging to these tiles are returned then
-        :param args:
-            arguments to pass to the constructor of the ``Mapper`` super class
-        :param kwargs:
-            key-word arguments to pass to the constructor of the ``Mapper``
-            super class
-        """
-        # initialize super-class
-        Mapper.__init__(self, *args, **kwargs)
-
-        object.__setattr__(self, "processing_level", processing_level)
-        object.__setattr__(self, "cloud_cover_threshold", cloud_cover_threshold)
-        object.__setattr__(self, "use_latest_pdgs_baseline", use_latest_pdgs_baseline)
-
-    def get_scenes(self) -> None:
-        """
-        Queries the Sentinel-2 metadata catalog for a selected time period and
-        feature collection.
-
-        NOTE:
-            By passing a list of Sentinel-2 tiles you can explicitly control
-            which Sentinel-2 tiles are considered. This might be useful for
-            mapper tasks where your feature collection lies completely within
-            a single Sentinel-2 tile but also overlaps with neighboring tiles.
-
-        The scene selection and processing workflow contains several steps:
-
-        1.  Query the metadata catalog for **ALL** available mapper that overlap
-            the bounding box of a given ``Polygon`` or ``MultiPolygon``
-            feature. **IMPORTANT**: By passing a list of Sentinel-2 tiles
-            to consider (``tile_ids``) you can explicitly control which
-            Sentinel-2 tiles are considered!
-        2.  Check if for a single sensing date several mapper are available
-        3.  If yes check if that's due to Sentinel-2 data take or tiling grid
-            design. If yes flag these mapper as potential merge candidates. A
-            second reason for multiple mapper are differences in PDGS baseline,
-            i.e., the dataset builds upon the **same** Sentinel-2 data but
-            was processed by different base-line version.
-        4.  If the mapper found have different spatial coordinate systems (CRS)
-            (usually different UTM zones) flag the data accordingly. The target
-            CRS is defined as that CRS the majority of mapper shares.
-        """
-        self._get_scenes(sensor='sentinel2')
-
-    def _resample_s2_scene(self, s2_scene) -> None:
-        """
-        Resamples the Sentinel-2 into a user-defined spatial resolution
-
-        :param s2_scene:
-            `~eodal.core.sentinel2.Sentinel2` object with loaded S2
-            data
-        """
-        # resample to target resolutionsorted_indices based on the MapperConfig settings
-        has_scl = False
-        band_selection = self.mapper_configs.band_names
-        if band_selection is None:
-            band_selection = s2_scene.band_names
-        # make sure SCL is always resampled to 10m using nearest neighbor
-        if "SCL" in band_selection:
-            band_selection.remove("SCL")
-            has_scl = True
-        s2_scene.resample(
-            band_selection=band_selection,
-            interpolation_method=self.mapper_configs.resampling_method,
-            target_resolution=self.mapper_configs.spatial_resolution,
-            inplace=True,
-        )
-        if has_scl:
-            s2_scene.resample(
-                band_selection=["SCL"],
-                interpolation_method=cv2.INTER_NEAREST_EXACT,
-                target_resolution=self.mapper_configs.spatial_resolution,
-                inplace=True,
-            )
-
-    def _read_multiple_scenes(
-        self, scenes_date: pd.DataFrame, feature_id: str, **kwargs
-    ) -> Union[gpd.GeoDataFrame, Sentinel2]:
-        """
-        Backend method for processing and reading scene data if more than one scene
-        is available for a given sensing date and feature (area of interest)
-
-        :param scenes_date:
-            `DataFrame` with all Sentinel-2 mapper of a single date
-        :param feature_id:
-            ID of the feature for which to extract data
-        :param kwargs:
-            optional key-word arguments to pass on to
-            `~eodal.core.sensors.Sentinel2.from_safe`
-        """
-        # check which baseline should be used
-        return_highest_baseline = kwargs.get("return_highest_baseline", True)
-        res = None
-
-        # get properties and geometry of the current feature from the collection
-        feature_dict = self.get_feature(feature_id)
-        feature_gdf = gpd.GeoDataFrame.from_features(feature_dict)
-        feature_gdf.crs = feature_dict["features"][0]["properties"]["epsg"]
-        # parse feature geometry in kwargs so that only a spatial subset is read
-        # in addition parse the S2 gain factor as "scale" argument
-        kwargs.update({"vector_features": feature_gdf})
-
-        # if the feature is a point we take the data set that is not blackfilled.
-        # If more than one data set is not blackfilled  we simply take the
-        # first data set
-        feature_geom = feature_gdf["geometry"].iloc[0].type
-        if feature_geom in ["Point", "MultiPoint"]:
-            for _, candidate_scene in scenes_date.iterrows():
-                if settings.USE_STAC:
-                    in_dir = candidate_scene["assets"]
-                else:
-                    in_dir = candidate_scene["real_path"]
-                feature_gdf = Sentinel2.read_pixels_from_safe(
-                    vector_features=feature_gdf,
-                    in_dir=in_dir,
-                    band_selection=self.mapper_configs.band_names,
-                )
-                # a empty data frame indicates black-fill
-                if feature_gdf.empty:
-                    continue
-                res = feature_gdf
-                # make sure the coordinates are reprojected if necessary
-                if res.crs.to_epsg() != scenes_date.target_crs.unique()[0]:
-                    res.to_crs(epsg=scenes_date.target_crs.unique()[0], inplace=True)
-
-                if isinstance(candidate_scene, (pd.Series, gpd.GeoSeries)):
-                    res["sensing_date"] = candidate_scene["sensing_date"]
-                    res['sensing_time'] = candidate_scene['sensing_time']
-                    res["scene_id"] = candidate_scene["scene_id"]
-                break
-        # in case of a (Multi-)Polygon: check if one of the candidate mapper complete
-        # contains the feature (i.e., its bounding box). If that's the case and the
-        # returned data is not black-filled, we can take that data set. If none of the
-        # candidate contains the scene complete, merging and (depending on the CRS)
-        # re-reprojection might be required. The result is then saved to disk in a temporary
-        # directory.
-        else:
-            # check processing baseline first (one dataset can appear in different processing
-            # baselines)
-            updated_scenes, old_scenes = identify_updated_scenes(
-                metadata_df=scenes_date, return_highest_baseline=return_highest_baseline
-            )
-            # only one scene left -> read the scene and return
-            if updated_scenes.shape[0] == 1:
-                if settings.USE_STAC:
-                    in_dir = updated_scenes["assets"].iloc[0]
-                else:
-                    in_dir = updated_scenes["real_path"].iloc[0]
-                # if there were only two input mapper we're done
-                # otherwise we have to check if we have to merge data
-                if scenes_date.shape[0] == 2:
-                    res = Sentinel2.from_safe(
-                        in_dir=in_dir,
-                        band_selection=self.mapper_configs.band_names,
-                        **kwargs,
-                    )
-                    self._resample_s2_scene(s2_scene=res)
-                    return res
-            # if updated mapper is not empty update the scenes_date DataFrame
-            if not updated_scenes.empty:
-                # drop "out-dated" mapper
-                appended = pd.concat([scenes_date, old_scenes])
-                appended.drop_duplicates(subset=['product_uri', 'tile_id'], keep=False, inplace=True)
-                scenes_date = appended.copy()
-                # if there is a single scene from a another tile in the
-                # "old" mapper append it to the scenes_date
-                old_scenes_grouped = old_scenes.groupby(by='tile_id')
-                for tile_scenes in old_scenes_grouped:
-                    if tile_scenes[0] not in scenes_date.tile_id.unique():
-                        scenes_date = pd.concat(
-                            [scenes_date, tile_scenes[1].copy()]
-                        )
-
-            # apply merge logic
-            tmp_fnames = []
-            scene_props = []
-            for _, candidate_scene in scenes_date.iterrows():
-                if settings.USE_STAC:
-                    in_dir = candidate_scene["assets"]
-                else:
-                    in_dir = candidate_scene["real_path"]
-                s2_scene = Sentinel2.from_safe(
-                    in_dir=in_dir,
-                    band_selection=self.mapper_configs.band_names,
-                    **kwargs,
-                )
-                scene_props.append(s2_scene.scene_properties)
-                self._resample_s2_scene(s2_scene=s2_scene)
-                # reproject the scene if its CRS is not the same as the target_crs
-                if (
-                    s2_scene[s2_scene.band_names[0]].geo_info.epsg
-                    != candidate_scene.target_crs
-                ):
-                    # make sure the pixel size remains the same after re-projection
-                    # to do so, construct an explicit affine transformation matrix.
-                    # Since all bands have the same spatial resolution this step can be
-                    # applied to all bands at the same time
-                    band = s2_scene.band_names[0]
-                    nodata = s2_scene[band].nodata
-                    geo_info_orig = s2_scene[band].geo_info
-                    pixres_x_dst = abs(geo_info_orig.pixres_x)
-                    pixres_y_dst = abs(geo_info_orig.pixres_y)
-                    s2_scene.reproject(
-                        band_selection=s2_scene.band_names,
-                        target_crs=candidate_scene.target_crs,
-                        dst_resolution=(pixres_x_dst, pixres_y_dst),
-                        dst_nodata=nodata,
-                        inplace=True,
-                    )
-                # write scene to temporary working directory and call rasterio.merge, returning
-                # the merged dataset (save as tif to avoid compression issues)
-                fname_scene = settings.TEMP_WORKING_DIR.joinpath(f"{uuid.uuid4()}.tif")
-                s2_scene.to_rasterio(fname_scene)
-                tmp_fnames.append(fname_scene)
-
-            # merge datasets
-            vector_features = kwargs.get("vector_features", None)
-            band_options = {
-                "band_names_dst": s2_scene.band_names,
-                "band_aliases": s2_scene.band_aliases,
-            }
-            # adopt scene properties
-            scene_dict = scene_props[0].__dict__
-            scene_dict_keys = [x[1::] for x in scene_dict.keys()]
-            scene_dict_vals = list(scene_dict.values())
-            scene_dict = dict(zip(scene_dict_keys, scene_dict_vals))
-            # combine product_uri's from the single datasets by '&'
-            new_product_uri = "".join([x.product_uri + "&" for x in scene_props])[:-1]
-            scene_dict.update({"product_uri": new_product_uri})
-            scene_properties = SceneProperties(**scene_dict)
-
-            try:
-                res = merge_datasets(
-                    datasets=tmp_fnames,
-                    vector_features=vector_features,
-                    band_options=band_options,
-                    scene_properties=scene_properties,
-                    sensor="sentinel2",
-                )
-            except Exception as e:
-                raise ValueError(f"Could not merge Sentinel-2 datasets: {e}")
-            # clean up working directory
-            for tmp_fname in tmp_fnames:
-                os.remove(tmp_fname)
-        return res
-
-    def get_observation(
-        self, feature_id: Any, sensing_date: date, **kwargs
-    ) -> Union[gpd.GeoDataFrame, Sentinel2, None]:
-        """
-        Returns the scene data (observations) for a selected feature and date.
-
-        If for the date provided no mapper are found, the data from the scene(s)
-        closest in time is returned
-
-        :param feature_id:
-            identifier of the feature for which to extract observations
-        :param sensing_date:
-            date for which to extract observations (or the closest date if
-            no observations are available for the given date)
-        :param kwargs:
-            optional key-word arguments to pass on to
-            `~eodal.core.sensors.Sentinel2.from_safe`
-        :returns:
-            depending on the geometry type of the feature either a
-            ``GeoDataFrame`` (geometry type: ``Point``) or ``Sentinel2Handler``
-            (geometry types ``Polygon`` or ``MultiPolygon``) is returned. if
-            the observation contains nodata, only, None is returned.
-        """
-        # call super class method for getting the observation
-        res = self._get_observation(feature_id=feature_id, sensing_date=sensing_date,
-                                    sensor='sentinel2', **kwargs)
-        # for multiple mapper a Sentinel-2 specific class must be called
-        if isinstance(res, tuple):
-            _, scenes_date, _ = res
-            res = self._read_multiple_scenes(
-                scenes_date=scenes_date, feature_id=feature_id, **kwargs
-            )
-        return res
-
-    def get_complete_timeseries(
-        self, feature_selection: Optional[List[Any]] = None,
-        drop_blackfilled_scenes: Optional[bool] = True, **kwargs
-    ) -> Dict[Any, Union[gpd.GeoDataFrame, List[Sentinel2]]]:
-        """
-        Extracts all observation with a time period for a feature collection.
-
-        This function takes the Sentinel-2 mapper retrieved from the metadata DB query
-        in `~Mapper.get_sentinel2_scenes` and extracts the Sentinel-2 data from the
-        original .SAFE archives for all available mapper.
-
-        :param feature_selection:
-            optional subset of features ids (you can only select features included
-            in the current feature collection)
-        :param drop_blackfilled_scenes:
-            drop mapper having no data values only (default)
-        :param kwargs:
-            optional key-word arguments to pass to `~eodal.core.band.Band.from_rasterio`
-        """
-        assets = {}
-        # check if band selection is passed in kwargs
-        band_selection_kwargs = kwargs.get('band_selection', None)
-        if band_selection_kwargs is not None:
-            object.__setattr__(self, 'band_names', band_selection_kwargs)
-        # loop over features (AOIs) in feature dict
-        for feature, scenes_df in self.observations.items():
-            # in case a feature selection is available check if the current
-            # feature is part of it
-            if feature_selection is not None:
-                if feature not in feature_selection:
-                    continue
-
-            # loop over mapper, they are already ordered by date (ascending)
-            # and check for each date which mapper are relevant and require
-            # potential reprojection or merging
-            sensing_dates = scenes_df.sensing_date.unique()
-            n_sensing_dates = len(sensing_dates)
-            feature_res = []
-            for idx, sensing_date in enumerate(sensing_dates):
-                try:
-                    res = self.get_observation(feature, sensing_date, **kwargs)
-                    if isinstance(res, gpd.GeoDataFrame):
-                        if res.empty: continue
-                    if drop_blackfilled_scenes:
-                        if hasattr(res, 'is_blackfilled'):
-                            if res.is_blackfilled:
-                                logger.info(
-                                    f"Feature {feature}: "
-                                    f"Skipped data due to blackfill from {sensing_date} "
-                                    f"({idx+1}/{n_sensing_dates})"
-                                )
-                                continue
-                    feature_res.append(res)
-                    logger.info(
-                        f"Feature {feature}: "
-                        f"Extracted data from {sensing_date} "
-                        f"({idx+1}/{n_sensing_dates})"
-                    )
-                except Exception as e:
-                    logger.error(
-                        f"Feature {feature}: "
-                        f"Extracting data from {sensing_date} "
-                        f"({idx+1}/{n_sensing_dates}) failed: {e}"
-                    )
-                    continue
-            # if res is a GeoDataFrame the list can be concated
-            if isinstance(res, gpd.GeoDataFrame):
-                assets[feature] = pd.concat(feature_res)
-            else:
-                # order mapper by acquisition time
-                timestamps = [x.scene_properties.acquisition_time for x in feature_res]
-                sorted_indices = np.argsort(np.array(timestamps))
-                feature_res_ordered = [feature_res[idx] for idx in sorted_indices]
-                assets[feature] = feature_res_ordered
-                
-        return assets
diff --git a/eodal/operational/mapping/utils.py b/eodal/operational/mapping/utils.py
deleted file mode 100644
index 8376bd12..00000000
--- a/eodal/operational/mapping/utils.py
+++ /dev/null
@@ -1,115 +0,0 @@
-'''
-Mapper class utility functions.
-
-Copyright (C) 2022 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/>.
-'''
-from __future__ import annotations
-
-import matplotlib.pyplot as plt
-import numpy as np
-
-from typing import Dict, List, Optional
-
-from eodal.core.raster import RasterCollection
-
-def plot_feature(feature_scenes: List[RasterCollection], band_selection: str | List[str],
-                 max_scenes_in_row: Optional[int] = 6,
-                 eodal_plot_kwargs: Optional[Dict] = {},
-                 feature_name: Optional[str] = '',
-                 **kwargs) -> plt.Figure:
-        """
-        Plots all mapper retrieved for a feature
-
-        :param band_selection:
-            selection of band(s) to use for plotting. Must be either a single
-            band or a set of three bands
-        :returns:
-            `Figure` object
-        """
-        # check number of passed bands
-        if isinstance(band_selection, str):
-            band_selection = [band_selection]
-        if not len(band_selection) == 1 and not len(band_selection) == 3:
-            raise ValueError('You must pass a single band name or three band names')
-
-        plot_multiple_bands = True
-        if len(band_selection) == 1:
-            plot_multiple_bands = False
-
-        # check number of mapper in feature_scenes and determine figure size
-        n_scenes = len(feature_scenes)
-        nrows = 1
-        ncols = 1
-        if n_scenes == 0:
-            raise ValueError('No mapper available for plotting')
-        elif n_scenes == 1:
-            f, ax = plt.subplots(**kwargs)
-            # cast to array to allow indexing
-            ax = np.array([ax]).reshape(1,1)
-        else:
-            if n_scenes <= max_scenes_in_row:
-                ncols = n_scenes
-                f, ax = plt.subplots(ncols=ncols, nrows=nrows, **kwargs)
-                # reshape to match the shape of ax array with >1 rows
-                ax = ax.reshape(1, ax.size)
-            else:
-                nrows = int(np.ceil(n_scenes / max_scenes_in_row))
-                ncols = max_scenes_in_row
-                f, ax = plt.subplots(ncols=ncols, nrows=nrows, **kwargs)
-
-        # get acquisition times of the mapper if available. If not label the
-        # plots by ascending numbers (Scene 1, Scene 2, Scene 3,...)
-        scene_labels = [
-            f'{x.scene_properties.acquisition_time} {x.scene_properties.platform}' \
-            if hasattr(x, 'scene_properties')  else f'Scene {idx+1}' for idx, x in \
-            enumerate(feature_scenes)
-        ]
-
-        row_idx, col_idx = 0, 0
-        for idx, feature_scene in enumerate(feature_scenes):
-            if plot_multiple_bands:
-                feature_scene.plot_multiple_bands(
-                    band_selection=band_selection,
-                    ax=ax[row_idx, col_idx],
-                    **eodal_plot_kwargs
-                )
-            else:
-                feature_scene[band_selection[0]].plot(
-                    ax=ax[row_idx, col_idx],
-                    **eodal_plot_kwargs
-                )
-            ax[row_idx, col_idx].set_title(scene_labels[idx])
-
-            # switch off axes labels if sharex == True and sharey=True
-            if kwargs.get('sharex', False) and kwargs.get('sharey', False) \
-            and n_scenes > 1:
-                if nrows > 1:
-                    if row_idx < (nrows - 1): ax[row_idx, col_idx].set_xlabel('')
-                if nrows > 1:
-                    if col_idx > 0: ax[row_idx, col_idx].set_ylabel('')
-
-            # increase column (and row) counter accordingly
-            col_idx += 1
-            # begin a new row when all columns are filled
-            if col_idx == max_scenes_in_row:
-                col_idx = 0
-                row_idx += 1
-
-        # make sure sub-plot labels do not overlap
-        f.tight_layout()
-        # figure title
-        f.suptitle(feature_name)
-        return f
diff --git a/eodal/operational/resampling/__init__.py b/eodal/operational/resampling/__init__.py
deleted file mode 100644
index e69de29b..00000000
diff --git a/eodal/operational/resampling/sentinel2/__init__.py b/eodal/operational/resampling/sentinel2/__init__.py
deleted file mode 100644
index bf05ff8f..00000000
--- a/eodal/operational/resampling/sentinel2/__init__.py
+++ /dev/null
@@ -1,275 +0,0 @@
-"""
-Module that could be used in (parallelized) way for pre-processing Sentinel-2 data
-as part of the ``eodal`` operational processing pipeline:
-             
-- resampling from 20 to 10m spatial resolution
-- generation of a RGB preview image per scene
-- merging of split mapper (due to data take issues)
-- resampling of SCL data (L2A processing level, only) and generation of a preview
-- generation of metadata of the processed data (links to the input datasets)
-
-The module relies on metadata from eodal's metadata DB.
-
-Depending on your hardware this might take a while!
-
-Copyright (C) 2022 Gregor Perich & 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 time
-import pandas as pd
-from datetime import date
-from pathlib import Path
-from joblib import Parallel, delayed
-from typing import Optional
-from typing import Tuple
-from typing import Dict
-from datetime import datetime
-
-from .resample_and_stack import resample_and_stack_s2
-from .merge_blackfill import merge_split_scenes
-
-from eodal.config import get_settings
-from eodal.metadata.sentinel2.database.querying import find_raw_data_by_tile
-from eodal.metadata.utils import reconstruct_path
-from eodal.operational.resampling.utils import identify_split_scenes
-from eodal.utils.constants import ProcessingLevels
-
-
-Settings = get_settings()
-logger = Settings.logger
-
-
-def do_parallel(
-    in_df: pd.DataFrame, loopcounter: int, out_dir: Path, **kwargs
-) -> Dict[str, Path]:
-    """
-    Wrapper function for (potential) parallel execution of S2 resampling & stacking,
-    and SCL resampling (L2A, only).
-
-    Returns a dict containing the file-paths to the generated datasets.
-
-    :param in_df:
-        dataframe containing metadata of S2 mapper (must follow AgripySat convention)
-    :param loopcounter:
-        Index to get actual S2 scene in loop.
-    :param out_dir:
-        path where the bandstacks (and the SCL files) get written to.
-    :param **kwargs:
-        Optional key-word arguments to pass to ``resample_and_stack_s2``
-    :return innerdict:
-        dictionary with file-paths to generated datasets returned from
-        ``resample_and_stack_s2``
-    """
-
-    innerdict = {}
-
-    try:
-        # reconstruct storage location based on DB records
-        in_dir = reconstruct_path(record=in_df.iloc[loopcounter])
-        innerdict = resample_and_stack_s2(in_dir=in_dir, out_dir=out_dir, **kwargs)
-        # preserve scene_id and product_uri
-        in_dir["scene_id"] = in_df.scene_id.iloc[loopcounter]
-        in_dir["product_uri"] = in_df.product_uri.iloc[loopcounter]
-
-    except Exception as e:
-        logger.error(e)
-        return innerdict
-
-    # write to successful_scenes.txt to indicate that the processing of
-    # the scene was accomplished without any errors
-    scenes_log_file = out_dir.joinpath("log").joinpath(
-        Settings.PROCESSING_CHECK_FILE_NO_BF
-    )
-    creation_time = datetime.now().strftime("%Y%m%d-%H%M%S")
-    with open(scenes_log_file, "a+") as src:
-        line = ""
-        for key in innerdict.keys():
-            line += str(innerdict[key]) + ","
-        line += creation_time
-        src.write(line + "\n")
-
-    return innerdict
-
-
-def exec_pipeline(
-    processed_data_archive: Path,
-    date_start: date,
-    date_end: date,
-    tile: str,
-    processing_level: Optional[ProcessingLevels] = ProcessingLevels.L2A,
-    n_threads: Optional[int] = 1,
-    **kwargs,
-) -> Tuple[pd.DataFrame, pd.DataFrame]:
-    """
-    Execution of the Sentinel-2 pre-processing pipeline, including:
-
-    -> resampling from 20 to 10m spatial resolution
-    -> generation of a RGB preview image per scene
-    -> merging of split mapper (due to data take issues)
-    -> resampling of SCL data (L2A processing level, only)
-
-    NOTE:
-        The function works on single tiles and supports L1C and L2A processing level.
-
-    NOTE:
-        The main part of the pipeline except the blackfill merging can be executed
-        using multiple threads if ``n_threads`` is set to a value larger 1. This feature
-        is experimental.
-
-    :param processed_data_archive:
-        target archive where the processed data should be stored.
-    :param date_start:
-        start date of the period to process. NOTE: The temporal
-        selection is bound by the available (i.e., downloaded) Sentinel-2 data!
-    :param date_end:
-        end_date of the period to process. NOTE: The temporal
-        selection is bound by the available (i.e., downloaded) Sentinel-2 data!
-    :param tile:
-        Sentinel-2 tile to process (e.g. 'T32TLT'). NOTE: The spatial (tile)
-        selection is bound by the available (i.e., downloaded) Sentinel-2 data!
-    :param processing_level:
-        Sentinel-2 processing level (either L1C or L2A, L2A by default).
-    :param n_threads:
-        number of threads to use for execution of the preprocessing depending on
-        your computer hardware. Per default set to 1 (no parallelization).
-    :param kwargs:
-        kwargs to pass to resample_and_stack_S2 and scl_10m_resampling (L2A, only)
-    :return:
-        tuple containing metadata of the processed datasets [0] and eventually
-        failed datasets in the second tuple item [1]
-    """
-
-    # make sub-directory for logging successfully processed mapper in out_dir
-    log_dir = processed_data_archive.joinpath("log")
-    if not log_dir.exists():
-        log_dir.mkdir(parents=True, exist_ok=True)
-
-    # query metadata DB for available Sentinel-2 mapper
-    try:
-        metadata = find_raw_data_by_tile(
-            date_start=date_start,
-            date_end=date_end,
-            processing_level=processing_level,
-            tile=tile,
-        )
-    except Exception as e:
-        logger.error(f"Could not find Sentinel-2 data: {e}")
-        return
-
-    # get "duplicates", i.e, scenes that have the same sensing date because of datastrip
-    # beginning/end issue
-    num_scenes = metadata.shape[0]
-    meta_blackfill = identify_split_scenes(metadata_df=metadata)
-
-    # exclude these duplicated mapper from the main (parallelized) workflow!
-    metadata = metadata[~metadata.product_uri.isin(meta_blackfill["product_uri"])]
-    if meta_blackfill.empty:
-        logger.info(f"Found {num_scenes} mapper out of which 0 must be merged")
-    else:
-        logger.info(
-            f"Found {num_scenes} mapper out of which {meta_blackfill.shape[0]} must be merged"
-        )
-
-    t = time.time()
-    result = Parallel(n_jobs=n_threads)(
-        delayed(do_parallel)(
-            in_df=metadata, loopcounter=idx, out_dir=processed_data_archive, **kwargs
-        )
-        for idx in range(metadata.shape[0])
-    )
-    logger.info(f"Execution time: {time.time()-t} seconds.")
-
-    # concatenate the metadata of the stacked image files into a pandas dataframe
-    bandstack_meta = pd.DataFrame(result)
-
-    # merge black-fill mapper (data take issue) if any
-    if not meta_blackfill.empty:
-        logger.info("Starting merging of blackfill mapper")
-
-        # after regular scene processsing, process the blackfill mapper single-threaded
-        for date in meta_blackfill.sensing_date.unique():
-            scenes = meta_blackfill[meta_blackfill.sensing_date == date]
-            product_id_1 = scenes.product_uri.iloc[0]
-            product_id_2 = scenes.product_uri.iloc[1]
-            # reconstruct storage location
-            raw_data_archive = reconstruct_path(record=scenes.iloc[0]).parent
-            scene_1 = raw_data_archive.joinpath(product_id_1)
-            scene_2 = raw_data_archive.joinpath(product_id_2)
-            logger.info(f"Working on {scene_1} and {scene_2}")
-
-            # the usual naming problem with the .SAFE directories
-            if not scene_1.exists():
-                # product_id_1 = Path(str(product_id_1.replace('.SAFE', '')))
-                scene_1 = Path(str(scene_1).replace(".SAFE", ""))
-            if not scene_2.exists():
-                # product_id_2 = Path(str(product_id_2.replace('.SAFE', '')))
-                scene_2 = Path(str(scene_2).replace(".SAFE", ""))
-
-            try:
-                res = merge_split_scenes(
-                    scene_1=scene_1,
-                    scene_2=scene_2,
-                    out_dir=processed_data_archive,
-                    **kwargs,
-                )
-                res.update(
-                    {
-                        "scene_was_merged": True,
-                        "spatial_resolution": kwargs.get("target_resolution", 10),
-                    }
-                )
-                # check if bandstack_meta is already populated, otherwise create it
-                if bandstack_meta.empty:
-                    bandstack_meta = pd.DataFrame.from_records([res])
-                else:
-                    bandstack_meta = bandstack_meta.append(res, ignore_index=True)
-
-                # write to log file
-                scenes_log_file = processed_data_archive.joinpath("log").joinpath(
-                    Settings.PROCESSING_CHECK_FILE_BF
-                )
-                creation_time = datetime.now().strftime("%Y%m%d-%H%M%S")
-                with open(scenes_log_file, "a+") as src:
-                    line = f"{str(scene_1)}, {str(scene_2)}, {creation_time}"
-                    src.write(line + "\n")
-
-            except Exception as e:
-                logger.error(f"Failed to merge {scene_1} and {scene_2}: {e}")
-
-        # also the storage location shall be inserted into the database later
-        # bandstack_meta['storage_share'] = target_s2_archive
-
-        logger.info("Finished merging of blackfill mapper")
-
-    # check for any empty bandstack paths (indicates that something went wrong
-    # during the processing
-    emtpy_bandstacks = bandstack_meta.loc[
-        (bandstack_meta.bandstack == None) | (bandstack_meta.bandstack == "")
-    ]
-
-    # in this remove the errored records from bandstack_meta and reported the
-    # failed records
-    if not emtpy_bandstacks.empty:
-        logger.error(f"Resampling failed for {emtpy_bandstacks.shape[0]} datasets")
-        errored_records = pd.merge(
-            bandstack_meta, emtpy_bandstacks, how="inner", on=["scene_id"]
-        )["scene_id"]
-        bandstack_meta = bandstack_meta[
-            ~bandstack_meta["scene_id"].isin(errored_records)
-        ]
-
-    # save metadata of all stacked files and return
-    return bandstack_meta, emtpy_bandstacks
diff --git a/eodal/operational/resampling/sentinel2/merge_blackfill.py b/eodal/operational/resampling/sentinel2/merge_blackfill.py
deleted file mode 100644
index 61017efe..00000000
--- a/eodal/operational/resampling/sentinel2/merge_blackfill.py
+++ /dev/null
@@ -1,232 +0,0 @@
-"""
-Sentinel-2 records data in so-called datatakes. When a datatake is over and a new
-begins the acquired image data is written to different files (based on the datatake
-time). Sometimes, this cause mapper of a single acquisition date to be split into
-two datasets which differ in their datatake. Thus, both datasets have a reasonable
-amount of blackfill in those areas not covered by the datatake they belong to. For
-users of satellite data, however, it is much more convenient to have those split
-mapper merged into one since the division into two mapper by the datatake has
-technical reasons only.
-
-Copyright (C) 2022 Gregor Perich & 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/>.
-"""
-
-from __future__ import annotations
-
-import shutil
-import rasterio as rio
-
-from pathlib import Path
-from typing import Dict
-
-from eodal.config import get_settings
-from eodal.core.sensors import Sentinel2
-from eodal.core.raster import RasterCollection
-from eodal.operational.mapping.merging import merge_datasets
-from eodal.operational.resampling.sentinel2.resample_and_stack import (
-    resample_and_stack_s2,
-)
-from eodal.operational.resampling.sentinel2.resample_and_stack import create_rgb_preview
-from eodal.operational.resampling.sentinel2.resample_and_stack import create_scl_preview
-from eodal.utils.constants.sentinel2 import s2_band_mapping
-
-Settings = get_settings()
-logger = Settings.logger
-
-
-def merge_split_scenes(
-    scene_1: Path, scene_2: Path, out_dir: Path, **kwargs
-) -> Dict[str, Path]:
-    """
-    merges two Sentinel-2 datasets in .SAFE formatof the same sensing date and tile
-    split by the datatake beginning/ end.
-
-    First, both mapper are resampled to 10m and stacked in a temporary working directory;
-    second, they are merged together so that the blackfill of the first mapper is replaced
-    by the values of the second one.
-    SCL (if available) and previews are managed accordingly.
-
-    :param scene_1:
-        .SAFE directory containing the first of two mapper split by the datatake beginning/
-        end of Sentinel-2
-    :param scene_2:
-        .SAFE directory containing the first of two mapper split by the datatake beginning/
-        end of Sentinel-2
-    :param out_dir:
-        directory where to save the final outputs to. In this directory a temporary
-        working directory for storing intermediate products is generated.
-    :param kwargs:
-        key word arguments to pass to resample_and_stack_S2.
-    :return:
-        dictionary with filepaths to bandstack, rgb_quicklook, and (L2A, only) SCL
-    """
-
-    # check if working directory exists
-    working_dir = out_dir.joinpath("temp_blackfill")
-    if not working_dir.exists():
-        working_dir.mkdir()
-
-    # save the outputs of the two mapper to different sub-directories within the working
-    # directory to avoid to override the output
-    out_dirs = [working_dir.joinpath("1"), working_dir.joinpath("2")]
-    for _out_dir in out_dirs:
-        if _out_dir.exists():
-            # for clean start of next loop iteration
-            shutil.rmtree(_out_dir)
-        _out_dir.mkdir()
-
-    # do the spatial resampling for the two mapper
-    # first scene
-    try:
-        scene_out_1 = resample_and_stack_s2(
-            in_dir=scene_1, out_dir=out_dirs[0], **kwargs
-        )
-    except Exception as e:
-        logger.error(f"Resampling of {scene_1} failed: {e}")
-        return {}
-
-    # second scene
-    try:
-        scene_out_2 = resample_and_stack_s2(
-            in_dir=scene_2, out_dir=out_dirs[1], **kwargs
-        )
-    except Exception as e:
-        logger.error(f"Resampling of {scene_1} failed: {e}")
-        return {}
-
-    # merge band stacks
-    # merged file has some name as the first scene and is directly written to target directory
-    out_file = out_dir.joinpath(scene_out_1["bandstack"].name)
-    # we create a temporary file first
-    out_file_temp = working_dir.joinpath(out_file.name.split(".")[0] + "_temp.jp2")
-    datasets = [scene_out_1["bandstack"], scene_out_2["bandstack"]]
-
-    logger.info(f"Starting merging files {datasets[0]} and {datasets[1]}")
-
-    # merge reflectance data
-    try:
-        merge_datasets(datasets=datasets, out_file=out_file_temp)
-    except Exception as e:
-        logger.error(f"Merging of {datasets[0]} and {datasets[1]} failed: {e}")
-        return {}
-
-    out_band_names = list(s2_band_mapping.keys())
-    out_band_names.remove("SCL")
-    if kwargs.get("skip_60m_bands", True):
-        out_band_names.remove("B01")
-        out_band_names.remove("B09")
-
-    # set correct band descriptions and write final image to target directory
-    meta = rio.open(out_file_temp, "r").meta
-    meta.update({"QUALITY": "100", "REVERSIBLE": "YES"})
-
-    with rio.open(out_file, "w", **meta) as dst:
-
-        # open target file and write band by band + setting band names
-        with rio.open(out_file_temp, "r") as src:
-            for idx in range(meta["count"]):
-                dst.set_band_description(idx + 1, out_band_names[idx])
-                band_data = src.read(idx + 1)
-                dst.write(band_data, idx + 1)
-
-    logger.info(f"Merged datasets into {out_file}")
-
-    # generate RGB preview image from merged dataset (and SCL if applicable)
-    vis_bands = ["B02", "B03", "B04"]
-    handler = RasterCollection.from_multi_band_raster(
-        fpath_raster=out_file,
-        band_names_src=vis_bands,
-        band_names_dst=["blue", "green", "red"],
-    )
-    fname_rgb_preview = create_rgb_preview(
-        out_dir=out_dir, reader=handler, out_filename=scene_out_1["rgb_preview"].name
-    )
-    handler = None
-
-    fnames_out = {
-        "bandstack": out_file,
-        "rgb_preview": fname_rgb_preview,
-        "scl": "",
-        "scl_preview": "",
-    }
-
-    # check if SCL must be merged as well (L2A)
-    if "scl" in scene_out_1.keys():
-
-        scl_datasets = [scene_out_1["scl"], scene_out_2["scl"]]
-        scl_out_dir = out_dir.joinpath(Path(Settings.SUBDIR_SCL_FILES))
-
-        # create SCL directory if not yet there
-        if not scl_out_dir.exists():
-            scl_out_dir.mkdir()
-
-        out_file_scl = scl_out_dir.joinpath(scene_out_1["scl"].name)
-
-        try:
-            merge_datasets(datasets=scl_datasets, out_file=out_file_scl)
-        except Exception as e:
-            logger.error(
-                f"Could not merge {scl_datasets[0]} and {scl_datasets[1]}: {e}"
-            )
-            return fnames_out
-
-        logger.info(f"Merged SCL files into {out_file_scl}")
-        fnames_out["scl"] = out_file_scl
-
-        # finally plot the scl file
-        try:
-            handler = Sentinel2().from_multi_band_raster(
-                fpath_raster=out_file_scl, band_idx=1, band_names_dst=["SCL"]
-            )
-
-            fname_scl_preview = fname_rgb_preview.parent.joinpath(
-                scene_out_1["scl_preview"].name
-            )
-            create_scl_preview(
-                out_dir=out_dir, reader=handler, out_filename=fname_scl_preview
-            )
-            handler = None
-        except Exception as e:
-            logger.error(f"Could not generate SCL preview from {out_file}: {e}")
-            return fnames_out
-
-        fnames_out["scl_preview"] = fname_scl_preview
-
-    # remove working directory
-    try:
-        shutil.rmtree(working_dir)
-    except Exception as e:
-        logger.error(f"Could not remove working directory {working_dir}: {e}")
-
-    return fnames_out
-
-
-if __name__ == "__main__":
-
-    sat_dir = Path(
-        "/home/graflu/public/Evaluation/Satellite_data/Sentinel-2/Rawdata/L2A/CH/2018"
-    )
-
-    scene_1 = sat_dir.joinpath(
-        "S2A_MSIL2A_20180518T104021_N0207_R008_T32TLT_20180518T124554"
-    )
-    scene_2 = sat_dir.joinpath(
-        "S2A_MSIL2A_20180518T104021_N0207_R008_T32TLT_20180518T125548"
-    )
-
-    out_dir = Path("/mnt/ides/Lukas/debug/Processed")
-
-    fnames_out = merge_split_scenes(scene_1, scene_2, out_dir)
diff --git a/eodal/operational/resampling/sentinel2/resample_and_stack.py b/eodal/operational/resampling/sentinel2/resample_and_stack.py
deleted file mode 100644
index 9b87cb9d..00000000
--- a/eodal/operational/resampling/sentinel2/resample_and_stack.py
+++ /dev/null
@@ -1,421 +0,0 @@
-"""
-This function is part of the eodal's Sentinel-2 default processing chain.
-
-It allows for spatial resampling of Sentinel-2 images to bring
-the 10 and 20m spectral bands into either 10m (default) or 20m spatial resolution.
-The function works on the full extent of a Sentinel-2 scene and writes
-the output bands into a single multi-band geoTiff file.
-
-For resampling spatial subsets of a scene which are considerably smaller than
-the spatial extent of Sentinel-2 scene (almost 110km by 110km) consider using
-``eodal.core.sentinel.S2_Band_Reader`` and its ``resample()`` method directly.
-
-Copyright (C) 2022 Gregor Perich & 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/>.
-"""
-
-from __future__ import annotations
-
-import cv2
-import rasterio as rio
-import matplotlib.pyplot as plt
-
-from pathlib import Path
-from typing import Dict
-from typing import List
-from typing import Optional
-from typing import Union
-
-from eodal.config import get_settings
-from eodal.config.sentinel2 import Sentinel2 as s2
-from eodal.core.raster import RasterCollection
-from eodal.core.sensors.sentinel2 import Sentinel2
-from eodal.core.utils.sentinel2 import read_s2_tcifile
-from eodal.core.utils.sentinel2 import read_s2_sclfile
-from eodal.utils.constants.sentinel2 import s2_band_mapping
-from eodal.utils.constants.sentinel2 import ProcessingLevels
-from eodal.utils.sentinel2 import get_S2_processing_level
-
-Settings = get_settings()
-logger = Settings.logger
-S2 = s2()
-
-
-def _get_output_file_names(
-    in_dir: Path, resampling_method: str, target_resolution: Union[int, float]
-) -> Dict[str, str]:
-    """
-    auxiliary method to get the output file names
-    for the band-stack, the quicklooks and (if applicable) the
-    SCL.
-
-    The file-naming convention for the band-stack and previews is
-    ``
-    <date>_<tile>_<processing_level>_<sensor>_<resampling_method>_<spatial_resolution>
-    ``
-
-    :param in_dir:
-        path of the .SAFE directory where the S2 data resides.
-    :param resampling_method:
-        name of the resampling method used
-    :param target_resolution:
-        spatial resolution of the output
-    :return:
-        dict with output file names
-    """
-    # get S2 UID
-    s2_uid = in_dir.name
-
-    splitted = s2_uid.split("_")
-    date = splitted[2].split("T")[0]
-    tile = splitted[-2]
-    level = splitted[1]
-    sensor = splitted[0]
-    resolution = f"{int(target_resolution)}m"
-
-    # define filenames
-    basename = (
-        date
-        + "_"
-        + tile
-        + "_"
-        + level
-        + "_"
-        + sensor
-        + "_"
-        + resampling_method
-        + "_"
-        + resolution
-    )
-
-    return {
-        "bandstack": f"{basename}.jp2",
-        "rgb_preview": f"{basename}.png",
-        "scl_preview": f"{basename}_SCL.png",
-        "scl": f"{basename}_SCL.tiff",
-    }
-
-
-def _get_resampling_name(resampling_method: int) -> str:
-    """
-    auxiliary method to map opencv's integer codes to meaningful resampling names
-    (unknown if the method is not known)
-
-    :param resampling_method:
-        integer code from opencv2 for one of its image resizing methods
-    :return:
-        resampling method name or 'unknown' if the integer code cannot be
-        translated
-    """
-    translator = {
-        0: "nearest",
-        1: "linear",
-        2: "cubic",
-        3: "area",
-        4: "lanczos",
-        5: "linear-exact",
-        6: "nearest-exact",
-    }
-
-    return translator.get(resampling_method, "unknown")
-
-
-def create_rgb_preview(
-    out_dir: Path, reader: RasterCollection, out_filename: str
-) -> Path:
-    """
-    Creates the RGB quicklook image (stored in a sub-directory).
-
-    :param out_dir:
-        directory where the band-stacked geoTiff are written to
-    :param s2_stack:
-        opened Sat_Data_Reader with 'tci' band
-    :param out_filename:
-        file name of the resulting RGB quicklook image (*.png)
-    """
-    # RGB previews are stored in their own sub-directory
-    rgb_subdir = out_dir.joinpath(Settings.SUBDIR_RGB_PREVIEWS)
-    if not rgb_subdir.exists():
-        rgb_subdir.mkdir()
-
-    out_file = rgb_subdir.joinpath(out_filename)
-
-    fig_rgb = reader.plot_multiple_bands(["red", "green", "blue"])
-    fig_rgb.savefig(fname=out_file, bbox_inches="tight", dpi=150)
-    plt.close(fig_rgb)
-
-    return out_file
-
-
-def create_scl_preview(out_dir: Path, reader: Sentinel2, out_filename: str) -> Path:
-    """
-    Creates the SCL quicklook image (stored in a sub-directory).
-
-    :param out_dir:
-        directory where the band-stacked geoTiff are written to
-    :param s2_stack:
-        opened S2_Band_Reader with 'scl' band
-    :param out_filename:
-        file name of the resulting SCL quicklook image (*.png)
-    :return:
-        path to output file
-    """
-    # SCL previews are stored in their own sub-directory alongside with the RGBs
-    rgb_subdir = out_dir.joinpath(Settings.SUBDIR_RGB_PREVIEWS)
-    if not rgb_subdir.exists():
-        rgb_subdir.mkdir()
-
-    out_file = rgb_subdir.joinpath(out_filename)
-
-    fig_scl = reader.plot_scl()
-    fig_scl.savefig(fname=out_file, bbox_inches="tight", dpi=150)
-    plt.close(fig_scl)
-
-    return out_file
-
-
-def create_scl(out_dir: Path, reader: Sentinel2, out_filename: str) -> Path:
-    """
-    Creates the SCL raster datasets (stored in a sub-directory).
-
-    :param out_dir:
-        directory where the band-stacked geoTiff are written to
-    :param reader:
-        opened ``Sentinel2Handler`` with 'scl' band
-    :param out_filename:
-        file name of the resulting SCL raster image (*.tiff)
-    :return:
-        file-path of output dataset
-    """
-    scl_subdir = out_dir.joinpath(Settings.SUBDIR_SCL_FILES)
-    if not scl_subdir.exists():
-        scl_subdir.mkdir()
-
-    out_file = scl_subdir.joinpath(out_filename)
-    reader.to_rasterio(fpath_raster=out_file, band_selection=["SCL"])
-    return out_file
-
-
-TARGET_RESOLUTIONS: List[int] = [10, 20]
-
-
-def resample_and_stack_s2(
-    in_dir: Path,
-    out_dir: Path,
-    target_resolution: Optional[Union[int, float]] = 10,
-    interpolation_method: Optional[int] = cv2.INTER_NEAREST_EXACT,
-    skip_60m_bands: Optional[bool] = True,
-) -> Dict[str, Union[Path, str, int, float]]:
-    """
-    Function to spatially resample a S2 scene in *.SAFE format and write it to a
-    single, stacked geoTiff. Creates also a RGB preview png-file of the scene and
-    stores the scene classification layer that comes with L2A products in 10m spatial
-    resolution. Possible spatial resolution values are 10 and 20m. The 60m bands
-    are not considered.
-
-    The function checks the processing level of the data (L1C or L2A) based on the
-    name of the .SAFE dataset.
-
-    Depending on the processing level the output will look a bit differently:
-
-    * in L1C level (top-of-atmosphere) the band-stack of the spectral bands and
-      the RGB quicklook is produced
-    * in L2A level (bottom-of-atmosphere) the same inputs as in the L1C case
-      are generated PLUS the scene classification layer (SCL) resampled to 10m
-      spatial resolution
-
-    IMPORTANT: If only a small area of interest shall be processed, also consider
-    ``eodal.core.sensors.sentinel2.Sat_Data_Reader.resample`` since this function works
-    on the **full** scene extent, only.
-
-    :param in_dir:
-        path of the .SAFE directory where the S2 data resides.
-    :param out_dir:
-        path where to save the resampled & stacked geoTiff files to.
-    :param target_resolution:
-        target spatial resolution you want to resample to. Must be one of
-        [10, 20]. The default is 10 (meters).
-    :param interpolation_method:
-        The interpolation algorithm you want to use for spatial resampling.
-        The default is opencv's ``cv2.INTER_NEAREST_EXACT``. See the opencv documentation
-        for other options such as ``cv2.INTER_LINEAR``.
-    :param skip_60m_bands:
-        if False (default) does not resample the 60m bands. If True, also includes the
-        60m bands (B01, B09). B10 is never processed.
-    :returns:
-        dictionary with filepaths to bandstack, rgb_quicklook, and (L2A, only) SCL
-        and related metadata
-    """
-    # check passed spatial resolution
-    if target_resolution not in TARGET_RESOLUTIONS:
-        raise ValueError(f"Spatial resolution must be one out of {TARGET_RESOLUTIONS}")
-
-    # determine name of the resampling method
-    resampling_method_str = _get_resampling_name(resampling_method=interpolation_method)
-
-    # get output filenames
-    out_file_names = _get_output_file_names(
-        in_dir=in_dir,
-        resampling_method=resampling_method_str,
-        target_resolution=target_resolution,
-    )
-    # save resampling method and spatial resolution
-    out_file_names["resampling_method"] = resampling_method_str
-    out_file_names["spatial_resolution"] = target_resolution
-
-    # get the TCI quicklook image first
-    try:
-        tci = read_s2_tcifile(in_dir)
-    except Exception as e:
-        logger.error(f"Could not read TCI file from {in_dir}: {e}")
-        return {}
-    # save to file in RGB sub-directory
-    try:
-        out_file_rgb_preview = create_rgb_preview(
-            out_dir=out_dir, reader=tci, out_filename=out_file_names["rgb_preview"]
-        )
-        out_file_names.update({"rgb_preview": out_file_rgb_preview})
-    except Exception as e:
-        logger.error(f"Generation of RGB preview from {in_dir} failed: {e}")
-        return {}
-    logger.info(f"Generated RGB preview image from {in_dir}")
-
-    # its meta information serves a blue-print for writing the output
-    meta = tci["blue"].meta
-    tci = None
-
-    # get scene classification layer if available (L2A processing level)
-    processing_level = get_S2_processing_level(dot_safe_name=in_dir)
-    if processing_level == ProcessingLevels.L2A:
-        try:
-            scl = read_s2_sclfile(in_dir)
-        except Exception as e:
-            logger.error(f"Could not read SCL file from {in_dir}: {e}")
-            return {}
-        # save as preview image and geoTiff dataset after resampling to 10m
-        out_file_scl_preview = create_scl_preview(
-            out_dir=out_dir, reader=scl, out_filename=out_file_names["scl_preview"]
-        )
-        out_file_names.update({"scl_preview": out_file_scl_preview})
-        # resample to 10m spatial resolution and save to geoTiff
-        # we always use pixel division for the SCL file
-        try:
-            scl.resample(target_resolution=target_resolution, inplace=True)
-        except Exception as e:
-            logger.error(f"Resampling of SCL file from {in_dir} failed: {e}")
-            return {}
-        logger.info(f"Generated SCL preview image from {in_dir}")
-        try:
-            out_file_scl = create_scl(
-                out_dir=out_dir, reader=scl, out_filename=out_file_names["scl"]
-            )
-            out_file_names.update({"scl": out_file_scl})
-        except Exception as e:
-            logger.error(f"Generation of SCL file from {in_dir} failed: {e}")
-            return {}
-        logger.info(f"Generated SCL geoTiff (10m) from {in_dir}")
-    # if scl is not available set it to empty string in the dictionary of filenames
-    else:
-        out_file_names["scl"] = ""
-        out_file_names["scl_preview"] = ""
-
-    # get driver for output file
-    fname_bandstack = out_dir.joinpath(out_file_names["bandstack"])
-    out_file_names.update({"bandstack": fname_bandstack})
-    driver = rio.drivers.driver_from_extension(fname_bandstack)
-
-    # check if the resolution fits
-    if abs(meta["transform"][0]) != target_resolution:
-        if processing_level == ProcessingLevels.L2A:
-            # take meta from SCL layer since this layer has been resampled, already
-            meta = scl["SCL"].meta
-        else:
-            raise NotImplementedError(
-                "Uups - we need to implement this functionality first!"
-            )
-
-    # update meta dictionary
-    meta.update(
-        {
-            "count": 10,  # Sentinel-2 has 10 bands relevant for agriculture
-            "dtype": "uint16",  # Sentinel-2 data is stored as unsigned integers (16bit)
-            "driver": driver,
-            "REVERSIBLE": "YES",
-            "QUALITY": 100,
-        }
-    )
-    # open output dataset
-    s2_bands = list(s2_band_mapping.keys())
-    # SCL is a separate file (see above)
-    s2_bands.remove("SCL")
-    # 60m bands ignored if selected
-    if skip_60m_bands:
-        if "B01" in s2_bands:
-            s2_bands.remove("B01")
-        if "B09" in s2_bands:
-            s2_bands.remove("B09")
-
-    # write bands to the output dataset iteratively
-    with rio.open(fname_bandstack, "w+", **meta) as dst:
-
-        logger.info(f"Opening output dataset {fname_bandstack}")
-        # loop over S2 bands: read from .SAFE, resample if required and write to output
-        for idx, s2_band in enumerate(s2_bands):
-            try:
-                src = Sentinel2().from_safe(
-                    in_dir=in_dir,
-                    band_selection=[s2_band],
-                    read_scl=False,
-                    apply_scaling=False,
-                )
-            except Exception as e:
-                logger.error(f"Could not read band {s2_band} from {in_dir}: {e}")
-                return {}
-            # get color name of the band (used in reader)
-            band_alias = s2_band_mapping[s2_band]
-            # check resolution of the band
-            band_res = abs(src[s2_band].geo_info.pixres_x)
-            # resample band if required
-            if band_res != target_resolution:
-                try:
-                    src.resample(
-                        target_resolution=target_resolution,
-                        interpolation_method=interpolation_method,
-                        band_selection=[s2_band],
-                        inplace=True,
-                    )
-                except Exception as e:
-                    logger.error(
-                        f"Resampling of band {s2_band}/{band_alias} "
-                        f"({idx+1}/{len(s2_bands)} failed: {e}"
-                    )
-                    return {}
-
-                logger.info(
-                    f'Resampled band {s2_band} - "{band_alias.upper()}" '
-                    f"({idx+1}/{len(s2_bands)}) to {int(target_resolution)}m "
-                    f"from {in_dir}"
-                )
-            # write band to dst
-            dst.set_band_description(idx + 1, s2_band)
-            dst.write(src[s2_band].values, idx + 1)
-            logger.info(
-                f'Wrote band {s2_band} - "{band_alias.upper()}" ({idx+1}/{len(s2_bands)}) '
-                f"from {in_dir} into {fname_bandstack}"
-            )
-            src = None
-        logger.info(f"Finished writing bands to {fname_bandstack}")
-
-    return out_file_names
diff --git a/eodal/operational/resampling/utils.py b/eodal/operational/resampling/utils.py
deleted file mode 100644
index ff01796b..00000000
--- a/eodal/operational/resampling/utils.py
+++ /dev/null
@@ -1,40 +0,0 @@
-"""
-Some auxiliary functions for spatial resampling of raster data
-
-Copyright (C) 2022 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/>.
-"""
-from __future__ import annotations
-
-import pandas as pd
-
-
-def identify_split_scenes(
-    metadata_df: pd.DataFrame,
-) -> pd.DataFrame:
-    """
-    Returns entries in a pandas ``DataFrame`` retrieved from a query in eodal's
-    metadata base that have the same sensing date. This could indicate, e.g.,
-    that mapper have been split because of data take changes which sometimes cause
-    Sentinel-2 mapper to be split into two separate .SAFE archives, each of them
-    with a large amount of blackfill.
-
-    :param metadata_df:
-        dataframe from metadata base query in which to search for mapper with
-        the same sensing_date
-    :return:
-        mapper with the same sensing date (might also be empty)
-    """
-    return metadata_df[metadata_df.sensing_date.duplicated(keep=False)]
diff --git a/eodal/utils/__init__.py b/eodal/utils/__init__.py
index 762c245e..5a439f74 100644
--- a/eodal/utils/__init__.py
+++ b/eodal/utils/__init__.py
@@ -1,6 +1,5 @@
 from .stacking import stack_dataframes
 
-from .sentinel2 import s2
 from .sentinel2 import get_S2_bandfiles
 from .sentinel2 import get_S2_sclfile
 from .sentinel2 import get_S2_bandfiles_with_res
diff --git a/eodal/utils/constants/__init__.py b/eodal/utils/constants/__init__.py
index 7ccdc632..6e810e69 100644
--- a/eodal/utils/constants/__init__.py
+++ b/eodal/utils/constants/__init__.py
@@ -6,11 +6,11 @@ class ProcessingLevels(Enum):
     possible processing levels of remotely sensed data
     """
 
-    L1A = "LEVEL1A"
-    L1B = "LEVEL1B"
-    L1C = "LEVEL1C"
-    L2A = "LEVEL2A"
-    L2B = "LEVEL2B"
-    L3 = "LEVEL3"
-    L4 = "LEVEL4"
+    L1A = "Level-1A"
+    L1B = "Level-1B"
+    L1C = "Level-1C"
+    L2A = "Level-2A"
+    L2B = "Level-2B"
+    L3 = "Level-3"
+    L4 = "Level-4"
     UNKNOWN = "UNKNOWN"
diff --git a/eodal/utils/decorators.py b/eodal/utils/decorators.py
index a6df2b20..e1f10ffa 100644
--- a/eodal/utils/decorators.py
+++ b/eodal/utils/decorators.py
@@ -29,7 +29,6 @@
 from eodal.core.utils.geometry import multi_to_single_points
 from eodal.utils.exceptions import UnknownProcessingLevel, BandNotFoundError
 from eodal.utils.geometry import box_to_geojson
-from eodal.core.utils.geometry import multi_to_single_points
 
 Settings = get_settings()
 
@@ -57,21 +56,20 @@ def prepare_point_features(f):
     """
     @wraps(f)
     def wrapper(*args, **kwargs):
-        if len(args) >= 2:
-            vector_features = args[1]
-        else:
-            vector_features = kwargs.get('vector_features')
+        vector_features = kwargs.get('vector_features')
+        if vector_features is None:
+            vector_features = args[2]
         # cast to single point geometries
         try:
             vector_features_updated = multi_to_single_points(vector_features)
         except Exception as e:
             print(e)
-        if len(args) >= 2:
+        if 'vector_features' in kwargs.keys():
+            kwargs.update({'vector_features': vector_features_updated})
+        else:
             arg_list = list(args)
-            arg_list[1] = vector_features_updated
+            arg_list[2] = vector_features_updated
             args = tuple(arg_list)
-        else:
-            kwargs.update({'vector_features': vector_features_updated})
         return f(*args, **kwargs)
     return wrapper
 
diff --git a/eodal/utils/sentinel2.py b/eodal/utils/sentinel2.py
index 108426df..e3351dee 100644
--- a/eodal/utils/sentinel2.py
+++ b/eodal/utils/sentinel2.py
@@ -20,7 +20,6 @@
 
 from __future__ import annotations
 
-import os
 import glob
 import pandas as pd
 import planetary_computer
@@ -32,15 +31,11 @@
 
 from eodal.config import get_settings
 from eodal.config.stac_providers import STAC_Providers
-from eodal.config import Sentinel2
 from eodal.utils.exceptions import BandNotFoundError, MetadataNotFoundError
 from eodal.utils.constants.sentinel2 import ProcessingLevels
 
-# global definition of spectral bands and their spatial resolution
-s2 = Sentinel2()
 Settings = get_settings()
 
-
 def _url_to_safe_name(stac_asset: Union[str, Dict[str, Any]]) -> str:
     """
     extracts the .SAFE name from the asset returned by the STAC query
diff --git a/examples/satellite_archive_query.py b/examples/satellite_archive_query.py
deleted file mode 100644
index 253d118e..00000000
--- a/examples/satellite_archive_query.py
+++ /dev/null
@@ -1,96 +0,0 @@
-"""
-sample script showing how to perform a simple metadata query to identify the
-number of available mapper for a Sentinel-2 tile below a user-defined cloud
-cover threshold on data already downloaded and ingested into the metadata base
-(offline mode)
-
-For calling the Copernicus archive (and optionally downloading data) refer
-to the script './copernicus_archive_query.py'
-
-The called function also plots the cloud cover over time.
-
-Copyright (C) 2022 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 matplotlib
-import matplotlib.pyplot as plt
-import pandas as pd
-
-from datetime import date, datetime
-from pathlib import Path
-from eodal.operational.cli import cli_s2_scene_selection
-from eodal.utils.constants import ProcessingLevels
-
-# settings for plotting
-plt.style.use('seaborn-darkgrid')
-matplotlib.rc('ytick', labelsize=16)
-matplotlib.rc('xtick', labelsize=16)
-matplotlib.rc('font', size=16)
-
-# user inputs
-tile = 'T32TLT'
-processing_level = ProcessingLevels.L2A
-out_dir = Path(f'../data')
-date_start = '2019-01-01'
-date_end = '2020-12-31'
-cc_threshold = 100.
-
-# date range
-date_start = datetime.strptime(date_start, '%Y-%m-%d')
-date_end = datetime.strptime(date_end, '%Y-%m-%d')
-
-# execute scene selection (plots cloud cover over time)
-metadata = cli_s2_scene_selection(
-    tile=tile,
-    processing_level=processing_level,
-    cloud_cover_threshold=cc_threshold,
-    date_start=date_start,
-    date_end=date_end,
-    out_dir=Path(out_dir)
-)
-
-# plot cloud cover by month
-date_start = date(2017,1,1)
-date_end = date(2021,12,31)
-
-metadata = cli_s2_scene_selection(
-    tile=tile,
-    processing_level=processing_level,
-    cloud_cover_threshold=cc_threshold,
-    date_start=date_start,
-    date_end=date_end,
-    out_dir=Path(out_dir)
-)
-
-# group by month and plot the average cloud cover
-metadata_monthly = metadata['cloudy_pixel_percentage'].groupby(
-    by=pd.to_datetime(metadata.sensing_date).dt.month).agg('median')
-metadata_monthly = metadata_monthly.reset_index()
-months = {1: 'Jan', 2: 'Feb', 3: 'Mar', 4: 'Apr', 5: 'May', 6: 'Jun', 7: 'Jul', 8: 'Aug', 9: 'Sep',
-          10: 'Oct', 11: 'Nov', 12: 'Dec'}
-metadata_monthly.sensing_date = metadata_monthly.sensing_date.map(months)
-
-f, ax = plt.subplots(figsize=(8,6))
-ax.plot(metadata_monthly.sensing_date, metadata_monthly.cloudy_pixel_percentage,
-        marker='o', label='Median')
-ax.set_xlabel('Month', fontsize=16)
-ax.set_ylabel('Cloudy Pixel Percentage [%]', fontsize=16)
-ax.set_title(f'Sentinel-2 Tile {tile[1::]} ({date_start} - {date_end})\nNumber of Scenes: {metadata.shape[0]}',
-             size=18)
-ax.set_ylim(0,100)
-ax.set_xlim(-1,12)
-ax.legend(fontsize=16)
-f.savefig(out_dir.joinpath(f'monthly_cloudy_pixel_percentage_{date_start}-{date_end}.png'), bbox_inches='tight')
diff --git a/examples/sentinel1_mapping_example.py b/examples/sentinel1_mapping_example.py
deleted file mode 100644
index d0d7468c..00000000
--- a/examples/sentinel1_mapping_example.py
+++ /dev/null
@@ -1,73 +0,0 @@
-"""
-Example script to extract a time series of Sentinel-1 mapper for a
-custom area of interest (AOI).
-
-The script shows how to use the Sentinel1Mapper class that takes over
-data handling such as
-
-    * querying of spatio-temporal metadata catalogs to identify
-      available Sentinel-1 mapper
-    * merging data from different Sentinel-1 tiles if required
-    * re-projection of imagery from one UTM zone into another
-      if required
-
-This script works either using local data sources or by retrieving Sentinel-1
-imagery from Microsoft Planetary Computer (https://planetarycomputer.microsoft.com).
-To use the RTC collection (sentinel-1-rtc, radiometrically terrain corrected) you
-have to specify a valid PL API key.
-
-To use Planetary Computer make sure to set the `USE_STAC` variable to True
-and specify Microsoft as STAC provider (default). If you want to use the RTC product
-make sure to provide a valid Planetary Computer SDK subscription key.
-
-.. code-block:: shell
-
-    export USE_STAC = "True"
-    export PC_SDK_SUBSCRIPTION_KEY = <your-api-key-(optional-for-RTC)>
-
-Copyright (C) 2022 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/>.
-"""
-
-from datetime import date
-from eodal.operational.mapping import Sentinel1Mapper
-from eodal.operational.mapping import plot_feature
-from pathlib import Path
-
-#%% user-inputs
-
-# ------------------------- Time Range ---------------------------------
-date_start: date = date(2022,4,1)          # year, month, day (incl.)
-date_end: date = date(2022,4,14)           # year, month, day (incl.)
-
-# ---------------------- Area of Interest ------------------------------
-aoi: Path = Path('../data/sample_polygons/lake_lucerne.gpkg')
-
-# ---------------------- Sentinel-1 collection --------------------------
-collection = 'sentinel-1-rtc'
-
-
-#%% executable part
-# get a new mapper instance
-mapper = Sentinel1Mapper(
-    date_start=date_start,
-    date_end=date_end,
-    feature_collection=aoi,
-    collection=collection
-)
-
-# retrieve metadata of mapper found (no reading)
-mapper.get_scenes()
-mapper.observations
diff --git a/examples/sentinel2_l2a_safe.py b/examples/sentinel2_l2a_safe.py
deleted file mode 100644
index c3a52669..00000000
--- a/examples/sentinel2_l2a_safe.py
+++ /dev/null
@@ -1,83 +0,0 @@
-"""
-Sample script showing how to load data from a Sentinel-2 scene
-organized in .SAFE (Standard Archive Format for Europe) into a
-`Sentinel-2` object.
-
-The Sentinel-2 scene can be downloaded here:
-https://data.mendeley.com/datasets/ckcxh6jskz/1
-
-Copyright (C) 2022 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 os
-import geopandas as gpd
-from pathlib import Path
-from shapely.geometry import Polygon
-
-from eodal.config import get_settings
-from eodal.core.sensors import Sentinel2
-
-# make EOdal use local data sources
-Settings = get_settings()
-Settings.USE_STAC = False
-
-# file-path to the .SAFE dataset
-dot_safe_dir = Path('../data/S2A_MSIL2A_20190524T101031_N0212_R022_T32UPU_20190524T130304.SAFE')
-
-# construct a bounding box for reading a spatial subset of the scene (geographic coordinates)
-ymin, ymax = 47.949, 48.027
-xmin, xmax = 11.295, 11.385
-bbox = Polygon([(xmin, ymin), (xmax, ymin), (xmax, ymax), (xmin, ymax)])
-
-# eodal expects a vector file or a GeoDataFrame for spatial sub-setting
-bbox_gdf = gpd.GeoDataFrame(geometry=[bbox], crs=4326)
-
-# read data from .SAFE (all 10 and 20m bands + scene classification layer)
-s2_ds = Sentinel2().from_safe(
-    in_dir=dot_safe_dir,
-    vector_features=bbox_gdf
-)
-
-# eodal support band aliasing. Thus, you can access the bands by their name ...
-s2_ds.band_names
-# >>> ['B02', 'B03', 'B04', 'B05', 'B06', 'B07', 'B08', 'B8A', 'B11', 'B12', 'SCL']
-# ... or by their alias, i.e., their color name
-s2_ds.band_aliases
-# >>> ['blue', 'green', 'red', 'red_edge_1', 'red_edge_2', 'red_edge_3', 'nir_1', 'nir_2', 'swir_1', 'swir_2', 'scl']
-
-# plot false-color infrared preview
-fig_nir = s2_ds.plot_multiple_bands(band_selection=['nir_1','red','green'])
-fig_nir.savefig('../img/eodal_Sentinel-2_NIR.png', dpi=150, bbox_inches='tight')
-
-# plot scene classification layer
-fig_scl = s2_ds.plot_scl()
-fig_scl.savefig('../img/eodal_Sentinel-2_SCL.png', dpi=150, bbox_inches='tight')
-
-# calculate the NDVI using 10m bands (no spatial resampling required)
-s2_ds.calc_si('NDVI', inplace=True)
-fig_ndvi = s2_ds.plot_band('NDVI', colormap='summer', vmin=-1, vmax=1)
-fig_ndvi.savefig('../img/eodal_Sentinel-2_NDVI.png', dpi=150, bbox_inches='tight')
-
-# mask the water (SCL class 6); requires resampling to 10m spatial resolution
-s2_ds.resample(target_resolution=10, inplace=True)
-s2_ds.mask(
-    mask='SCL',
-    mask_values=[6],
-    bands_to_mask=['NDVI'],
-    inplace=True
-)
-fig_ndvi = s2_ds.plot_band('NDVI', colormap='summer', vmin=-1, vmax=1)
-fig_ndvi.savefig('../img/eodal_Sentinel-2_NDVI_masked.png', dpi=150, bbox_inches='tight')
diff --git a/examples/sentinel2_mapping_example.py b/examples/sentinel2_mapping_example.py
deleted file mode 100644
index c0c5fe31..00000000
--- a/examples/sentinel2_mapping_example.py
+++ /dev/null
@@ -1,104 +0,0 @@
-"""
-Example script to extract a time series of Sentinel-2 mapper for a
-custom area of interest (AOI).
-
-The script shows how to use the Sentinel2Mapper class that takes over
-data handling such as
-
-	* querying of spatio-temporal metadata catalogs to identify
-	  available Sentinel-2 mapper
-	* merging data from different Sentinel-2 tiles if required
-	* re-projection of imagery from one UTM zone into another
-	  if required
-	* removal of black-filled mapper
-
-This script works either using local data sources or by retrieving Sentinel-2
-imagery from Microsoft Planetary Computer (https://planetarycomputer.microsoft.com).
-No authentication required.
-
-To use Planetary Computer make sure to set the `USE_STAC` variable to True
-and specify Microsoft as STAC provider (default).
-
-.. code-block:: shell
-
-    export USE_STAC = "True"
-
-Copyright (C) 2022 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 cv2
-
-from datetime import date
-from eodal.operational.mapping import MapperConfigs, Sentinel2Mapper
-from eodal.operational.mapping import plot_feature
-from eodal.utils.sentinel2 import ProcessingLevels
-from pathlib import Path
-
-#%% user-inputs
-# ---------------------- Spatial Resolution ----------------------------
-spatial_resolution: int = 10 			# meters
-resampling_method: int = cv2.INTER_NEAREST_EXACT
-
-# ----------------------- Processing Level -----------------------------
-processing_level: ProcessingLevels = ProcessingLevels.L2A 	# BOA
-
-# ------------------------ Cloud Coverage ------------------------------
-scene_cloud_cover_threshold: int = 50 	# percent
-
-# ------------------------- Time Range ---------------------------------
-date_start: date = date(2022,4,1)  		# year, month, day (incl.)
-date_end: date = date(2022,4,14)   		# year, month, day (incl.)
-
-# ---------------------- Area of Interest ------------------------------
-aoi: Path = Path('../data/sample_polygons/lake_lucerne.gpkg')
-
-#%% executable part
-# Sentinel-2 mapper configuration
-mapper_configs = MapperConfigs(
-	spatial_resolution=spatial_resolution,
-	resampling_method=resampling_method,
-)
-
-# get a new mapper instance
-mapper = Sentinel2Mapper(
-	date_start=date_start,
-	date_end=date_end,
-	processing_level=processing_level,
-	cloud_cover_threshold=scene_cloud_cover_threshold,
-	mapper_configs=mapper_configs,
-	feature_collection=aoi
-)
-
-# retrieve metadata of mapper found (no reading)
-mapper.get_scenes()
-mapper.observations # displays mapper found
-
-# read data into eodal's RasterCollection objects
-s2_data = mapper.get_complete_timeseries()
-features = mapper.get_feature_ids()
-
-# loop features (in this case it's just a single one) and plot them
-for feature in features:
-	fig_feature = plot_feature(
-		feature_scenes=s2_data[feature],
-		band_selection=['nir_1', 'red', 'green'],
-		figsize=(16,14),
-		max_scenes_in_row=2,
-		sharex=True,
-		sharey=True
-	)
-	fig_feature.savefig('../img/s2_lake_lucerne.png', bbox_inches='tight')
-
-# further program logic ...
diff --git a/examples/setup_db.py b/examples/setup_db.py
deleted file mode 100644
index 81b9b148..00000000
--- a/examples/setup_db.py
+++ /dev/null
@@ -1,9 +0,0 @@
-"""
-Creates tables in a PostgreSQL/PostGIS DBMS instance if not available yet
-for offline mode
-"""
-
-from eodal.metadata.database.db_model import create_tables
-
-if __name__ == '__main__':
-    create_tables()
diff --git a/img/EOdal_Concept.png b/img/EOdal_Concept.png
deleted file mode 100644
index 2e4363a5..00000000
Binary files a/img/EOdal_Concept.png and /dev/null differ
diff --git a/img/eodal_data-providers.png b/img/eodal_data-providers.png
deleted file mode 100644
index ef0d3a2d..00000000
Binary files a/img/eodal_data-providers.png and /dev/null differ
diff --git a/requirements.txt b/requirements.txt
index 152d504a..08b591c7 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -23,3 +23,4 @@ requests
 planetary_computer
 alive_progress
 rasterstats
+pyyaml
diff --git a/scripts/README.md b/scripts/README.md
new file mode 100644
index 00000000..cd55a748
--- /dev/null
+++ b/scripts/README.md
@@ -0,0 +1,7 @@
+# README
+
+This directory contains a more or less consolidated collection of scripts useful for our day-to-day work.
+The scripts fuflil two purposes:
+
+	a) they show how to E:earth_africa:dal
+	b) they allow you to automate workflows such as keeping a local Sentinel-1 or Sentinel-2 archive up-to-date or how to download PlanetScope data.
diff --git a/examples/digital_terrain_model.py b/scripts/cloud_optimized_geoTiff.py
similarity index 100%
rename from examples/digital_terrain_model.py
rename to scripts/cloud_optimized_geoTiff.py
diff --git a/examples/planet_download.py b/scripts/planet_download.py
similarity index 100%
rename from examples/planet_download.py
rename to scripts/planet_download.py
diff --git a/scripts/sentinel2_mapper.py b/scripts/sentinel2_mapper.py
new file mode 100644
index 00000000..d1bf6c97
--- /dev/null
+++ b/scripts/sentinel2_mapper.py
@@ -0,0 +1,137 @@
+"""
+Script to extract a collection of Sentinel-2 scenes for a
+custom area of interest (AOI).
+
+The script shows how to use the EOdal Mapper class that takes over
+data handling such as
+
+	* querying of spatio-temporal metadata catalogs to identify
+	  available Sentinel-2 scenes
+	* merging data from different Sentinel-2 tiles if required
+	* re-projection of imagery from one UTM zone into another
+	  if required
+	* removal of black-filled (i.e., no-data) scenes
+
+This script works by retrieving Sentinel-2 scenes from
+Microsoft Planetary Computer (https://planetarycomputer.microsoft.com).
+This requires no authentication required. Alternatively, the same code
+can be used to read data from a local EOdal Sentinel archive or a different
+STAC provider.
+
+Copyright (C) 2022/23 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
+
+from datetime import datetime
+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 archove
+Settings.USE_STAC = True
+
+def preprocess_sentinel2_scenes(
+		ds: Sentinel2,
+		target_resolution: int,
+	) -> Sentinel2:
+	"""
+	Resample Sentinel-2 scenes and mask clouds, shadows, and snow
+	based on the Scene Classification Layer (SCL).
+
+	NOTE:
+		Depending on your needs, the pre-processing function can be
+		fully customized using the full power of EOdal and its
+		interfacing libraries!
+
+	:param target_resolution:
+		spatial target resolution to resample all bands to.
+	:returns:
+		resampled, cloud-masked Sentinel-2 scene.
+	"""
+	# resample scene
+	ds.resample(inplace=True, target_resolution=target_resolution)
+	# mask clouds, shadows, and snow
+	ds.mask_clouds_and_shadows(inplace=True)
+	return ds
+
+if __name__ == '__main__':
+	#%% user-inputs
+	# -------------------------- Collection -------------------------------
+	collection: str = 'sentinel2-msi'
+	
+	# ------------------------- Time Range ---------------------------------
+	time_start: datetime = datetime(2022,3,1)  		# year, month, day (incl.)
+	time_end: datetime = datetime(2022,6,30)   		# year, month, day (incl.)
+	
+	# ---------------------- Spatial Feature  ------------------------------
+	geom: Path = Path('../data/sample_polygons/lake_lucerne.gpkg')
+	
+	# ------------------------- Metadata Filters ---------------------------
+	metadata_filters: List[Filter] = [
+		Filter('cloudy_pixel_percentage','<', 80),
+		Filter('processing_level', '==', 'Level-2A')
+	]
+	
+	#%% 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,
+	    time_end=time_end,
+	    feature=feature,
+	    metadata_filters=metadata_filters
+	)
+	# to enhance reproducibility and provide proper documentation, the MapperConfigs
+	# can be saved as yaml (and also then be loaded again from yaml)
+	mapper_configs.to_yaml('../data/sample_mapper_call.yaml')
+	
+	# now, a new Mapper instance is created
+	mapper = Mapper(mapper_configs)
+	mapper.query_scenes()
+	# the metadata is loaded into a GeoPandas GeoDataFrame
+	mapper.metadata
+	
+	#%% load the scenes available from STAC
+	scene_kwargs = {
+	    'scene_constructor': Sentinel2.from_safe,
+	    'scene_constructor_kwargs': {'band_selection': ['B02', 'B03', 'B04']},
+	    'scene_modifier': preprocess_sentinel2_scenes,
+	    'scene_modifier_kwargs': {'target_resolution': 10}
+	}
+	mapper.load_scenes(scene_kwargs=scene_kwargs)
+	# the data loaded into `mapper.data` as a EOdal SceneCollection
+	mapper.data
+
+	# plot scenes in collection
+	f_scenes = mapper.data.plot(band_selection=['red', 'green', 'blue'])
+
+	# EOdal SceneCollections can be made persistent by storing them as serialized pickled
+	# objects on disk (and can be loaded from there)
+	fpath = Path('../data/sample_mapper_data.pkl')
+	with open(fpath, 'wb+') as dst:
+		dst.write(mapper.data.to_pickle())
+
+	# read data from pickled file object into SceneCollectio
+	scoll = SceneCollection.from_pickle(stream=fpath)
diff --git a/eodal/operational/archive/planet_scope.py b/scripts/setup_db.py
similarity index 61%
rename from eodal/operational/archive/planet_scope.py
rename to scripts/setup_db.py
index 480cb5c4..0b157f82 100644
--- a/eodal/operational/archive/planet_scope.py
+++ b/scripts/setup_db.py
@@ -1,7 +1,10 @@
-'''
-Function to keep a local Planet archive up-to-date
+"""
+This scripts shows how to create a local metadata DBMS for storing satellite
+data locally and make it searchable.
+Creates tables in a PostgreSQL/PostGIS DBMS instance if not available yet.
+Existing tables remain untouched!
 
-Copyright (C) 2022 Samuel Wildhaber
+Copyright (C) 2022 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
@@ -15,4 +18,9 @@
 
 You should have received a copy of the GNU General Public License
 along with this program.  If not, see <http://www.gnu.org/licenses/>.
-'''
+"""
+
+from eodal.metadata.database.db_model import create_tables
+
+if __name__ == '__main__':
+    create_tables()
diff --git a/scripts/update_sentinel_archive.py b/scripts/update_sentinel_archive.py
new file mode 100644
index 00000000..0bc7867f
--- /dev/null
+++ b/scripts/update_sentinel_archive.py
@@ -0,0 +1,403 @@
+"""
+Function to keep a local Sentinel (1/2) archive up-to-date by comparing data available
+in your local archive for a given geographic region, time period, product type and
+sensor mode.
+
+The function not only checks if CREODIAS has new datasets available, it also automatically
+downloads them into a user-defined location.
+
+IMPORTANT:
+    In order to receive results the region (i.e., geographic extent) for which
+    to download data must be defined in the metadata DB.
+
+IMPORTANT:
+    CREODIAS does not allow more than 2000 records (each Sentinel scene is a record)
+    to be queried at once. If you might exceed this threshold (e.g., your region is large and/or
+    your time period is long) split your query into smaller chunks by using, e.g., shorter time
+    periods for querying.
+
+Copyright (C) 2022/23 Gregor Perich & 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/>.
+"""
+
+from __future__ import annotations
+
+import pandas as pd
+import numpy as np
+import shutil
+
+from datetime import date, datetime
+from pathlib import Path
+from typing import Dict, Optional
+
+from eodal.config import get_settings
+from eodal.downloader.sentinel1.creodias import query_creodias as s1_creodias_query
+from eodal.downloader.sentinel2.creodias import query_creodias as s2_creodias_query
+from eodal.downloader.utils.creodias import download_datasets
+from eodal.downloader.utils import unzip_datasets
+from eodal.metadata.database.querying import get_region
+from eodal.metadata.sentinel1.database.ingestion import meta_df_to_database as s1_meta_to_db
+from eodal.metadata.sentinel1.database.querying import find_raw_data_by_bbox as s1_db_query
+from eodal.metadata.sentinel1.parsing import parse_s1_metadata
+from eodal.metadata.sentinel2.database.ingestion import meta_df_to_database as s2_meta_to_db
+from eodal.metadata.sentinel2.database.querying import find_raw_data_by_bbox as s2_db_query
+from eodal.metadata.sentinel2.parsing import parse_s2_scene_metadata
+from eodal.utils.constants import ProcessingLevels
+
+settings = get_settings()
+logger = settings.logger
+max_records = settings.CREODIAS_MAX_RECORDS
+
+def pull_from_creodias(
+    date_start: date,
+    date_end: date,
+    path_out: Path,
+    region: str,
+    sensor: str,
+    unzip: Optional[bool] = True,
+    overwrite_existing_zips: Optional[bool] = False,
+    **kwargs
+) -> pd.DataFrame:
+    """
+    Checks if CREODIAS has Sentinet datasets not yet available locally
+    and downloads these datasets from CREODIAS.
+
+    IMPORTANT:
+        Unless specified otherwise, default products of the Sentinel
+        platforms are used (e.g., GRD-IW for Sentinel-1, L2A for Sentinel-2)
+
+    NOTE:
+        CREODIAS limits the maximum amount of datasets to download within a
+        single query to somewhat around 2000. We therefore recommend to split
+        the query into smaller ones if the query is likely to hit this limit.
+
+    :param date start:
+        Start date of the database & creodias query
+    :param date end:
+        End date of the database & creodias query
+    :param path_out:
+        Out directory where the additional data from CREODIAS should be
+        downloaded to
+    :param region:
+        Region identifier of the Sentinel-2 archive. By region we mean a
+        geographic extent (bounding box) in which the data is organized. The bounding
+        box extent is taken from the metadata DB based on the region identifier.
+    :param sensor:
+        name of the sensor to download data. Currently, `sentinel1` and `sentinel2`
+        are supported.
+    :param unzip:
+        if True (default) datasets are unzipped and zip archives are deleted
+    :param overwrite_existing_zips:
+        if False (default) overwrites eventually existing zip files. If the download
+        process was interrupted (e.g., due to a connection timeout) setting the flag
+        to True can save time because datasets already downloaded are ignored. NOTE:
+        The function does **not** check if a dataset was downloaded completely!
+    :param kwargs:
+        optional key-word arguments to pass to 
+        `~eodal.downloader.sentinel1.creodias.query_creodias` such sensor_mode and
+        product_type
+    :return:
+        dataframe with references to downloaded datasets
+    """
+
+    # query database to get the bounding box of the selected region
+    try:
+        region_gdf = get_region(region)
+    except Exception as e:
+        logger.error(f"Failed to query region: {e}")
+        return
+
+    # parse the region's geometry as extended well-known-text
+    bounding_box = region_gdf.geometry.iloc[0]
+    bounding_box_ewkt = f"SRID=4326;{bounding_box.wkt}"
+
+    # local database query to check what is already available locally
+    if sensor == 'sentinel1':
+        try:
+            product_type = kwargs.get('product_type', 'GRD')
+            sensor_mode = kwargs.get('sensor_mode', 'IW')
+            meta_db_df = s1_db_query(
+                date_start=date_start,
+                date_end=date_end,
+                product_type=product_type,
+                sensor_mode=sensor_mode,
+                bounding_box=bounding_box_ewkt,
+            )
+
+            # check for available datasets on CREODIAS
+            datasets = s1_creodias_query(
+                start_date=date_start,
+                end_date=date_end,
+                max_records=max_records,
+                bounding_box=bounding_box,
+                product_type=product_type,
+                sensor_mode=sensor_mode
+            )
+        except Exception as e:
+            logger.error(f'Failed to update Sentinel1 archive: {e}')
+            return pd.DataFrame([])
+
+    elif sensor == 'sentinel2':
+        try:
+            processing_level = kwargs.get('processing_level', ProcessingLevels.L2A)
+            cloud_cover_threshold = kwargs.get('cloud_cover_threshold', 100)
+            meta_db_df = s2_db_query(
+                date_start=date_start,
+                date_end=date_end,
+                processing_level=processing_level,
+                bounding_box=bounding_box_ewkt,
+            )
+
+            # check for available datasets
+            datasets = s2_creodias_query(
+                start_date=date_start,
+                end_date=date_end,
+                max_records=max_records,
+                processing_level=processing_level,
+                bounding_box=bounding_box,
+                cloud_cover_threshold=cloud_cover_threshold,
+            )
+        except Exception as e:
+            logger.error(f'Failed to update Sentinel2 archive: {e}')
+            return pd.DataFrame([])
+
+    else:
+        raise ValueError(f'Unknown sensor: {sensor}')
+
+    # get .SAFE datasets from CREODIAS
+    datasets["product_uri"] = datasets.properties.apply(
+        lambda x: Path(x["productIdentifier"]).name
+    )
+
+    # compare with records from local metadata DB and keep those records
+    # not available locally
+    missing_datasets = np.setdiff1d(
+        datasets["product_uri"].values,
+        meta_db_df["product_uri"].values
+    )
+    datasets_filtered = datasets[datasets.product_uri.isin(missing_datasets)].copy()
+
+    # download those mapper not available in the local database from CREODIAS
+    download_datasets(
+        datasets=datasets_filtered,
+        download_dir=path_out,
+        overwrite_existing_zips=overwrite_existing_zips,
+    )
+
+    # unzip datasets
+    if unzip:
+        unzip_datasets(download_dir=path_out)
+
+    return datasets_filtered
+
+def sentinel_creodias_update(
+    sentinel_raw_data_archive: Path,
+    region: str,
+    path_options: Optional[Dict[str, str]] = None,
+    overwrite_existing_zips: Optional[bool] = False,
+    **kwargs
+) -> None:
+    """
+    Loops over an existing Sentinel rawdata (i.e., *.SAFE datasets) archive
+    and checks the datasets available locally with those datasets available at
+    CREODIAS for a defined Area of Interest (AOI) and time period (we store the
+    data year by year, therefore, we always check an entire year).
+
+    The missing data (if any) is downloaded into a temporary download directory
+    `temp_dl` and unzip. The data is then copied into the actual Sentinel
+    archive and the metadata is extracted and ingested into the metadata base.
+
+    IMPORTANT:
+        Requires a CREODIAS user account (user name and password).
+
+    IMPORTANT:
+        For database ingestion it is important to map the paths correctly
+        and store them in the database in a way that allows accessing the data from
+        all your system components (``file_system_options``).
+
+    * In the easiest case, this the absolute path to the datasets (or URI)
+    * If your data is stored on a NAS, you might specify the address of the NAS
+      in the variable `storage_device_ip` and provide a `mount_point`, i.e.,
+      the path (or drive on Windows) where the NAS is mounted into your local
+      file system. Also aliasing of the is supported (`storage_device_ip_alias`),
+      however, if possible try to avoid it.
+
+
+    :param sentinel_raw_data_archive:
+        Sentinel (1 or 2) raw data archive (containing *.SAFE datasets) to monitor.
+        Existing datasets must have been already ingested into the metadata
+        base!
+    :param region:
+        eodal's archive philosophy organizes datasets by geographic regions.
+        Each region is identified by a unique region identifier (e.g., we use
+        `CH` for Switzerland) and has a geographic extent described by a polygon
+        geometry (bounding box) in geographic coordinates (WGS84). The geometry
+        also defines the geographic dimension of the CREODIAS query. It is
+        stored as a entry in the metadata base.
+    :param processing_level:
+        Sentinel-2 processing level (L1C or L2A) to check.
+    :param cloud_cover_threshold:
+        optional cloud cover threshold to filter out to cloudy mapper as integer
+        between 0 and 100%.
+    :param path_options:
+        optional dictionary specifying storage_device_ip, storage_device_ip_alias
+        (if applicable) and mount point in case the data is stored on a NAS
+        and should be accessible from different operating systems or file systems
+        with different mount points of the NAS share. If not provided, the absolute
+        path of the dataset is used in the database.
+    :param overwrite_existing_zips:
+        if False (default) overwrites eventually existing zip files. If the download
+        process was interrupted (e.g., due to a connection timeout) setting the flag
+        to True can save time because datasets already downloaded are ignored. NOTE:
+        The function does **not** check if a dataset was downloaded completely!
+
+    Example
+    -------
+
+    .. code-block:: python
+
+    from pathlib import Path
+    from eodal.utils.constants import ProcessingLevels
+
+    # this example shows the workflow for Sentinel-2 (Sentinel-1 would also be possible
+    # but requires slightly different inputs)
+    sensor = 'sentinel2'
+    
+    # define processing level (usually L2A but also L1C works)
+    processing_level = ProcessingLevels.L2A
+
+    # specifiy region, we use Switzerland (the extent of Switzerland is defined in the database)
+    region = 'CH'
+
+    # the archive should be always mounted in the same way for each user
+    user_name = '<your_username>'
+    s2_raw_data_archive = Path(f'/home/{user_name}/public/Evaluation/Satellite_data/Sentinel-2/Rawdata')
+
+    # file-system specific handling: this allows to store the paths in the database
+    # in such way that the dataset paths can be found from Linux and Windows machines
+    file_system_options = {
+        'storage_device_ip': '<your_nas_ip>',
+        'storage_device_ip_alias': '<alternative_nas_ip_if_any',
+        'mount_point': f'/home/{user_name}/public/'
+    }
+
+    sentinel_creodias_update(
+        s2_raw_data_archive=s2_raw_data_archive,
+        region=region,
+        path_options=file_system_options,
+        sensor='sentinel2',
+        processing_level=processing_level
+    )
+
+    """
+    # since the data is stored by year (each year is a single sub-directory) we
+    # can simple loop over the sub-directories and do the check
+    for path in Path(sentinel_raw_data_archive).iterdir():
+        if path.is_dir():
+            # get year automatically
+            try:
+                year = int(path.name)
+            except ValueError:
+                logger.info(f'{path.name} is not provided in YYYY format - skipping')
+                continue
+
+            # create temporary download directory
+            path_out = path.joinpath(f"temp_dl_{year}")
+            path_out.mkdir(exist_ok=True)
+
+            # download data from CREODIAS
+            downloaded_ds = pull_from_creodias(
+                date_start=date(year,1,1),
+                date_end=date(year,12,31),
+                path_out=path_out,
+                region=region,
+                overwrite_existing_zips=overwrite_existing_zips,
+                **kwargs
+            )
+
+            if downloaded_ds.empty:
+                logger.info(f"No new datasets found for year {year} on CREODIAS")
+                continue
+
+            # move the datasets into the actual SAT archive (on level up)
+            error_happened = False
+            errored_datasets = []
+            error_msgs = []
+            parent_dir = path_out.parent
+            for _, record in downloaded_ds.iterrows():
+                try:
+                    shutil.move(record.dataset_name, "..")
+                    logger.info(f"Moved {record.dataset_name} to {parent_dir}")
+                    # once the dataset is moved successfully parse its metadata and
+                    # ingest it into the database
+                    in_dir = path.joinpath(record.dataset_name)
+                    sensor = kwargs.get('sensor')
+                    if sensor == 'sentinel1':
+                        scene_metadata = parse_s1_metadata(in_dir=in_dir)
+                    elif sensor == 'sentinel2':
+                        scene_metadata, _ = parse_s2_scene_metadata(in_dir=in_dir)
+
+                    # some path handling if required
+                    if path_options != {}:
+                        scene_metadata["storage_device_ip"] = path_options.get(
+                            "storage_device_ip", ""
+                        )
+                        scene_metadata["storage_device_ip_alias"] = path_options.get(
+                            "storage_device_ip_alias", ""
+                        )
+                        mount_point = path_options.get("mount_point", "")
+                        mount_point_replacement = path_options.get(
+                            "mount_point_replacement", ""
+                        )
+                        scene_metadata["storage_share"] = scene_metadata[
+                            "storage_share"
+                        ].replace(mount_point, mount_point_replacement)
+
+                    # database insert
+                    if sensor == 'sentinel1':
+                        s1_meta_to_db(meta_df=scene_metadata)
+                    elif sensor == 'sentinel2':
+                        s2_meta_to_db(meta_df=scene_metadata)
+
+                    logger.info(
+                        f"Ingested scene metadata for {record.dataset_name} into DB"
+                    )
+
+                except Exception as e:
+                    error_happened = True
+                    logger.error(f"{record.dataset_name} produced an error: {e}")
+                    errored_datasets.append(record.dataset_name)
+                    error_msgs.append(e)
+
+            # delete the temp_dl directory
+            shutil.rmtree(path_out)
+
+            # log the errored datasets if any error happened
+            if error_happened:
+                # get the timestamp (date) of the current run
+                processing_date = datetime.now().date().strftime("%Y-%m-%d")
+                # write to log directory
+                log_dir = path_out.parent.joinpath("logs")
+                if not log_dir.exists():
+                    log_dir.mkdir()
+                errored_logfile = log_dir.joinpath("datasets_errored.csv")
+                if not errored_logfile.exists():
+                    with open(errored_logfile, "w") as src:
+                        line = "download_date,dataset_name,error_message\n"
+                        src.writelines(line)
+                with open(errored_logfile, "a") as src:
+                    for error in list(zip(errored_datasets, error_msgs)):
+                        line = processing_date + "," + error[0] + "," + str(error[1])
+                        src.writelines(line + "\n")


Test patch:
diff --git a/tests/mapper/test_feature.py b/tests/mapper/test_feature.py
new file mode 100644
index 00000000..db45febc
--- /dev/null
+++ b/tests/mapper/test_feature.py
@@ -0,0 +1,65 @@
+'''
+Tests for the Feature class (spatial filters for mapper)
+
+.. versionadded:: 0.1.1
+
+Copyright (C) 2022 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 pandas as pd
+import pytest
+
+from shapely.geometry import Point
+from eodal.mapper.feature import Feature
+
+def test_feature():
+
+    # working constructor calls
+    geom = Point([49,11])
+    epsg = 4326
+    name = 'Test Point'
+    feature = Feature(name, geom, epsg)
+
+    assert feature.geometry == geom, 'geometry differs'
+    assert feature.epsg == epsg, 'EPSG code differs'
+    assert feature.name == name, 'name differs'
+    assert feature.attributes == {}, 'attributes must be empty'
+
+    attributes = {'key': 'value'}
+    feature = Feature(name, geom, epsg, attributes)
+    assert feature.attributes == attributes, 'attributes differ'
+
+    attributes = pd.Series({'key1': 'value1', 'key2': 'value2'})
+    feature = Feature(name, geom, epsg, attributes)
+    assert feature.attributes == attributes.to_dict(), 'attributes differ'
+
+    gds = feature.to_geoseries()
+    assert gds.name == feature.name, 'name differs'
+    assert gds.crs.to_epsg() == feature.epsg, 'EPSG differs'
+    assert gds.attrs == feature.attributes, 'attributes differ'
+
+    # from_geoseries class method
+    gds.attrs = {}
+    feature = Feature.from_geoseries(gds)
+    assert gds.name == feature.name, 'name differs'
+    assert gds.crs.to_epsg() == feature.epsg, 'EPSG differs'
+    assert gds.attrs == feature.attributes, 'attributes differ'
+
+    # project into another spatial reference system
+    feature_utm = feature.to_epsg(epsg=32632)
+    assert feature_utm.epsg == 32632, 'projection had no effect'
+    assert feature_utm.name == feature.name, 'name got lost'
+    assert feature_utm.attributes == feature.attributes, 'attributes got lost'
diff --git a/tests/mapper/test_filter.py b/tests/mapper/test_filter.py
new file mode 100644
index 00000000..4dd0e3ff
--- /dev/null
+++ b/tests/mapper/test_filter.py
@@ -0,0 +1,42 @@
+'''
+Tests for the filter class.
+
+.. versionadded:: 0.1.1
+
+Copyright (C) 2022 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 pytest
+
+from eodal.mapper.filter import Filter
+
+def test_filter():
+    # filter by cloud cover
+    cc_filter = Filter(entity='cloudy_pixel_percentage', operator='<', value=30)
+    assert cc_filter.expression == 'cloudy_pixel_percentage < 30'
+    assert cc_filter.entity == 'cloudy_pixel_percentage'
+    assert cc_filter.operator == '<'
+    assert cc_filter.value == 30
+
+    # wrong data types
+    with pytest.raises(ValueError):
+        cc_filter = Filter(entity='cloudy_pixel_percentage', operator='lt', value=30)
+    with pytest.raises(TypeError):
+        cc_filter = Filter(entity=4, operator='<', value=30)
+
+    # passing None
+    with pytest.raises(ValueError):
+        cc_filter = Filter(entity='cloudy_pixel_percentage', operator='<', value=None)
diff --git a/tests/mapper/test_mapper.py b/tests/mapper/test_mapper.py
new file mode 100644
index 00000000..855f91dc
--- /dev/null
+++ b/tests/mapper/test_mapper.py
@@ -0,0 +1,264 @@
+'''
+Created on Dec 14, 2022
+
+@author: graflu
+'''
+
+import pytest
+
+import geopandas as gpd
+
+from datetime import datetime
+from shapely.geometry import Point, Polygon
+
+from eodal.config import get_settings
+from eodal.core.scene import SceneCollection
+from eodal.core.sensors import Sentinel1, Sentinel2
+from eodal.mapper.feature import Feature
+from eodal.mapper.filter import Filter
+from eodal.mapper.mapper import MapperConfigs, Mapper
+
+Settings = get_settings()
+
+def test_mapper_configs(tmppath):
+
+    # construct a Feature
+    collection = 'sentinel2-msi'
+    geom = Point([49,11])
+    epsg = 4326
+    name = 'Test Point'
+    attributes = {'this': 'is a test', 'a': 123}
+    feature = Feature(name, geom, epsg, attributes)
+
+    # define further search criteria
+    time_start = datetime(2022,12,1)
+    time_end = datetime.now()
+    metadata_filters = [
+        Filter('cloudy_pixel_percentage', '<', 30), Filter('snow_ice_percentage', '<', 10)
+    ]
+
+    mapper_configs = MapperConfigs(collection, feature, time_start, time_end, metadata_filters)
+    fpath = tmppath.joinpath('test.yml')
+    mapper_configs.to_yaml(fpath)
+
+    mapper_configs_rl = MapperConfigs.from_yaml(fpath)
+    assert mapper_configs.feature.name == mapper_configs_rl.feature.name, 'wrong feature name'
+    assert mapper_configs.feature.epsg == mapper_configs_rl.feature.epsg, 'wrong feature EPSG'
+    assert mapper_configs.feature.geometry == mapper_configs_rl.feature.geometry, 'wrong feature geometry'
+    assert set(mapper_configs.feature.attributes) == set(mapper_configs_rl.feature.attributes), \
+        'wrong feature attributes'
+    assert mapper_configs.time_start == mapper_configs_rl.time_start, 'wrong start time'
+    assert mapper_configs.time_end == mapper_configs_rl.time_end, 'wrong end time'
+    assert mapper_configs.collection == mapper_configs_rl.collection, 'wrong collection'
+    for idx, _filter in enumerate(mapper_configs.metadata_filters):
+        assert _filter.expression == mapper_configs_rl.metadata_filters[idx].expression, \
+            'wrong filtering expression'
+        assert _filter.entity == mapper_configs_rl.metadata_filters[idx].entity, \
+            'wrong filtering entity'
+
+@pytest.mark.parametrize(
+    'collection,time_start,time_end,geom,metadata_filters',
+    [(
+        'sentinel2-msi',
+        datetime(2022,7,1),
+        datetime(2022,7,15),
+        Polygon(
+            [[493504.953633058525156, 5258840.576098721474409],
+            [493511.206339373020455, 5258839.601945200935006],
+            [493510.605988947849255, 5258835.524093257263303],
+            [493504.296645800874103, 5258836.554883609525859],
+            [493504.953633058525156, 5258840.576098721474409]]
+        ),
+        [Filter('cloudy_pixel_percentage','<', 100), Filter('processing_level', '==', 'Level-2A')],
+    ),
+    (
+        'sentinel1-grd',
+        datetime(2020,7,1),
+        datetime(2020,7,15),
+        Polygon(
+            [[493504.953633058525156, 5258840.576098721474409],
+            [493511.206339373020455, 5258839.601945200935006],
+            [493510.605988947849255, 5258835.524093257263303],
+            [493504.296645800874103, 5258836.554883609525859],
+            [493504.953633058525156, 5258840.576098721474409]]
+        ),
+        [Filter('product_type','==', 'GRD')],
+    )]
+)
+def test_mapper_get_scenes_db(collection, time_start, time_end, geom, metadata_filters):
+    """
+    testing the get_scenes() method of the mapper class
+
+    IMPORTANT:
+        This tests only works within the ETH EOdal environment as it expects
+        certain entries (scenes) in the database (see TODO-statement)
+
+    TODO: Create a test database instance
+    """
+    Settings.USE_STAC = False
+    feature = Feature(
+        name='Test Area',
+        geometry=geom,
+        epsg=32632,
+        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, 'expected some items to be returned'
+
+
+@pytest.mark.parametrize(
+    'collection,time_start,time_end,geom,metadata_filters',
+    [(
+        'sentinel2-msi',
+        datetime(2022,7,1),
+        datetime(2022,7,15),
+        Polygon(
+            [[493504.953633058525156, 5258840.576098721474409],
+            [493511.206339373020455, 5258839.601945200935006],
+            [493510.605988947849255, 5258835.524093257263303],
+            [493504.296645800874103, 5258836.554883609525859],
+            [493504.953633058525156, 5258840.576098721474409]]
+        ),
+        [Filter('cloudy_pixel_percentage','<', 100), Filter('processing_level', '==', 'Level-2A')],
+    ),
+    (
+        'sentinel1-grd',
+        datetime(2020,7,1),
+        datetime(2020,7,15),
+        Polygon(
+            [[493504.953633058525156, 5258840.576098721474409],
+            [493511.206339373020455, 5258839.601945200935006],
+            [493510.605988947849255, 5258835.524093257263303],
+            [493504.296645800874103, 5258836.554883609525859],
+            [493504.953633058525156, 5258840.576098721474409]]
+        ),
+        [Filter('product_type','==', 'GRD')],
+    )]
+)
+def test_mapper_get_scenes_stac(collection, time_start, time_end, geom, metadata_filters):
+    """"""
+    Settings.USE_STAC = True
+    feature = Feature(
+        name='Test Area',
+        geometry=geom,
+        epsg=32632,
+        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, 'expected some items to be returned'
+
+@pytest.fixture
+def get_mapper():
+    def _get_mapper():
+
+        Settings.USE_STAC = False
+
+        collection = 'sentinel2-msi'
+        time_start = datetime(2022,7,1)
+        time_end = datetime(2022,7,15)
+        geometry = Polygon(
+            [[493504.953633058525156, 5258840.576098721474409],
+            [493511.206339373020455, 5258839.601945200935006],
+            [493510.605988947849255, 5258835.524093257263303],
+            [493504.296645800874103, 5258836.554883609525859],
+            [493504.953633058525156, 5258840.576098721474409]]
+        )
+        filters = [Filter('cloudy_pixel_percentage','<', 100), Filter('processing_level', '==', 'Level-2A')]
+        feature = Feature(
+            name='Test Area',
+            geometry=geometry,
+            epsg=32632,
+            attributes={'id': 1}
+        )
+        mapper_configs = MapperConfigs(
+            collection=collection,
+            time_start=time_start,
+            time_end=time_end,
+            feature=feature,
+            metadata_filters=filters
+        )
+        mapper = Mapper(mapper_configs)
+        mapper.query_scenes()
+
+        return mapper
+    return _get_mapper
+
+def test_mapper_load_scenes(get_mapper):
+    """
+    test loading of Sentinel-2 scenes into a SceneCollection instance
+    """
+    mapper = get_mapper()
+    scene_constructor = Sentinel2.from_safe
+    scene_constructor_kwargs = {'band_selection': ['B02', 'B03', 'B04']}
+    load_scenes_kwargs = {
+        'scene_constructor': scene_constructor,
+        'scene_constructor_kwargs': scene_constructor_kwargs
+    }
+    mapper.load_scenes(scene_kwargs=load_scenes_kwargs)
+
+    assert isinstance(mapper.data, SceneCollection), 'expected a SceneCollection'
+    assert isinstance(mapper.data[mapper.data.identifiers[0]], Sentinel2), 'expected a Sentinel 2 scene'
+
+@pytest.mark.parametrize(
+    'collection,time_start,time_end,geom,metadata_filters',
+    [(
+        'sentinel2-msi',
+        datetime(2022,7,1),
+        datetime(2022,7,15),
+        Point(493504.953633058525156, 5258840.576098721474409),
+        [Filter('cloudy_pixel_percentage','<', 100), Filter('processing_level', '==', 'Level-2A')],
+    )]
+)
+def test_mapper_get_pixels_stac(collection, time_start, time_end, geom, metadata_filters):
+    """
+    Test extraction of single pixel values for a single point location over time.
+    """
+    Settings.USE_STAC = True
+    feature = Feature(
+        name='Test Area',
+        geometry=geom,
+        epsg=32632,
+        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, time_column='sensing_date')
+    mapper.query_scenes()
+
+    if collection == 'sentinel2-msi':
+        pixel_kwargs = {'pixel_reader': Sentinel2.read_pixels_from_safe}
+    # elif collection == 'sentinel1-grd':
+    #     pixel_kwargs = {'pixel_reader': Sentinel1.read_pixels_from_safe}
+    mapper.load_scenes(pixel_kwargs=pixel_kwargs)
+
+    assert isinstance(mapper.metadata, gpd.GeoDataFrame), 'expected a GeoDataFrame'
+    assert not mapper.metadata.empty, 'expected some items to be returned'
+    assert isinstance(mapper.data, gpd.GeoDataFrame), 'expected a GeoDataFrame'
+    assert mapper.metadata.shape[0] == len(mapper.data), \
+        'mis-match between length of metadata and data'
+    
\ No newline at end of file
diff --git a/tests/mapper/test_sentinel2_mapper.py b/tests/mapper/test_sentinel2_mapper.py
new file mode 100644
index 00000000..7dd1a13f
--- /dev/null
+++ b/tests/mapper/test_sentinel2_mapper.py
@@ -0,0 +1,77 @@
+'''
+Created on Jan 11, 2023
+
+@author: graflu
+'''
+
+import pytest
+
+from datetime import datetime
+from shapely.geometry import Polygon
+
+from eodal.config import get_settings
+from eodal.core.raster import RasterCollection
+from eodal.core.scene import SceneCollection
+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()
+
+@pytest.mark.parametrize(
+    'collection,time_start,time_end,geom,metadata_filters',
+    [(
+        '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')],
+    ),]
+)
+def test_sentinel2_mapper(collection, time_start, time_end, geom, metadata_filters):
+    """
+    Test the mapper class for handling Sentinel-2 data including mosaicing data
+    from two different Sentinel-2 tiles
+    """
+    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, sensor='sentinel2')
+    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': ['B04', 'B08']},
+        'scene_modifier': resample,
+        'scene_modifier_kwargs': {'target_resolution': 10}
+    }
+    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'
diff --git a/tests/metadata/test_stac_client.py b/tests/metadata/test_stac_client.py
index 35989ee6..bde5aa22 100644
--- a/tests/metadata/test_stac_client.py
+++ b/tests/metadata/test_stac_client.py
@@ -6,7 +6,6 @@
 
 from datetime import date
 
-from eodal.config import get_settings, STAC_Providers
 from eodal.metadata.stac import sentinel1, sentinel2
 from eodal.utils.sentinel1 import _url_to_safe_name
 from eodal.utils.sentinel2 import ProcessingLevels