diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index a7ef147f..635956fd 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -3,12 +3,24 @@ Changelog
 #########
 
 
+2025.8.1 - 2025-08-12
+----------------------
+
+* core
+
+  * remove `blend_modes` (https://github.com/flrs/blend_modes) as external dependency as the package is not being maintained
+  * former `blend_modes` into `image_operations.blend_modes.blending_functions`, the original repository is under MIT LICENSE thus we will use it under those terms (https://github.com/flrs/blend_modes/blob/master/LICENSE.txt)
+  * added the `blend_modes` tests and images (`tests/testdata/blend_modes/{blend_mode}.png`)
+  * removed some unused mapchete test configs
+
+
 2025.8.0 - 2025-08-07
 ----------------------
 
 * add rudimentary example
 * add init docs, TODO: decide on CI for docs build and publish
 
+
 2025.7.0 - 2025-07-30
 ----------------------
 
diff --git a/mapchete_eo/__init__.py b/mapchete_eo/__init__.py
index 117f7b0d..17521c0a 100644
--- a/mapchete_eo/__init__.py
+++ b/mapchete_eo/__init__.py
@@ -1 +1 @@
-__version__ = "2025.8.0"
+__version__ = "2025.8.1"
diff --git a/mapchete_eo/image_operations/blend_modes/blending_functions.py b/mapchete_eo/image_operations/blend_modes/blending_functions.py
new file mode 100644
index 00000000..d4cbbb35
--- /dev/null
+++ b/mapchete_eo/image_operations/blend_modes/blending_functions.py
@@ -0,0 +1,198 @@
+"""
+
+Original LICENSE:
+
+MIT License
+
+Copyright (c) 2016 Florian Roscheck
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated
+documentation files (the "Software"), to deal in the Software without restriction, including without limitation the
+rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit
+persons to whom the Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the
+Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
+INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE
+OR OTHER DEALINGS IN THE SOFTWARE.
+
+Overview
+--------
+
+.. currentmodule:: blend_modes.blending_functions
+
+.. autosummary::
+    :nosignatures:
+
+    addition
+    darken_only
+    difference
+    divide
+    dodge
+    grain_extract
+    grain_merge
+    hard_light
+    lighten_only
+    multiply
+    normal
+    overlay
+    screen
+    soft_light
+    subtract
+"""
+
+import numpy as np
+from typing import Callable
+from mapchete_eo.image_operations.blend_modes.type_checks import (
+    assert_image_format,
+    assert_opacity,
+)
+
+
+class BlendBase:
+    def __init__(
+        self,
+        opacity=1.0,
+        disable_type_checks=False,
+        dtype=np.float16,
+        fcn_name="BlendBase",
+    ):
+        self.opacity = opacity
+        self.disable_type_checks = disable_type_checks
+        self.dtype = dtype
+        self.fcn_name = fcn_name
+
+    def _prepare(self, src: np.ndarray, dst: np.ndarray):
+        if not self.disable_type_checks:
+            assert_image_format(src, fcn_name=self.fcn_name, arg_name="src")
+            assert_image_format(dst, fcn_name=self.fcn_name, arg_name="dst")
+            assert_opacity(self.opacity, fcn_name=self.fcn_name)
+        if src.dtype != self.dtype:
+            src = src.astype(self.dtype)
+        if dst.dtype != self.dtype:
+            dst = dst.astype(self.dtype)
+        return src, dst
+
+    def blend(self, src: np.ndarray, dst: np.ndarray, blend_func: Callable):
+        src, dst = self._prepare(src, dst)
+        blended = blend_func(src, dst)
+        result = (blended * self.opacity) + (dst * (1 - self.opacity))
+        return np.clip(result, 0, 1).astype(self.dtype)
+
+
+def make_blend_function(blend_func: Callable):
+    # This function returns a wrapper that uses a shared BlendBase instance
+    base = BlendBase()
+
+    def func(
+        src: np.ndarray,
+        dst: np.ndarray,
+        opacity: float = 1.0,
+        disable_type_checks: bool = False,
+        dtype: np.dtype = np.float16,
+    ) -> np.ndarray:
+        # If parameters differ from base, create new BlendBase (rare)
+        if (
+            opacity != base.opacity
+            or disable_type_checks != base.disable_type_checks
+            or dtype != base.dtype
+        ):
+            base_local = BlendBase(opacity, disable_type_checks, dtype)
+            return base_local.blend(src, dst, blend_func)
+        return base.blend(src, dst, blend_func)
+
+    return func
+
+
+normal = make_blend_function(lambda s, d: s)
+multiply = make_blend_function(lambda s, d: s * d)
+screen = make_blend_function(lambda s, d: 1 - (1 - s) * (1 - d))
+darken_only = make_blend_function(lambda s, d: np.minimum(s, d))
+lighten_only = make_blend_function(lambda s, d: np.maximum(s, d))
+difference = make_blend_function(lambda s, d: np.abs(d - s))
+subtract = make_blend_function(lambda s, d: np.clip(d - s, 0, 1))
+
+
+def divide_blend(s: np.ndarray, d: np.ndarray) -> np.ndarray:
+    with np.errstate(divide="ignore", invalid="ignore"):
+        res = np.true_divide(d, s)
+        res[~np.isfinite(res)] = 0
+        return np.clip(res, 0, 1)
+
+
+divide = make_blend_function(divide_blend)
+
+
+def grain_extract_blend(s: np.ndarray, d: np.ndarray) -> np.ndarray:
+    return np.clip(d - s + 0.5, 0, 1)
+
+
+grain_extract = make_blend_function(grain_extract_blend)
+
+
+def grain_merge_blend(s: np.ndarray, d: np.ndarray) -> np.ndarray:
+    return np.clip(d + s - 0.5, 0, 1)
+
+
+grain_merge = make_blend_function(grain_merge_blend)
+
+
+def overlay_blend(s: np.ndarray, d: np.ndarray) -> np.ndarray:
+    mask = d <= 0.5
+    result = np.empty_like(d)
+    result[mask] = 2 * s[mask] * d[mask]
+    result[~mask] = 1 - 2 * (1 - s[~mask]) * (1 - d[~mask])
+    return np.clip(result, 0, 1)
+
+
+overlay = make_blend_function(overlay_blend)
+
+
+def hard_light_blend(s: np.ndarray, d: np.ndarray) -> np.ndarray:
+    mask = s <= 0.5
+    result = np.empty_like(d)
+    result[mask] = 2 * s[mask] * d[mask]
+    result[~mask] = 1 - 2 * (1 - s[~mask]) * (1 - d[~mask])
+    return np.clip(result, 0, 1)
+
+
+hard_light = make_blend_function(hard_light_blend)
+
+
+def soft_light_blend(s: np.ndarray, d: np.ndarray) -> np.ndarray:
+    result = (1 - 2 * s) * d**2 + 2 * s * d
+    return np.clip(result, 0, 1)
+
+
+soft_light = make_blend_function(soft_light_blend)
+
+
+def dodge_blend(s: np.ndarray, d: np.ndarray) -> np.ndarray:
+    with np.errstate(divide="ignore", invalid="ignore"):
+        res = np.true_divide(d, 1 - s)
+        res[~np.isfinite(res)] = 1
+        return np.clip(res, 0, 1)
+
+
+dodge = make_blend_function(dodge_blend)
+
+
+def burn_blend(s: np.ndarray, d: np.ndarray) -> np.ndarray:
+    with np.errstate(divide="ignore", invalid="ignore"):
+        res = 1 - np.true_divide(1 - d, s)
+        res[~np.isfinite(res)] = 0
+        return np.clip(res, 0, 1)
+
+
+burn = make_blend_function(burn_blend)
+
+
+def addition_blend(s: np.ndarray, d: np.ndarray) -> np.ndarray:
+    return np.clip(s + d, 0, 1)
+
+
+addition = make_blend_function(addition_blend)
diff --git a/mapchete_eo/image_operations/blend_modes/type_checks.py b/mapchete_eo/image_operations/blend_modes/type_checks.py
new file mode 100644
index 00000000..1b91e5bb
--- /dev/null
+++ b/mapchete_eo/image_operations/blend_modes/type_checks.py
@@ -0,0 +1,99 @@
+"""This module includes functions to check if variable types match expected formats."""
+
+import numpy as np
+
+
+def assert_image_format(image, fcn_name: str, arg_name: str, force_alpha: bool = True):
+    """Assert if image arguments have the expected format.
+
+    Checks:
+        - Image is a numpy array
+        - Array is of floating-point type
+        - Array is 3D (height x width x channels)
+        - Array has the correct number of layers (3 or 4)
+
+    Args:
+        image: The image to be checked.
+        fcn_name (str): Calling function name for display in error messages.
+        arg_name (str): Relevant argument name for display in error messages.
+        force_alpha (bool): Whether the image must include an alpha layer.
+
+    Raises:
+        TypeError: If type or shape are incorrect.
+    """
+    if not isinstance(image, np.ndarray):
+        raise TypeError(
+            f"\n[Invalid Type]\n"
+            f"Function: {fcn_name}\n"
+            f"Argument: {arg_name}\n"
+            f"Expected: numpy.ndarray\n"
+            f"Got: {type(image).__name__}\n"
+            f'Hint: Pass a numpy.ndarray for "{arg_name}".'
+        )
+
+    if image.dtype.kind != "f":
+        raise TypeError(
+            f"\n[Invalid Data Type]\n"
+            f"Function: {fcn_name}\n"
+            f"Argument: {arg_name}\n"
+            f'Expected dtype kind: "f" (floating-point)\n'
+            f'Got dtype kind: "{image.dtype.kind}"\n'
+            f"Hint: Convert the array to float, e.g., image.astype(float)."
+        )
+
+    if len(image.shape) != 3:
+        raise TypeError(
+            f"\n[Invalid Dimensions]\n"
+            f"Function: {fcn_name}\n"
+            f"Argument: {arg_name}\n"
+            f"Expected: 3D array (height x width x channels)\n"
+            f"Got: {len(image.shape)}D array\n"
+            f"Hint: Ensure the array has three dimensions."
+        )
+
+    if force_alpha and image.shape[2] != 4:
+        raise TypeError(
+            f"\n[Invalid Channel Count]\n"
+            f"Function: {fcn_name}\n"
+            f"Argument: {arg_name}\n"
+            f"Expected: 4 layers (R, G, B, Alpha)\n"
+            f"Got: {image.shape[2]} layers\n"
+            f"Hint: Include all four channels if force_alpha=True."
+        )
+
+
+def assert_opacity(opacity, fcn_name: str, arg_name: str = "opacity"):
+    """Assert if opacity has the expected format.
+
+    Checks:
+        - Opacity is float or int
+        - Opacity is within 0.0 <= x <= 1.0
+
+    Args:
+        opacity: The opacity value to be checked.
+        fcn_name (str): Calling function name for display in error messages.
+        arg_name (str): Argument name for display in error messages.
+
+    Raises:
+        TypeError: If type is not float or int.
+        ValueError: If opacity is out of range.
+    """
+    if not isinstance(opacity, (float, int)):
+        raise TypeError(
+            f"\n[Invalid Type]\n"
+            f"Function: {fcn_name}\n"
+            f"Argument: {arg_name}\n"
+            f"Expected: float (or int)\n"
+            f"Got: {type(opacity).__name__}\n"
+            f"Hint: Pass a float between 0.0 and 1.0."
+        )
+
+    if not 0.0 <= opacity <= 1.0:
+        raise ValueError(
+            f"\n[Out of Range]\n"
+            f"Function: {fcn_name}\n"
+            f"Argument: {arg_name}\n"
+            f"Expected: value in range 0.0 <= x <= 1.0\n"
+            f"Got: {opacity}\n"
+            f"Hint: Clamp or normalize the value to the valid range."
+        )
diff --git a/mapchete_eo/image_operations/compositing.py b/mapchete_eo/image_operations/compositing.py
index 8e7d080b..0abe1f26 100644
--- a/mapchete_eo/image_operations/compositing.py
+++ b/mapchete_eo/image_operations/compositing.py
@@ -2,13 +2,15 @@
 from enum import Enum
 from typing import Callable, Optional
 
-import blend_modes
 import cv2
 import numpy as np
 import numpy.ma as ma
 from mapchete import Timer
 from rasterio.plot import reshape_as_image, reshape_as_raster
 
+from mapchete_eo.image_operations.blend_modes import blending_functions
+
+
 logger = logging.getLogger(__name__)
 
 
@@ -68,63 +70,63 @@ def _blend_base(
 
 
 def normal(bg: np.ndarray, fg: np.ndarray, opacity: float = 1) -> ma.MaskedArray:
-    return _blend_base(bg, fg, opacity, blend_modes.normal)
+    return _blend_base(bg, fg, opacity, blending_functions.normal)
 
 
 def soft_light(bg: np.ndarray, fg: np.ndarray, opacity: float = 1) -> ma.MaskedArray:
-    return _blend_base(bg, fg, opacity, blend_modes.soft_light)
+    return _blend_base(bg, fg, opacity, blending_functions.soft_light)
 
 
 def lighten_only(bg: np.ndarray, fg: np.ndarray, opacity: float = 1) -> ma.MaskedArray:
-    return _blend_base(bg, fg, opacity, blend_modes.lighten_only)
+    return _blend_base(bg, fg, opacity, blending_functions.lighten_only)
 
 
 def screen(bg: np.ndarray, fg: np.ndarray, opacity: float = 1) -> ma.MaskedArray:
-    return _blend_base(bg, fg, opacity, blend_modes.screen)
+    return _blend_base(bg, fg, opacity, blending_functions.screen)
 
 
 def dodge(bg: np.ndarray, fg: np.ndarray, opacity: float = 1) -> ma.MaskedArray:
-    return _blend_base(bg, fg, opacity, blend_modes.dodge)
+    return _blend_base(bg, fg, opacity, blending_functions.dodge)
 
 
 def addition(bg: np.ndarray, fg: np.ndarray, opacity: float = 1) -> ma.MaskedArray:
-    return _blend_base(bg, fg, opacity, blend_modes.addition)
+    return _blend_base(bg, fg, opacity, blending_functions.addition)
 
 
 def darken_only(bg: np.ndarray, fg: np.ndarray, opacity: float = 1) -> ma.MaskedArray:
-    return _blend_base(bg, fg, opacity, blend_modes.darken_only)
+    return _blend_base(bg, fg, opacity, blending_functions.darken_only)
 
 
 def multiply(bg: np.ndarray, fg: np.ndarray, opacity: float = 1) -> ma.MaskedArray:
-    return _blend_base(bg, fg, opacity, blend_modes.multiply)
+    return _blend_base(bg, fg, opacity, blending_functions.multiply)
 
 
 def hard_light(bg: np.ndarray, fg: np.ndarray, opacity: float = 1) -> ma.MaskedArray:
-    return _blend_base(bg, fg, opacity, blend_modes.hard_light)
+    return _blend_base(bg, fg, opacity, blending_functions.hard_light)
 
 
 def difference(bg: np.ndarray, fg: np.ndarray, opacity: float = 1) -> ma.MaskedArray:
-    return _blend_base(bg, fg, opacity, blend_modes.difference)
+    return _blend_base(bg, fg, opacity, blending_functions.difference)
 
 
 def subtract(bg: np.ndarray, fg: np.ndarray, opacity: float = 1) -> ma.MaskedArray:
-    return _blend_base(bg, fg, opacity, blend_modes.subtract)
+    return _blend_base(bg, fg, opacity, blending_functions.subtract)
 
 
 def grain_extract(bg: np.ndarray, fg: np.ndarray, opacity: float = 1) -> ma.MaskedArray:
-    return _blend_base(bg, fg, opacity, blend_modes.grain_extract)
+    return _blend_base(bg, fg, opacity, blending_functions.grain_extract)
 
 
 def grain_merge(bg: np.ndarray, fg: np.ndarray, opacity: float = 1) -> ma.MaskedArray:
-    return _blend_base(bg, fg, opacity, blend_modes.grain_merge)
+    return _blend_base(bg, fg, opacity, blending_functions.grain_merge)
 
 
 def divide(bg: np.ndarray, fg: np.ndarray, opacity: float = 1) -> ma.MaskedArray:
-    return _blend_base(bg, fg, opacity, blend_modes.divide)
+    return _blend_base(bg, fg, opacity, blending_functions.divide)
 
 
 def overlay(bg: np.ndarray, fg: np.ndarray, opacity: float = 1) -> ma.MaskedArray:
-    return _blend_base(bg, fg, opacity, blend_modes.overlay)
+    return _blend_base(bg, fg, opacity, blending_functions.overlay)
 
 
 METHODS = {
diff --git a/pyproject.toml b/pyproject.toml
index a24eefb6..6f525c1c 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -22,7 +22,6 @@ classifiers = [
     "Topic :: Scientific/Engineering :: GIS",
 ]
 dependencies = [
-    "blend-modes",
     "click",
     "croniter",
     "lxml",
diff --git a/tests/conftest.py b/tests/conftest.py
index aa95096a..3ad94dad 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -7,6 +7,7 @@
 from mapchete.path import MPath
 from mapchete.testing import ProcessFixture
 from mapchete.tile import BufferedTilePyramid
+from PIL import Image
 from pystac_client import Client
 from rasterio import Affine
 from shapely import wkt
@@ -177,33 +178,6 @@ def eoxcloudless_8bit_dtype_scale_mapchete(tmp_path, testdata_dir):
         yield example
 
 
-@pytest.fixture
-def eoxcloudless_sentinel2_color_correction_mapchete(tmp_path, testdata_dir):
-    with ProcessFixture(
-        testdata_dir / "eoxcloudless_sentinel2_color_correction.mapchete",
-        output_tempdir=tmp_path,
-    ) as example:
-        yield example
-
-
-@pytest.fixture
-def eoxcloudless_rgb_map_mapchete(tmp_path, testdata_dir):
-    with ProcessFixture(
-        testdata_dir / "eoxcloudless_sentinel2_rgb_map.mapchete",
-        output_tempdir=tmp_path,
-    ) as example:
-        yield example
-
-
-@pytest.fixture
-def eoxcloudless_mosaic_mapchete(tmp_path, testdata_dir):
-    with ProcessFixture(
-        testdata_dir / "eoxcloudless_mosaic.mapchete",
-        output_tempdir=tmp_path,
-    ) as example:
-        yield example
-
-
 @pytest.fixture
 def sentinel2_antimeridian_east_mapchete(tmp_path, testdata_dir):
     with ProcessFixture(
@@ -222,15 +196,6 @@ def sentinel2_antimeridian_west_mapchete(tmp_path, testdata_dir):
         yield example
 
 
-@pytest.fixture
-def eoxcloudless_mosaic_regions_merge_mapchete(tmp_path, testdata_dir):
-    with ProcessFixture(
-        testdata_dir / "eoxcloudless_mosaic_regions_merge.mapchete",
-        output_tempdir=tmp_path,
-    ) as example:
-        yield example
-
-
 @pytest.fixture
 def sentinel2_mapchete(tmp_path, testdata_dir):
     with ProcessFixture(
@@ -800,3 +765,52 @@ def set_cdse_test_env(monkeypatch, request):
         monkeypatch.delenv("AWS_REQUEST_PAYER", raising=False)
     else:
         pytest.fail("CDSE AWS credentials not found in environment")
+
+
+@pytest.fixture(scope="session")
+def testdata_blend_modes_dir():
+    return (
+        MPath(os.path.dirname(os.path.realpath(__file__))) / "testdata" / "blend_modes"
+    )
+
+
+@pytest.fixture(scope="session")
+def src_image(testdata_blend_modes_dir):
+    img_path = testdata_blend_modes_dir / "orig.png"
+    img = Image.open(img_path).convert("RGBA")
+    arr = np.array(img) / 255.0
+    return arr.astype(np.float16)
+
+
+@pytest.fixture(scope="session")
+def dst_images(testdata_blend_modes_dir):
+    blend_names = [
+        "normal",
+        "multiply",
+        "screen",
+        "darken_only",
+        "lighten_only",
+        "difference",
+        "subtract",
+        "divide",
+        "grain_extract",
+        "grain_merge",
+        "overlay",
+        "hard_light",
+        "soft_light",
+        "dodge",
+        "burn",
+        "addition",
+    ]
+
+    dst_imgs = {}
+    for fname in os.listdir(testdata_blend_modes_dir):
+        if fname.endswith(".png"):
+            key = fname[:-4]  # remove .png
+            # Use only if file matches blend_names exactly (ignore *_50p etc)
+            if key in blend_names:
+                img = Image.open(testdata_blend_modes_dir / fname).convert("RGBA")
+                arr = np.array(img) / 255.0
+                dst_imgs[key] = arr.astype(np.float16)
+
+    return dst_imgs
diff --git a/tests/image_operations/test_blend_modes.py b/tests/image_operations/test_blend_modes.py
new file mode 100644
index 00000000..ec76e4bd
--- /dev/null
+++ b/tests/image_operations/test_blend_modes.py
@@ -0,0 +1,302 @@
+import pytest
+import numpy as np
+import mapchete_eo.image_operations.blend_modes.blending_functions as bf
+from mapchete_eo.image_operations.blend_modes.blending_functions import BlendBase
+
+
+# List of blend modes (replace with actual blend function names)
+blend_names = [
+    "normal",
+    "multiply",
+    "screen",
+    "overlay",
+    "lighten_only",
+    "darken_only",
+    # add your modes here...
+]
+
+# Blend modes that allow opacity parameter
+blend_modes_with_opacity = {
+    "normal",
+    "multiply",
+    "screen",
+    "overlay",
+    # add other modes supporting opacity here
+}
+
+# Blend modes where result can equal dst_image at opacity=1.0 (no difference expected)
+blend_modes_allowing_equal = {"lighten_only", "some_other_modes_if_any"}
+
+
+def test_func_behavior_with_blendbase():
+    base = BlendBase()
+
+    def test_blend_func(s, d):
+        return s  # just src, opacity applied outside in BlendBase.blend()
+
+    src = np.ones((2, 2, 4), dtype=np.float16)
+    dst = np.zeros((2, 2, 4), dtype=np.float16)
+
+    def func(
+        src: np.ndarray,
+        dst: np.ndarray,
+        opacity: float = 1.0,
+        disable_type_checks: bool = False,
+        dtype: np.dtype = np.float16,
+    ) -> np.ndarray:
+        if (
+            opacity != base.opacity
+            or disable_type_checks != base.disable_type_checks
+            or dtype != base.dtype
+        ):
+            base_local = BlendBase(opacity, disable_type_checks, dtype)
+
+            def local_blend_func(s, d):
+                return s  # raw src only
+
+            return base_local.blend(src, dst, local_blend_func)
+        return base.blend(src, dst, test_blend_func)
+
+    # opacity=1.0 => output == src
+    result_default = func(src, dst)
+    assert np.allclose(result_default, src)
+    assert result_default.dtype == np.float16
+    assert result_default.shape == src.shape
+
+    # opacity=0.3 => output == 0.3*src + 0.7*dst = 0.3 * 1 + 0.7 * 0 = 0.3
+    result_opacity = func(src, dst, opacity=0.3)
+    expected = 0.3 * src + 0.7 * dst
+    assert np.allclose(result_opacity, expected)
+    assert result_opacity.dtype == np.float16
+    assert result_opacity.shape == src.shape
+
+
+def test_make_blend_function_default_path():
+    # Create a dummy blend function, e.g. returns src
+    def dummy_blend(s, d):
+        return s
+
+    func = bf.make_blend_function(dummy_blend)
+
+    src = np.ones((2, 2, 4), dtype=np.float16)
+    dst = np.zeros((2, 2, 4), dtype=np.float16)
+
+    # Call with default params to hit the return base.blend() line
+    result = func(src, dst)
+
+    assert np.allclose(result, src)
+    assert result.dtype == np.float16
+    assert result.shape == src.shape
+
+
+def test_prepare_type_checks_enabled_casting():
+    blend = BlendBase()
+    blend.disable_type_checks = False
+    blend.fcn_name = "test_blend"
+    blend.dtype = np.float16
+    blend.opacity = 1.0
+
+    # Valid 3D float64 arrays with 4 channels -> should cast to float16
+    src = np.ones((2, 2, 4), dtype=np.float64)
+    dst = np.zeros((2, 2, 4), dtype=np.float64)
+
+    src_out, dst_out = blend._prepare(src, dst)
+
+    assert src_out.dtype == np.float16
+    assert dst_out.dtype == np.float16
+    assert src_out.shape == src.shape
+    assert dst_out.shape == dst.shape
+
+
+def test_prepare_type_checks_enabled_no_cast():
+    blend = BlendBase()
+    blend.disable_type_checks = False
+    blend.fcn_name = "test_blend"
+    blend.dtype = np.float16
+    blend.opacity = 1.0
+
+    # Correct dtype and shape: no cast, outputs share memory
+    src = np.ones((2, 2, 4), dtype=np.float16)
+    dst = np.zeros((2, 2, 4), dtype=np.float16)
+
+    src_out, dst_out = blend._prepare(src, dst)
+
+    assert src_out.dtype == np.float16
+    assert dst_out.dtype == np.float16
+    assert np.shares_memory(src, src_out)
+    assert np.shares_memory(dst, dst_out)
+
+
+def test_prepare_type_checks_disabled_casting():
+    blend = BlendBase()
+    blend.disable_type_checks = True  # disables type and shape checks
+    blend.fcn_name = "test_blend"
+    blend.dtype = np.float16
+    blend.opacity = 1.0
+
+    # Invalid shape (2D), but no error because checks disabled; dtype cast applied
+    src = np.ones((2, 2), dtype=np.float64)
+    dst = np.zeros((2, 2), dtype=np.float64)
+
+    src_out, dst_out = blend._prepare(src, dst)
+
+    assert src_out.dtype == np.float16
+    assert dst_out.dtype == np.float16
+    assert src_out.shape == src.shape
+    assert dst_out.shape == dst.shape
+
+
+def test_prepare_type_checks_enabled_invalid_shape_raises():
+    blend = BlendBase()
+    blend.disable_type_checks = False
+    blend.fcn_name = "test_blend"
+    blend.dtype = np.float16
+    blend.opacity = 1.0
+
+    # Invalid 2D shape arrays; should raise TypeError due to assert_image_format
+    src = np.ones((2, 2), dtype=np.float16)
+    dst = np.zeros((2, 2), dtype=np.float16)
+
+    with pytest.raises(TypeError, match="Expected: 3D array"):
+        blend._prepare(src, dst)
+
+
+def test_prepare_type_checks_enabled_invalid_channels_raises():
+    blend = BlendBase()
+    blend.disable_type_checks = False
+    blend.fcn_name = "test_blend"
+    blend.dtype = np.float16
+    blend.opacity = 1.0
+
+    # 3D shape but with 3 channels instead of 4, should raise on channel count
+    src = np.ones((2, 2, 3), dtype=np.float16)
+    dst = np.zeros((2, 2, 3), dtype=np.float16)
+
+    with pytest.raises(TypeError, match="Expected: 4 layers"):
+        blend._prepare(src, dst)
+
+
+@pytest.mark.parametrize("blend_name", blend_names)
+def test_blend_functions_opacity_zero(src_image, dst_images, blend_name):
+    if blend_name not in dst_images:
+        pytest.skip(f"No destination image for {blend_name}")
+
+    if blend_name not in blend_modes_with_opacity:
+        pytest.skip(f"Blend mode {blend_name} does not support opacity parameter")
+
+    dst_image = dst_images[blend_name]
+    blend_func = getattr(bf, blend_name)
+
+    opacity = 0.0
+    result = blend_func(
+        src_image,
+        dst_image,
+        opacity=opacity,
+        disable_type_checks=True,
+        dtype=np.float16,
+    )
+
+    assert result.shape == src_image.shape
+    assert result.dtype == np.float16, f"{blend_name} output dtype is not float16"
+    assert (result >= 0).all() and (result <= 1).all()
+
+    # Fully transparent blend: result should equal destination image
+    np.testing.assert_allclose(result, dst_image.astype(np.float16), rtol=1e-3)
+
+
+@pytest.mark.parametrize("blend_name", blend_names)
+def test_blend_functions_opacity_one(src_image, dst_images, blend_name):
+    if blend_name not in dst_images:
+        pytest.skip(f"No destination image for {blend_name}")
+
+    if blend_name not in blend_modes_with_opacity:
+        pytest.skip(f"Blend mode {blend_name} does not support opacity parameter")
+
+    dst_image = dst_images[blend_name]
+    blend_func = getattr(bf, blend_name)
+
+    opacity = 1.0
+    result = blend_func(
+        src_image,
+        dst_image,
+        opacity=opacity,
+        disable_type_checks=True,
+        dtype=np.float16,
+    )
+
+    assert result.shape == src_image.shape
+    assert result.dtype == np.float16, f"{blend_name} output dtype is not float16"
+    assert (result >= 0).all() and (result <= 1).all()
+
+    # Only assert difference if blend mode expected to differ from dst_image at full opacity
+    if blend_name not in blend_modes_allowing_equal:
+        assert not np.allclose(
+            result, dst_image.astype(np.float16)
+        ), f"Blend mode {blend_name} with opacity=1.0 result equals dst_image, which is unexpected"
+
+
+@pytest.mark.parametrize("blend_name", blend_names)
+@pytest.mark.parametrize("opacity", [0.25, 0.5, 0.75])
+def test_blend_functions_opacity_mid(src_image, dst_images, blend_name, opacity):
+    if blend_name not in dst_images:
+        pytest.skip(f"No destination image for {blend_name}")
+
+    if blend_name not in blend_modes_with_opacity:
+        pytest.skip(f"Blend mode {blend_name} does not support opacity parameter")
+
+    dst_image = dst_images[blend_name]
+    blend_func = getattr(bf, blend_name)
+
+    result = blend_func(
+        src_image,
+        dst_image,
+        opacity=opacity,
+        disable_type_checks=True,
+        dtype=np.float16,
+    )
+
+    assert result.shape == src_image.shape
+    assert result.dtype == np.float16, f"{blend_name} output dtype is not float16"
+    assert (result >= 0).all() and (result <= 1).all()
+
+
+def test_burn_blend_basic():
+    s = np.array([[0.5, 1.0], [0.2, 0.0]], dtype=np.float32)
+    d = np.array([[0.2, 0.3], [0.9, 0.5]], dtype=np.float32)
+    result = bf.burn_blend(s, d)
+
+    # output shape and dtype checks
+    assert result.shape == s.shape
+    assert result.dtype == s.dtype
+
+    # output range check
+    assert (result >= 0).all() and (result <= 1).all()
+
+    # no NaNs or Infs in output
+    assert np.isfinite(result).all()
+
+    # Check known values manually:
+    # When s=0.5, d=0.2 => 1 - (1 - 0.2)/0.5 = 1 - 0.8/0.5 = 1 - 1.6 = -0.6 clipped to 0
+    assert result[0, 0] == 0
+
+    # When s=1.0, d=0.3 => 1 - (1 - 0.3)/1.0 = 1 - 0.7 = 0.3
+    assert np.isclose(result[0, 1], 0.3)
+
+    # When s=0 (division by zero), output should be 0, no NaNs
+    assert result[1, 1] == 0
+
+
+def test_burn_blend_all_ones():
+    s = np.ones((3, 3), dtype=np.float32)
+    d = np.ones((3, 3), dtype=np.float32)
+    result = bf.burn_blend(s, d)
+    # burn_blend(1,1) == 1 - (1-1)/1 = 1
+    assert np.allclose(result, 1)
+
+
+def test_burn_blend_zero_s_and_d():
+    s = np.zeros((2, 2), dtype=np.float32)
+    d = np.zeros((2, 2), dtype=np.float32)
+    result = bf.burn_blend(s, d)
+    # division by zero case, all values set to 0 safely
+    assert np.all(result == 0)
diff --git a/tests/image_operations/test_blend_modes_type_checks.py b/tests/image_operations/test_blend_modes_type_checks.py
new file mode 100644
index 00000000..da2821d4
--- /dev/null
+++ b/tests/image_operations/test_blend_modes_type_checks.py
@@ -0,0 +1,65 @@
+import numpy as np
+import pytest
+
+from mapchete_eo.image_operations.blend_modes.type_checks import (
+    assert_image_format,
+    assert_opacity,
+)  # replace with actual module name
+
+
+def test_valid_image_with_alpha():
+    img = np.zeros((10, 10, 4), dtype=float)
+    assert_image_format(img, "blend_func", "image")  # Should not raise
+
+
+def test_valid_image_without_alpha_when_not_forced():
+    img = np.zeros((10, 10, 3), dtype=float)
+    assert_image_format(
+        img, "blend_func", "image", force_alpha=False
+    )  # Should not raise
+
+
+def test_image_not_numpy_array():
+    with pytest.raises(TypeError, match=r"\[Invalid Type\]"):
+        assert_image_format("not an array", "blend_func", "image")
+
+
+def test_image_wrong_dtype():
+    img = np.zeros((10, 10, 4), dtype=np.uint8)
+    with pytest.raises(TypeError, match=r"\[Invalid Data Type\]"):
+        assert_image_format(img, "blend_func", "image")
+
+
+def test_image_wrong_dimensions():
+    img = np.zeros((10, 10), dtype=float)  # 2D
+    with pytest.raises(TypeError, match=r"\[Invalid Dimensions\]"):
+        assert_image_format(img, "blend_func", "image")
+
+
+def test_image_wrong_channel_count_with_force_alpha():
+    img = np.zeros((10, 10, 3), dtype=float)  # No alpha channel
+    with pytest.raises(TypeError, match=r"\[Invalid Channel Count\]"):
+        assert_image_format(img, "blend_func", "image", force_alpha=True)
+
+
+def test_valid_opacity_float():
+    assert_opacity(0.5, "blend_func")  # Should not raise
+
+
+def test_valid_opacity_int():
+    assert_opacity(1, "blend_func")  # Should not raise
+
+
+def test_opacity_wrong_type():
+    with pytest.raises(TypeError, match=r"\[Invalid Type\]"):
+        assert_opacity("not a number", "blend_func")
+
+
+def test_opacity_below_range():
+    with pytest.raises(ValueError, match=r"\[Out of Range\]"):
+        assert_opacity(-0.1, "blend_func")
+
+
+def test_opacity_above_range():
+    with pytest.raises(ValueError, match=r"\[Out of Range\]"):
+        assert_opacity(1.1, "blend_func")
diff --git a/tests/testdata/blend_modes/addition.png b/tests/testdata/blend_modes/addition.png
new file mode 100644
index 00000000..3e12e706
Binary files /dev/null and b/tests/testdata/blend_modes/addition.png differ
diff --git a/tests/testdata/blend_modes/darken_only.png b/tests/testdata/blend_modes/darken_only.png
new file mode 100644
index 00000000..b375c032
Binary files /dev/null and b/tests/testdata/blend_modes/darken_only.png differ
diff --git a/tests/testdata/blend_modes/difference.png b/tests/testdata/blend_modes/difference.png
new file mode 100644
index 00000000..0f18c30a
Binary files /dev/null and b/tests/testdata/blend_modes/difference.png differ
diff --git a/tests/testdata/blend_modes/divide.png b/tests/testdata/blend_modes/divide.png
new file mode 100644
index 00000000..87f14daf
Binary files /dev/null and b/tests/testdata/blend_modes/divide.png differ
diff --git a/tests/testdata/blend_modes/dodge.png b/tests/testdata/blend_modes/dodge.png
new file mode 100644
index 00000000..d793c162
Binary files /dev/null and b/tests/testdata/blend_modes/dodge.png differ
diff --git a/tests/testdata/blend_modes/grain_extract.png b/tests/testdata/blend_modes/grain_extract.png
new file mode 100644
index 00000000..451be7c3
Binary files /dev/null and b/tests/testdata/blend_modes/grain_extract.png differ
diff --git a/tests/testdata/blend_modes/grain_merge.png b/tests/testdata/blend_modes/grain_merge.png
new file mode 100644
index 00000000..e0919373
Binary files /dev/null and b/tests/testdata/blend_modes/grain_merge.png differ
diff --git a/tests/testdata/blend_modes/hard_light.png b/tests/testdata/blend_modes/hard_light.png
new file mode 100644
index 00000000..128975b9
Binary files /dev/null and b/tests/testdata/blend_modes/hard_light.png differ
diff --git a/tests/testdata/blend_modes/layer.png b/tests/testdata/blend_modes/layer.png
new file mode 100644
index 00000000..fe124863
Binary files /dev/null and b/tests/testdata/blend_modes/layer.png differ
diff --git a/tests/testdata/blend_modes/layer_50p.png b/tests/testdata/blend_modes/layer_50p.png
new file mode 100644
index 00000000..2ba54c56
Binary files /dev/null and b/tests/testdata/blend_modes/layer_50p.png differ
diff --git a/tests/testdata/blend_modes/lighten_only.png b/tests/testdata/blend_modes/lighten_only.png
new file mode 100644
index 00000000..22fbebf1
Binary files /dev/null and b/tests/testdata/blend_modes/lighten_only.png differ
diff --git a/tests/testdata/blend_modes/multiply.png b/tests/testdata/blend_modes/multiply.png
new file mode 100644
index 00000000..0f840218
Binary files /dev/null and b/tests/testdata/blend_modes/multiply.png differ
diff --git a/tests/testdata/blend_modes/normal_100p.png b/tests/testdata/blend_modes/normal_100p.png
new file mode 100644
index 00000000..98f3ea18
Binary files /dev/null and b/tests/testdata/blend_modes/normal_100p.png differ
diff --git a/tests/testdata/blend_modes/normal_50p.png b/tests/testdata/blend_modes/normal_50p.png
new file mode 100644
index 00000000..8ff8e34f
Binary files /dev/null and b/tests/testdata/blend_modes/normal_50p.png differ
diff --git a/tests/testdata/blend_modes/orig.png b/tests/testdata/blend_modes/orig.png
new file mode 100644
index 00000000..26b09410
Binary files /dev/null and b/tests/testdata/blend_modes/orig.png differ
diff --git a/tests/testdata/blend_modes/overlay.png b/tests/testdata/blend_modes/overlay.png
new file mode 100644
index 00000000..42835fca
Binary files /dev/null and b/tests/testdata/blend_modes/overlay.png differ
diff --git a/tests/testdata/blend_modes/screen.png b/tests/testdata/blend_modes/screen.png
new file mode 100644
index 00000000..957d1234
Binary files /dev/null and b/tests/testdata/blend_modes/screen.png differ
diff --git a/tests/testdata/blend_modes/soft_light.png b/tests/testdata/blend_modes/soft_light.png
new file mode 100644
index 00000000..a7129366
Binary files /dev/null and b/tests/testdata/blend_modes/soft_light.png differ
diff --git a/tests/testdata/blend_modes/soft_light_50p.png b/tests/testdata/blend_modes/soft_light_50p.png
new file mode 100644
index 00000000..12aef6ad
Binary files /dev/null and b/tests/testdata/blend_modes/soft_light_50p.png differ
diff --git a/tests/testdata/blend_modes/subtract.png b/tests/testdata/blend_modes/subtract.png
new file mode 100644
index 00000000..0a3dcb8a
Binary files /dev/null and b/tests/testdata/blend_modes/subtract.png differ
diff --git a/tests/testdata/eoxcloudless_mosaic.mapchete b/tests/testdata/eoxcloudless_mosaic.mapchete
deleted file mode 100644
index 850053c2..00000000
--- a/tests/testdata/eoxcloudless_mosaic.mapchete
+++ /dev/null
@@ -1,72 +0,0 @@
-process: mapchete_eo.processes.eoxcloudless_mosaic
-input:
-  sentinel2:
-    format: Sentinel-2
-    level: L2A
-    time:
-      start: 2023-08-10
-      end: 2023-08-13
-    cat_baseurl: sentinel2/full_products/catalog.json
-    # DEBUG: Nice weather time frame below 
-    # time:
-    #  start: 2023-07-10
-    #  end: 2023-07-15
-
-    # DEBUG: Cache for local here below, switch off local catalogue baseurl
-    # cache:
-    #  keep: true
-    #  path: sentinel2/full_products
-    #  product_path_generation_method: date_year_first
-    #  assets: ["red", "green", "blue", "nir"]
-    #  max_disk_usage: 100.0
-    #  check_cached_files_exist: true
-    #  brdf:
-    #    bands: ["red", "green", "blue", "nir"]
-output:
-  format: GTiff
-  bands: 3
-  path: eoxcloudless_tmp
-  dtype: uint16
-pyramid:
-  grid: geodetic
-# DEBUG: use zoom 13 with the fancy date range
-# zoom_levels: 13
-zoom_levels: 9
-bounds: [16, 46, 16.1, 46.1]
-
-
-process_parameters:
-  resampling: bilinear
-  target_height: 10
-  read_masks: true
-  mask_config:
-    # !! dicts with one entry cause a bug in mapchete
-    footprint: true
-    buffer: 16
-    cloud_probability_threshold: 20
-    cloud_probability_resolution: 20
-    scl_classes: ["vegetation"]
-    # SCL classes and its corresponding mapchete_eo names
-    # mask using one or more of the SCL classes
-        # nodata = 0
-        # saturated_or_defected = 1
-        # dark_area_pixels = 2
-        # cloud_shadows = 3
-        # vegetation = 4
-        # not_vegetated = 5
-        # water = 6
-        # unclassified = 7
-        # cloud_medium_probability = 8
-        # cloud_high_probability = 9
-        # thin_cirrus = 10
-        # snow = 11
-  custom_mask_config:
-    footprint: true
-    l1c_cloud_type: all
-    buffer: 0
-    cloud_probability_threshold: 60
-    cloud_probability_resolution: 60     
-  # DEBUG: RGB + NIR for fun
-  # assets: ["red", "green", "blue", "nir"]
-  assets: ["red", "green", "blue"]
-  from_brightness_extract_method: third_quartile
\ No newline at end of file
diff --git a/tests/testdata/eoxcloudless_mosaic_regions_merge.mapchete b/tests/testdata/eoxcloudless_mosaic_regions_merge.mapchete
deleted file mode 100644
index ca099946..00000000
--- a/tests/testdata/eoxcloudless_mosaic_regions_merge.mapchete
+++ /dev/null
@@ -1,56 +0,0 @@
-process: mapchete_eo.processes.eoxcloudless_mosaic_merge
-input:
-  sentinel2:
-    region1:
-      format: Sentinel-2
-      level: L2A
-      time:
-        start: 2023-08-10
-        end: 2023-08-13
-      cat_baseurl: sentinel2/full_products/catalog.json
-      # upper half of tile 13 2001 8910
-      area: "POLYGON ((15.79833984375 46.021728515625, 15.79833984375 46.03271484375, 15.7763671875 46.03271484375, 15.7763671875 46.021728515625, 15.79833984375 46.021728515625))"
-    region2:
-      format: Sentinel-2
-      level: L2A
-      time:
-        start: 2023-08-15
-        end: 2023-08-18
-      cat_baseurl: sentinel2/full_products/catalog.json
-      # lower half of tile 13 2001 8910
-      area: "POLYGON ((15.79833984375 46.0107421875, 15.79833984375 46.021728515625, 15.7763671875 46.021728515625, 15.7763671875 46.0107421875, 15.79833984375 46.0107421875))"
-output:
-  format: GTiff
-  bands: 3
-  path: eoxcloudless_tmp
-  dtype: uint16
-pyramid:
-  grid: geodetic
-  # this has to be set to match gradient_buffer
-  pixelbuffer: 10
-zoom_levels: 13
-# tile 13 2001 8910
-bounds: [15.7763671875, 46.0107421875, 15.79833984375, 46.03271484375]
-
-process_parameters:
-  # mosaic settings
-  resampling: bilinear
-  target_height: 10
-  read_masks: true
-  mask_config:
-    footprint: true
-    buffer: 16
-    cloud_probability_threshold: 20
-    cloud_probability_resolution: 20
-    scl_classes: ["vegetation"]
-  custom_mask_config:
-    footprint: true
-    l1c_cloud_type: all
-    buffer: 0
-    cloud_probability_threshold: 60
-    cloud_probability_resolution: 60     
-  assets: ["red", "green", "blue"]
-  from_brightness_extract_method: third_quartile
-  # merge settings:
-  gradient_buffer: 10
-  merge_method: footprint_gradient
\ No newline at end of file
diff --git a/tests/testdata/eoxcloudless_sentinel2_color_correction.mapchete b/tests/testdata/eoxcloudless_sentinel2_color_correction.mapchete
deleted file mode 100644
index 1720a94b..00000000
--- a/tests/testdata/eoxcloudless_sentinel2_color_correction.mapchete
+++ /dev/null
@@ -1,36 +0,0 @@
-process: mapchete_eo.processes.eoxcloudless_sentinel2_color_correction
-input:
-  mosaic: eoxcloudless/2022_rgbnir/
-  desert_mask: 14-3758-17876.geojson
-output:
-  format: GTiff
-  bands: 3
-  path: eoxcloudless/color_correction/
-  dtype: uint8
-pyramid:
-  grid: geodetic
-zoom_levels: 13
-# tmx bounds 13 1879 8938
-bounds: [16.3916015625, 48.69140625, 16.41357421875, 48.71337890625]
-
-process_parameters:
-  bands: [1, 2, 3, 4]
-  resampling: nearest
-  matching_method: min
-  matching_max_zoom: 13
-  matching_precision: 12
-  fallback_to_higher_zoom: false
-  out_dtype: uint8
-  rgb_composite:
-    red: [0, 2400]
-    green: [5, 2400]
-    blue: [0, 2400]
-    gamma: 1.15
-    saturation: 1.6
-    clahe_clip_limit: 3.5
-    clahe_flag: true
-    sigmoidal_flag: true
-    sigmoidal_contrast: 2
-    sigmoidal_bias: 0.3
-    sharpen: true
-    smooth_water: true
\ No newline at end of file
diff --git a/tests/testdata/eoxcloudless_sentinel2_rgb_map.mapchete b/tests/testdata/eoxcloudless_sentinel2_rgb_map.mapchete
deleted file mode 100644
index b3d202c9..00000000
--- a/tests/testdata/eoxcloudless_sentinel2_rgb_map.mapchete
+++ /dev/null
@@ -1,43 +0,0 @@
-process: mapchete_eo.processes.eoxcloudless_rgb_map
-input:
-  mosaic: eoxcloudless/color_correction/
-  ocean_mask: eox-data/naturalearth/ne_10m_ocean.fgb
-  land_mask: eox-data/naturalearth/ne_10m_land.fgb
-  fuzzy_ocean_mask:
-    zoom<=9: eox-data/fuzzy_ocean_mask_v5/
-  ocean_depth:
-    zoom<=9: eox-data/gebco/
-  bathymetry:
-    zoom<=9: eox-data/bathymetry_v2/
-  mosaic_mask: eox-data/eoxcloudless_masks/antarctica.fgb
-
-output:
-  format: GTiff
-  bands: 3
-  path: eoxcloudless/rgb_map/
-  dtype: uint8
-pyramid:
-  grid: geodetic
-  pixelbuffer: 4
-zoom_levels: 6
-# tmx bounds 6 17 68
-bounds: [11.25, 39.375, 14.0625, 42.1875]
-
-process_parameters:
-  resampling: bilinear
-  matching_method: min
-  matching_max_zoom: 13
-  matching_precision: 12
-  fallback_to_higher_zoom: false
-
-  ocean_color: "#182c4e"
-  land_color: "#ffffff"
-  ocean_depth_opacity: 0.8
-  bathymetry_opacity: 0.8
-
-  fillnodata: true
-  fillnodata_method: nodata_neighbors  
-  fillnodata_max_patch_size: 9
-  fillnodata_max_nodata_neighbors: 9
-  fillnodata_max_search_distance: 0.25
-  fillnodata_smoothing_iterations: 1