diff --git a/eodal/__meta__.py b/eodal/__meta__.py
index 2377e1cd..5ea3e413 100644
--- a/eodal/__meta__.py
+++ b/eodal/__meta__.py
@@ -15,4 +15,4 @@
 description = "Earth Observation Data Analysis Library"  # One-liner
 url = "https://github.com/EOA-team/eodal"  # your project home-page
 license = "GNU General Public License version 3"  # See https://choosealicense.com
-version = "0.2.2"
+version = "0.2.3"
diff --git a/eodal/core/algorithms.py b/eodal/core/algorithms.py
index e088c4fc..0bf2887f 100644
--- a/eodal/core/algorithms.py
+++ b/eodal/core/algorithms.py
@@ -22,6 +22,7 @@
 import eodal
 import os
 import geopandas as gpd
+import numpy as np
 import uuid
 
 from pathlib import Path
@@ -38,9 +39,10 @@
 
 def _get_crs_and_attribs(
     in_file: Path, **kwargs
-) -> Tuple[GeoInfo, List[Dict[str, Any]]]:
+) -> Tuple[GeoInfo, List[Dict[str, Any]], str]:
     """
-    Returns the ``GeoInfo`` from a multi-band raster dataset
+    Returns the ``GeoInfo``, attributes and data type from
+    a multi-band raster dataset.
 
     :param in_file:
         raster datasets from which to extract the ``GeoInfo`` and
@@ -55,7 +57,8 @@ def _get_crs_and_attribs(
     ds = RasterCollection.from_multi_band_raster(fpath_raster=in_file, **kwargs)
     geo_info = ds[ds.band_names[0]].geo_info
     attrs = [ds[x].get_attributes() for x in ds.band_names]
-    return geo_info, attrs
+    dtype = ds[ds.band_names[0]].values.dtype
+    return geo_info, attrs, dtype
 
 
 def merge_datasets(
@@ -110,7 +113,7 @@ def merge_datasets(
     crs_list = []
     attrs_list = []
     for dataset in datasets:
-        geo_info, attrs = _get_crs_and_attribs(in_file=dataset)
+        geo_info, attrs, dtype = _get_crs_and_attribs(in_file=dataset)
         crs_list.append(geo_info.epsg)
         attrs_list.append(attrs)
 
@@ -130,7 +133,9 @@ def merge_datasets(
     # use rasterio merge to get a new raster dataset
     dst_kwds = {"QUALITY": "100", "REVERSIBLE": "YES"}
     try:
-        res = merge(datasets=datasets, dst_path=out_file, dst_kwds=dst_kwds, **kwargs)
+        res = merge(
+            datasets=datasets, dst_path=out_file, dst_kwds=dst_kwds,
+            dtype=dtype, **kwargs)
         if res is not None:
             out_ds, out_transform = res[0], res[1]
     except Exception as e:
@@ -188,16 +193,24 @@ def merge_datasets(
             else:
                 band_alias = f"B{idx+1}"
 
+        # get the mask from nodata values
+        nodata_mask = out_ds[idx, :, :] == nodata
+        band_arr = np.ma.masked_array(
+            data=out_ds[idx, :, :],
+            mask=nodata_mask,
+            fill_value=nodata
+        )
+
         raster.add_band(
             band_constructor=Band,
             band_name=band_name,
-            values=out_ds[idx, :, :],
+            values=band_arr,
             geo_info=geo_info,
             is_tiled=is_tiled,
             scale=scale,
             offset=offset,
             band_alias=band_alias,
-            unit=unit,
+            unit=unit
         )
 
     # clip raster collection if required to vector_features to keep consistency
diff --git a/eodal/core/band.py b/eodal/core/band.py
index af9874d8..cf8556c4 100644
--- a/eodal/core/band.py
+++ b/eodal/core/band.py
@@ -861,6 +861,13 @@ def from_rasterio(
             nodata, nodata_vals = None, attrs.get("nodatavals", None)
             if nodata_vals is not None:
                 nodata = nodata_vals[band_idx - 1]
+        # make sure the nodata type matches the datatype of the
+        # band values
+        if band_data.dtype.kind not in 'fc' and np.isnan(nodata):
+            raise TypeError(
+                f"The datatype of the band data is {band_data.dtype} " +
+                f"while the nodata value ({nodata}) is float.\n" +
+                "Please provide an appropriate nodata value")
 
         if masking:
             # make sure to set the EPSG code
@@ -1372,6 +1379,9 @@ def get_meta(self, driver: Optional[str] = "gTiff", **kwargs) -> Dict[str, Any]:
         meta["transform"] = self.transform
         meta["is_tile"] = self.is_tiled
         meta["driver"] = driver
+        # "compress" as suggested here:
+        # https://github.com/rasterio/rasterio/discussions/2933#discussioncomment-7208578
+        meta["compress"] = "DEFLATE"
         meta.update(kwargs)
 
         return meta
@@ -1699,10 +1709,18 @@ def mask(self, mask: np.ndarray, inplace: Optional[bool] = False):
                     # ignore pixels already masked
                     if not orig_mask[row, col]:
                         orig_mask[row, col] = mask[row, col]
+            # re-use original fill value
+            fill_value = self.values.fill_value
             # update band data array
-            masked_array = np.ma.MaskedArray(data=self.values.data, mask=orig_mask)
+            masked_array = np.ma.MaskedArray(
+                data=self.values.data,
+                mask=orig_mask,
+                fill_value=fill_value)
         elif self.is_ndarray:
-            masked_array = np.ma.MaskedArray(data=self.values, mask=mask)
+            # determine fill value from nodata value
+            fill_value = self.nodata
+            masked_array = np.ma.MaskedArray(
+                data=self.values, mask=mask, fill_value=fill_value)
         elif self.is_zarr:
             raise NotImplementedError()
 
@@ -1890,8 +1908,10 @@ def resample(
             out_mask = cv2.resize(in_mask, dim_resampled, cv2.INTER_NEAREST_EXACT)
             # convert mask back to boolean array
             out_mask = out_mask.astype(bool)
+            # re-use fill value of the original array
+            fill_value = self.values.fill_value
             # save as masked array
-            res = np.ma.masked_array(data=res, mask=out_mask)
+            res = np.ma.masked_array(data=res, mask=out_mask, fill_value=fill_value)
 
         # update the geo_info with new pixel resolution. The upper left x and y
         # coordinate must be changed if the pixel coordinates refer to the center
@@ -1962,6 +1982,7 @@ def reproject(
         if self.is_masked_array:
             band_data = deepcopy(self.values.data).astype(float)
             band_mask = deepcopy(self.values.mask).astype(float)
+            fill_value = self.values.fill_value
         elif self.is_ndarray:
             band_data = deepcopy(self.values).astype(float)
         elif self.is_zarr:
@@ -2004,7 +2025,8 @@ def reproject(
             # due to the raster alignment
             nodata = reprojection_options.get("src_nodata", 0)
             out_mask[out_data == nodata] = True
-            out_data = np.ma.MaskedArray(data=out_data, mask=out_mask)
+            out_data = np.ma.MaskedArray(
+                data=out_data, mask=out_mask, fill_value=fill_value)
 
         new_geo_info = GeoInfo.from_affine(affine=out_transform, epsg=target_crs)
         if inplace:
@@ -2192,23 +2214,24 @@ def scale_data(
         scale, offset = self.scale, self.offset
         if self.is_masked_array:
             if pixel_values_to_ignore is None:
-                scaled_array = scale * (self.values.data + offset)
+                scaled_array = scale * self.values.data - offset
             else:
                 scaled_array = self.values.data.copy().astype(float)
-                scaled_array[~np.isin(scaled_array, pixel_values_to_ignore)] = scale * (
-                    scaled_array[~np.isin(scaled_array, pixel_values_to_ignore)]
-                    + offset
-                )
-            scaled_array = np.ma.MaskedArray(data=scaled_array, mask=self.values.mask)
+                scaled_array[~np.isin(scaled_array, pixel_values_to_ignore)] = scale * \
+                    scaled_array[~np.isin(scaled_array, pixel_values_to_ignore)] \
+                    - offset
+            # reuse fill value
+            fill_value = self.values.fill_value
+            scaled_array = np.ma.MaskedArray(
+                data=scaled_array, mask=self.values.mask, fill_value=fill_value)
         elif self.is_ndarray:
             if pixel_values_to_ignore is None:
-                scaled_array = scale * (self.values + offset)
+                scaled_array = scale * self.values - offset
             else:
                 scaled_array = self.values.copy().astype(float)
-                scaled_array[~np.isin(scaled_array, pixel_values_to_ignore)] = scale * (
-                    scaled_array[~np.isin(scaled_array, pixel_values_to_ignore)]
-                    + offset
-                )
+                scaled_array[~np.isin(scaled_array, pixel_values_to_ignore)] = scale * \
+                    scaled_array[~np.isin(scaled_array, pixel_values_to_ignore)] \
+                    - offset
         elif self.is_zarr:
             raise NotImplementedError()
 
@@ -2356,6 +2379,9 @@ def to_rasterio(self, fpath_raster: Path, **kwargs) -> None:
         with rio.open(fpath_raster, "w+", **meta) as dst:
             # set band name
             dst.set_band_description(1, self.band_name)
+            # set scale and offset
+            dst._set_all_scales([self.scale])
+            dst._set_all_offsets([self.offset])
             # write band data
             if self.is_masked_array:
                 vals = self.values.data
diff --git a/eodal/core/raster.py b/eodal/core/raster.py
index de196338..a7b7023d 100644
--- a/eodal/core/raster.py
+++ b/eodal/core/raster.py
@@ -90,6 +90,7 @@ class allows thereby to handle ``Band`` instances with different spatial referen
 from numbers import Number
 from pathlib import Path
 from rasterio.drivers import driver_from_extension
+from rasterio.io import MemoryFile
 from typing import Any
 from typing import Callable
 from typing import Dict
@@ -101,6 +102,7 @@ class allows thereby to handle ``Band`` instances with different spatial referen
 from eodal.core.band import Band
 from eodal.core.operators import Operator
 from eodal.core.spectral_indices import SpectralIndices
+from eodal.core.utils import get_highest_dtype
 from eodal.utils.constants import ProcessingLevels
 from eodal.utils.decorators import check_band_names
 from eodal.utils.exceptions import BandNotFoundError
@@ -1675,6 +1677,7 @@ def to_rasterio(
         fpath_raster: Path,
         band_selection: Optional[List[str]] = None,
         use_band_aliases: Optional[bool] = False,
+        as_cog: Optional[bool] = False
     ) -> None:
         """
         Writes bands in collection to a raster dataset on disk using
@@ -1689,14 +1692,30 @@ def to_rasterio(
         :param use_band_aliases:
             use band aliases instead of band names for setting raster
             band descriptions to the output dataset
+        :param as_cog:
+            write the raster dataset as cloud-optimized GeoTIFF. This
+            requires the ``rio-cogeo`` package to be installed. Disabled
+            by default.
         """
-        # check output file naming and driver
-        try:
-            driver = driver_from_extension(fpath_raster)
-        except Exception as e:
-            raise ValueError(
-                f"Could not determine GDAL driver for " f"{fpath_raster.name}: {e}"
-            )
+        # check if COG output is enabled
+        if as_cog:
+            try:
+                from rio_cogeo.cogeo import cog_translate
+                from rio_cogeo.profiles import cog_profiles
+            except ModuleNotFoundError:
+                raise ModuleNotFoundError(
+                    "rio-cogeo is required for writing cloud-optimized GeoTIFFs\n" +
+                    "Install it with `pip install rio-cogeo`"
+                )
+            driver = "GTiff"
+        else:
+            # check output file naming and driver
+            try:
+                driver = driver_from_extension(fpath_raster)
+            except Exception as e:
+                raise ValueError(
+                    f"Could not determine GDAL driver for " f"{fpath_raster.name}: {e}"
+                )
 
         # check band_selection, if not provided use all available bands
         if band_selection is None:
@@ -1721,46 +1740,67 @@ def to_rasterio(
         # check meta and update it with the selected driver for writing the result
         meta = deepcopy(self[band_selection[0]].meta)
         dtypes = [self[x].values.dtype for x in band_selection]
-        if len(set(dtypes)) != 1:
-            UserWarning(
-                f"Multiple data types found in arrays to write ({set(dtypes)}). "
-                f"Casting to highest data type"
-            )
 
-        if len(set(dtypes)) == 1:
-            dtype_str = str(dtypes[0])
-        else:
-            # TODO: determine highest dtype
-            dtype_str = "float32"
+        # data type checking. We need to get the highest data type
+        highest_dtype = get_highest_dtype(dtypes)
 
         # update driver, the number of bands and the metadata value
         meta.update(
             {
                 "driver": driver,
                 "count": len(band_selection),
-                "dtype": dtype_str,
+                "dtype": str(highest_dtype),
                 "nodata": self[band_selection[0]].nodata,
+                "compression": "DEFLATE"
             }
         )
 
         # open the result dataset and try to write the bands
-        with rio.open(fpath_raster, "w+", **meta) as dst:
-            for idx, band_name in enumerate(band_selection):
-                # check with band name to set
-                dst.set_band_description(idx + 1, band_name)
-                # write band data
-                band_data = self.get_band(band_name).values.astype(dtype_str)
-                # set masked pixels to nodata
-                if self[band_name].is_masked_array:
-                    vals = band_data.data
-                    mask = band_data.mask
-                    vals[mask] = self[band_name].nodata
-                dst.write(band_data, idx + 1)
+        if as_cog:
+            with MemoryFile() as memfile:
+                with memfile.open(**meta) as mem:
+                    # set scales and offsets
+                    scales = [self[band_name].scale for band_name in band_selection]
+                    offsets = [self[band_name].offset for band_name in band_selection]
+                    mem._set_all_scales(scales)
+                    mem._set_all_offsets(offsets)
+                    # polulate data
+                    for idx, band_name in enumerate(band_selection):
+                        # check with band name to set
+                        mem.set_band_description(idx + 1, band_name)
+                        # write band data. Cast to highest data type if necessary.
+                        band_data = self.get_band(band_name).values.astype(highest_dtype)
+                        mem.write(band_data, idx + 1)
+
+                    # write the COG
+                    dst_profile = cog_profiles.get("deflate")
+                    cog_translate(
+                        mem,
+                        fpath_raster,
+                        dst_profile,
+                        in_memory=True,
+                        quiet=True
+                    )
+
+        else:
+            with rio.open(fpath_raster, "w+", **meta) as dst:
+                # set scales and offsets
+                scales = [self[band_name].scale for band_name in band_selection]
+                offsets = [self[band_name].offset for band_name in band_selection]
+                dst._set_all_scales(scales)
+                dst._set_all_offsets(offsets)
+                # polulate data
+                for idx, band_name in enumerate(band_selection):
+                    # check with band name to set
+                    dst.set_band_description(idx + 1, band_name)
+                    # write band data. Cast to highest data type if necessary.
+                    band_data = self.get_band(band_name).values.astype(highest_dtype)
+                    dst.write(band_data, idx + 1)
 
     @check_band_names
     def to_xarray(self, band_selection: Optional[List[str]] = None) -> xr.DataArray:
         """
-        Converts bands in collection a ``xarray.DataArray``
+        Converts bands in collection into a ``xarray.DataArray``
 
         :param band_selection:
             selection of bands to process. If not provided uses all
diff --git a/eodal/core/sensors/sentinel2.py b/eodal/core/sensors/sentinel2.py
index 542b1791..fdfe4974 100644
--- a/eodal/core/sensors/sentinel2.py
+++ b/eodal/core/sensors/sentinel2.py
@@ -112,10 +112,11 @@ def _get_gain_and_offset(in_dir: Union[str, Path]) -> Tuple[Number, Number]:
         baseline = get_S2_processing_baseline_from_safe(dot_safe_name=in_dir)
         # starting with baseline N0400 (400) S2 reflectances have an offset
         # value of -1000, i.e., the values reported in the .jp2 files must
-        # be subtracted by 1000 to obtain the actual reflectance factor values
+        # be subtracted by 1000 (0.1 in scaled representation) to obtain the
+        # actual reflectance factor values.
         s2_offset = 0
-        if baseline == 400:
-            s2_offset = -1000
+        if baseline >= 400:
+            s2_offset = 0.1
         return (s2_gain_factor, s2_offset)
 
     @staticmethod
@@ -293,7 +294,7 @@ def from_safe(
         masking_after_read_required = False
 
         if kwargs.get("vector_features") is not None:
-            bounds_df, shape_mask, lowest_resolution = adopt_vector_features_to_mask(
+            bounds_df, shape_mask, _ = adopt_vector_features_to_mask(
                 band_df=band_df_safe,
                 vector_features=kwargs.get("vector_features")
             )
@@ -382,9 +383,6 @@ def from_safe(
                 )
             # 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")
@@ -394,7 +392,11 @@ def from_safe(
                 )
                 # cast back to boolean
                 mask = res.astype("bool")
-                sentinel2.mask(mask=mask, bands_to_mask=[band_name], inplace=True)
+                sentinel2.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 = sentinel2.band_names
@@ -405,6 +407,10 @@ def from_safe(
                 band_selection=sel_bands,
                 pixel_values_to_ignore=[sentinel2[sentinel2.band_names[0]].nodata],
             )
+        # set SCL fill value to 0
+        if "SCL" in sentinel2.band_names:
+            if sentinel2["SCL"].is_masked_array:
+                sentinel2["SCL"].values.fill_value = 0
         return sentinel2
 
     @classmethod
diff --git a/eodal/core/utils/__init__.py b/eodal/core/utils/__init__.py
index e69de29b..f5d44c43 100644
--- a/eodal/core/utils/__init__.py
+++ b/eodal/core/utils/__init__.py
@@ -0,0 +1,38 @@
+
+import numpy as np
+
+
+# ranking of Python data types according to numpy
+DATATYPES = [
+    np.dtype('int'), np.dtype('int8'), np.dtype('int16'),
+    np.dtype('int32'),np.dtype('int64'), np.dtype('uint'),
+    np.dtype('uint8'), np.dtype('uint16'), np.dtype('uint32'),
+    np.dtype('uint64'), np.dtype('float16'), np.dtype('float32'),
+    np.dtype('float64'), np.dtype('complex64'),
+    np.dtype('complex128')
+]
+DTYPE_RANKS = {x: x.num for x in DATATYPES}
+
+
+def get_rank(dtype: np.dtype | str) -> int:
+    """
+    Get the rank of a data type
+
+    :param dtype: data type for which to get the rank.
+    :return: rank of `dtype` as defined by `numpy`.
+    """
+    if dtype in DTYPE_RANKS:
+        return DTYPE_RANKS[dtype]
+    else:
+        raise ValueError(f"Unknown data type: {dtype}")
+
+
+def get_highest_dtype(dtype_list: list[np.dtype | str]) -> np.dtype | str:
+    """
+    Get the highest data type of a list of data types based
+    on the data types ranks defined by `numpy`
+
+    :param dtype_list: list of data types
+    :return: higest data type in its `numpy.dtype` or string representation
+    """
+    return max(dtype_list, key=get_rank)
diff --git a/eodal/mapper/mapper.py b/eodal/mapper/mapper.py
index f78d2bbf..f3cd1bc0 100644
--- a/eodal/mapper/mapper.py
+++ b/eodal/mapper/mapper.py
@@ -393,8 +393,9 @@ def query_scenes(self) -> None:
 
         # make sure the time column is handled as pandas datetime
         # objects
-        scenes_df[self.time_column] = pd.to_datetime(
-            scenes_df[self.time_column], utc=True)
+        if not scenes_df.empty:
+            scenes_df[self.time_column] = pd.to_datetime(
+                scenes_df[self.time_column], utc=True)
 
         # populate the metadata attribute
         self.metadata = scenes_df
diff --git a/eodal/utils/geometry.py b/eodal/utils/geometry.py
index c57a5d1c..61e419fb 100644
--- a/eodal/utils/geometry.py
+++ b/eodal/utils/geometry.py
@@ -96,6 +96,9 @@ def prepare_gdf(
     :returns:
         resulting GeoDataFrame sorted by sensing time.
     """
+    # return an empty GeoDataFrame if there are no entries
+    if len(metadata_list) == 0:
+        return gpd.GeoDataFrame()
     df = pd.DataFrame(metadata_list)
     df.sort_values(by="sensing_time", inplace=True)
     return gpd.GeoDataFrame(df, geometry="geom", crs=df.epsg.unique()[0])
@@ -112,8 +115,6 @@ def adopt_vector_features_to_mask(
     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:
@@ -155,7 +156,7 @@ def adopt_vector_features_to_mask(
             dataset=src,
             shapes=vector_features_geom,
             all_touched=True,
-            crop=True,
+            crop=True
         )
     # get upper left coordinates rasterio takes for the band
     # with the coarsest spatial resolution
diff --git a/examples/sentinel2_mapper.py b/examples/sentinel2_mapper.py
index 4e115557..abd634e7 100644
--- a/examples/sentinel2_mapper.py
+++ b/examples/sentinel2_mapper.py
@@ -126,9 +126,9 @@ def preprocess_sentinel2_scenes(
     # 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()
+    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.
@@ -139,6 +139,7 @@ def preprocess_sentinel2_scenes(
         'scene_constructor': Sentinel2.from_safe,
         'scene_constructor_kwargs': {'band_selection':
                                      ['B02', 'B03', 'B04', 'B08'],
+                                     'apply_scaling': False,
                                      'read_scl': True},
         'scene_modifier': preprocess_sentinel2_scenes,
         'scene_modifier_kwargs': {'target_resolution': 10}
@@ -147,7 +148,14 @@ def preprocess_sentinel2_scenes(
     # 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
+    scoll = mapper.data
+
+    # save scenes as cloud-optimized GeoTiff
+    for timestamp, scene in scoll:
+        scene.to_rasterio(
+            f'data/{timestamp}.tiff',
+            band_selection=['red', 'green', 'blue', 'nir_1'],
+            as_cog=True)
 
     # plot scenes in collection
     f_scenes = mapper.data.plot(band_selection=['blue', 'green', 'red'])
diff --git a/setup.py b/setup.py
index 8150ad56..199bdd51 100755
--- a/setup.py
+++ b/setup.py
@@ -140,7 +140,7 @@ def run(self):
     name='eodal',
     # setup_requires=['setuptools_scm'],
     # use_scm_version={'version_scheme': 'python-simplified-semver'},
-    version='0.2.2',
+    version='0.2.3',
     description='The Earth Observation Data Analysis Library EOdal',
     long_description=long_description,
     long_description_content_type='text/markdown',
diff --git a/tests/conftest.py b/tests/conftest.py
index 520a4a59..5b5edb36 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -264,7 +264,8 @@ def _get_scene_collection():
         # open three scenes
         scene_list = []
         for i in range(3):
-            ds = RasterCollection.from_multi_band_raster(fpath_raster=fpath_raster)
+            ds = RasterCollection.from_multi_band_raster(
+                fpath_raster=fpath_raster, nodata=0)
             ds.scene_properties.acquisition_time = 1000 * (i+1)
             scene_list.append(ds)
         scoll = SceneCollection.from_raster_collections(
diff --git a/tests/core/test_band.py b/tests/core/test_band.py
index 96e4716d..1ab7580a 100644
--- a/tests/core/test_band.py
+++ b/tests/core/test_band.py
@@ -249,6 +249,19 @@ def test_masking(datadir, get_test_band, get_bandstack, get_points3):
         'band data not the same after writing'
 
 
+def test_wrong_nodata_type(get_bandstack):
+    # usage of an incorrect no-data type (nan) for INT array
+    fpath_raster = get_bandstack()
+    with pytest.raises(TypeError):
+        band = Band.from_rasterio(
+            fpath_raster=fpath_raster,
+            band_idx=1,
+            band_name_dst='B02',
+            full_bounding_box_only=True,
+            nodata=np.nan
+        )
+
+
 def test_read_pixels(get_bandstack, get_test_band, get_polygons, get_points3):
     # read single pixels from raster dataset
     fpath_raster = get_bandstack()
@@ -285,7 +298,8 @@ def test_read_pixels(get_bandstack, get_test_band, get_polygons, get_points3):
         band_idx=1,
         band_name_dst='B02',
         vector_features=vector_features,
-        full_bounding_box_only=True
+        full_bounding_box_only=True,
+        nodata=0
     )
 
     assert not band.is_masked_array, 'data should not be masked'
diff --git a/tests/core/test_raster_algebra.py b/tests/core/test_raster_algebra.py
index d74f2052..d29d3da5 100644
--- a/tests/core/test_raster_algebra.py
+++ b/tests/core/test_raster_algebra.py
@@ -7,7 +7,7 @@
 def test_raster_algebra_scalar(get_bandstack):
     """test algebraic operations using scalar values on RasterCollections"""
     fpath = get_bandstack()
-    rcoll = RasterCollection.from_multi_band_raster(fpath)
+    rcoll = RasterCollection.from_multi_band_raster(fpath, nodata=0)
     scalar = 2
 
     # arithmetic operators
@@ -78,7 +78,7 @@ def test_raster_algebra_scalar(get_bandstack):
 def test_raster_algebra_band_and_raster(get_bandstack):
     """test algebraic operations using Bands and Rasters on RasterCollections"""
     fpath = get_bandstack()
-    rcoll = RasterCollection.from_multi_band_raster(fpath)
+    rcoll = RasterCollection.from_multi_band_raster(fpath, nodata=0)
     band = rcoll['B02'].copy()
 
     # RasterCollection <-> Band
@@ -104,7 +104,7 @@ def test_raster_algebra_band_and_raster(get_bandstack):
     assert rcoll_le.get_values().any(), 'wrong results'
 
     # RasterCollection <-> RasterCollection
-    other = RasterCollection.from_multi_band_raster(fpath)
+    other = RasterCollection.from_multi_band_raster(fpath, nodata=0)
     rcoll_add = rcoll + other
     assert (rcoll_add.get_values() == rcoll.get_values() + other.get_values()).all(), 'wrong result'
     rcoll_sub = rcoll - other
diff --git a/tests/core/test_raster_apply.py b/tests/core/test_raster_apply.py
index c3a55a49..ab78f265 100644
--- a/tests/core/test_raster_apply.py
+++ b/tests/core/test_raster_apply.py
@@ -57,7 +57,7 @@ def test_apply_custom_function(get_bandstack):
     """
     fpath_raster = get_bandstack()
     gTiff_collection = RasterCollection.from_multi_band_raster(
-        fpath_raster=fpath_raster
+        fpath_raster=fpath_raster, nodata=0
     )
 
     # define a custom function for calculation the square root of values
diff --git a/tests/core/test_raster_collection.py b/tests/core/test_raster_collection.py
index 7ada4202..060615e8 100644
--- a/tests/core/test_raster_collection.py
+++ b/tests/core/test_raster_collection.py
@@ -47,7 +47,8 @@ def test_ndarray_constructor():
         band_name=band_name,
         values=values,
         band_alias=color_name,
-        geo_info=geo_info
+        geo_info=geo_info,
+        nodata=np.nan
     )
 
     assert len(handler) == 1, 'wrong number of bands in collection'
@@ -105,7 +106,12 @@ def test_rasterio_constructor_single_band(get_bandstack):
     # add a band from rasterio
     fpath_raster = get_bandstack()
     band_idx = 1
-    handler.add_band(Band.from_rasterio, fpath_raster=fpath_raster, band_idx=band_idx)
+    handler.add_band(
+        Band.from_rasterio,
+        fpath_raster=fpath_raster,
+        band_idx=band_idx,
+        nodata=0
+    )
 
     assert set(handler.band_names) == {'B1'}, \
         'band names not set properly in collection'
@@ -137,14 +143,15 @@ def test_scale():
         values=ones,
         geo_info=geo_info,
         scale=scale,
-        offset=offset
+        offset=offset,
+        nodata=0
     )
 
     before_scaling = handler.get_values()
     handler.scale(inplace=True)
     after_scaling = handler.get_values()
-    assert (after_scaling == (scale * (before_scaling + offset))).all(), 'values not scaled correctly'
-    assert (after_scaling == 300).all(), 'wrong value after scaling'
+    assert (after_scaling == scale * before_scaling - offset).all(), 'values not scaled correctly'
+    assert (after_scaling == 98).all(), 'wrong value after scaling'
 
 
 def test_rasterio_constructor_multi_band(get_bandstack):    
@@ -154,13 +161,13 @@ def test_rasterio_constructor_multi_band(get_bandstack):
     # read multi-band geoTiff into new handler
     fpath_raster = get_bandstack()
     gTiff_collection = RasterCollection.from_multi_band_raster(
-        fpath_raster=fpath_raster
+        fpath_raster=fpath_raster, nodata=0
     )
 
     assert gTiff_collection.band_names == \
         ['B02', 'B03', 'B04', 'B05', 'B06', 'B07', 'B08', 'B8A', 'B11', 'B12'], \
         'wrong list of band names in collection'
-    assert gTiff_collection.is_bandstack(), 'collection must be bandstacked'
+    assert gTiff_collection.is_bandstack(), 'collection must be band-stacked'
     gTiff_collection['B02'].crs == 32632, 'wrong CRS'
 
     # read multi-band geoTiff into new handler with custom destination names
@@ -168,7 +175,8 @@ def test_rasterio_constructor_multi_band(get_bandstack):
               'nir_1', 'nir_2', 'swir_1', 'swir_2']
     gTiff_collection = RasterCollection.from_multi_band_raster(
         fpath_raster=fpath_raster,
-        band_names_dst=colors
+        band_names_dst=colors,
+        nodata=0
     )
     assert gTiff_collection.band_names == colors, 'band names not set properly'
     gTiff_collection.calc_si('NDVI', inplace=True)
@@ -189,7 +197,8 @@ def test_rasterio_constructor_multi_band(get_bandstack):
     # with band aliases
     gTiff_collection = RasterCollection.from_multi_band_raster(
         fpath_raster=fpath_raster,
-        band_aliases=colors
+        band_aliases=colors,
+        nodata=0
     )
     assert gTiff_collection.has_band_aliases, 'band aliases must exist'
 
@@ -202,7 +211,8 @@ def test_reprojection(get_bandstack):
     fpath_raster = get_bandstack()
     gTiff_collection = RasterCollection.from_multi_band_raster(
         fpath_raster=fpath_raster,
-        band_names_dst=colors
+        band_names_dst=colors,
+        nodata=0
     )
     # try reprojection of raster bands to geographic coordinates
     reprojected = gTiff_collection.reproject(target_crs=4326)
@@ -227,7 +237,8 @@ def test_resampling(datadir,get_bandstack):
     # resample all bands to 5m
     gTiff_collection = RasterCollection.from_multi_band_raster(
         fpath_raster=fpath_raster,
-        band_aliases=colors
+        band_aliases=colors,
+        nodata=0
     )
     resampled = gTiff_collection.resample(
         target_resolution=5
@@ -240,11 +251,13 @@ def test_resampling(datadir,get_bandstack):
     resampled.to_rasterio(fpath_out)
     assert fpath_out.exists(), 'output-file not created'
 
+
 def test_clipping(get_bandstack):
     """Spatial clipping (subsetting) of RasterCollections"""
     fpath_raster = get_bandstack()
     rcoll = RasterCollection.from_multi_band_raster(
-        fpath_raster=fpath_raster
+        fpath_raster=fpath_raster,
+        nodata=0
     )
     # clip the collection to a Polygon buffered 100m inwards of the
     # bounding box of the first band
@@ -262,11 +275,13 @@ def test_clipping(get_bandstack):
     with pytest.raises(ValueError):
         rcoll.clip_bands(clipping_bounds=clipping_bounds, inplace=True)
 
+
 def test_band_summaries(get_bandstack, get_polygons):
     """test band summary statistics"""
     fpath_raster = get_bandstack()
     rcoll = RasterCollection.from_multi_band_raster(
-        fpath_raster=fpath_raster
+        fpath_raster=fpath_raster,
+        nodata=0
     )
     # try band summary statistics for polygons
     polys = get_polygons()
diff --git a/tests/core/test_raster_copy.py b/tests/core/test_raster_copy.py
index ce134f46..6e5f546e 100644
--- a/tests/core/test_raster_copy.py
+++ b/tests/core/test_raster_copy.py
@@ -13,7 +13,8 @@ def test_raster_copy(get_bandstack):
     fpath_raster = get_bandstack()
     rcoll = RasterCollection.from_multi_band_raster(
         fpath_raster=fpath_raster,
-        band_aliases=['a','b','c','d','e','f','g','h','i','j']
+        band_aliases=['a','b','c','d','e','f','g','h','i','j'],
+        nodata=0
     )
     rcoll_copy = rcoll.copy()
 
diff --git a/tests/core/test_raster_iterator.py b/tests/core/test_raster_iterator.py
index 25ce191c..1af7e646 100644
--- a/tests/core/test_raster_iterator.py
+++ b/tests/core/test_raster_iterator.py
@@ -15,7 +15,7 @@ def test_raster_iterator(get_bandstack):
 
     fpath_raster = get_bandstack()
     ds = RasterCollection.from_multi_band_raster(
-        fpath_raster=fpath_raster
+        fpath_raster=fpath_raster, nodata=0
     )
     band_names = ds.band_names
 
diff --git a/tests/core/test_raster_slicing.py b/tests/core/test_raster_slicing.py
index 9615263f..f6a45240 100644
--- a/tests/core/test_raster_slicing.py
+++ b/tests/core/test_raster_slicing.py
@@ -16,7 +16,8 @@ def test_raster_slice(get_bandstack):
                    'nir1', 'nir2', 'swir1', 'swir2']
     ds = RasterCollection.from_multi_band_raster(
         fpath_raster=fpath_raster,
-        band_aliases=color_names
+        band_aliases=color_names,
+        nodata=0
     )
 
     # slicing using band names
diff --git a/tests/core/test_scene_collection.py b/tests/core/test_scene_collection.py
index 4efcbda6..821a9935 100644
--- a/tests/core/test_scene_collection.py
+++ b/tests/core/test_scene_collection.py
@@ -51,7 +51,7 @@ def test_raster_is_scene(get_bandstack):
 
     fpath_raster = get_bandstack()
     ds = RasterCollection.from_multi_band_raster(
-        fpath_raster=fpath_raster
+        fpath_raster=fpath_raster, nodata=0
     )
     assert not ds.is_scene, 'scene metadata have not been set, so it is not a scene'
 
@@ -105,7 +105,8 @@ def test_scene_collection(get_s2_safe_l2a, get_polygons_2, get_bandstack):
     with pytest.raises(ValueError):
         scoll = SceneCollection(
             scene_constructor=RasterCollection.from_multi_band_raster,
-            fpath_raster=fpath_no_scene
+            fpath_raster=fpath_no_scene,
+            nodata=0
         )
 
     # add another scene