diff --git a/.github/ci-environment.yml b/.github/ci-environment.yml
new file mode 100644
index 00000000..a859c5a0
--- /dev/null
+++ b/.github/ci-environment.yml
@@ -0,0 +1,5 @@
+name: ci
+channels:
+  - conda-forge
+dependencies:
+  - gdal
diff --git a/.github/workflows/install-methods.yml b/.github/workflows/install-methods.yml
new file mode 100644
index 00000000..ae3f720c
--- /dev/null
+++ b/.github/workflows/install-methods.yml
@@ -0,0 +1,57 @@
+name: Install Methods
+
+on:
+  schedule:
+    - cron: '0 6 * * 1'
+  push:
+    branches: [main]
+  pull_request:
+    branches: [main]
+  workflow_dispatch:
+
+jobs:
+  conda-pip:
+    name: conda + pip
+    runs-on: ubuntu-24.04
+    defaults:
+      run:
+        shell: bash -l {0}
+    steps:
+      - uses: actions/checkout@v6
+      - uses: conda-incubator/setup-miniconda@v4
+        with:
+          python-version: "3.12"
+          activate-environment: install-test
+          channels: conda-forge
+          channel-priority: strict
+          conda-remove-defaults: true
+      - run: conda install -y gdal
+      - run: pip install .[implementations]
+      - run: python -c "import openeo_processes_dask; print('OK')"
+
+  micromamba-pip:
+    name: micromamba + pip
+    runs-on: ubuntu-24.04
+    defaults:
+      run:
+        shell: bash -l {0}
+    steps:
+      - uses: actions/checkout@v6
+      - uses: mamba-org/setup-micromamba@v2
+        with:
+          environment-name: install-test
+          create-args: python=3.12 gdal
+      - run: pip install .[implementations]
+      - run: python -c "import openeo_processes_dask; print('OK')"
+
+  system-pip:
+    name: system packages + pip
+    runs-on: ubuntu-24.04
+    steps:
+      - uses: actions/checkout@v6
+      - run: sudo apt-get update && sudo apt-get install -y gdal-bin libgdal-dev
+      - uses: actions/setup-python@v6
+        with:
+          python-version: "3.12"
+      - run: pip install "gdal==$(gdal-config --version)" .[implementations]
+      - run: python -c "import openeo_processes_dask; print('OK')"
diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index cf96d3bb..10ef4b24 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -14,8 +14,7 @@ on:
   workflow_dispatch:
 
 env:
-  POETRY_VERSION: 1.5.1
-
+  POETRY_VERSION: 2.3.3
 
 jobs:
   tests:
@@ -24,77 +23,83 @@ jobs:
     strategy:
       matrix:
         os: [Ubuntu]
-        python-version: ["3.10", "3.11", "3.12"]
+        python-version: ["3.10", "3.11", "3.12", "3.13"]
         include:
           - os: Ubuntu
-            image: ubuntu-24.04   # switch to Noble (24.04)
+            image: ubuntu-24.04
       fail-fast: false
     defaults:
       run:
-        shell: bash
+        shell: bash -l {0}
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v6
         with:
           submodules: recursive
 
-      - name: Set up Python ${{ matrix.python-version }}
-        uses: actions/setup-python@v4
+      # --- Conda env with GDAL from conda-forge ---
+      - name: Set up Miniforge + Python ${{ matrix.python-version }}
+        uses: conda-incubator/setup-miniconda@v4
         with:
+          miniforge-version: latest
           python-version: ${{ matrix.python-version }}
+          activate-environment: ci
+          channels: conda-forge
+          channel-priority: strict
+          environment-file: .github/ci-environment.yml
+          use-mamba: true
+          conda-remove-defaults: true
+
+      # Cache downloaded conda tarballs (NOT the resolved environment)
+      - name: Cache conda package downloads
+        uses: actions/cache@v5
+        with:
+          path: ~/conda_pkgs_dir
+          key: conda-pkgs-${{ runner.os }}-py${{ matrix.python-version }}-${{ hashFiles('.github/ci-environment.yml') }}
+          restore-keys: |
+            conda-pkgs-${{ runner.os }}-py${{ matrix.python-version }}-
 
-      - name: Install system GDAL 3.11.3 (UbuntuGIS)
-        run: |
-          sudo apt-get update
-          sudo apt-get install -y software-properties-common
-          sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable -y
-          sudo apt-get update
-          # Install GDAL from the PPA (brings matching libgdal)
-          sudo apt-get install -y libgdal-dev gdal-bin
-          gdalinfo --version
-          # Show where libgdal comes from (should be /usr/lib/* from UbuntuGIS)
-          apt-cache policy libgdal-dev | sed -n '1,20p'
+      - name: Configure conda pkg cache directory
+        run: conda config --add pkgs_dirs ~/conda_pkgs_dir
 
-      - name: Get full Python version
-        id: full-python-version
-        run: echo version=$(python -c "import sys; print('-'.join(str(v) for v in sys.version_info))") >> $GITHUB_OUTPUT
+      # Always do a full solve + install — never skip, so dependency breakage is caught
+      - name: Install GDAL via conda-forge
+        run: |
+          mamba install -y gdal
+          conda clean --all -f -y
 
-      - name: Bootstrap poetry
-        run: curl -sSL https://install.python-poetry.org | python3 - -y --version $POETRY_VERSION
+      - name: Verify GDAL
+        run: |
+          gdalinfo --version
+          python -c "from osgeo import gdal; print('GDAL Python:', gdal.__version__)"
 
-      - name: Update PATH
-        run: echo "$HOME/.local/bin" >> $GITHUB_PATH
+      # --- Poetry ---
+      - name: Install Poetry
+        run: |
+          curl -sSL https://install.python-poetry.org | python3 - -y --version $POETRY_VERSION
+          echo "$HOME/.local/bin" >> $GITHUB_PATH
 
-      - name: Configure poetry
-        run: poetry config virtualenvs.in-project true
+      - name: Configure Poetry to use conda env
+        run: poetry config virtualenvs.create false
 
-      - name: Set up cache
-        uses: actions/cache@v3
-        id: cache
+      # Cache downloaded pip wheels (NOT installed packages)
+      - name: Cache pip downloads
+        uses: actions/cache@v5
         with:
-          path: .venv
-          key: venv-${{ runner.os }}-${{ steps.full-python-version.outputs.version }}-${{ hashFiles('**/poetry.lock') }}
+          path: ~/.cache/pip
+          key: pip-${{ runner.os }}-py${{ matrix.python-version }}-${{ hashFiles('**/poetry.lock') }}
+          restore-keys: |
+            pip-${{ runner.os }}-py${{ matrix.python-version }}-
 
-      - name: Ensure cache is healthy
-        if: steps.cache.outputs.cache-hit == 'true'
-        run: |
-          [ "$(command -v timeout)" ] || function timeout() { perl -e 'alarm shift; exec @ARGV' "$@"; }
-          timeout 10s poetry run pip --version || rm -rf .venv
-
-      # Important: make Poetry resolve GDAL against system libgdal=3.11
-      # If you explicitly pin GDAL in pyproject, set it to "==3.11.3".
       - name: Install dependencies
-        run: |
-          ver="$(gdal-config --version)"
-          poetry add "gdal==${ver}"          # adds to main group
-          poetry install --with dev --all-extras
+        run: poetry install --with dev --all-extras
 
       - name: Pre-commit hooks
-        run: poetry run pre-commit run --all-files
+        run: pre-commit run --all-files
 
       - name: Run pytest
-        run: poetry run pytest --cov=openeo_processes_dask --cov-report=xml
+        run: pytest --cov=openeo_processes_dask --cov-report=xml
 
       - name: Upload coverage to Codecov
-        uses: codecov/codecov-action@v3
+        uses: codecov/codecov-action@v6
         with:
           files: ./coverage.xml
diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml
index 05106e61..6d412bef 100644
--- a/.github/workflows/release.yml
+++ b/.github/workflows/release.yml
@@ -6,29 +6,53 @@ on:
       - '*.*.*'
 
 env:
-  POETRY_VERSION: 1.5.1
+  POETRY_VERSION: 2.3.3
 
 jobs:
   release:
     name: Release
-    runs-on: ubuntu-latest
+    runs-on: ubuntu-24.04
+    defaults:
+      run:
+        shell: bash -l {0}
     steps:
       - name: Checkout code
-        uses: actions/checkout@v3
+        uses: actions/checkout@v6
         with:
           submodules: recursive
 
-      - name: Set up Python 3.11
-        uses: actions/setup-python@v4
+      - name: Set up Miniforge + Python 3.12
+        uses: conda-incubator/setup-miniconda@v4
         with:
-          python-version: "3.11"
+          miniforge-version: latest
+          python-version: "3.12"
+          activate-environment: release
+          channels: conda-forge
+          channel-priority: strict
+          environment-file: .github/ci-environment.yml
+          use-mamba: true
+          conda-remove-defaults: true
 
-      - name: Bootstrap poetry
+      - name: Install GDAL via conda-forge
         run: |
-          curl -sL https://install.python-poetry.org | python - -y
+          mamba install -y gdal
+          conda clean --all -f -y
 
-      - name: Update PATH
-        run: echo "$HOME/.local/bin" >> $GITHUB_PATH
+      - name: Verify GDAL
+        run: |
+          gdalinfo --version
+          python -c "from osgeo import gdal; print('GDAL Python:', gdal.__version__)"
+
+      - name: Install Poetry
+        run: |
+          curl -sSL https://install.python-poetry.org | python3 - -y --version $POETRY_VERSION
+          echo "$HOME/.local/bin" >> $GITHUB_PATH
+
+      - name: Configure Poetry to use conda env
+        run: poetry config virtualenvs.create false
+
+      - name: Install dependencies
+        run: poetry install
 
       - name: Build project for distribution
         run: poetry build
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 52e68ec4..2e9e6532 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,38 +1,39 @@
 # See https://pre-commit.com for more information
 # See https://pre-commit.com/hooks.html for more hooks
 default_language_version:
-    python: python3
+  python: python3
+
 repos:
--   repo: https://github.com/asottile/pyupgrade
-    rev: v3.10.1
+  - repo: https://github.com/asottile/pyupgrade
+    rev: v3.21.2 # upgrading to 3.21.2 from 3.10.1 to accomodate py3.14 in test_matrix
     hooks:
       - id: pyupgrade
         args:
           - --keep-runtime-typing
           - --py39-plus
 
--   repo: https://github.com/python-poetry/poetry
-    rev: 1.5.1
+  - repo: https://github.com/python-poetry/poetry
+    rev: 2.3.3
     hooks:
       - id: poetry-check
 
--   repo: https://github.com/pycqa/isort
+  - repo: https://github.com/pycqa/isort
     rev: 5.12.0
     hooks:
       - id: isort
         name: isort (python)
         args: ["--profile", "black"]
 
--   repo: https://github.com/pre-commit/pre-commit-hooks
+  - repo: https://github.com/pre-commit/pre-commit-hooks
     rev: v4.4.0
     hooks:
-    -   id: trailing-whitespace
-    -   id: end-of-file-fixer
-    -   id: check-yaml
-    -   id: check-added-large-files
+      - id: trailing-whitespace
+      - id: end-of-file-fixer
+      - id: check-yaml
+      - id: check-added-large-files
 
--   repo: https://github.com/psf/black
+  - repo: https://github.com/psf/black
     rev: 23.7.0
     hooks:
-    -   id: black
+      - id: black
         exclude: ^specs/openeo-processes/
diff --git a/README.md b/README.md
index fbd20ad8..b9fe96d2 100644
--- a/README.md
+++ b/README.md
@@ -3,61 +3,46 @@
 [![Poetry](https://img.shields.io/endpoint?url=https://python-poetry.org/badge/v0.json)](https://python-poetry.org/)
 ![PyPI - Status](https://img.shields.io/pypi/status/openeo-processes-dask)
 ![PyPI](https://img.shields.io/pypi/v/openeo-processes-dask)
-![PyPI - Python Version](https://img.shields.io/pypi/pyversions/openeo-processes-dask)
-[![codecov](https://codecov.io/github/Open-EO/openeo-processes-dask/branch/main/graph/badge.svg?token=RA82MUN9RZ)](https://codecov.io/github/Open-EO/openeo-processes-dask)
+![Python Version](https://img.shields.io/badge/python-3.10%20|%203.11%20|%203.12%20|%203.13-blue)
+[![codecov](https://codecov.io/github/Eurac-Research-Institute-for-EO/openeo-processes-dask/branch/fix/dependencies_update/graph/badge.svg?token=RA82MUN9RZ)](https://codecov.io/github/Eurac-Research-Institute-for-EO/openeo-processes-dask)
 
 `openeo-processes-dask` is a collection of Python implementations of [OpenEO processes](https://processes.openeo.org/) based on the [xarray](https://github.com/pydata/xarray)/[dask](https://github.com/dask/dask) ecosystem. It is intended to be used alongside with [openeo-pg-parser-networkx](https://github.com/Open-EO/openeo-pg-parser-networkx), which handles the parsing and execution of [OpenEO process graphs](https://openeo.org/documentation/1.0/developers/api/reference.html#section/Processes/Process-Graphs). There you'll also find a tutorial on how to register process implementations from an arbitrary source (e.g. this repo) to the registry of available processes.
 
 ## Installation
 
-### Recommended installation (with system GDAL)
+### Conda (recommended — mirrors CI)
 
-If you already have GDAL installed on your system (or in a conda/micromamba env), always install the matching Python bindings first, then this package:
-
-#### Installation
-
-Install this project via pip:
+Conda-forge provides GDAL system libraries and Python bindings as a single coherent package. Once installed, pip does not need to touch GDAL at all.
 
 ```bash
-pip install "gdal==$(gdal-config --version)" openeo-processes-dask
+conda create -n openeo_processes_dask -c conda-forge python=3.12 gdal
+conda activate openeo_processes_dask
+pip install openeo-processes-dask[implementations]
 ```
 
-Note that by default this only installs the JSON process specs.
-In order to install the actual implementations, add the `implementations` extra:
-
+**Micromamba (lightweight alternative):**
 
 ```bash
-pip install "gdal==$(gdal-config --version)" openeo-processes-dask[implementations]
-````
-
-This ensures that the Python bindings link against your system GDAL libraries and avoids pip pulling a mismatched GDAL wheel.
+micromamba create -n openeo_processes_dask -c conda-forge python=3.12 gdal
+micromamba activate openeo_processes_dask
+pip install openeo-processes-dask[implementations]
+```
 
 ---
 
-### Installing GDAL if you don’t have it yet
+### System packages (Ubuntu/Debian)
 
-**System packages (Ubuntu/Debian):**
+If you already have GDAL system libraries, pin the pip `gdal` wheel to the matching version:
 
 ```bash
 sudo apt-get install gdal-bin libgdal-dev python3-gdal
 pip install "gdal==$(gdal-config --version)" openeo-processes-dask[implementations]
 ```
 
-**Conda (recommended for most users):**
-
-```bash
-conda create -n openeo_processes_dask -c conda-forge python=3.12 gdal
-conda activate openeo_processes_dask
-pip install openeo-processes-dask[implementations]
-```
-
-**Micromamba (lightweight alternative to conda):**
+> This pin avoids the common mismatch between a pip `gdal` wheel and your system `libgdal`. When using conda (above) this is unnecessary because conda-forge provides both the library and the Python binding together.
 
-```bash
-micromamba create -n openeo_processes_dask -c conda-forge python=3.12 gdal
-micromamba activate openeo_processes_dask
-pip install openeo-processes-dask[implementations]
-```
+Note that by default `pip install openeo-processes-dask` only installs the JSON process specs.
+In order to install the actual implementations, add the `implementations` extra as shown in the examples above.
 
 ---
 
@@ -70,15 +55,13 @@ A subset of process implementations with heavy or unstable dependencies are hidd
   ```bash
   pip install openeo-processes-dask[ml]
   ```
-* **Experimental processes:**
 
-  ```bash
-  pip install openeo-processes-dask[experimental]
-  ```
 
 ⚠️ **Note on GDAL:**
-Some extras (e.g. `implementations`, `ml`, `experimental`) may trigger installation of packages that depend on GDAL.
-To avoid version conflicts, make sure you have installed GDAL first (via conda/micromamba or system packages) and then install the extras as shown above.
+The `implementations` extra depends on GDAL transitively via `rasterio`, `rioxarray`, `odc-stac`, and `geopandas`.
+Always install GDAL **first** (via conda-forge or system packages) before pip-installing extras.
+The `ml` (`xgboost`) and `deforestation` (`rqadeforestation`) extras do not directly depend on GDAL.
+This project requires **GDAL >=3.9** and is CI-tested against conda-forge GDAL on Python 3.10–3.13.
 
 ---
 
@@ -92,22 +75,24 @@ Clone the repository with `--recurse-submodules` to also fetch the process specs
 git clone --recurse-submodules git@github.com:Open-EO/openeo-processes-dask.git
 ```
 
-To setup the python venv and install this project into it run:
-
-```bash
-poetry install --all-extras
-```
-
-⚠️ **Note on GDAL for development:**
-When using `poetry install --all-extras`, Poetry will attempt to install GDAL via pip, which may pull the latest GDAL wheels and cause conflicts with system libraries.
-It is strongly recommended to create a conda/micromamba environment with GDAL preinstalled before running `poetry install`. For example:
+**Development setup (CI pattern — conda-forge GDAL + Poetry):**
 
 ```bash
+# 1. Create conda env with GDAL (mirrors `.github/ci-environment.yml`)
 conda create -n openeo_processes_dask_dev -c conda-forge python=3.12 gdal
 conda activate openeo_processes_dask_dev
+
+# 2. Install Poetry deps into the conda env
+poetry config virtualenvs.create false
 poetry install --all-extras
+
+# 3. Verify GDAL
+gdalinfo --version
+python -c "from osgeo import gdal; print('GDAL Python:', gdal.__version__)"
 ```
 
+> `poetry config virtualenvs.create false` ensures GDAL from conda-forge is visible to pip-installed geospatial packages (`rasterio`, `rioxarray`, etc.). Without this, Poetry's isolated venv may not find the conda-forge GDAL libraries.
+
 ---
 
 To add a new core dependency run:
diff --git a/openeo_processes_dask/process_implementations/arrays.py b/openeo_processes_dask/process_implementations/arrays.py
index 5e5e0bbd..19fa8905 100644
--- a/openeo_processes_dask/process_implementations/arrays.py
+++ b/openeo_processes_dask/process_implementations/arrays.py
@@ -296,28 +296,35 @@ def array_find(
     if reverse:
         data = np.flip(data, axis=axis)
 
-    idxs = (data == value).argmax(axis=axis)
+    eq = data == value
+    idxs = eq.argmax(axis=axis)
+    mask = ~np.array(eq.any(axis=axis))
+
+    # np.isnan is only defined for numeric-like values; non-numeric search
+    # values are not NaN and should keep the normal not-found mask.
+    try:
+        if bool(np.isnan(value)):
+            mask = True
+    except (TypeError, ValueError):
+        pass
 
-    mask = ~np.array((data == value).any(axis=axis))
-    if not isinstance(value, str) and np.isnan(value):
-        mask = True
     if reverse:
-        if axis is None:
-            size = data.size
-        else:
-            size = data.shape[axis]
+        size = data.size if axis is None else data.shape[axis]
         idxs = size - 1 - idxs
 
-    logger.warning(
-        "array_find: numpy has no sentinel value for missing data in integer arrays, therefore np.masked_array is used to return the indices of found elements. Further operations might fail if not defined for masked arrays."
-    )
     if isinstance(idxs, da.Array):
         idxs = idxs.compute_chunk_sizes()
-        masked_idxs = np.atleast_1d(da.ma.masked_array(idxs, mask=mask))
-    else:
-        masked_idxs = np.atleast_1d(np.ma.masked_array(idxs, mask=mask))
+        masked_idxs = da.ma.masked_array(idxs, mask=mask)
+        filled_idxs = da.ma.filled(masked_idxs)
+        return da.atleast_1d(filled_idxs)
 
-    return masked_idxs
+    masked_idxs = np.ma.masked_array(idxs, mask=mask)
+    filled_idxs = np.ma.filled(masked_idxs)
+    filled_idxs = np.atleast_1d(filled_idxs)
+
+    if axis is None and filled_idxs.size == 1:
+        return filled_idxs[0]
+    return filled_idxs
 
 
 def array_find_label(data: ArrayLike, label: Union[str, int, float], dim_labels=None):
@@ -421,15 +428,15 @@ def array_interpolate_linear(data: ArrayLike, axis=None, dim_labels=None):
             x = np.arange(len(data))
 
     def interp(data):
+        data = np.asarray(data)
         valid = np.isfinite(data)
-        if (valid == 1).all():
+        if valid.all():
             return data
         if len(x[valid]) < 2:
             return data
         data[~valid] = np.interp(
             x[~valid], x[valid], data[valid], left=np.nan, right=np.nan
         )
-
         return data
 
     if axis:
@@ -438,12 +445,17 @@ def interp(data):
                 raise Exception(
                     f"Cannot load data of shape: {data.shape} into memory. "
                 )
-            # currently, there seems to be no way around loading the values,
-            # apply_along_axis cannot handle dask arrays
+            # np.apply_along_axis cannot handle dask arrays
             data = data.compute()
         data = np.apply_along_axis(interp, axis=axis, arr=data)
     else:
-        data = interp(data)
+        if isinstance(data, da.Array):
+            # Interpolation needs the full 1D context. Rechunk to a single block
+            # to preserve the eager implementation's result while staying lazy.
+            data = data.rechunk(data.shape)
+            data = da.map_blocks(interp, data, dtype=data.dtype)
+        else:
+            data = interp(data)
     if return_label:
         return array_create_labeled(data=data, labels=dim_labels)
     return data
diff --git a/openeo_processes_dask/process_implementations/comparison.py b/openeo_processes_dask/process_implementations/comparison.py
index f5e688b1..61247223 100644
--- a/openeo_processes_dask/process_implementations/comparison.py
+++ b/openeo_processes_dask/process_implementations/comparison.py
@@ -26,7 +26,7 @@
 
 
 def is_infinite(x: ArrayLike):
-    if np.issubsctype(get_scalar_type(x), np.number):
+    if np.issubdtype(get_scalar_type(x), np.number):
         return np.isinf(x)
     else:
         return False
@@ -56,26 +56,31 @@ def eq(
     x_dtype = get_scalar_type(x)
     y_dtype = get_scalar_type(y)
 
-    if np.issubsctype(x_dtype, np.number) and np.issubsctype(y_dtype, np.number):
+    if np.issubdtype(x_dtype, np.number) and np.issubdtype(y_dtype, np.number):
         if delta:
             ar_eq = np.isclose(x, y, atol=delta)
         else:
             ar_eq = x == y
 
-    elif np.issubsctype(x_dtype, np.bool_) and np.issubsctype(y_dtype, np.bool_):
+    elif np.issubdtype(x_dtype, np.bool_) and np.issubdtype(y_dtype, np.bool_):
         ar_eq = x == y
 
-    elif np.issubsctype(x_dtype, np.flexible) and np.issubsctype(y_dtype, np.flexible):
+    elif np.issubdtype(x_dtype, np.flexible) and np.issubdtype(y_dtype, np.flexible):
         if not case_sensitive:
-            if np.issubsctype(get_scalar_type(x), np.character):
+            if np.issubdtype(get_scalar_type(x), np.character):
                 x = np.char.lower(x)
-            if np.issubsctype(get_scalar_type(y), np.character):
+            if np.issubdtype(get_scalar_type(y), np.character):
                 y = np.char.lower(y)
         ar_eq = x == y
 
-    elif np.issubsctype(x_dtype, np.flexible) and np.issubsctype(y_dtype, np.flexible):
+    elif np.issubdtype(x_dtype, np.flexible) and np.issubdtype(y_dtype, np.flexible):
         ar_eq = x == y
     else:
+        # Types are incompatible, return False (but preserve array shape)
+        if hasattr(x, "shape"):
+            return np.zeros_like(x, dtype=bool)
+        elif hasattr(y, "shape"):
+            return np.zeros_like(y, dtype=bool)
         return False
 
     null_mask = np.logical_and(notnull(x), notnull(y))
diff --git a/openeo_processes_dask/process_implementations/cubes/apply_neighborhood_intertwin.py b/openeo_processes_dask/process_implementations/cubes/apply_neighborhood_intertwin.py
index 324acfcc..4b227b9e 100644
--- a/openeo_processes_dask/process_implementations/cubes/apply_neighborhood_intertwin.py
+++ b/openeo_processes_dask/process_implementations/cubes/apply_neighborhood_intertwin.py
@@ -16,6 +16,9 @@ def apply_neighborhood_intertwin(
     overlap: Optional[dict[int]] = None,
     context: Optional[dict] = None,
 ) -> RasterCube:
+    for dim, dim_size in zip(data.dims, data.shape):
+        if dim_size == 0:
+            raise ValueError(f"Data has a zero-size dimension: {dim}")
     positional_parameters = {"data": 0}
     named_parameters = {"context": context}
 
diff --git a/openeo_processes_dask/process_implementations/cubes/load.py b/openeo_processes_dask/process_implementations/cubes/load.py
index 400328a9..19dd77c2 100644
--- a/openeo_processes_dask/process_implementations/cubes/load.py
+++ b/openeo_processes_dask/process_implementations/cubes/load.py
@@ -5,7 +5,7 @@
 from datetime import datetime
 from pathlib import PurePosixPath
 from typing import Any, Callable, Dict, List, Optional, Tuple, Union
-from urllib.parse import unquote, urljoin, urlparse
+from urllib.parse import unquote, urlparse
 
 import numpy as np
 import odc.stac
@@ -83,116 +83,6 @@ def _search_for_parent_catalog(url):
     return catalog_url, collection_id
 
 
-def _load_with_xcube_eopf(
-    data_id: str,
-    spatial_extent: Optional[BoundingBox] = None,
-    temporal_extent: Optional[TemporalInterval] = None,
-    bands: Optional[list[str]] = None,
-    resolution: Optional[float] = None,
-    projection: Optional[Union[int, str]] = None,
-) -> RasterCube:
-    """Load data using xcube-eopf package for EOPF STAC endpoints."""
-    try:
-        from xcube.core.store import new_data_store
-    except ImportError:
-        raise ImportError(
-            "xcube-eopf package is required for loading EOPF STAC data. "
-            "Please install it with: pip install xcube-eopf"
-        )
-
-    from pyproj import CRS
-
-    # Convert spatial extent to bbox
-    bbox = None
-    if spatial_extent is not None:
-        try:
-            crs_obj = (
-                CRS(spatial_extent.crs) if spatial_extent.crs else CRS("EPSG:4326")
-            )
-            epsg = crs_obj.to_epsg()
-
-            if epsg == 4326:
-                projection = "EPSG:4326"
-                resolution = 10 / 111320
-
-            elif epsg and ((32601 <= epsg <= 32660) or (32701 <= epsg <= 32760)):
-                # Any UTM zone (north/south)
-                projection = f"EPSG:{epsg}"
-                resolution = 10
-
-            else:
-                raise ValueError(
-                    f"Unsupported CRS: {crs_obj.to_string()} "
-                    "(only EPSG:4326 or UTM EPSG:32601–32660 / 32701–32760 are supported)"
-                )
-
-            spatial_extent_used = spatial_extent
-            bbox = [
-                spatial_extent_used.west,
-                spatial_extent_used.south,
-                spatial_extent_used.east,
-                spatial_extent_used.north,
-            ]
-
-        except Exception as e:
-            raise Exception(f"Unable to parse the provided spatial extent: {e}")
-
-    # Convert temporal extent to time_range
-    time_range = None
-    if temporal_extent is not None:
-        start_date = (
-            str(temporal_extent[0].to_numpy().astype("datetime64[D]"))
-            if temporal_extent[0] is not None
-            else None
-        )
-        end_date = (
-            str(temporal_extent[1].to_numpy().astype("datetime64[D]"))
-            if temporal_extent[1] is not None
-            else None
-        )
-        # print(start_date, end_date)
-        time_range = [start_date, end_date]
-
-    # Set CRS
-    crs = projection if projection else "EPSG:4326"
-
-    # Convert resolution from degrees to meters if needed
-    spatial_res = resolution if resolution else 10 / 111320
-
-    # Create store and open data
-    store = new_data_store("eopf-zarr")
-    ds = store.open_data(
-        data_id=data_id,
-        bbox=bbox,
-        time_range=time_range,
-        spatial_res=spatial_res,
-        crs=crs,
-        variables=bands,
-    )
-
-    # Convert to dataarray if it's a dataset
-    if isinstance(ds, xr.Dataset):
-        # Find the appropriate dimension name for bands
-        band_dim = None
-        for dim in ds.dims:
-            if dim.lower() in ["band", "bands", "variable", "variables"]:
-                band_dim = dim
-                break
-
-        if band_dim is None:
-            # If no band dimension found, try to create one
-            data_vars = list(ds.data_vars.keys())
-            if len(data_vars) > 1:
-                ds = ds.to_array(dim="bands")
-            else:
-                # Single variable, convert to dataarray
-                ds = ds[data_vars[0]]
-        else:
-            ds = ds.to_array(dim=band_dim)
-
-    return ds
-
-
 def load_stac(
     url: str,
     spatial_extent: Optional[BoundingBox] = None,
@@ -203,26 +93,6 @@ def load_stac(
     projection: Optional[Union[int, str]] = None,
     resampling: Optional[str] = None,
 ) -> RasterCube:
-    # Check if this is an EOPF STAC URL
-    if "stac.core.eopf.eodc.eu" in url:
-        logger.info(f"Detected EOPF STAC URL: {url}, using xcube-eopf backend")
-
-        # Extract data_id from URL
-        parsed_url = urlparse(url)
-        path_parts = PurePosixPath(unquote(parsed_url.path)).parts
-        data_id = path_parts[-1] if path_parts else "sentinel-2-l2a"  # default fallback
-
-        # print(properties)
-
-        # Use xcube-eopf for loading
-        return _load_with_xcube_eopf(
-            data_id=data_id,
-            spatial_extent=spatial_extent,
-            temporal_extent=temporal_extent,
-            bands=bands,
-        )
-
-    # Original implementation for non-EOPF STAC URLs
     stac_type = _validate_stac(url)
 
     # TODO: load_stac should have a parameter to enable scale and offset?
diff --git a/openeo_processes_dask/process_implementations/cubes/resample.py b/openeo_processes_dask/process_implementations/cubes/resample.py
index 1657d131..ffffd737 100644
--- a/openeo_processes_dask/process_implementations/cubes/resample.py
+++ b/openeo_processes_dask/process_implementations/cubes/resample.py
@@ -196,12 +196,9 @@ def resample_cube_temporal(data, target, dimension=None, valid_within=None):
     for d in target[dimension].values:
         difference = np.abs(d - data[dimension].values)
         nearest = np.argwhere(difference == np.min(difference))
-        # The rare case of ties is resolved by choosing the earlier timestamps. (index 0)
-        if np.shape(nearest) == (2, 1):
-            nearest = nearest[0]
-        if np.shape(nearest) == (1, 2):
-            nearest = nearest[:, 0]
-        index.append(int(nearest))
+        # Flatten to handle all shapes consistently
+        nearest = nearest.flatten()
+        index.append(int(nearest[0]))  # always take first (earliest timestamp)
     times_at_target_time = data[dimension].values[index]
     new_data = data.loc[{dimension: times_at_target_time}]
     filter_values = new_data[dimension].values
diff --git a/openeo_processes_dask/process_implementations/ml/random_forest.py b/openeo_processes_dask/process_implementations/ml/random_forest.py
index de248b64..2bcda62e 100644
--- a/openeo_processes_dask/process_implementations/ml/random_forest.py
+++ b/openeo_processes_dask/process_implementations/ml/random_forest.py
@@ -28,7 +28,12 @@ def fit_regr_random_forest(
     target_var: str = None,
     **kwargs,
 ) -> Booster:
-    import xgboost as xgb
+    try:
+        from xgboost import dask as dxgb
+    except ImportError:
+        raise ImportError(
+            "xgboost[dask] is required for fit_regr_random_forest."
+        ) from None
 
     def load_geometries(geometries):
         if isinstance(geometries, str):
@@ -114,8 +119,8 @@ def drop_col(df, keep_var):
     y = drop_col(target, target_var)
 
     client = dask.distributed.default_client()
-    dtrain = xgb.dask.DaskDMatrix(client, X, y)
-    output = xgb.dask.train(client, params, dtrain, num_boost_round=1)
+    dtrain = dxgb.DaskDMatrix(client, X, y)
+    output = dxgb.train(client, params, dtrain, num_boost_round=1)
 
     return output["booster"]
 
@@ -126,7 +131,12 @@ def predict_random_forest(
     axis: int = -1,
     context: dict = None,
 ) -> RasterCube:
-    import xgboost as xgb
+    try:
+        from xgboost import dask as dxgb
+    except ImportError:
+        raise ImportError(
+            "xgboost[dask] is required for predict_random_forest."
+        ) from None
 
     if not model:
         if isinstance(context, dict) and "model" in context:
@@ -142,7 +152,7 @@ def predict_random_forest(
 
     # Run prediction
     client = dask.distributed.default_client()
-    preds_flat = xgb.dask.inplace_predict(client, model, X)
+    preds_flat = dxgb.inplace_predict(client, model, X)
 
     output_shape = list(data.shape)
     output_shape[axis] = 1
diff --git a/openeo_processes_dask/process_implementations/utils.py b/openeo_processes_dask/process_implementations/utils.py
index d6070d1d..be763c0e 100644
--- a/openeo_processes_dask/process_implementations/utils.py
+++ b/openeo_processes_dask/process_implementations/utils.py
@@ -5,7 +5,7 @@
 
 def get_scalar_type(obj):
     if np.isscalar(obj):
-        return np.obj2sctype(type(obj))
+        return np.dtype(type(obj)).type
     if hasattr(obj, "dtype"):
-        return obj.dtype
+        return np.dtype(obj.dtype).type
     return np.object_
diff --git a/pyproject.toml b/pyproject.toml
index b9ac2482..33cde78e 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -16,6 +16,7 @@ classifiers = [
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
+    "Programming Language :: Python :: 3.13",
 ]
 
 packages = [
@@ -23,14 +24,15 @@ packages = [
 ]
 
 [tool.poetry.dependencies]
-python = ">=3.10,<3.13"
-geopandas = { version = ">=0.11.1,<1", optional = true }
+python = ">=3.10,<3.14"
+gdal = { version = ">=3.9", optional = true }
+geopandas = { version = ">=0.11.1", optional = true }
 pandas = { version = ">=2.0.0", optional = true }
-xarray = { version = ">=2022.11.0,<2025.08.01", optional = true }
-dask = {extras = ["array", "dataframe", "distributed"], version = ">=2023.4.0,<2025.2.0", optional = true}
+xarray = { version = ">=2022.11.0", optional = true }
+dask = {extras = ["array", "dataframe", "distributed"], version = ">=2023.4.0", optional = true}
 rasterio = { version = "^1.3.4", optional = true }
-dask-geopandas = { version = "0.4.3", optional = true }
-xgboost = { version = ">=1.5.1,<2.1.4", optional = true }
+dask-geopandas = { version = ">=0.4.3", optional = true }
+xgboost = { version = ">=1.5.1", extras = ["dask"], optional = true }
 rioxarray = { version = ">=0.12.0,<1", optional = true }
 openeo-pg-parser-networkx = { version = ">=2025.10", optional = true }
 rqadeforestation = { version = ">=0.1", optional = true }
@@ -43,12 +45,11 @@ scipy = "^1.11.3"
 xvec = { version = "0.2.0", optional = true }
 joblib = { version = ">=1.3.2", optional = true }
 geoparquet = "^0.0.3"
-pyarrow = "^15.0.2"
+pyarrow = ">=15.0.2"
 openeo = ">=0.36.0"
-numpy = { version = "<2.0.0", optional = false }
-pystac = { version = "<1.12.0", optional = false }
+numpy = { version = ">=1.26.3,<3", optional = false }
+pystac = { version = ">=1.8.0", optional = false }
 zarr = "<=2.18.7"
-xcube-eopf = ">=0.2.0"
 
 
 [tool.poetry.group.dev.dependencies]
diff --git a/tests/conftest.py b/tests/conftest.py
index 51c41583..618e0f5c 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -53,7 +53,7 @@ def bounding_box(
         "north": north,
         "crs": crs,
     }
-    return BoundingBox.parse_obj(spatial_extent)
+    return BoundingBox.model_validate(spatial_extent)
 
 
 @pytest.fixture
@@ -67,7 +67,7 @@ def bounding_box_small(
         "north": north,
         "crs": crs,
     }
-    return BoundingBox.parse_obj(spatial_extent)
+    return BoundingBox.model_validate(spatial_extent)
 
 
 @pytest.fixture
@@ -95,7 +95,7 @@ def polygon_geometry_small(
 
 @pytest.fixture
 def temporal_interval(interval=["2018-05-01", "2018-06-01"]) -> TemporalInterval:
-    return TemporalInterval.parse_obj(interval)
+    return TemporalInterval.model_validate(interval)
 
 
 @pytest.fixture
diff --git a/tests/test_aggregate.py b/tests/test_aggregate.py
index f0a7d155..905f1a12 100644
--- a/tests/test_aggregate.py
+++ b/tests/test_aggregate.py
@@ -125,7 +125,7 @@ def test_aggregate_temporal_period(
     input_cube = create_fake_rastercube(
         data=random_raster_data,
         spatial_extent=bounding_box,
-        temporal_extent=TemporalInterval.parse_obj(temporal_extent),
+        temporal_extent=TemporalInterval.model_validate(temporal_extent),
         bands=["B02", "B03", "B04", "B08"],
     )
 
diff --git a/tests/test_arrays.py b/tests/test_arrays.py
index 6ab55c09..305738fb 100644
--- a/tests/test_arrays.py
+++ b/tests/test_arrays.py
@@ -258,7 +258,7 @@ def test_array_contains_object_dtype():
         ([1, 0, 3, 0, 2], 0, 3, None, True),
         ([[1, 0, 3, 2], [5, 3, 6, 8]], 3, [999999, 1, 0, 999999], 0, True),
         ([[1, 0, 3, 2], [5, 3, 6, 8]], 3, [2, 1], 1, True),
-        (["A", "B", "C"], "b", 99999, None, False),
+        (["A", "B", "C"], "b", 999999, None, False),
     ],
 )
 def test_array_find(data, value, expected, axis, reverse):
@@ -340,6 +340,15 @@ def test_array_interpolate_linear(data, expected):
     )
 
 
+def test_array_interpolate_linear_dask_multichunk_uses_global_context():
+    data = da.from_array(np.array([0.0, np.nan, np.nan, 3.0]), chunks=(2,))
+
+    result = array_interpolate_linear(data)
+
+    assert isinstance(result, da.Array)
+    np.testing.assert_allclose(result.compute(), [0.0, 1.0, 2.0, 3.0])
+
+
 @pytest.mark.parametrize(
     "data, expected, labels",
     [
@@ -648,3 +657,28 @@ def test_count(temporal_interval, bounding_box, random_raster_data, process_regi
     )
     assert output_cube.dims == ("x", "y", "t")
     xr.testing.assert_equal(output_cube, xr.zeros_like(output_cube))
+
+
+def _masked_fill_scalar():
+    return np.ma.array([0], mask=[True]).filled()[0]
+
+
+def test_array_find_non_numeric_string_not_found_returns_scalar_fill_value():
+    data = np.array(["a", "b", "c"], dtype=object)
+    result = array_find(data=data, value="z")
+    assert np.isscalar(result)
+    assert result == _masked_fill_scalar()
+
+
+def test_array_find_nan_value_returns_scalar_fill_value():
+    data = np.array([1.0, np.nan, 3.0])
+    result = array_find(data=data, value=np.nan)
+    assert np.isscalar(result)
+    assert result == _masked_fill_scalar()
+
+
+def test_array_find_axis_none_not_found_returns_scalar_fill_value():
+    data = np.array([1, 2, 3])
+    result = array_find(data=data, value=99, axis=None)
+    assert np.isscalar(result)
+    assert result == _masked_fill_scalar()
diff --git a/tests/test_comparison.py b/tests/test_comparison.py
index d1ae9930..8bf6a339 100644
--- a/tests/test_comparison.py
+++ b/tests/test_comparison.py
@@ -12,6 +12,7 @@
 from openeo_processes_dask.process_implementations.comparison import *
 from openeo_processes_dask.process_implementations.cubes.apply import apply
 from openeo_processes_dask.process_implementations.cubes.reduce import reduce_dimension
+from openeo_processes_dask.process_implementations.utils import get_scalar_type
 from tests.general_checks import assert_numpy_equals_dask_numpy, general_output_checks
 from tests.mockdata import create_fake_rastercube
 
@@ -67,6 +68,20 @@ def test_is_inf(value, expected, is_dask):
         assert hasattr(output, "dask")
 
 
+@pytest.mark.parametrize(
+    "value, expected",
+    [
+        (1, np.int64),
+        ("test", np.str_),
+        (None, np.object_),
+        (np.array([1, 2]), np.int64),
+        (da.from_array(np.array([1, 2])), np.int64),
+    ],
+)
+def test_get_scalar_type(value, expected):
+    assert get_scalar_type(value) is expected
+
+
 @pytest.mark.parametrize(
     "value,expected",
     [
@@ -185,7 +200,7 @@ def test_eq_mask():
     data = np.array([[10, 10], [10, 0]])
     data = da.from_array(data)
     m = eq(data, 10)
-    assert (m == data / 10).all()
+    assert (m == data / 10).all().compute()  # add .compute()
 
 
 @pytest.mark.parametrize(
diff --git a/tests/test_filter.py b/tests/test_filter.py
index 9b29906d..0cc655fb 100644
--- a/tests/test_filter.py
+++ b/tests/test_filter.py
@@ -29,7 +29,7 @@ def test_filter_temporal(temporal_interval, bounding_box, random_raster_data):
         backend="dask",
     )
 
-    temporal_interval_part = TemporalInterval.parse_obj(
+    temporal_interval_part = TemporalInterval.model_validate(
         ["2018-05-15T00:00:00", "2018-06-01T00:00:00"]
     )
     output_cube = filter_temporal(data=input_cube, extent=temporal_interval_part)
@@ -57,7 +57,9 @@ def test_filter_temporal(temporal_interval, bounding_box, random_raster_data):
             extent=["2018-05-31T23:59:59", "2018-05-15T00:00:00"],
         )
 
-    temporal_interval_open = TemporalInterval.parse_obj([None, "2018-05-03T00:00:00"])
+    temporal_interval_open = TemporalInterval.model_validate(
+        [None, "2018-05-03T00:00:00"]
+    )
     output_cube = filter_temporal(data=input_cube, extent=temporal_interval_open)
 
     xr.testing.assert_equal(
diff --git a/tests/test_ml.py b/tests/test_ml.py
index 6e757a39..97ab307a 100644
--- a/tests/test_ml.py
+++ b/tests/test_ml.py
@@ -1,3 +1,4 @@
+import builtins
 from functools import partial
 
 import dask
@@ -15,6 +16,7 @@
     fit_curve,
     fit_regr_random_forest,
     predict_curve,
+    random_forest,
 )
 from tests.mockdata import create_fake_rastercube
 
@@ -49,6 +51,40 @@ def test_fit_regr_random_forest_inline_geojson(
     assert isinstance(model, xgb.core.Booster)
 
 
+@pytest.mark.parametrize(
+    "func, args, message",
+    [
+        (
+            random_forest.fit_regr_random_forest,
+            (None, None),
+            "xgboost[dask] is required for fit_regr_random_forest.",
+        ),
+        (
+            random_forest.predict_random_forest,
+            (None, None),
+            "xgboost[dask] is required for predict_random_forest.",
+        ),
+    ],
+)
+def test_random_forest_missing_xgboost_dask_suppresses_import_context(
+    monkeypatch, func, args, message
+):
+    real_import = builtins.__import__
+
+    def fake_import(name, globals=None, locals=None, fromlist=(), level=0):
+        if name == "xgboost" and "dask" in fromlist:
+            raise ImportError("missing xgboost.dask")
+        return real_import(name, globals, locals, fromlist, level)
+
+    monkeypatch.setattr(builtins, "__import__", fake_import)
+
+    with pytest.raises(ImportError) as exc:
+        func(*args)
+
+    assert str(exc.value) == message
+    assert exc.value.__suppress_context__ is True
+
+
 @pytest.mark.parametrize("size", [(6, 5, 4, 3)])
 @pytest.mark.parametrize("dtype", [np.float64])
 def test_curve_fitting(temporal_interval, bounding_box, random_raster_data):
diff --git a/tests/test_resample.py b/tests/test_resample.py
index 7af4e5be..0d5dc361 100644
--- a/tests/test_resample.py
+++ b/tests/test_resample.py
@@ -161,14 +161,14 @@ def test_aggregate_temporal_period(
     input_cube = create_fake_rastercube(
         data=random_raster_data,
         spatial_extent=bounding_box,
-        temporal_extent=TemporalInterval.parse_obj(temporal_extent_1),
+        temporal_extent=TemporalInterval.model_validate(temporal_extent_1),
         bands=["B02", "B03", "B04", "B08"],
     )
 
     target_cube = create_fake_rastercube(
         data=random_raster_data,
         spatial_extent=bounding_box,
-        temporal_extent=TemporalInterval.parse_obj(temporal_extent_2),
+        temporal_extent=TemporalInterval.model_validate(temporal_extent_2),
         bands=["B02", "B03", "B04", "B08"],
     )