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

Issue:
Landsat Support
This issue is a tracker for the progress on the full integration of Landsat into EOdal. The goal is to provide users with a similar smooth experience as with Sentinel-2 including the usage of the `Mapper` class. We focus on Collection-2 products (1982 to present) but also an integration of Collection-1 is feasible.

Thanks to the modular and sensor-agnostic design of the `core` and `mapper` submodule in EOdal, only few coding efforts are necessary.

- [x] define a class `Landsat` in `eodal.core.sensors` that inherits from the `RasterCollection` class
- [x] add a class-method `from_usgs` that allows to read Landsat data from GeoTiff files issued by [USGS](https://earthexplorer.usgs.gov/)
- [x] add a method to conveniently mask out clouds and other artifacts
- [x] ensure Landsat data can be fetched from STAC from Microsoft Planetary Computer
- [x] test the integration with the `Mapper` class


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index c0fde734..8298ab9b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -242,3 +242,4 @@ eodal.code-workspace
 .vscode/
 
 *.xml
+.vscode/launch.json
diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index f40b1fa3..b15d6dd6 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -9,7 +9,7 @@ version: 2
 build:
   os: ubuntu-22.04
   tools:
-    python: "3.10"
+    python: "3.11"
 
 # Build documentation in the docs/ directory with Sphinx
 sphinx:
@@ -21,7 +21,6 @@ sphinx:
 
 # Optionally declare the Python requirements required to build your docs
 python:
-   version: 3.10
    install:
       - requirements: requirements.txt
       - method: pip
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 635a97c7..db8f6e59 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -10,6 +10,16 @@ The format is based on `Keep a Changelog`_, and this project adheres to `Semanti
 
 Categories for changes are: Added, Changed, Deprecated, Removed, Fixed, Security.
 
+Version `0.2.2 < https://github.com/EOA-team/eodal/releases/tag/v0.2.2>`__
+--------------------------------------------------------------------------------
+
+Release date: 2023-08-XX
+
+- Added: Full support for USGS Landsat Collection-2 data given access to data >50 years of data starting with Landsat-1 (#66).
+- Changed: The spectral index module has undergone some flexibilization efforts. Now, custom bands can be passed for index calculation.
+- Fixed: Calls to new `pydantic_settings` have been added where necessary to provide compatability with `pydantic` version 2.+ (#73).
+
+
 Version `0.2.1 < https://github.com/EOA-team/eodal/releases/tag/v0.2.1>`__
 --------------------------------------------------------------------------------
 
diff --git a/eodal/core/band.py b/eodal/core/band.py
index 8473a45f..e4450e46 100644
--- a/eodal/core/band.py
+++ b/eodal/core/band.py
@@ -680,7 +680,7 @@ def _get_pixel_geometries(
 
         # convert to centroids if the geometries are not of type Point
         gdf.geometry = gdf.geometry.apply(
-            lambda x: x.centroid if x.type in ["Polygon", "MultiPolygon"] else x
+            lambda x: x.centroid if x.geom_type in ["Polygon", "MultiPolygon"] else x
         )
 
         return gdf
diff --git a/eodal/core/raster.py b/eodal/core/raster.py
index 560217d0..de196338 100644
--- a/eodal/core/raster.py
+++ b/eodal/core/raster.py
@@ -1538,7 +1538,7 @@ def mask(
             if inplace:
                 self[band_name].mask(mask=_mask, inplace=inplace)
             else:
-                band = self.get_band(band_name)
+                band = deepcopy(self.get_band(band_name))
                 collection.add_band(
                     band_constructor=band.mask, mask=_mask, inplace=inplace
                 )
@@ -1581,7 +1581,6 @@ def scale(
             if inplace:
                 self.collection[band_name].scale_data(inplace=inplace, **kwargs)
             else:
-                # TODO: there seems to be a bug here
                 band = self.get_band(band_name)
                 collection.add_band(
                     band_constructor=band.scale_data,
@@ -1599,18 +1598,32 @@ def join(self, other):
         pass
 
     def calc_si(
-        self, si_name: str, inplace: Optional[bool] = False
+        self, si_name: str, inplace: Optional[bool] = False,
+        band_mapping: Optional[Dict[str, str]] = {}
     ) -> Union[None, np.ndarray, np.ma.MaskedArray]:
         """
         Calculates a spectral index based on color-names (set as band aliases)
 
         :param si_name:
-            name of the spectral index to calculate (e.g., 'NDVI')
+            name of the spectral index to calculate (e.g., 'NDVI').
+        :param inplace:
+            if False, returns the calculated SI as numpy array (default), else
+            adds the calculated SI as a new band to the RasterCollection
+            object.
+        :param band_mapping:
+            ..versionadd 0.2.2::
+            optional band mapping to use different color names instead
+            of the default ones (e.g., nir_2 instead of nir_1 as default
+            near-infrared band, which might be helpful for some sensor
+            with more than one NIR band such as Sentinel-2 MSI).
+            See ~`eodal.core.spectral_indices.SpectralIndices` how
+            `band_mapping` must look like.
         :returns:
             ``np.ndarray`` or ``np.ma.MaskedArray`` if inplace is False, None
             otherwise (is added as band to the collection)
         """
-        si_values = SpectralIndices.calc_si(si_name, self)
+        spectral_indices = SpectralIndices(band_mapping)
+        si_values = spectral_indices.calc_si(si_name, self)
         # since SIs are floats by nature set the nodata value to np.nan
         nodata = np.nan
         if inplace:
diff --git a/eodal/core/sensors/__init__.py b/eodal/core/sensors/__init__.py
index feb38224..e685dfd9 100644
--- a/eodal/core/sensors/__init__.py
+++ b/eodal/core/sensors/__init__.py
@@ -1,3 +1,4 @@
 from .sentinel1 import Sentinel1
 from .sentinel2 import Sentinel2
 from .planet_scope import PlanetScope, SuperDove
+from .landsat import Landsat
diff --git a/eodal/core/sensors/landsat.py b/eodal/core/sensors/landsat.py
new file mode 100644
index 00000000..2582eaa9
--- /dev/null
+++ b/eodal/core/sensors/landsat.py
@@ -0,0 +1,704 @@
+"""
+This module contains the ``Landsat`` class that inherits from
+eodal's core ``RasterCollection`` class.
+
+The ``Landsat`` class enables reading one or more spectral bands from Landsat 5 to 9
+science products in L2 processing level.
+
+The class handles data in L1C and L2A processing level.
+
+Copyright (C) 2023 Lukas Valentin Graf
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from __future__ import annotations
+
+import cv2
+import io
+import json
+import numpy as np
+import pandas as pd
+import planetary_computer
+import requests
+import warnings
+
+from copy import deepcopy
+from datetime import datetime
+from pathlib import Path
+from typing import Dict, List, Optional
+from xml.etree import ElementTree
+
+from eodal.config import get_settings, STAC_Providers
+from eodal.core.band import Band
+from eodal.core.raster import RasterCollection, SceneProperties
+from eodal.utils.geometry import adopt_vector_features_to_mask
+from eodal.utils.constants import ProcessingLevels
+from eodal.utils.constants.landsat import (
+    band_resolution,
+    landsat_band_mapping,
+    platform_sensor_mapping)
+from eodal.utils.types import xml_to_dict_recursive
+
+Settings = get_settings()
+
+
+class Landsat(RasterCollection):
+    """
+    Class for working with Landsat science products (Level-2).
+    """
+
+    @staticmethod
+    def _get_available_bands(
+            in_dir: Path | Dict[str, str],
+    ) -> list[str]:
+        """
+        Check the available bands in the input directory and return the
+        bands found in a list.
+
+        :param in_dir:
+            Path to the directory containing the Landsat data or collection of
+            hyper-links in the case of STAC. If a dictionary is passed, we assume
+            that all files are stored in the same directory.
+        :return:
+            List of available bands.
+        """
+        if isinstance(in_dir, dict):
+            available_bands = list(in_dir.keys())
+            # remove non-band files
+            available_bands = [
+                x for x in available_bands if Path(x).suffix
+                not in ['.txt', '.xml', '.json']
+                and x != 'tilejson'
+                and x != 'rendered_preview']
+        elif isinstance(in_dir, Path):
+            available_bands = [f.stem.split('_')[-1] for f in in_dir.glob("*.TIF")]
+        else:
+            raise TypeError("in_dir must be a Path or a dictionary.")
+        return available_bands
+
+    @staticmethod
+    def _check_platform(
+            in_dir: Path | Dict[str, str]
+    ) -> str:
+        """
+        Check which Landsat platform (spacecraft) produced the data.
+        """
+        if isinstance(in_dir, dict):
+            # the green band should be always available
+            test_band = in_dir['green']['href'].split('/')[-1]
+        elif isinstance(in_dir, Path):
+            test_band = in_dir.glob("*_B*.TIF")[0]
+        # get the platform from the band name
+        platform = test_band.split('_')[0]
+        return platform
+
+    @staticmethod
+    def _get_metadata_json(
+            in_dir: Path | Dict[str, str]
+    ) -> dict:
+        """
+        Get the Landsat metadata JSON file from the input directory.
+        """
+        if isinstance(in_dir, dict):
+            metadata_url = planetary_computer.sign_url(in_dir['mtl.json']['href'])
+            response_code = requests.get(metadata_url).status_code
+            # check for JSON or XML
+            if response_code == 200:
+                metadata = json.loads(requests.get(metadata_url).content)
+            elif response_code == 404:
+                metadata_url = planetary_computer.sign_url(in_dir['mtl.xml']['href'])
+                metadata_xml_tree = ElementTree.parse(
+                    io.BytesIO(requests.get(metadata_url).content))
+                metadata_xml_root = metadata_xml_tree.getroot()
+                # convert to dictionary
+                metadata = xml_to_dict_recursive(metadata_xml_root)
+        elif isinstance(in_dir, Path):
+            fpath_json = in_dir.joinpath('MTL.json')
+            if fpath_json.exists():
+                metadata = json.loads(metadata_json.read_text())
+            else:
+                fpath_xml = in_dir.joinpath('MTL.xml')
+                metadata_xml_tree = ElementTree.parse(fpath_xml)
+                metadata_xml_root = metadata_xml_tree.getroot()
+                # convert to dictionary
+                metadata = xml_to_dict_recursive(metadata_xml_root)
+        return metadata
+
+    @staticmethod
+    def _product_uri_from_filename(
+            in_dir: Path | Dict[str, str]
+    ) -> str:
+        """
+        Get the product URI from the file name.
+
+        :param in_dir:
+            Path to the directory containing the Landsat data or collection of
+            hyper-links in the case of STAC. If a dictionary is passed, we assume
+            that all files are stored in the same directory.
+        :return:
+            Product URI.
+        """
+        if isinstance(in_dir, dict):
+            # the green band should be always available
+            test_band = in_dir['green']['href'].split('/')[-1]
+        elif isinstance(in_dir, Path):
+            test_band = in_dir.glob("*.TIF")[0]
+        # get the product URI from the band name
+        product_uri = '_'.join(test_band.split('_')[:-2])
+        return product_uri
+
+    @staticmethod
+    def _sensing_time_from_filename(
+            in_dir: Path | Dict[str, str]
+    ) -> datetime:
+        """
+        Get the sensing time from the file name.
+
+        :param in_dir:
+            Path to the directory containing the Landsat data or collection of
+            hyper-links in the case of STAC. If a dictionary is passed, we assume
+            that all files are stored in the same directory.
+        :return:
+            Sensing time.
+        """
+        if isinstance(in_dir, dict):
+            # the green band should be always available
+            test_band = in_dir['red']['href'].split('/')[-1]
+        elif isinstance(in_dir, Path):
+            test_band = in_dir.glob("*.TIF")[0]
+        # get the sensing time from the metadata file
+        metadata = Landsat._get_metadata_json(in_dir)['LANDSAT_METADATA_FILE']
+
+        # unpack dictionaries from lists
+        if isinstance(metadata, list):
+            metadata_as_dict = {}
+            for entry in metadata:
+                for k, v in entry.items():
+                    metadata_as_dict.update({k: {}})
+                    for inner_item in v:
+                        metadata_as_dict[k].update(inner_item)
+            metadata = metadata_as_dict
+
+        date = metadata['IMAGE_ATTRIBUTES']['DATE_ACQUIRED']
+        scene_center_time = metadata['IMAGE_ATTRIBUTES']['SCENE_CENTER_TIME']
+        combined = date + ' ' + scene_center_time.split('.')[0]
+        return datetime.strptime(combined, '%Y-%m-%d %H:%M:%S')
+
+    @staticmethod
+    def _sensor_from_platform(
+            platform: str
+    ) -> str:
+        """
+        Get the sensor name from the abbrevated platform name.
+
+        :param platform:
+            Landsat platform name (e.g., LC09).
+        :return:
+            Sensor name (e.g., Operational_Land_Imager).
+        """
+        landsat_platform_number = int(platform[2:4])
+        landsat_platform_full_name = f'LANDSAT_{landsat_platform_number}'
+        return platform_sensor_mapping[landsat_platform_full_name]
+
+    def _preprocess_band_selection(
+            self,
+            in_dir: Path | Dict[str, str],
+            band_selection: Optional[List[str]] = None,
+            read_qa: Optional[bool] = True,
+            read_atcor: Optional[bool] = False
+    ):
+        """
+        Process the user-defined band selection if any and return the available
+        bands in a pandas DataFrame.
+
+        :param in_dir:
+            Path to the directory containing the Landsat data or collection of
+            hyper-links in the case of STAC. If a dictionary is passed, we assume
+            that all files are stored in the same directory.
+        :param band_selection:
+            List of bands to be loaded. If None, all bands will be loaded.
+            The bands loaded depend on the sensor, e.g., TM, ETM, OLI, TIRS.
+            The sensor will be auto-detected. To make sure, the band selection
+            works for multiple sensors, use the aliases defined in
+            `~eodal.utils.constants.landsat_mapper.landsat_band_mapping`.
+            For instance, to load the VIS bands pass `["blue", "green", "red"]`.
+        :param read_qa:
+            If True, the QA bands will be loaded as well. Default is True.
+        :param read_atcor:
+            If True, the bands stemming from the atmospheric correction are read
+            as well. Default is False.
+        :return:
+            Pandas DataFrame with the available bands.
+        """
+        # get available bands
+        available_bands = self._get_available_bands(in_dir)
+
+        # check which platform and hence sensor is used
+        platform = self._check_platform(in_dir)
+
+        # get the sensor name from the platform
+        sensor = self._sensor_from_platform(platform)
+        sensor_bands = landsat_band_mapping[sensor]
+
+        # compare against the band selection. In case the user passed the
+        # sensor specific names (e.g,. "TM_B1") we need to convert them to
+        # the generic names (e.g., "blue")
+        # convert available bands to their generic names first
+        for available_band in available_bands:
+            if available_band in sensor_bands.keys():
+                available_bands[available_bands.index(available_band)] = \
+                    sensor_bands[available_band]
+
+        # then check the user-defined selection (if any)
+        if band_selection is not None:
+            for band_name in band_selection:
+                if band_name not in available_bands:
+                    if band_name in sensor_bands.keys():
+                        band_selection[band_selection.index(band_name)] = \
+                            sensor_bands[band_name]
+                    else:
+                        warnings.warn(
+                            f'Band {band_name} not found - ' +
+                            'remove it from selection')
+                        band_selection.remove(band_name)
+        else:
+            band_selection = available_bands
+
+        # check if the user wants to read the atmospheric correction bands
+        atcor_bands = landsat_band_mapping['atmospheric_correction']
+        for atcor_band in atcor_bands:
+            if atcor_band in available_bands:
+                # read_atcor is True:
+                # add those bands to the band selection that are available
+                if read_atcor:
+                    if not atcor_band in band_selection:
+                        band_selection.append(atcor_band)
+                # read_atcor is False:
+                # remove atmospheric correction bands
+                else:
+                    if atcor_band in band_selection:
+                        band_selection.remove(atcor_band)
+            else:
+                # band_selection.remove(atcor_band)
+                pass
+
+        # check if the user wants to read the QA bands
+        qa_bands = landsat_band_mapping['quality_flags']
+        for qa_band in qa_bands:
+            if qa_band in available_bands:
+                # read_qa is True:
+                # add those bands to the band selection that are available
+                if read_qa:
+                    if not qa_band in band_selection:
+                        band_selection.append(qa_band)
+                # read_atcor is False:
+                # remove atmospheric correction bands
+                else:
+                    if qa_band in band_selection:
+                        band_selection.remove(qa_band)
+            else:
+                # band_selection.remove(qa_band)
+                pass
+ 
+        # next, we need the paths to the bands
+        band_list = []
+        for band_name in band_selection:
+            if isinstance(in_dir, dict):
+                if Settings.STAC_BACKEND == STAC_Providers.MSPC:
+                    band_fpath = planetary_computer.sign(in_dir[band_name]["href"])
+                else:
+                    band_fpath = in_dir[band_name]["href"]
+            elif isinstance(in_dir, Path):
+                # for local resources, we have to use the sensor specific names
+                # to find the bands
+                sensor_specific_band_name = [
+                    k for k, v in sensor_bands.items() if v == 'red'][0]
+                band_fpath = in_dir.glob(f"*_{sensor_specific_band_name}.TIF")[0]
+            band_res = None
+            if band_name in sensor_bands.values():
+                band_res = band_resolution[sensor][band_name]
+            elif read_qa or band_name in qa_bands:
+                band_res = band_resolution['quality_flags'][band_name]
+            elif read_atcor or band_name in atcor_bands:
+                band_res = band_resolution['atmospheric_correction'][band_name]
+
+            item = {
+                'band_name': band_name,
+                'band_path': band_fpath,
+                'band_resolution': band_res}
+            band_list.append(item)
+
+        # construct pandas DataFrame with all band entries and return
+        band_df = pd.DataFrame(band_list)
+        band_df.dropna(inplace=True)
+
+        return band_df
+
+    @classmethod
+    def from_usgs(
+        cls,
+        in_dir: Path | Dict[str, str],
+        band_selection: Optional[List[str]] = None,
+        read_qa: Optional[bool] = True,
+        read_atcor: Optional[bool] = False,
+        apply_scaling: Optional[bool] = True,
+        **kwargs
+    ):
+        """
+        Load Landsat Level-2 science products (spectral bands, thermal bands,
+        QA bands) into a EOdal `RasterCollection.
+
+        :param in_dir:
+            Path to the directory containing the Landsat data or collection of
+            hyper-links in the case of STAC. If a dictionary is passed, we assume
+            that all files are stored in the same directory.
+        :param band_selection:
+            List of bands to be loaded. If None, all bands will be loaded.
+            The bands loaded depend on the sensor, e.g., TM, ETM, OLI, TIRS.
+            The sensor will be auto-detected. To make sure, the band selection
+            works for multiple sensors, use the aliases defined in
+            `~eodal.utils.constants.landsat_band_mapping`. For instance, to load
+            the VIS bands pass `["blue", "green", "red"]`.
+        :param read_qa:
+            If True, the QA bands will be loaded as well. Default is True.
+        :param read_atcor:
+            If True, the bands stemming from the atmospheric correction are read
+            as well. Default is False.
+        :param apply_scaling:
+            If True, scales reflectance values to the range [0, 1] from original
+            unsigned 16-bit integers. Default is True.
+        :param kwargs:
+            Additional keyword arguments passed to the `RasterCollection` constructor.
+        :return:
+            RasterCollection containing the Landsat bands.
+        """
+        # check band selection and determine the platform and sensor
+        band_df = cls._preprocess_band_selection(
+            cls,
+            in_dir=in_dir,
+            band_selection=band_selection,
+            read_qa=read_qa,
+            read_atcor=read_atcor)
+
+        # check the clipping extent of the raster with the lowest (coarsest) spatial
+        # resolution and remember it for all other bands with higher spatial
+        # resolutions.
+        # By doing so, it is ensured that all bands will be clipped to the same spatial
+        # extent regardless of their pixel size. This works since all Landsat bands
+        # share the same coordinate origin.
+        # get lowest spatial resolution (maximum pixel size) band
+        masking_after_read_required = False
+
+        if kwargs.get("vector_features") is not None:
+            bounds_df, shape_mask, lowest_resolution = adopt_vector_features_to_mask(
+                band_df=band_df,
+                vector_features=kwargs.get("vector_features")
+            )
+            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})
+
+        # get platform, sensor and sensing time
+        platform = cls._check_platform(in_dir=in_dir)
+        sensor = cls._sensor_from_platform(platform=platform)
+        sensing_time = cls._sensing_time_from_filename(in_dir=in_dir)
+        # get the product uri
+        product_uri = cls._product_uri_from_filename(in_dir=in_dir)
+        # determine processing level
+        if sensor in \
+            ['Multispectral_Scanner_System_L1-3',
+             'Multispectral_Scanner_System_L4-5']:
+            processing_level = ProcessingLevels.L1C
+        else:
+            processing_level = ProcessingLevels.L2A
+
+        # construct the SceneProperties object
+        scene_properties = SceneProperties(
+            platform=f'LANDSAT_{platform[-1]}',
+            sensor=sensor,
+            acquisition_time=sensing_time,
+            product_uri=product_uri,
+            processing_level=processing_level
+        )
+
+        # set proper scaling factors to allow for conversion to
+        # reflectance [0, 1]
+        gain, offset = 1, 0
+        if processing_level == ProcessingLevels.L2A:
+            gain, offset = 0.00001, 0.0
+
+        # loop over bands and add them to the collection of bands
+        landsat = cls(scene_properties=scene_properties)
+        for band_name in list(band_df.band_name):
+            # get entry from dataframe with file-path of band
+            band_safe = band_df[band_df.band_name == band_name]
+            band_fpath = band_safe.band_path.values[0]
+
+            # get alias name of band
+            band_alias = None
+            if band_name in landsat_band_mapping[sensor].values():
+                band_alias = [
+                    k for k, v in landsat_band_mapping[sensor].items()
+                    if v == band_name][0]
+            elif band_name in landsat_band_mapping['quality_flags'].keys() \
+                    or band_name in landsat_band_mapping['quality_flags'].values():
+                band_alias = band_name.lower()
+            elif band_name in landsat_band_mapping['atmospheric_correction'].keys() \
+                    or band_name in landsat_band_mapping['atmospheric_correction'
+                                                         ].values():
+                band_alias = band_name.lower()
+
+            # read band
+            try:
+                if band_name in landsat_band_mapping[sensor].values():
+                    kwargs.update({'scale': gain, 'offset': offset})
+                landsat.add_band(
+                    Band.from_rasterio,
+                    fpath_raster=band_fpath,
+                    band_idx=1,
+                    band_name_dst=band_name,
+                    band_alias=band_alias,
+                    **kwargs,
+                )
+            except Exception as e:
+                raise Exception(
+                    f"Could not add band {band_name} " +
+                    f"from {scene_properties.product_uri}: {e}"
+                )
+            # apply actual vector features if masking is required
+            if masking_after_read_required:
+                # nothing to do when the lowest resolution is passed
+                if band_safe.band_resolution.values == lowest_resolution:
+                    continue
+                # otherwise resample the mask of the lowest resolution to the
+                # current resolution using nearest neighbor interpolation
+                tmp = shape_mask.astype("uint8")
+                dim_resampled = (landsat[band_name].ncols, landsat[band_name].nrows)
+                res = cv2.resize(
+                    tmp, dim_resampled, interpolation=cv2.INTER_NEAREST_EXACT
+                )
+                # cast back to boolean
+                mask = res.astype("bool")
+                landsat.mask(mask=mask, bands_to_mask=[band_name], inplace=True)
+
+        # scaling of reflectance values (i.e., do not scale SCL)
+        if apply_scaling:
+            sel_bands = landsat.band_names
+            # only scale the reflectance bands
+            sel_bands = [x for x in sel_bands if x in
+                         landsat_band_mapping[sensor].values()]
+            landsat.scale(
+                inplace=True,
+                band_selection=sel_bands,
+                pixel_values_to_ignore=[landsat[landsat.band_names[0]].nodata],
+            )
+
+        return landsat
+
+    def calc_si(
+        self,
+        si_name: str,
+        inplace: Optional[bool] = False,
+        band_mapping: Optional[Dict[str, str]] = {}
+    ) -> None | np.ndarray | np.ma.MaskedArray:
+        """
+        Calculates a spectral index based on color-names (set as band aliases)
+
+        :param si_name:
+            name of the spectral index to calculate (e.g., 'NDVI').
+        :param inplace:
+            if False, returns the calculated SI as numpy array (default), else
+            adds the calculated SI as a new band to the RasterCollection
+            object.
+        :param band_mapping:
+            ..versionadd 0.2.2::
+            optional band mapping to use different color names instead
+            of the default ones (e.g., nir_2 instead of nir_1 as default
+            near-infrared band, which might be helpful for some sensor
+            with more than one NIR band such as Sentinel-2 MSI).
+            See ~`eodal.core.spectral_indices.SpectralIndices` how
+            `band_mapping` must look like.
+        """
+
+        # set color names properly before calling SI class
+        # (background: EOdal internally uses slightly different color names
+        # compared to Landsat collection 2 products. The differences
+        # only affect the NIR and SWIR bands)
+        if len(band_mapping) == 0:
+            band_mapping = {
+                'nir_1': 'nir08',
+                'swir_1': 'swir16',
+                'swir_2': 'swir22'
+            }
+        return RasterCollection.calc_si(
+            self, si_name=si_name, inplace=inplace,
+            band_mapping=band_mapping)
+
+    def get_water_mask(
+            self,
+            mask_band: Optional[str] = 'qa_pixel',
+            inplace: Optional[bool] = False,
+            name_water_mask: Optional[str] = 'water_mask'
+    ) -> Band | None:
+        """
+        Get a water mask from the pixel quality band (bit 7).
+
+        NOTE:
+            Since the `mask_band` can be set to *any* `Band` it is also
+            possible to use a different water mask, e.g., from
+            a custom classifier as long as the bit map used for the classes
+            is the same.
+            See also
+            https://www.usgs.gov/media/images/landsat-collection-2-pixel-quality-assessment-bit-index
+            for more details.
+
+        :param mask_band:
+            The name of the band to use for masking. Default is 'qa_pixel'.
+        :param inplace:
+            Whether to return a new `Band` instance or add a `Band` to the
+            current `Landsat` instance. Default is False.
+        :param name_water_mask:
+            The name of the water mask band. Default is 'water_mask'.
+        :return:
+            A `Band` instance containing the water mask or `None` if
+            `inplace` is True.
+        """
+        water_mask = self.mask_from_qa_bits(
+            bit_range=(7, 7),
+            band_name=mask_band
+        )
+        water_mask_band = Band(
+                band_name=name_water_mask,
+                values=water_mask,
+                geo_info=self[mask_band].geo_info,
+                nodata=0,
+            )
+        if inplace:
+            self.add_band(water_mask_band)
+            return None
+        else:
+            return water_mask_band
+
+    def mask_from_qa_bits(
+            self,
+            bit_range: tuple,
+            band_name: Optional[str] = 'qa_pixel'
+    ) -> np.ndarray | np.ma.MaskedArray:
+        """
+        Extract Quality Accurance (QA) bits from the QA band of
+        Landsat to obtain a binary mask.
+
+        This function can be used to obtain a binary cloud or
+        cloud shadow mask, for example.
+
+        :param bit_range:
+            A tuple of two integers indicating the start and end bit
+            of the bit range to extract.
+        :param band_name:
+            The name of the QA band. Default is 'qa_pixel'.
+        :return:
+            A 2-dimensional numpy array containing the obtained mask.
+        """
+        # get the band
+        band = deepcopy(self[band_name].values)
+        band_mask = None
+        if isinstance(band, np.ndarray):
+            band_data = band
+        elif isinstance(band, np.ma.MaskedArray):
+            band_data = band.data
+            band_mask = band.mask
+        # compute the bits to extract
+        pattern = 0
+        for i in range(bit_range[0], bit_range[1] + 1):
+            pattern += 2 ** i
+        # assign the extracted bits to the band
+        band_data = np.bitwise_and(band_data, pattern)
+        band_data = np.right_shift(band_data, bit_range[0])
+
+        # create a masked array if the input was also masked
+        if band_mask is not None:
+            band_data = np.ma.MaskedArray(
+                data=band_data,
+                mask=band_mask
+            )
+
+        return band_data
+
+    def mask_clouds_and_shadows(
+            self,
+            bands_to_mask: Optional[list[str]] = None,
+            cloud_classes: Optional[list[int]] = [1, 2, 3, 5],
+            mask_band: Optional[str] = 'qa_pixel',
+            **kwargs
+    ) -> None | Landsat:
+        """
+        A wrapper around the inherited ``mask`` method to mask cirrus,
+        cloud and cloud shadow pixels in Landsat scenes based on the
+        pixel quality band.
+
+        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 as long as the bit map used for the classes
+            is the same.
+            See also
+            https://www.usgs.gov/media/images/landsat-collection-2-pixel-quality-assessment-bit-index
+            for more details.
+
+        NOTE:
+            You might also use the mask function from
+            `eodal.core.raster.RasterCollection` directly.
+
+        :param bands_to_mask:
+            A list of bands to mask. Default is all bands.
+        :param cloud_classes:
+            A list of integers indicating the cloud classes to mask.
+            Default is [1, 2, 3, 5].
+        :param mask_band:
+            The name of the band to use for masking. Default is 'qa_pixel'.
+        :param kwargs:
+            optional kwargs to pass to `~eodal.core.raster.RasterCollection.mask`
+        :returns:
+            depending on `inplace` (passed in the kwargs) a new `Sentinel2` instance
+            or None
+        """
+        if bands_to_mask is None:
+            bands_to_mask = self.band_names
+        # the mask band should never be masked
+        if mask_band in bands_to_mask:
+            bands_to_mask.remove(mask_band)
+        try:
+            mask_list = []
+            for cloud_class in cloud_classes:
+                # construct the bit range
+                bit_range = (cloud_class, cloud_class)
+                # get the binary mask as boolean array
+                mask = self.mask_from_qa_bits(
+                    bit_range=bit_range,
+                    band_name=mask_band).astype('bool')
+                # add the mask to the list
+                mask_list.append(mask)
+
+            # combine the masks
+            cloud_mask = np.any(mask_list, axis=0)
+            # mask the bands
+            return self.mask(
+                mask=cloud_mask,
+                bands_to_mask=bands_to_mask,
+                **kwargs
+            )
+        except Exception as e:
+            raise Exception(f"Could not mask clouds and shadows: {e}")
diff --git a/eodal/core/sensors/sentinel2.py b/eodal/core/sensors/sentinel2.py
index 6a19c928..5d1593dc 100644
--- a/eodal/core/sensors/sentinel2.py
+++ b/eodal/core/sensors/sentinel2.py
@@ -35,7 +35,6 @@
 from matplotlib import colors
 from numbers import Number
 from pathlib import Path
-from rasterio.mask import raster_geometry_mask
 from shapely.geometry import box
 from typing import Any, Dict, Optional, List, Tuple, Union
 
@@ -51,6 +50,7 @@
     SCL_Classes,
 )
 from eodal.utils.decorators import prepare_point_features
+from eodal.utils.geometry import adopt_vector_features_to_mask
 from eodal.utils.exceptions import BandNotFoundError
 from eodal.utils.sentinel2 import (
     get_S2_bandfiles_with_res,
@@ -59,7 +59,6 @@
     get_S2_acquistion_time_from_safe,
     get_S2_processing_baseline_from_safe,
 )
-from eodal.core.utils.geometry import convert_3D_2D
 from eodal.config import get_settings
 from eodal.utils.sentinel2 import _url_to_safe_name
 
@@ -294,73 +293,14 @@ def from_safe(
         masking_after_read_required = False
 
         if kwargs.get("vector_features") is not None:
-            lowest_resolution = band_df_safe["band_resolution"].max()
-            if band_df_safe["band_resolution"].unique().shape[0] > 1:
-                if kwargs.get("vector_features") is not None:
-                    low_res_band = band_df_safe[
-                        band_df_safe["band_resolution"] == lowest_resolution
-                    ].iloc[0]
-                    # get vector feature(s) for spatial subsetting
-                    vector_features = kwargs.get("vector_features")
-                    if isinstance(vector_features, Path):
-                        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(
-                            geometry=vector_features.copy()
-                        )
-                    else:
-                        raise TypeError(
-                            "Geometry must be vector file, GeoSeries or GeoDataFrame"
-                        )
-
-                    # drop Nones in geometry column
-                    none_idx = vector_features_df[
-                        vector_features_df.geometry == None  # noqa: E711
-                    ].index
-                    vector_features_df.drop(index=none_idx, inplace=True)
-
-                    with rio.open(low_res_band.band_path, "r") as src:
-                        # convert to raster CRS
-                        raster_crs = src.crs
-                        vector_features_df.to_crs(crs=raster_crs, inplace=True)
-                        # check if the geometry contains the z (3rd) dimension. If yes
-                        # convert it to 2d to avoid an error poping up from rasterio
-                        vector_features_geom = convert_3D_2D(
-                            vector_features_df.geometry
-                        )
-                        shape_mask, transform, window = raster_geometry_mask(
-                            dataset=src,
-                            shapes=vector_features_geom,
-                            all_touched=True,
-                            crop=True,
-                        )
-                    # get upper left coordinates rasterio takes for the band
-                    # with the coarsest spatial resolution
-                    ulx_low_res, uly_low_res = transform.c, transform.f
-                    # reconstruct the lower right corner
-                    llx_low_res = ulx_low_res + window.width * transform.a
-                    lly_low_res = uly_low_res + window.height * transform.e
-
-                    # overwrite original vector features' bounds in the S2 scene
-                    # geometry of the lowest spatial resolution
-                    low_res_feature_bounds_s2_grid = box(
-                        minx=ulx_low_res,
-                        miny=lly_low_res,
-                        maxx=llx_low_res,
-                        maxy=uly_low_res,
-                    )
-                    # update bounds and pass them on to the kwargs
-                    bounds_df = gpd.GeoDataFrame(
-                        geometry=[low_res_feature_bounds_s2_grid],
-                    )
-                    bounds_df.set_crs(crs=raster_crs, inplace=True)
-                    # 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})
+            bounds_df, shape_mask, lowest_resolution = adopt_vector_features_to_mask(
+                band_df=band_df_safe,
+                vector_features=kwargs.get("vector_features")
+            )
+            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:
@@ -409,8 +349,7 @@ def from_safe(
 
             # get color name and set it as alias
             color_name = s2_band_mapping[band_name]
-            kwargs.update({"scale": 1})
-            kwargs.update({"offset": 0})
+            kwargs.update({"scale": 1, "offset": 0})
 
             # store wavelength information per spectral band
             if band_name != "SCL":
@@ -425,8 +364,7 @@ def from_safe(
                 kwargs.update({"wavelength_info": wvl_info})
                 # do not apply the gain and offset factors from the spectral bands
                 # to the SCL file
-                kwargs.update({"scale": gain})
-                kwargs.update({"offset": offset})
+                kwargs.update({"scale": gain, "offset": offset})
 
             # read band
             try:
@@ -440,7 +378,7 @@ def from_safe(
                 )
             except Exception as e:
                 raise Exception(
-                    f"Could not add band {band_name} from {in_dir.name}: {e}"
+                    f"Could not add band {band_name} from {in_dir}: {e}"
                 )
             # apply actual vector features if masking is required
             if masking_after_read_required:
@@ -645,7 +583,7 @@ def mask_clouds_and_shadows(
         cloud_classes: Optional[List[int]] = [1, 2, 3, 7, 8, 9, 10, 11],
         mask_band: Optional[str] = "SCL",
         **kwargs,
-    ) -> Sentinel2:
+    ) -> None | Sentinel2:
         """
         A Wrapper around the inherited ``mask`` method to mask clouds,
         shadows, water and snow based on (by default) the SCL band.
@@ -678,16 +616,16 @@ def mask_clouds_and_shadows(
         # 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")
+            bands_to_mask.remove(mask_band)
         try:
             return self.mask(
                 mask=mask_band,
                 mask_values=cloud_classes,
                 bands_to_mask=bands_to_mask,
-                **kwargs,
+                **kwargs
             )
         except Exception as e:
-            raise Exception(f"Could not apply cloud mask: {e}")
+            raise Exception(f"Could not mask clouds and shadows: {e}")
 
     def get_scl_stats(self) -> pd.DataFrame:
         """
diff --git a/eodal/core/spectral_indices.py b/eodal/core/spectral_indices.py
index 782551b8..e7ee5a2a 100644
--- a/eodal/core/spectral_indices.py
+++ b/eodal/core/spectral_indices.py
@@ -24,31 +24,56 @@
 import numpy as np
 
 from numpy import inf
-from typing import List, Union
+from typing import Dict, List, Optional, Union
 
 
 class SpectralIndices(object):
     """generic spectral indices"""
 
     # define color names (for optical sensors)
-    blue = "blue"
-    green = "green"
-    red = "red"
-    red_edge_1 = "red_edge_1"
-    red_edge_2 = "red_edge_2"
-    red_edge_3 = "red_edge_3"
-    nir_1 = "nir_1"
-    nir_2 = "nir_2"
-    swir_1 = "swir_1"
-    swir_2 = "swir_2"
-
-    # polarizations (for SAR)
-    vv = "VV"
-    vh = "VH"
-    hh = "HH"
-
-    @classmethod
-    def get_si_list(cls) -> List[str]:
+    def __init__(self, band_mapping: Optional[Dict[str, str]] = {}):
+        """
+        Class constructor. To override the default
+        bands pass custom band mappings as a dictionary:
+        
+        Example
+        -------
+        To use Sentinel-2 band nir_2 (Band 8A) instead of
+        Sentinel-2 band nir_1 (Band 8) for, e.g., NDVI
+        calculation pass
+
+        .. highlight:: python
+        .. code-block:: python
+
+            # use nir_2 instead of nir_1
+            # the other bands remain unaffected
+            band_mapping = {'nir_1': 'nir_2'}
+            spectral_indices = SpectralIndices(band_mapping)
+
+        :param band_mapping:
+            optional band mapping to override default band
+            setting for SI calculation (see example above)
+        """
+
+        # common optical bands
+        self.blue = band_mapping.get("blue", "blue")
+        self.green = band_mapping.get("green", "green")
+        self.red = band_mapping.get("red", "red")
+        self.red_edge_1 = band_mapping.get("red_edge_1", "red_edge_1")
+        self.red_edge_2 = band_mapping.get("red_edge_2", "red_edge_2")
+        self.red_edge_3 = band_mapping.get("red_edge_3", "red_edge_3")
+        self.nir_1 = band_mapping.get("nir_1", "nir_1")
+        self.nir_2 = band_mapping.get("nir_2", "nir_2")
+        self.swir_1 = band_mapping.get("swir_1", "swir_1")
+        self.swir_2 = band_mapping.get("swir_2", "swir_2")
+    
+        # polarizations (for SAR)
+        self.vv = band_mapping.get("VV", "VV")
+        self.vh = band_mapping.get("VH", "VH")
+        self.hh = band_mapping.get("HH", "HH")
+
+    
+    def get_si_list(self) -> List[str]:
         """
         Returns a list of implemented Spectral Indices (SIs)
 
@@ -57,12 +82,11 @@ def get_si_list(cls) -> List[str]:
         """
         return [
             x
-            for x in dir(cls)
+            for x in dir(self)
             if not x.startswith("__") and not x.endswith("__") and not x.islower()
         ]
 
-    @classmethod
-    def calc_si(cls, si: str, collection: dict) -> Union[np.ndarray, np.ma.MaskedArray]:
+    def calc_si(self, si: str, collection) -> Union[np.ndarray, np.ma.MaskedArray]:
         """
         Calculates the selected spectral index (SI) for
         spectral band data derived from `~eodal.core.RasterCollection`.
@@ -71,24 +95,33 @@ def calc_si(cls, si: str, collection: dict) -> Union[np.ndarray, np.ma.MaskedArr
         :param si:
             name of the selected vegetation index (e.g., NDVI). Raises
             an error if the vegetation index is not implemented/ found.
+        :param collection:
+            ~`eodal.core.raster.RasterCollection` with spectral bands.
         :returns:
             2d ``numpy.ndarray`` or ``np.ma.MaskedArray`` with VI values
             (depends on array-type of the input)
         """
+        # SI calculation
+        if not hasattr(self, si.upper()):
+            raise NotImplementedError(
+                f'SI {si} not found.\nSee `SpectralIndices()' +
+                '.get_si_list()` for a list of currently implemented SIs.')
+
+        si_fun = getattr(self, si.upper())
+        si_data = si_fun.__call__(collection)
+
+        # post-processing
         try:
-            si_fun = eval(f"cls.{si.upper()}")
-            si_data = si_fun.__call__(collection)
-            # replace infinity values with nan
+            # replace infinity values with NaN
             si_data[si_data == inf] = np.nan
-            # replace masked values with nodata (might look weird otherwise)
+            # replace masked values with no-data (might look weird otherwise)
             if isinstance(si_data, np.ma.MaskedArray):
                 si_data.data[si_data.mask] = np.nan
         except Exception as e:
-            raise NotImplementedError(e)
+            raise ValueError from e
         return si_data
 
-    @classmethod
-    def NDVI(cls, collection) -> np.array:
+    def NDVI(self, collection) -> np.array:
         """
         Calculates the Normalized Difference Vegetation Index
         (NDVI) using the red and the near-infrared (NIR) channel.
@@ -99,13 +132,12 @@ def NDVI(cls, collection) -> np.array:
             NDVI values
         """
 
-        nir = collection.get(cls.nir_1).values.astype("float")
-        red = collection.get(cls.red).values.astype("float")
+        nir = collection.get(self.nir_1).values.astype("float")
+        red = collection.get(self.red).values.astype("float")
         ndvi = (nir - red) / (nir + red)
         return ndvi
 
-    @classmethod
-    def EVI(cls, collection):
+    def EVI(self, collection):
         """
         Calculates the Enhanced Vegetation Index (EVI) following the formula
         provided by Huete et al. (2002)
@@ -115,9 +147,9 @@ def EVI(cls, collection):
         :returns:
             EVI values
         """
-        blue = collection.get(cls.blue).values.astype("float")
-        nir = collection.get(cls.nir_1).values.astype("float")
-        red = collection.get(cls.red).values.astype("float")
+        blue = collection.get(self.blue).values.astype("float")
+        nir = collection.get(self.nir_1).values.astype("float")
+        red = collection.get(self.red).values.astype("float")
         numerator = nir - red
         denominator = nir + 6 * red - 7.5 * blue + 1
         # values larger 1 and smaller -1 might occur (e.g., on artificial
@@ -127,8 +159,7 @@ def EVI(cls, collection):
         evi[evi < -1.0] = -1.0
         return evi
 
-    @classmethod
-    def MSAVI(cls, collection) -> np.array:
+    def MSAVI(self, collection) -> np.array:
         """
         Calculates the Modified Soil-Adjusted Vegetation Index
         (MSAVI). MSAVI is sensitive to the green leaf area index
@@ -139,12 +170,11 @@ def MSAVI(cls, collection) -> np.array:
         :returns:
             MSAVI values
         """
-        nir = collection.get(cls.nir_1).values.astype("float")
-        red = collection.get(cls.red).values.astype("float")
+        nir = collection.get(self.nir_1).values.astype("float")
+        red = collection.get(self.red).values.astype("float")
         return 0.5 * (2 * nir + 1 - np.sqrt((2 * nir + 1) ** 2 - 8 * (nir - red)))
 
-    @classmethod
-    def CI_GREEN(cls, collection) -> np.array:
+    def CI_GREEN(self, collection) -> np.array:
         """
         Calculates the green chlorophyll index (CI_green).
         It is sensitive to canopy chlorophyll concentration (CCC).
@@ -155,12 +185,11 @@ def CI_GREEN(cls, collection) -> np.array:
             CI-green values
         """
 
-        nir = collection.get(cls.nir_1).values.astype("float")
-        green = collection.get(cls.green).values.astype("float")
+        nir = collection.get(self.nir_1).values.astype("float")
+        green = collection.get(self.green).values.astype("float")
         return (nir / green) - 1
 
-    @classmethod
-    def MTCARI_OSAVI(cls, collection) -> np.array:
+    def MTCARI_OSAVI(self, collection) -> np.array:
         """
         Calculates the ratio of the Transformed Chlorophyll Index (TCARI)
         and the Optimized Soil-Adjusted Vegetation Index (OSAVI). It is sensitive
@@ -171,10 +200,10 @@ def MTCARI_OSAVI(cls, collection) -> np.array:
         :returns:
             TCARI/OSAVI values
         """
-        green = collection.get(cls.green).values.astype("float")
-        red = collection.get(cls.red).values.astype("float")
-        red_edge_1 = collection.get(cls.red_edge_1).values.astype("float")
-        red_edge_3 = collection.get(cls.red_edge_3).values.astype("float")
+        green = collection.get(self.green).values.astype("float")
+        red = collection.get(self.red).values.astype("float")
+        red_edge_1 = collection.get(self.red_edge_1).values.astype("float")
+        red_edge_3 = collection.get(self.red_edge_3).values.astype("float")
 
         TCARI = 3 * (
             (red_edge_1 - red) - 0.2 * (red_edge_1 - green) * (red_edge_1 / red)
@@ -183,8 +212,7 @@ def MTCARI_OSAVI(cls, collection) -> np.array:
         tcari_osavi = TCARI / OSAVI
         return tcari_osavi
 
-    @classmethod
-    def NDRE(cls, collection) -> np.array:
+    def NDRE(self, collection) -> np.array:
         """
         Calculates the Normalized Difference Red Edge (NDRE). It extends
         the capabilities of the NDVI for middle and late season crops.
@@ -194,12 +222,11 @@ def NDRE(cls, collection) -> np.array:
         :returns:
             NDRE values
         """
-        red_edge_1 = collection.get(cls.red_edge_1).values.astype("float")
-        red_edge_3 = collection.get(cls.red_edge_3).values.astype("float")
+        red_edge_1 = collection.get(self.red_edge_1).values.astype("float")
+        red_edge_3 = collection.get(self.red_edge_3).values.astype("float")
         return (red_edge_3 - red_edge_1) / (red_edge_3 + red_edge_1)
 
-    @classmethod
-    def MCARI(cls, collection):
+    def MCARI(self, collection):
         """
         Calculates the Modified Chlorophyll Absorption Ratio Index (MCARI).
         It is sensitive to leaf chlorophyll concentration (LCC).
@@ -209,13 +236,12 @@ def MCARI(cls, collection):
         :returns:
             MCARI values
         """
-        green = collection.get(cls.green).values.astype("float")
-        red = collection.get(cls.red).values.astype("float")
-        red_edge_1 = collection.get(cls.red_edge_1).values.astype("float")
+        green = collection.get(self.green).values.astype("float")
+        red = collection.get(self.red).values.astype("float")
+        red_edge_1 = collection.get(self.red_edge_1).values.astype("float")
         return ((red_edge_1 - red) - 0.2 * (red_edge_1 - green)) * (red_edge_1 / red)
 
-    @classmethod
-    def BSI(cls, collection):
+    def BSI(self, collection):
         """
         Calculates the Bare Soil Index (BSI).
 
@@ -224,14 +250,13 @@ def BSI(cls, collection):
         :returns:
             BSI values
         """
-        blue = collection.get(cls.blue).values.astype("float")
-        red = collection.get(cls.red).values.astype("float")
-        nir = collection.get(cls.nir_1).values.astype("float")
-        swir_1 = collection.get(cls.swir_1).values.astype("float")
+        blue = collection.get(self.blue).values.astype("float")
+        red = collection.get(self.red).values.astype("float")
+        nir = collection.get(self.nir_1).values.astype("float")
+        swir_1 = collection.get(self.swir_1).values.astype("float")
         return ((swir_1 + red) - (nir + blue)) / ((swir_1 + red) + (nir + blue))
 
-    @classmethod
-    def VARI(cls, collection):
+    def VARI(self, collection):
         """
         Calculates the Visible Atmospherically Resistant Index (VARI)
 
@@ -240,13 +265,12 @@ def VARI(cls, collection):
         :returns:
             BSI values
         """
-        blue = collection.get(cls.blue).values.astype("float")
-        green = collection.get(cls.green).values.astype("float")
-        red = collection.get(cls.red).values.astype("float")
+        blue = collection.get(self.blue).values.astype("float")
+        green = collection.get(self.green).values.astype("float")
+        red = collection.get(self.red).values.astype("float")
         return (green - red) / (green + red - blue)
 
-    @classmethod
-    def NDYI(cls, collection):
+    def NDYI(self, collection):
         """
         Calculates the Normalized Difference Yellowness Index
 
@@ -255,12 +279,11 @@ def NDYI(cls, collection):
         :returns:
             NDYI values
         """
-        blue = collection.get(cls.blue).values.astype("float")
-        green = collection.get(cls.green).values.astype("float")
+        blue = collection.get(self.blue).values.astype("float")
+        green = collection.get(self.green).values.astype("float")
         return (green - blue) / (green + blue)
 
-    @classmethod
-    def NDWI(cls, collection):
+    def NDWI(self, collection):
         """
         Calculates the Normalized Difference Yellowness Index
 
@@ -269,35 +292,32 @@ def NDWI(cls, collection):
         :returns:
             NDYI values
         """
-        green = collection.get(cls.green).values.astype("float")
-        nir = collection.get(cls.nir_1).values.astype("float")
+        green = collection.get(self.green).values.astype("float")
+        nir = collection.get(self.nir_1).values.astype("float")
         return (green - nir) / (green + nir)
 
-    @classmethod
-    def CR(cls, collection):
+    def CR(self, collection):
         """
         Calculates the Cross Ratio  VH/VV
         """
-        vh = collection.get(cls.vh).values.astype("float")
-        vv = collection.get(cls.vv).values.astype("float")
+        vh = collection.get(self.vh).values.astype("float")
+        vv = collection.get(self.vv).values.astype("float")
         return vh / vv
 
-    @classmethod
-    def NHI(cls, collection):
+    def NHI(self, collection):
         """
-        Calculates the Normalized Humidity Index. First described for detecting ponds with vegetation in and on top.
+        Calculates the Normalized Humidity Index. First described for detecting
+        ponds with vegetation in and on top.
         :param collection:
             reflectance in the 'green' and 'swir_1' channel
         :returns:
             NHI values
         """
-        swir_1 = collection.get(cls.swir_1).values.astype("float")
-        green = collection.get(cls.green).values.astype("float")
+        swir_1 = collection.get(self.swir_1).values.astype("float")
+        green = collection.get(self.green).values.astype("float")
         return (swir_1 - green)/(swir_1 + green)
-        
-        
-    @classmethod
-    def NDTI(cls, collection):
+
+    def NDTI(self, collection):
         """
         Calculates the Normalized Tillage Index. Useful for tillage and crop residues.
         :param collection:
@@ -305,28 +325,25 @@ def NDTI(cls, collection):
         :returns:
             NDTI values
         """
-        swir_1 = collection.get(cls.swir_1).values.astype("float")
-        swir_2 = collection.get(cls.swir_2).values.astype("float")
+        swir_1 = collection.get(self.swir_1).values.astype("float")
+        swir_2 = collection.get(self.swir_2).values.astype("float")
         return (swir_1 - swir_2)/(swir_1 + swir_2)
-        
-        
-    @classmethod
-    def NDRI(cls, collection):
+
+    def NDRI(self, collection):
         """
-        Calculates the Normalized Difference Residue Index. Useful for detecting crop residues.
+        Calculates the Normalized Difference Residue Index. Useful for detecting
+        crop residues.
+
         :param collection:
             reflectance in the 'red' and 'swir_2' channel
         :returns:
             NDRI values
         """
-        red = collection.get(cls.red).values.astype("float")
-        swir_2 = collection.get(cls.swir_2).values.astype("float")
+        red = collection.get(self.red).values.astype("float")
+        swir_2 = collection.get(self.swir_2).values.astype("float")
         return (red - swir_2)/(red + swir_2)
-        
-        
-         
-    @classmethod
-    def GNDVI(cls, collection):
+
+    def GNDVI(self, collection):
         """
         Calculates the Green Normalized Difference Vegetation Index 
         :param collection:
@@ -334,6 +351,6 @@ def GNDVI(cls, collection):
         :returns:
             GNDVI values
         """
-        nir = collection.get(cls.nir_1).values.astype("float")
-        green = collection.get(cls.green).values.astype("float")
+        nir = collection.get(self.nir_1).values.astype("float")
+        green = collection.get(self.green).values.astype("float")
         return (nir - green)/(nir + green)
diff --git a/eodal/mapper/mapper.py b/eodal/mapper/mapper.py
index 76245968..72fcd811 100644
--- a/eodal/mapper/mapper.py
+++ b/eodal/mapper/mapper.py
@@ -356,12 +356,13 @@ def query_scenes(self) -> None:
         feature_wgs84 = self.mapper_configs.feature.to_epsg(4326)
         bbox = box(*feature_wgs84.geometry.bounds)
 
-        # determine platform
+        # determine platform and collection
         platform = self.mapper_configs.platform
+        collection = self.mapper_configs.collection
 
         # put kwargs together
         kwargs = {
-            "platform": platform,
+            "collection": collection,
             "time_start": self.mapper_configs.time_start,
             "time_end": self.mapper_configs.time_end,
             "bounding_box": bbox,
@@ -629,7 +630,7 @@ def _load_scenes_collection(
                     scoll.add_scene(scene)
                 except KeyError:
                     logger.warn(
-                        f"Scene with ID {merged_scene_properties.product_uri} "
+                        f"Scene with ID {scene.scene_properties.product_uri} "
                         "already added to SceneCollection - continue"
                     )
                     continue
diff --git a/eodal/metadata/stac/__init__.py b/eodal/metadata/stac/__init__.py
index 8919dbf3..c4dce3e5 100644
--- a/eodal/metadata/stac/__init__.py
+++ b/eodal/metadata/stac/__init__.py
@@ -1,1 +1,1 @@
-from .client import query_stac, sentinel1, sentinel2
+from .client import query_stac, sentinel1, sentinel2, landsat
diff --git a/eodal/metadata/stac/client.py b/eodal/metadata/stac/client.py
index 23f37a74..183b594a 100644
--- a/eodal/metadata/stac/client.py
+++ b/eodal/metadata/stac/client.py
@@ -18,7 +18,9 @@
 from eodal.config import get_settings, STAC_Providers
 from eodal.mapper.filter import Filter
 from eodal.utils.decorators import prepare_bbox
+from eodal.utils.geometry import box_from_transform, prepare_gdf
 from eodal.utils.reprojection import infer_utm_zone
+from eodal.utils.timestamps import datetime_to_date
 
 Settings = get_settings()
 
@@ -79,7 +81,7 @@ def query_stac(
     )
     # fetch items and convert them to GeoDataFrame, drop all records with
     # too many clouds
-    items = search.get_all_items()
+    items = search.item_collection()
     item_json = items.to_dict()
     scenes = item_json["features"]
     return scenes
@@ -99,22 +101,41 @@ def _filter_criteria_fulfilled(
         `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"]:
+    criteria_fulfilled = []
+    for _filter in metadata_filters:
+        if _filter.entity in ["processing_level", "product_type"]:
             continue
-        if filter.entity not in metadata_dict.keys():
+        if _filter.entity not in metadata_dict.keys():
             warnings.warn(
-                f"{filter.entity} could not be retrieved from STAC -> skipping filter"
+                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
+        # check if the metadata item to filter is a list or single value
+        if not isinstance(metadata_dict[_filter.entity], list):
+            if isinstance(metadata_dict[_filter.entity], str):
+                eval_str = \
+                    f'metadata_dict["{_filter.entity}"] ' + \
+                    f'{_filter.operator} "{_filter.value}"'
+            else:
+                eval_str = \
+                    f'metadata_dict["{_filter.entity}"] ' + \
+                    f'{_filter.operator} {_filter.value}'
+            condition_met = eval(eval_str)
+        else:
+            # TODO: this is not really elegant and might not deliver always the
+            # results we are looking for ...
+            # convert _filter.value to list if it is not yet
+            value = _filter.value
+            if not isinstance(value, list):
+                value = [value]
+            # redefine the operators to work with sets
+            if _filter.operator == '==':
+                condition_met = set(value).issubset(metadata_dict[_filter.entity])
+            elif _filter.operator == '!=':
+                condition_met = not set(value).issubset(metadata_dict[_filter.entity])
+        criteria_fulfilled.append(condition_met)
+
+    return all(criteria_fulfilled)
 
 
 @prepare_bbox
@@ -142,8 +163,8 @@ def sentinel2(metadata_filters: List[Filter], **kwargs) -> gpd.GeoDataFrame:
 
     # query STAC catalog
     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 scenes found and apply the Filters provided
@@ -178,9 +199,7 @@ def sentinel2(metadata_filters: List[Filter], **kwargs) -> gpd.GeoDataFrame:
             "scene_id": scene_id,
             "spacecraft_name": props[s2.platform],
             "tile_id": tile_id,
-            "sensing_date": datetime.strptime(
-                props[s2.sensing_time].split("T")[0], "%Y-%m-%d"
-            ).date(),
+            "sensing_date": datetime_to_date(props[s2.sensing_time]),
             "cloudy_pixel_percentage": props[s2.cloud_cover],
             "epsg": props[s2.epsg],
             "sensing_time": datetime.strptime(
@@ -199,9 +218,7 @@ def sentinel2(metadata_filters: List[Filter], **kwargs) -> gpd.GeoDataFrame:
             metadata_list.append(meta_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)
+    return prepare_gdf(metadata_list)
 
 
 @prepare_bbox
@@ -237,7 +254,6 @@ def sentinel1(metadata_filters: List[Filter], **kwargs) -> gpd.GeoDataFrame:
 
     # query the catalog
     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 = []
@@ -245,9 +261,7 @@ def sentinel1(metadata_filters: List[Filter], **kwargs) -> gpd.GeoDataFrame:
         metadata_dict = scene["properties"]
         metadata_dict["assets"] = scene["assets"]
         metadata_dict["sensing_time"] = metadata_dict["datetime"]
-        metadata_dict["sensing_date"] = datetime.strptime(
-            metadata_dict["sensing_time"].split("T")[0], "%Y-%m-%d"
-        ).date()
+        metadata_dict["sensing_date"] = datetime_to_date(metadata_dict['sensing_time'])
         del metadata_dict["datetime"]
         # infer EPSG code of the scene in UTM coordinates from its bounding box
         bbox = box(*scene["bbox"])
@@ -259,6 +273,60 @@ def sentinel1(metadata_filters: List[Filter], **kwargs) -> gpd.GeoDataFrame:
             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)
+    return prepare_gdf(metadata_list)
+
+
+@prepare_bbox
+def landsat(metadata_filters: List[Filter], **kwargs) -> gpd.GeoDataFrame:
+    """
+    Landsat specific STAC query function to retrieve mapper from MSPC.
+
+    :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 Landsat scenes
+    """
+
+    if Settings.STAC_BACKEND != STAC_Providers.MSPC:
+        raise ValueError("This method currently requires Microsoft Planetary Computer")
+
+    # query the catalog
+    stac_kwargs = kwargs.copy()
+    scenes = query_stac(**stac_kwargs)
+
+    if len(scenes) == 0:
+        raise ValueError(f'STAC query returned no results! {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_to_date(metadata_dict['sensing_time'])
+        del metadata_dict["datetime"]
+        # apply filters
+        append_scene = _filter_criteria_fulfilled(metadata_dict, metadata_filters)
+        if append_scene:
+            metadata_list.append(metadata_dict)
+        # reconstruct the scene footprint
+        transform = metadata_dict['proj:transform']
+        shape = metadata_dict['proj:shape']
+        metadata_dict['geom'] = box_from_transform(transform, shape)
+        metadata_dict['epsg'] = metadata_dict['proj:epsg']
+
+        # extract angles
+        metadata_dict['sun_azimuth_angle'] = metadata_dict['view:sun_azimuth']
+        del metadata_dict['view:sun_azimuth']
+        metadata_dict['sun_zenith_angle'] = 90. - metadata_dict['view:sun_elevation']
+        del metadata_dict['view:sun_elevation']
+        metadata_dict['sensor_zenith_angle'] = metadata_dict['view:off_nadir']
+        del metadata_dict['view:off_nadir']
+
+    if len(metadata_list) == 0:
+        raise ValueError(
+            f'No scenes fulfilling filter criteria: {metadata_filters}')
+    return prepare_gdf(metadata_list)
+    
\ No newline at end of file
diff --git a/eodal/utils/constants/landsat.py b/eodal/utils/constants/landsat.py
new file mode 100644
index 00000000..fdf5821e
--- /dev/null
+++ b/eodal/utils/constants/landsat.py
@@ -0,0 +1,154 @@
+"""
+Defines some static attributes of Landsat Collection 2 products.
+
+Copyright (C) 2023 Lukas Valentin Graf
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from enum import Enum
+
+# available processing levels
+class ProcessingLevels(Enum):  # noqa: F811
+    L1 = "LEVEL1"
+    L2 = "LEVEL2"
+
+
+# Landsat processing levels as defined in the metadatabase
+ProcessingLevelsDB = {"L1": "Level-1", "L2": "Level-2"}
+
+# band mapping organized by sensor type
+landsat_band_mapping = {
+    "Multispectral_Scanner_System_L1-3": {
+        "B4": "green",
+        "B5": "red",
+        "B6": "nir08",
+        "B7": "nir09"},
+    "Multispectral_Scanner_System_L4-5": {
+        "B1": "green",
+        "B2": "red",
+        "B3": "nir08",
+        "B4": "nir09"
+    },
+    "Thematic_Mapper": {
+        "B1": "blue",
+        "B2": "green",
+        "B3": "red",
+        "B4": "nir08",
+        "B5": "swir16",
+        "B6": "lwir",
+        "B7": "swir22"},
+    "Enhanced_Thematic_Mapper_Plus": {
+        "B1": "blue",
+        "B2": "green",
+        "B3": "red",
+        "B4": "nir08",
+        "B5": "swir16",
+        "B6": "lwir",
+        "B7": "swir22"},
+    "Operational_Land_Imager": {
+        "B1": "coastal",
+        "B2": "blue",
+        "B3": "green",
+        "B4": "red",
+        "B5": "nir08",
+        "B6": "swir16",
+        "B7": "swir22",
+        "B10": "lwir11"},
+    "quality_flags": {
+        "qa_pixel": "PIXEL",
+        "qa_aerosol": "AEROSOL",
+        "qa_radsat": "RADSAT",
+        "qa": "QA",
+        "cloud_qa": "cloud_qa"},
+    "atmospheric_correction": {
+        "cdist": "CDIST",
+        "drad": "DRAD",
+        "emis": "EMIS",
+        "emsd": "EMSD",
+        "lwir": "LWIR",
+        "trad": "TRAD",
+        "urad": "URAD",
+        "atran": "ATRAN",
+        "atmos_opacity": "OPACITY",
+        "ang": "ANG"}
+}
+
+# TODO: L4 and L5 actually have two instruments (TM and MSS)
+platform_sensor_mapping = {
+    "LANDSAT_1": "Multispectral_Scanner_System_L1-3",
+    "LANDSAT_2": "Multispectral_Scanner_System_L1-3",
+    "LANDSAT_3": "Multispectral_Scanner_System_L1-3",
+    "LANDSAT_4": "Thematic_Mapper",
+    "LANDSAT_5": "Thematic_Mapper",
+    "LANDSAT_7": "Enhanced_Thematic_Mapper_Plus",
+    "LANDSAT_8": "Operational_Land_Imager",
+    "LANDSAT_9": "Operational_Land_Imager"}
+
+# spatial resolutions of the Landsat bands organized by sensor and product
+# in meters
+band_resolution = {
+    "Multispectral_Scanner_System_L1-3": {
+        "green": 80,
+        "red": 80,
+        "nir08": 80,
+        "nir09": 80},
+    "Multispectral_Scanner_System_L4-5": {
+        "green": 80,
+        "red": 80,
+        "nir08": 80,
+        "nir09": 80},
+    "Thematic_Mapper": {
+        "blue": 30,
+        "green": 30,
+        "red": 30,
+        "nir08": 30,
+        "swir16": 30,
+        "lwir": 120,
+        "swir22": 30},
+    "Enhanced_Thematic_Mapper_Plus": {
+        "blue": 30,
+        "green": 30,
+        "red": 30,
+        "nir08": 30,
+        "swir16": 30,
+        "lwir": 60,
+        "swir22": 30},
+    "Operational_Land_Imager": {
+        "coastal": 30,
+        "blue": 30,
+        "green": 30,
+        "red": 30,
+        "nir08": 30,
+        "swir16": 30,
+        "swir22": 30,
+        "lwir11": 100},
+    "quality_flags": {
+        "qa_pixel": 30,
+        "qa_aerosol": 30,
+        "qa_radsat": 30,
+        "qa": 30,
+        "cloud_qa": 30},
+    "atmospheric_correction": {
+        "cdist": 30,
+        "drad": 30,
+        "emis": 30,
+        "emsd": 30,
+        "lwir": 30,
+        "trad": 30,
+        "urad": 30,
+        "atran": 30,
+        "atmos_opacity": 30,
+        "ang": 30}
+}
diff --git a/eodal/utils/constants/sentinel2.py b/eodal/utils/constants/sentinel2.py
index 6bbb0f23..7ae4cee6 100644
--- a/eodal/utils/constants/sentinel2.py
+++ b/eodal/utils/constants/sentinel2.py
@@ -31,6 +31,7 @@ class ProcessingLevels(Enum):  # noqa: F811
 ProcessingLevelsDB = {"L1C": "Level-1C", "L2A": "Level-2A"}
 
 # native spatial resolution of the S2 bands per processing level
+# (in meters)
 band_resolution = {
     ProcessingLevels.L1C: {
         "B01": 60,
diff --git a/eodal/utils/geometry.py b/eodal/utils/geometry.py
index f137aec3..c57a5d1c 100644
--- a/eodal/utils/geometry.py
+++ b/eodal/utils/geometry.py
@@ -1,7 +1,7 @@
 """
 Vector geometry operations.
 
-Copyright (C) 2022 Lukas Valentin Graf
+Copyright (C) 2022/2023 Lukas Valentin Graf
 
 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -20,10 +20,16 @@
 
 import geopandas as gpd
 import json
+import pandas as pd
+import rasterio as rio
 
 from copy import deepcopy
+from pathlib import Path
+from rasterio.mask import raster_geometry_mask
 from shapely.geometry import box, Point, Polygon
 
+from eodal.core.utils.geometry import convert_3D_2D
+
 
 def box_to_geojson(gdf: gpd.GeoDataFrame | Polygon) -> str:
     """
@@ -45,3 +51,130 @@ def box_to_geojson(gdf: gpd.GeoDataFrame | Polygon) -> str:
         bbox_poly = deepcopy(gdf)
     bbox_json = gpd.GeoSeries([bbox_poly]).to_json()
     return json.loads(bbox_json)["features"][0]["geometry"]
+
+
+def box_from_transform(
+        transform: list[float] | tuple,
+        shape: list[int] | tuple[int, int]
+) -> Polygon:
+    """
+    Get a `shapely` box geometry from a transform vector
+    and raster shape.
+
+    ..versionadd:: 0.3.0
+
+    :param transform:
+        affine-like transform parameters.
+    :param shape:
+        raster shape (nrows, ncols).
+    :returns:
+        resulting geometry as `shapely.geometry.box` object.
+    """
+    # get pixel resolution
+    pixres_x = transform[0]
+    pixres_y = transform[4]
+    # get origin of raster in spatial reference system coordinates
+    minx = transform[2]
+    maxy = transform[5]
+    # get extent of raster in spatial reference system coordinates
+    maxx = minx + (pixres_x * shape[1])
+    miny = maxy + (pixres_y * shape[0])
+    # construct the box geometry
+    return box(minx, miny, maxx, maxy)
+
+
+def prepare_gdf(
+        metadata_list: list[dict]
+) -> gpd.GeoDataFrame:
+    """
+    Convert list of metadata dictionaries to a GeoDataFrame.
+
+    ..versionadd:: 0.3.0
+
+    :param metadata_list:
+        list of metadata entries as dictionaries.
+    :returns:
+        resulting GeoDataFrame sorted by sensing time.
+    """
+    df = pd.DataFrame(metadata_list)
+    df.sort_values(by="sensing_time", inplace=True)
+    return gpd.GeoDataFrame(df, geometry="geom", crs=df.epsg.unique()[0])
+
+
+def adopt_vector_features_to_mask(
+        band_df: pd.DataFrame,
+        vector_features: gpd.GeoDataFrame | gpd.GeoSeries | Path
+) -> tuple[gpd.GeoDataFrame, tuple[int, int], int | float]:
+    """
+    Adopt the vector features used for clipping and/or masking data
+    to the spatial resolution of the band with the coarsest spatial
+    resolution. This is necessary to avoid artifacts (namely different
+    spatial extents) in the data caused by the spatial subsetting with
+    different pixel sizes.
+
+    ..versionadd:: 0.3.0
+
+    :param band_df:
+        DataFrame containing the band metadata.
+    :param vector_features:
+        vector features to be used for masking.
+    :returns:
+        Updated vector features, shape of the resulting mask and
+        spatial resolution of the band with the coarsest spatial
+        resolution.
+    """
+    # get lowest spatial resolution
+    lowest_resolution = band_df["band_resolution"].max()
+    # get band with lowest spatial resolution
+    low_res_band = band_df[band_df["band_resolution"] == lowest_resolution].iloc[0]
+    # get vector feature(s) for spatial subsetting
+    if isinstance(vector_features, Path):
+        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(geometry=vector_features.copy())
+    else:
+        raise TypeError(
+            "Geometry must be vector file, GeoSeries or GeoDataFrame"
+        )
+
+    # drop Nones in geometry column
+    none_idx = vector_features_df[
+        vector_features_df.geometry == None].index  # noqa: E711
+    vector_features_df.drop(index=none_idx, inplace=True)
+
+    with rio.open(low_res_band.band_path, "r") as src:
+        # convert to raster CRS
+        raster_crs = src.crs
+        vector_features_df.to_crs(crs=raster_crs, inplace=True)
+        # check if the geometry contains the z (3rd) dimension. If yes
+        # convert it to 2d to avoid an error poping up from rasterio
+        vector_features_geom = convert_3D_2D(vector_features_df.geometry)
+        shape_mask, transform, window = raster_geometry_mask(
+            dataset=src,
+            shapes=vector_features_geom,
+            all_touched=True,
+            crop=True,
+        )
+    # get upper left coordinates rasterio takes for the band
+    # with the coarsest spatial resolution
+    ulx_low_res, uly_low_res = transform.c, transform.f
+    # reconstruct the lower right corner
+    llx_low_res = ulx_low_res + window.width * transform.a
+    lly_low_res = uly_low_res + window.height * transform.e
+
+    # overwrite original vector features' bounds in the S2 scene
+    # geometry of the lowest spatial resolution
+    low_res_feature_bounds_s2_grid = box(
+        minx=ulx_low_res,
+        miny=lly_low_res,
+        maxx=llx_low_res,
+        maxy=uly_low_res,
+    )
+    # update bounds and pass them on to the kwargs
+    bounds_df = gpd.GeoDataFrame(
+        geometry=[low_res_feature_bounds_s2_grid],
+    )
+    bounds_df.set_crs(crs=raster_crs, inplace=True)
+    return bounds_df, shape_mask, lowest_resolution
diff --git a/eodal/utils/timestamps.py b/eodal/utils/timestamps.py
new file mode 100644
index 00000000..8d3b7a85
--- /dev/null
+++ b/eodal/utils/timestamps.py
@@ -0,0 +1,21 @@
+"""
+Utility functions for working with timestamps.
+"""
+
+from datetime import datetime, date
+
+
+def datetime_to_date(
+        timestamp: str,
+) -> date:
+    """
+    Convert `timestamp` received from STAC to `date` objects.
+
+    :param timestamp:
+        Timestamp received from STAC (e.g., '2023-05-16T10:15:46.380697Z').
+    :return:
+        Date object (e.g., 2023-05-16).
+    """
+    return datetime.strptime(
+        timestamp.split("T")[0], "%Y-%m-%d"
+    ).date()
diff --git a/eodal/utils/types.py b/eodal/utils/types.py
index edeb362c..33ede7d3 100644
--- a/eodal/utils/types.py
+++ b/eodal/utils/types.py
@@ -17,10 +17,21 @@
 along with this program.  If not, see <http://www.gnu.org/licenses/>.
 """
 
-from pandas import DataFrame
-from typing import Any
-from typing import Dict
-from typing import NewType
+from xml.etree.ElementTree import Element
 
 
-S2Scenes = NewType("S2Scenes", Dict[Any, DataFrame])
+def xml_to_dict_recursive(root: Element) -> dict:
+    """
+    Convert a XML Element to a dictionary recursively.
+
+    Taken from:
+    https://stackoverflow.com/questions/2148119/how-to-convert-an-xml-string-to-a-dictionary
+
+    :param root: XML Element to convert to dictionary.
+    :return: resulting dictionary
+    """
+
+    if len(list(root)) == 0:
+        return {root.tag:root.text}
+    else:
+        return {root.tag:list(map(xml_to_dict_recursive, list(root)))}
diff --git a/scripts/landsat_mapper.py b/scripts/landsat_mapper.py
new file mode 100644
index 00000000..7252c65b
--- /dev/null
+++ b/scripts/landsat_mapper.py
@@ -0,0 +1,139 @@
+"""
+Script to extract a collection of Landsat 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 Landsat scenes
+* merging data from different 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 Landsat Collection-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) 2023 Lukas Valentin Graf
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from datetime import datetime
+from pathlib import Path
+from shapely.geometry import box
+
+from eodal.config import get_settings
+from eodal.core.sensors import Landsat
+from eodal.mapper.feature import Feature
+from eodal.mapper.filter import Filter
+from eodal.mapper.mapper import Mapper, MapperConfigs
+
+Settings = get_settings()
+# we use STAC, i.e., Microsoft Planetary Computer
+Settings.USE_STAC = True
+
+
+def prepocess_landsat_scene(
+        ds: Landsat
+) -> Landsat:
+    """
+    Mask clouds and cloud shadows in a Landsat scene based
+    on the 'qa_pixel' band.
+
+    NOTE:
+        Depending on your needs, the pre-processing function can be
+        fully customized using the full power of EOdal and its
+        interfacing libraries!
+
+    :param ds:
+        Landsat scene before cloud mask applied.
+    :return:
+        Landsat scene with clouds and cloud shadows masked.
+    """
+    ds.mask_clouds_and_shadows(inplace=True)
+    return ds
+
+
+if __name__ == "__main__":
+
+    # user-inputs
+    # -------------------------- Collection -------------------------------
+    collection = 'landsat-c2-l2'
+
+    # ---------------------- Spatial Feature  ------------------------------
+    bbox = box(*[7.0, 47.0, 8.0, 48.0])  # can be also shp, gpkg, etc.
+    feature = Feature(
+        name='landsat-test',
+        geometry=bbox,
+        epsg=4326,
+        attributes={})
+
+    # ------------------------- Time Range ---------------------------------
+    time_start = datetime(2023, 5, 1)
+    time_end = datetime(2023, 5, 31)
+
+    # ------------------------- Metadata Filters ---------------------------
+    metadata_filters = [
+        Filter('eo:cloud_cover', '<', 70),
+        Filter('landsat:wrs_row', '==', '027')
+    ]
+
+    # set up the Mapper configuration
+    mapper_configs = MapperConfigs(
+        metadata_filters=metadata_filters,
+        collection=collection,
+        feature=feature,
+        time_start=time_start,
+        time_end=time_end)
+    # get a new mapper instance
+    mapper = Mapper(mapper_configs)
+
+    # query the scenes available (no I/O of scenes, this only fetches metadata)
+    mapper.query_scenes()
+    # the metadata is stored as a GeoDataFrame
+    mapper.metadata
+
+    # we tell EOdal how to load the Landsat scenes using `Landsat.from_usgs`
+    # and pass on some kwargs, e.g., the selection of bands we want to read.
+    # in addition, we tell EOdal to mask out clouds and shadows and the fly
+    # while reading the data using the qa_pixel band (therefore, we set the
+    # `read_qa` flag to True.
+    scene_kwargs = {
+        'scene_constructor': Landsat.from_usgs,
+        'scene_constructor_kwargs': {
+            'band_selection': ["blue", "green", "red", "nir08"],
+            'read_qa': True},
+        'scene_modifier': prepocess_landsat_scene,
+        'scene_modifier_kwargs': {}}
+
+    # now we load the scenes
+    mapper.load_scenes(scene_kwargs=scene_kwargs)
+    # the scenes are loaded into a EOdal SceneCollection object
+    mapper.data
+
+    # test scene for cloud masking
+    scene = mapper.data[mapper.data.timestamps[-1]]
+
+    # the scenes can be plotted
+    f_scenes = mapper.data.plot(['red', 'green', 'blue'], figsize=(15,5))
+
+    # make the SceneCollection obtained persistent so that we do not have to re-run
+    # the STAC query all the time we use the data.
+    fpath = Path('../data/sample_mapper_data.pkl')
+    with open(fpath, 'wb+') as dst:
+        dst.write(mapper.data.to_pickle())
diff --git a/scripts/landsat_mapper_level1.py b/scripts/landsat_mapper_level1.py
new file mode 100644
index 00000000..73fe5a23
--- /dev/null
+++ b/scripts/landsat_mapper_level1.py
@@ -0,0 +1,114 @@
+"""
+Script to extract a collection of Landsat 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 Landsat scenes
+* merging data from different 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 Landsat Collection-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) 2023 Lukas Valentin Graf
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from datetime import datetime
+from pathlib import Path
+from shapely.geometry import box
+
+from eodal.config import get_settings
+from eodal.core.sensors import Landsat
+from eodal.mapper.feature import Feature
+from eodal.mapper.filter import Filter
+from eodal.mapper.mapper import Mapper, MapperConfigs
+
+Settings = get_settings()
+# we use STAC, i.e., Microsoft Planetary Computer
+Settings.USE_STAC = True
+
+
+if __name__ == "__main__":
+
+    # user-inputs
+    # -------------------------- Collection -------------------------------
+    collection = 'landsat-c2-l1'
+
+    # ---------------------- Spatial Feature  ------------------------------
+    bbox = box(*[8.4183, 47.2544, 8.7639, 47.5176])  # can be also shp, gpkg, etc.
+    feature = Feature(
+        name='zurich area',
+        geometry=bbox,
+        epsg=4326,
+        attributes={})
+
+    # ------------------------- Time Range ---------------------------------
+    time_start = datetime(1972, 9, 1)
+    time_end = datetime(1972, 10, 31)
+
+    # ------------------------- Metadata Filters ---------------------------
+    metadata_filters = [
+        Filter('eo:cloud_cover', '<', 70)
+    ]
+
+    # set up the Mapper configuration
+    mapper_configs = MapperConfigs(
+        metadata_filters=metadata_filters,
+        collection=collection,
+        feature=feature,
+        time_start=time_start,
+        time_end=time_end)
+    # get a new mapper instance
+    mapper = Mapper(mapper_configs)
+
+    # query the scenes available (no I/O of scenes, this only fetches metadata)
+    mapper.query_scenes()
+    # the metadata is stored as a GeoDataFrame
+    mapper.metadata
+
+    # we tell EOdal how to load the Landsat scenes using `Landsat.from_usgs`
+    # and pass on some kwargs, e.g., the selection of bands we want to read.
+    # in addition, we tell EOdal to mask out clouds and shadows and the fly
+    # while reading the data using the qa_pixel band (therefore, we set the
+    # `read_qa` flag to True.
+    scene_kwargs = {
+        'scene_constructor': Landsat.from_usgs,
+        'scene_constructor_kwargs': {
+            'band_selection': ["green", "red", "nir08"],
+            'read_qa': True, 'apply_scaling': False}
+    }
+
+    # now we load the scenes
+    mapper.load_scenes(scene_kwargs=scene_kwargs)
+    # the scenes are loaded into a EOdal SceneCollection object
+    mapper.data
+
+    # the scenes can be plotted
+    f_scenes = mapper.data.plot(['nir08', 'red', 'green'], figsize=(15,15))
+    f_scenes.savefig('img/landsat1_Zurich.png', dpi=300)
+
+    # make the SceneCollection obtained persistent so that we do not have to re-run
+    # the STAC query all the time we use the data.
+    fpath = Path('landsat1_mapper_data.pkl')
+    with open(fpath, 'wb+') as dst:
+        dst.write(mapper.data.to_pickle())
diff --git a/scripts/sentinel2_mapper.py b/scripts/sentinel2_mapper.py
index 89cfa4ca..18ba1733 100644
--- a/scripts/sentinel2_mapper.py
+++ b/scripts/sentinel2_mapper.py
@@ -33,8 +33,6 @@
 along with this program.  If not, see <http://www.gnu.org/licenses/>.
 """
 
-import geopandas as gpd
-
 from datetime import datetime
 from pathlib import Path
 from shapely.geometry import box
@@ -62,9 +60,9 @@ def preprocess_sentinel2_scenes(
     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!
+        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.
@@ -80,10 +78,6 @@ def preprocess_sentinel2_scenes(
 
 if __name__ == '__main__':
 
-    import os
-    cwd = Path(__file__).absolute().parent.parent
-    os.chdir(cwd)
-
     # user-inputs
     # -------------------------- Collection -------------------------------
     collection: str = 'sentinel2-msi'
@@ -125,12 +119,28 @@ def preprocess_sentinel2_scenes(
     # the metadata is loaded into a GeoPandas GeoDataFrame
     mapper.metadata
 
-    # load the least cloudy scene available from STAC
+    # optional: get the least cloudy scene
+    # to apply it uncomment the statement below. This
+    # will return just a single scene no matter how long the time period chosen.
+    # mapper.metadata = mapper.metadata[
+    #       mapper.metadata.cloudy_pixel_percentage ==
+    #       mapper.metadata.cloudy_pixel_percentage.min()].copy()
+
+    # we tell EOdal how to load the Sentinel-2 scenes using `Sentinel2.from_safe`
+    # and pass on some kwargs, e.g., the selection of bands we want to read.
+    # Moreover, we tell EOdal to mask out clouds and shadows using the
+    # Scene Classification Layer (SCL) which is part of the L2A product.
+    # Therefore, we set the `read_scl` flag to True.
     scene_kwargs = {
         'scene_constructor': Sentinel2.from_safe,
-        'scene_constructor_kwargs': {'band_selection': [
-            "B02", "B03", "B04"], 'read_scl': False}}
-
+        'scene_constructor_kwargs': {'band_selection':
+                                     ['B02', 'B03', 'B04', 'B08'],
+                                     'read_scl': True},
+        'scene_modifier': preprocess_sentinel2_scenes,
+        'scene_modifier_kwargs': {'target_resolution': 10}
+    }
+
+    # load the scenes available from STAC
     mapper.load_scenes(scene_kwargs=scene_kwargs)
     # the data loaded into `mapper.data` as a EOdal SceneCollection
     mapper.data


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index 85241aa6..520a4a59 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -5,8 +5,8 @@
 import os
 import pytest
 import requests
+import shutil
 
-from distutils import dir_util
 from pathlib import Path
 
 from eodal.config import get_settings
@@ -40,7 +40,7 @@ def datadir(tmppath, request):
     test_dir, _ = os.path.splitext(filename)
 
     if os.path.isdir(test_dir):
-        dir_util.copy_tree(test_dir, str(tmppath))
+        shutil.copytree(test_dir, str(tmppath))
 
     return tmppath
 
diff --git a/tests/core/test_landsat.py b/tests/core/test_landsat.py
new file mode 100644
index 00000000..f09c4d35
--- /dev/null
+++ b/tests/core/test_landsat.py
@@ -0,0 +1,55 @@
+'''
+Tests for the Landsat class.
+'''
+
+import geopandas as gpd
+import pytest
+
+from datetime import datetime
+from eodal.core.sensors import Landsat
+from eodal.mapper.filter import Filter
+from eodal.metadata.stac.client import landsat as landsat_stac
+from shapely.geometry import box
+
+
+def test_landsat_from_usgs():
+
+    # query STAC for a custom region
+    collection = 'landsat_mapper-c2-l2'
+    bbox = box(*[7.0, 47.0, 8.0, 48.0])
+
+    years = [1995, 2023]
+
+    for year in years:
+        time_start = datetime(year, 5, 1)
+        time_end = datetime(year, 5, 30)
+    
+        metadata_filters = [
+            Filter('eo:cloud_cover', '<', 70),
+            Filter('landsat_mapper:wrs_row', '==', '028')
+        ]
+
+        landsat_items = landsat_stac(
+            metadata_filters=metadata_filters,
+            collection=collection,
+            bounding_box=bbox,
+            time_start=time_start,
+            time_end=time_end)
+    
+        # read only a part of the test scene
+        landsat_scene_item = landsat_items.iloc[0]
+        gdf = gpd.GeoSeries([bbox], crs=4326)
+        landsat = Landsat.from_usgs(
+            in_dir=landsat_scene_item['assets'],
+            vector_features=gdf,
+            read_qa=False,
+            read_atcor=False,
+            band_selection=['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'qa']
+        )
+        assert landsat.band_names == \
+            ['coastal', 'blue', 'green', 'red', 'nir08', 'swir16', 'qa'], 'wrong band names'
+        assert landsat.band_aliases == \
+            ['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'qa'], 'wrong band aliases'
+        for band_name in landsat.band_names[:-1]:
+            assert 0 < landsat[band_name].values.min() <= 1, 'wrong value'
+            assert 0 < landsat[band_name].values.max() <= 1, 'wrong value'
diff --git a/tests/core/test_raster_apply.py b/tests/core/test_raster_apply.py
index eeefe0f8..c3a55a49 100644
--- a/tests/core/test_raster_apply.py
+++ b/tests/core/test_raster_apply.py
@@ -1,7 +1,5 @@
 '''
-Created on Nov 20, 2022
-
-@author: graflu
+Tests for RasterCollection.apply() method.
 '''
 
 import numpy as np
diff --git a/tests/core/test_si_calculation.py b/tests/core/test_si_calculation.py
index 3dece2ae..d2f85a6a 100644
--- a/tests/core/test_si_calculation.py
+++ b/tests/core/test_si_calculation.py
@@ -1,9 +1,54 @@
 import pytest
 
+import numpy as np
+
 from eodal.config import get_settings
+from eodal.core.spectral_indices import SpectralIndices
 from eodal.core.sensors import Sentinel2
 
 
+def test_si_list():
+
+    si_list = SpectralIndices().get_si_list()
+    assert len(si_list) > 0, 'too few SIs found'
+    assert 'NDVI' in si_list, 'expected at least to find the NDVI'
+
+
+def test_si_bandmapping(get_s2_safe_l2a, get_polygons_2):
+    settings = get_settings()
+    settings.USE_STAC = False
+
+    s2_dir = get_s2_safe_l2a()
+    vector_features = get_polygons_2()
+
+    s2 = Sentinel2().from_safe(
+        in_dir=s2_dir,
+        vector_features=vector_features
+    )
+
+    # calculate default NDVI using band 8 (nir_1)
+    ndvi_b08 = s2.calc_si('ndvi')
+
+    # calculate the NDVI using band 8A instead
+    # this should fail because Band 8A is still in 20 m spatial resolution
+    # while the red band is 10 m.
+    with pytest.raises(ValueError):
+        s2.calc_si('ndvi', band_mapping={'nir_1': 'nir_2'})
+    # but it should succeed if we resample the bands to 10 m
+    s2.resample(target_resolution=10, inplace=True)
+    ndvi_b8a = s2.calc_si('ndvi', band_mapping={'nir_1': 'nir_2'})
+
+    assert (ndvi_b08 != ndvi_b8a).any(), \
+        'B08 and B8A should not yield the exactly same result'
+    assert not (ndvi_b08 == ndvi_b8a).all(), \
+        'B08 and B8A should not yield the exactly same result'
+
+    assert np.min(ndvi_b08) >= -1 and np.min(ndvi_b8a) >= -1, \
+        'implausible NDVI values encountered'
+    assert np.max(ndvi_b08) <= 1 and np.max(ndvi_b8a) <= 1, \
+        'implausible NDVI values encountered'
+
+
 @pytest.mark.parametrize('apply_scaling', [False, True])
 def test_sentinel2_vi(get_s2_safe_l2a, get_polygons_2, apply_scaling):