Repository: ungarj/mapchete (language: python)

Issue:
enable caching of remote vector data
e.g. a GML or GeoPackage on S3 would have to be fully read per process tile. There should be an option to either download it locally into a temporary directory or create a remote tile directory


Discussion / hints:

Fix patch:
diff --git a/doc/source/process_input.rst b/doc/source/process_input.rst
index 444bb55c..c0594e07 100644
--- a/doc/source/process_input.rst
+++ b/doc/source/process_input.rst
@@ -30,24 +30,48 @@ a file path (either remote or local) is sufficient.
         dem: path/to/elevation_data.tif
         land_polygons: path/to/land_polygons.shp
 
-from_command_line
------------------
 
-If the process is designed for a single file, you can also use the
-``--input_file`` parameter of the :doc:`cli`. In this case, the ``input``
-parameter has to be set to ``from_command_line``.
+Processing performance is best if the input files support quick access to subsets of the data.
+For rasters, converting them to COG (Cloud Optimized GeoTIFF) can sometimes do wonders especially
+if the data sits on an object storage. For vector data a format which support spatial indexing (e.g.
+GeoPackage or FlatGeoBuf) can do similar tricks although this works just in a limited way if the data
+is on an object storage.
 
 **Example:**
 
 .. code-block:: yaml
 
-    input: from_command_line
+    input:
+        dem:
+            format: raster_file
+            path: path/to/elevation_data.tif
+            cache:
+                path: tmp/path/to/elevation_data.tif
+                format: COG
+        land_polygons:
+            format: vector_file
+            path: path/to/land_polygons.shp
+            cache:
+                path: tmp/path/to/land_polygons.fgb
+                format: FlatGeobuf
+                keep: false  # if set to false cached data will not be deleted after the process finishes
+
+In some cases it can be beneficial if the data is not cached on disk but directly in memory. This feature
+should be used with care and only for smaller datasets.
 
-Run the process afterwards like this:
+**Example:**
 
-.. code-block:: shell
+.. code-block:: yaml
 
-    mapchete execute my_process.mapchete --input_file path/to/elevation_data.tif
+    input:
+        dem:
+            format: raster_file
+            path: path/to/elevation_data.tif
+            cache: memory
+        land_polygons:
+            format: vector_file
+            path: path/to/land_polygons.shp
+            cache: memory
 
 It is also possible to define input data groups e.g. for extracted Sentinel-2
 granules, where bands are stored in separate files:
diff --git a/mapchete/__init__.py b/mapchete/__init__.py
index 92622cf3..8b3400c8 100644
--- a/mapchete/__init__.py
+++ b/mapchete/__init__.py
@@ -2,7 +2,8 @@
 
 from mapchete._core import open, Mapchete
 from mapchete._executor import Executor, FakeFuture, SkippedFuture
-from mapchete._processing import Job, MapcheteProcess, ProcessInfo
+from mapchete._processing import Job, ProcessInfo
+from mapchete._user_process import MapcheteProcess
 from mapchete.tile import count_tiles
 from mapchete._timer import Timer
 
diff --git a/mapchete/_core.py b/mapchete/_core.py
index f7e8fca6..049686ca 100644
--- a/mapchete/_core.py
+++ b/mapchete/_core.py
@@ -12,12 +12,14 @@
 from mapchete.formats import read_output_metadata
 from mapchete.io import fs_from_path, tiles_exist
 from mapchete._processing import (
+    compute,
     _run_on_single_tile,
     _run_area,
     _preprocess,
     ProcessInfo,
-    TileProcess,
+    task_batches,
 )
+from mapchete._tasks import TileTask
 from mapchete.tile import count_tiles
 from mapchete._timer import Timer
 from mapchete.validate import validate_tile, validate_zooms
@@ -199,6 +201,18 @@ def skip_tiles(self, tiles=None, tiles_batches=None):
                 for tile in tiles:
                     yield (tile, False)
 
+    def task_batches(self, zoom=None, tile=None, skip_output_check=False, **kwargs):
+        """
+        Generate task batches from preprocessing tasks and tile tasks.
+        """
+        yield from task_batches(
+            self, zoom=zoom, tile=tile, skip_output_check=skip_output_check, **kwargs
+        )
+
+    def compute(self, **kwargs):
+        """Compute preprocessing tasks and tile tasks in one go."""
+        yield from compute(self, **kwargs)
+
     def batch_preprocessor(
         self,
         dask_scheduler=None,
@@ -399,10 +413,10 @@ def batch_processor(
                 dask_scheduler=dask_scheduler,
                 process=self,
                 tile=self.config.process_pyramid.tile(*tuple(tile)),
-            )
+            ).result()
         # run area
         else:
-            for process_info in _run_area(
+            for future in _run_area(
                 process=self,
                 zoom_levels=list(_get_zoom_level(zoom, self)),
                 dask_scheduler=dask_scheduler,
@@ -414,7 +428,7 @@ def batch_processor(
                 skip_output_check=skip_output_check,
                 executor=executor,
             ):
-                yield process_info
+                yield future.result()
 
     def count_tasks(self, minzoom=None, maxzoom=None, init_zoom=0):
         """
@@ -490,7 +504,7 @@ def execute(self, process_tile, raise_nodata=False):
         self.batch_preprocess()
         try:
             return self.config.output.streamline_output(
-                TileProcess(tile=process_tile, config=self.config).execute()
+                TileTask(tile=process_tile, config=self.config).execute()
             )
         except MapcheteNodataTile:
             if raise_nodata:  # pragma: no cover
@@ -708,6 +722,9 @@ def __exit__(self, exc_type, exc_value, exc_traceback):
             self.current_processes = None
             self.process_lock = None
 
+    def __repr__(self):  # pragma: no cover
+        return f"Mapchete <process_name={self.process_name}>"
+
 
 def _get_zoom_level(zoom, process):
     """Determine zoom levels."""
diff --git a/mapchete/_executor.py b/mapchete/_executor.py
index cfeb4e94..0b9991d1 100644
--- a/mapchete/_executor.py
+++ b/mapchete/_executor.py
@@ -181,6 +181,7 @@ def _finished_future(self, future, result=None, _dask=False):
         fut_exception = future.exception(timeout=FUTURE_TIMEOUT)
         if fut_exception:  # pragma: no cover
             logger.error("exception caught in future %s", future)
+            logger.exception(fut_exception)
             raise fut_exception
         result = result or future.result(timeout=FUTURE_TIMEOUT)
         if isinstance(result, CancelledError):  # pragma: no cover
@@ -204,12 +205,12 @@ def __enter__(self):
 
     def __exit__(self, *args):
         """Exit context manager."""
-        logger.debug("closing executor %s...", self._executor)
+        logger.info("closing executor %s...", self._executor)
         try:
             self._executor.close()
         except Exception:
             self._executor.__exit__(*args)
-        logger.debug("closed executor %s", self._executor)
+        logger.info("closed executor %s", self._executor)
 
     def __repr__(self):  # pragma: no cover
         return f"<Executor ({self._executor_cls})>"
@@ -230,6 +231,7 @@ def __init__(
         from dask.distributed import as_completed, Client, LocalCluster
 
         self._executor_client = dask_client
+        self._local_cluster = None
         if self._executor_client:  # pragma: no cover
             logger.info("using existing dask client: %s", dask_client)
         else:
@@ -237,9 +239,9 @@ def __init__(
                 n_workers=max_workers or os.cpu_count(), threads_per_worker=1
             )
             self._executor_cls = Client
-            self._executor_kwargs = dict(
-                address=dask_scheduler or LocalCluster(**local_cluster_kwargs),
-            )
+            if dask_scheduler is None:
+                self._local_cluster = LocalCluster(**local_cluster_kwargs)
+            self._executor_kwargs = dict(address=dask_scheduler or self._local_cluster)
             logger.info(
                 "starting dask.distributed.Client with kwargs %s", self._executor_kwargs
             )
@@ -458,14 +460,17 @@ def _executor(self):
     def __exit__(self, *args):
         """Exit context manager."""
         if self._executor_client:  # pragma: no cover
-            logger.debug("client not closing as it was passed on as kwarg")
+            logger.info("client not closing as it was passed on as kwarg")
         else:
-            logger.debug("closing executor %s...", self._executor)
+            logger.info("closing executor %s...", self._executor)
             try:
                 self._executor.close()
             except Exception:  # pragma: no cover
                 self._executor.__exit__(*args)
-            logger.debug("closed executor %s", self._executor)
+            logger.info("closed executor %s", self._executor)
+        if self._local_cluster:
+            logger.info("closing %s", self._local_cluster)
+            self._local_cluster.close()
 
 
 class ConcurrentFuturesExecutor(_ExecutorBase):
@@ -615,7 +620,7 @@ def cancelled(self):  # pragma: no cover
 
     def __repr__(self):  # pragma: no cover
         """Return string representation."""
-        return f"FakeFuture(result={self._result}, exception={self._exception})"
+        return f"<FakeFuture: type: {type(self._result)}, exception: {type(self._exception)})"
 
 
 class SkippedFuture:
@@ -637,11 +642,15 @@ def cancelled(self):  # pragma: no cover
         """Nothing to cancel here."""
         return False
 
+    def __repr__(self):  # pragma: no cover
+        """Return string representation."""
+        return f"<SkippedFuture: type: {type(self._result)}, exception: {type(self._exception)})"
+
 
 class FinishedFuture:
     """Wrapper class to mimick future interface."""
 
-    def __init__(self, future, result=None):
+    def __init__(self, future=None, result=None):
         """Set attributes."""
         try:
             self._result, self._exception = (
@@ -668,4 +677,4 @@ def cancelled(self):  # pragma: no cover
 
     def __repr__(self):  # pragma: no cover
         """Return string representation."""
-        return f"FinishedFuture(result={self._result}, exception={self._exception})"
+        return f"<FinishedFuture: type: {type(self._result)}, exception: {type(self._exception)})"
diff --git a/mapchete/_processing.py b/mapchete/_processing.py
index 09c9b9ec..c085b492 100644
--- a/mapchete/_processing.py
+++ b/mapchete/_processing.py
@@ -1,27 +1,40 @@
 """Internal processing classes and functions."""
 
 from collections import namedtuple
-
+from contextlib import ExitStack
 from itertools import chain
 import logging
 import multiprocessing
 from traceback import format_exc
 from typing import Generator
 
-from mapchete.commons import clip as commons_clip
-from mapchete.commons import contours as commons_contours
-from mapchete.commons import hillshade as commons_hillshade
 from mapchete.config import get_process_func
-from mapchete._executor import Executor, SkippedFuture
-from mapchete.errors import MapcheteNodataTile, MapcheteProcessException
-from mapchete.io import raster
+from mapchete._executor import DaskExecutor, Executor, SkippedFuture, FinishedFuture
+from mapchete.errors import MapcheteNodataTile
+from mapchete._tasks import to_dask_collection, TileTaskBatch, TileTask, TaskBatch, Task
 from mapchete._timer import Timer
-from mapchete.validate import deprecated_kwargs
+from mapchete.validate import validate_zooms
 
 logger = logging.getLogger(__name__)
 
 
-ProcessInfo = namedtuple("ProcessInfo", "tile processed process_msg written write_msg")
+ProcessInfo = namedtuple(
+    "ProcessInfo",
+    "tile processed process_msg written write_msg data",
+    defaults=(None, None, None, None, None, None),
+)
+
+TileProcessInfo = namedtuple(
+    "TileProcessInfo",
+    "tile processed process_msg written write_msg data",
+    defaults=(None, None, None, None, None, None),
+)
+
+PreprocessingProcessInfo = namedtuple(
+    "PreprocessingProcessInfo",
+    "task_key processed process_msg written write_msg data",
+    defaults=(None, None, None, None, None, None),
+)
 
 
 class Job:
@@ -70,6 +83,22 @@ def _run(self):
             yield from self.func(*self.fargs, executor=self.executor, **self.fkwargs)
             self.status = "finished"
 
+    def set_executor_kwargs(self, executor_kwargs):
+        """
+        Overwrite default or previously set Executor creation kwargs.
+
+        This only has an effect if Job was initialized with 'as_iterator' and has not yet run.
+        """
+        self.executor_kwargs = executor_kwargs
+
+    def set_executor_concurrency(self, executor_concurrency):
+        """
+        Overwrite default or previously set Executor concurrency.
+
+        This only has an effect if Job was initialized with 'as_iterator' and has not yet run.
+        """
+        self.executor_concurrency = executor_concurrency
+
     def cancel(self):
         """Cancel all running and pending Job tasks."""
         if self._as_iterator:
@@ -92,328 +121,148 @@ def __repr__(self):  # pragma: no cover
         return f"<{self.status} Job with {self._total} tasks.>"
 
 
-class TileProcess:
-    """
-    Class to process on a specific process tile.
-
-    If skip is set to True, all attributes will be set to None.
-    """
-
-    def __init__(self, tile=None, config=None, skip=False):
-        """Set attributes depending on baselevels or not."""
-        self.tile = (
-            config.process_pyramid.tile(*tile) if isinstance(tile, tuple) else tile
-        )
-        self.skip = skip
-        self.config_zoom_levels = None if skip else config.zoom_levels
-        self.config_baselevels = None if skip else config.baselevels
-        self.process = None if skip else config.process
-        self.config_dir = None if skip else config.config_dir
-        if (
-            skip
-            or self.tile.zoom not in self.config_zoom_levels
-            or self.tile.zoom in self.config_baselevels
-        ):
-            self.input, self.process_func_params, self.output_params = {}, {}, {}
-        else:
-            self.input = config.get_inputs_for_tile(tile)
-            self.process_func_params = config.get_process_func_params(tile.zoom)
-            self.output_params = config.output_reader.output_params
-        self.mode = None if skip else config.mode
-        self.output_reader = (
-            None if skip or not config.baselevels else config.output_reader
+def task_batches(process, zoom=None, tile=None, skip_output_check=False):
+    with Timer() as duration:
+        # preprocessing tasks
+        yield TaskBatch(
+            id="preprocessing_tasks",
+            tasks=process.config.preprocessing_tasks().values(),
+            func=_preprocess_task_wrapper,
         )
+    logger.debug("preprocessing tasks batch generated in %s", duration)
 
-    def execute(self):
-        """
-        Run the Mapchete process and return the result.
-
-        If baselevels are defined it will generate the result from the other zoom levels
-        accordingly.
-
-        Parameters
-        ----------
-        process_tile : Tile or tile index tuple
-            Member of the process tile pyramid (not necessarily the output
-            pyramid, if output has a different metatiling setting)
-
-        Returns
-        -------
-        data : NumPy array or features
-            process output
-        """
-        if self.mode not in ["memory", "continue", "overwrite"]:  # pragma: no cover
-            raise ValueError("process mode must be memory, continue or overwrite")
-
-        if self.tile.zoom not in self.config_zoom_levels:
-            raise MapcheteNodataTile
-
-        return self._execute()
-
-    def _execute(self):
-        # If baselevel is active and zoom is outside of baselevel,
-        # interpolate from other zoom levels.
-        if self.config_baselevels:
-            if self.tile.zoom < min(self.config_baselevels["zooms"]):
-                return self._interpolate_from_baselevel("lower")
-            elif self.tile.zoom > max(self.config_baselevels["zooms"]):
-                return self._interpolate_from_baselevel("higher")
-        # Otherwise, execute from process file.
-        process_func = get_process_func(
-            process=self.process, config_dir=self.config_dir
-        )
-        try:
-            with Timer() as duration:
-                # Actually run process.
-                process_data = process_func(
-                    MapcheteProcess(
-                        tile=self.tile,
-                        params=self.process_func_params,
-                        input=self.input,
-                        output_params=self.output_params,
-                    ),
-                    **self.process_func_params,
-                )
-        except MapcheteNodataTile:
-            raise
-        except Exception as e:
-            # Log process time
-            logger.exception(
-                (self.tile.id, "exception in user process", e, str(duration))
+    with Timer() as duration:
+        if tile:
+            zoom_levels = [tile.zoom]
+            tiles = {tile.zoom: [tile]}
+        else:
+            zoom_levels = list(
+                reversed(process.config.zoom_levels)
+                if zoom is None
+                else validate_zooms(zoom)
             )
-            new = MapcheteProcessException(format_exc())
-            new.old = e
-            raise new
-
-        return process_data
-
-    def _interpolate_from_baselevel(self, baselevel=None):
-        # This is a special tile derived from a pyramid which has the pixelbuffer setting
-        # from the output pyramid but metatiling from the process pyramid. This is due to
-        # performance reasons as for the usual case overview tiles do not need the
-        # process pyramid pixelbuffers.
-        tile = self.config_baselevels["tile_pyramid"].tile(*self.tile.id)
-
-        # get output_tiles that intersect with process tile
-        output_tiles = (
-            list(self.output_reader.pyramid.tiles_from_bounds(tile.bounds, tile.zoom))
-            if tile.pixelbuffer > self.output_reader.pyramid.pixelbuffer
-            else self.output_reader.pyramid.intersecting(tile)
-        )
-
-        with Timer() as duration:
-            # resample from parent tile
-            if baselevel == "higher":
-                parent_tile = self.tile.get_parent()
-                process_data = raster.resample_from_array(
-                    self.output_reader.read(parent_tile),
-                    in_affine=parent_tile.affine,
-                    out_tile=self.tile,
-                    resampling=self.config_baselevels["higher"],
-                    nodata=self.output_reader.output_params["nodata"],
+            tiles = {
+                zoom: (
+                    tile
+                    for tile, skip, process_msg in _filter_skipable(
+                        process=process,
+                        tiles_batches=process.get_process_tiles(zoom, batch_by="row"),
+                        target_set=None,
+                        skip_output_check=skip_output_check,
+                    )
                 )
-            # resample from children tiles
-            elif baselevel == "lower":
-                if self.output_reader.pyramid.pixelbuffer:  # pragma: no cover
-                    lower_tiles = set(
-                        y
-                        for y in chain(
-                            *[
-                                self.output_reader.pyramid.tiles_from_bounds(
-                                    x.bounds, x.zoom + 1
-                                )
-                                for x in output_tiles
-                            ]
+                for zoom in zoom_levels
+            }
+        if process.config.output.write_in_parent_process:
+            func = _execute
+            fkwargs = dict(append_data=True)
+        else:
+            func = _execute_and_write
+            fkwargs = dict(append_data=True, output_writer=process.config.output)
+
+        # tile tasks
+        for zoom in zoom_levels:
+            yield TileTaskBatch(
+                id=f"zoom_{zoom}",
+                tasks=(
+                    TileTask(
+                        tile=tile,
+                        config=process.config,
+                        skip=(
+                            process.mode == "continue"
+                            and process.config.output_reader.tiles_exist(tile)
                         )
+                        if skip_output_check
+                        else False,
                     )
-                else:
-                    lower_tiles = list(chain(*[x.get_children() for x in output_tiles]))
-                mosaic = raster.create_mosaic(
-                    [
-                        (lower_tile, self.output_reader.read(lower_tile))
-                        for lower_tile in lower_tiles
-                    ],
-                    nodata=self.output_reader.output_params["nodata"],
-                )
-                process_data = raster.resample_from_array(
-                    in_raster=mosaic.data,
-                    in_affine=mosaic.affine,
-                    out_tile=self.tile,
-                    resampling=self.config_baselevels["lower"],
-                    nodata=self.output_reader.output_params["nodata"],
-                )
-        logger.debug((self.tile.id, "generated from baselevel", str(duration)))
-        return process_data
-
-
-class MapcheteProcess(object):
-    """
-    Process class inherited by user process script.
-
-    Its attributes and methods can be accessed via "self" from within a
-    Mapchete process Python file.
-
-    Parameters
-    ----------
-    tile : BufferedTile
-        Tile process should be run on
-    config : MapcheteConfig
-        process configuration
-    params : dictionary
-        process parameters
-
-    Attributes
-    ----------
-    identifier : string
-        process identifier
-    title : string
-        process title
-    version : string
-        process version string
-    abstract : string
-        short text describing process purpose
-    tile : BufferedTile
-        Tile process should be run on
-    tile_pyramid : BufferedTilePyramid
-        process tile pyramid
-    output_pyramid : BufferedTilePyramid
-        output tile pyramid
-    params : dictionary
-        process parameters
-    """
-
-    def __init__(
-        self, tile=None, params=None, input=None, output_params=None, config=None
-    ):
-        """Initialize Mapchete process."""
-        self.identifier = ""
-        self.title = ""
-        self.version = ""
-        self.abstract = ""
-
-        self.tile = tile
-        self.tile_pyramid = tile.tile_pyramid
-        if config is not None:
-            input = config.get_inputs_for_tile(tile)
-            params = config.params_at_zoom(tile.zoom)
-        self.params = dict(params, input=input)
-        self.input = input
-        self.output_params = output_params
-
-    def write(self, data, **kwargs):
-        """Deprecated."""
-        raise DeprecationWarning(
-            "Please return process output data instead of using self.write()."
-        )
-
-    def read(self, **kwargs):
-        """
-        Read existing output data from a previous run.
-
-        Returns
-        -------
-        process output : NumPy array (raster) or feature iterator (vector)
-        """
-        raise DeprecationWarning(
-            "Read existing output from within a process is deprecated"
-        )
+                    for tile in tiles[zoom]
+                ),
+                func=func,
+                fkwargs=fkwargs,
+            )
+    logger.debug("tile task batches generated in %s", duration)
 
-    @deprecated_kwargs
-    def open(self, input_id, **kwargs):
-        """
-        Open input data.
-
-        Parameters
-        ----------
-        input_id : string
-            input identifier from configuration file or file path
-        kwargs : driver specific parameters (e.g. resampling)
-
-        Returns
-        -------
-        tiled input data : InputTile
-            reprojected input data within tile
-        """
-        if input_id not in self.input:
-            raise ValueError("%s not found in config as input" % input_id)
-        return self.input[input_id]
 
-    def hillshade(self, elevation, azimuth=315.0, altitude=45.0, z=1.0, scale=1.0):
-        """
-        Calculate hillshading from elevation data.
-
-        Parameters
-        ----------
-        elevation : array
-            input elevation data
-        azimuth : float
-            horizontal angle of light source (315: North-West)
-        altitude : float
-            vertical angle of light source (90 would result in slope shading)
-        z : float
-            vertical exaggeration factor
-        scale : float
-            scale factor of pixel size units versus height units (insert 112000
-            when having elevation values in meters in a geodetic projection)
-
-        Returns
-        -------
-        hillshade : array
-        """
-        return commons_hillshade.hillshade(
-            elevation, self.tile, azimuth, altitude, z, scale
-        )
+def compute(
+    process,
+    zoom=None,
+    tile=None,
+    executor=None,
+    concurrency="processes",
+    workers=None,
+    dask_scheduler=None,
+    multiprocessing_start_method=None,
+    multiprocessing_module=None,
+    skip_output_check=False,
+    dask_compute_graph=True,
+    raise_errors=True,
+    with_results=False,
+    **kwargs,
+):
+    """Computes all tasks and yields progress."""
+    concurrency = "dask" if dask_scheduler else concurrency
+    num_processed = 0
+    with ExitStack() as exit_stack:
+        # create fresh Executor if it was not passed on
+        if executor is None:
+            executor = exit_stack.enter_context(
+                Executor(
+                    max_workers=workers,
+                    concurrency=concurrency,
+                    start_method=multiprocessing_start_method,
+                    multiprocessing_module=multiprocessing_module,
+                    dask_scheduler=dask_scheduler,
+                )
+            )
 
-    def contours(self, elevation, interval=100, field="elev", base=0):
-        """
-        Extract contour lines from elevation data.
-
-        Parameters
-        ----------
-        elevation : array
-            input elevation data
-        interval : integer
-            elevation value interval when drawing contour lines
-        field : string
-            output field name containing elevation value
-        base : integer
-            elevation base value the intervals are computed from
-
-        Returns
-        -------
-        contours : iterable
-            contours as GeoJSON-like pairs of properties and geometry
-        """
-        return commons_contours.extract_contours(
-            elevation, self.tile, interval=interval, field=field, base=base
-        )
+        logger.info("run process on area")
+        duration = exit_stack.enter_context(Timer())
+        if tile:
+            tile = process.config.process_pyramid.tile(*tile)
+            zoom_levels = [tile.zoom]
+        else:
+            zoom_levels = list(
+                reversed(process.config.zoom_levels)
+                if zoom is None
+                else validate_zooms(zoom)
+            )
+        if dask_compute_graph and isinstance(executor, DaskExecutor):
+            for num_processed, future in enumerate(
+                _compute_task_graph(
+                    executor=executor,
+                    process=process,
+                    zoom_levels=zoom_levels,
+                    tile=tile,
+                    skip_output_check=skip_output_check,
+                    with_results=with_results,
+                    raise_errors=raise_errors,
+                ),
+                1,
+            ):
+                if raise_errors and future.exception():  # pragma: no cover
+                    logger.exception(future.exception())
+                    raise future.exception()
+                logger.debug("task %s finished: %s", num_processed, future)
+                yield future
+        else:
+            for num_processed, future in enumerate(
+                _compute_tasks(
+                    executor=executor,
+                    process=process,
+                    zoom_levels=zoom_levels,
+                    tile=tile,
+                    skip_output_check=skip_output_check,
+                    **kwargs,
+                ),
+                1,
+            ):
+                if raise_errors and future.exception():  # pragma: no cover
+                    logger.exception(future.exception())
+                    raise future.exception()
+                logger.debug("task %s finished: %s", num_processed, future)
+                yield future
 
-    def clip(self, array, geometries, inverted=False, clip_buffer=0):
-        """
-        Clip array by geometry.
-
-        Parameters
-        ----------
-        array : array
-            raster data to be clipped
-        geometries : iterable
-            geometries used to clip source array
-        inverted : bool
-            invert clipping (default: False)
-        clip_buffer : int
-            buffer (in pixels) geometries before applying clip
-
-        Returns
-        -------
-        clipped array : array
-        """
-        return commons_clip.clip_array_with_vector(
-            array,
-            self.tile.affine,
-            geometries,
-            inverted=inverted,
-            clip_buffer=clip_buffer * self.tile.pixel_x_size,
-        )
+    logger.info("computed %s tasks in %s", num_processed, duration)
 
 
 #######################
@@ -421,9 +270,11 @@ def clip(self, array, geometries, inverted=False, clip_buffer=0):
 #######################
 
 
-def _preprocess_task_wrapper(task_tuple):
-    task_key, (func, fargs, fkwargs) = task_tuple
-    return task_key, func(*fargs, **fkwargs)
+def _preprocess_task_wrapper(task, append_data=True, **kwargs):
+    data = task.execute(**kwargs)
+    return PreprocessingProcessInfo(
+        task_key=task.id, data=data if append_data else None
+    )
 
 
 def _preprocess(
@@ -440,46 +291,45 @@ def _preprocess(
     # If preprocessing tasks already finished, don't run them again.
     if process.config.preprocessing_tasks_finished:  # pragma: no cover
         return
+    num_processed = 0
 
     # If an Executor is passed on, don't close after processing. If no Executor is passed on,
     # create one and properly close it afterwards.
-    create_executor = executor is None
-    executor = executor or Executor(
-        max_workers=workers,
-        concurrency="dask" if dask_scheduler else "processes",
-        start_method=multiprocessing_start_method,
-        multiprocessing_module=multiprocessing_module,
-        dask_scheduler=dask_scheduler,
-    )
-    try:
-        with Timer() as duration:
-            logger.info(
-                "run preprocessing on %s tasks using %s workers", len(tasks), workers
-            )
-
-            # process all remaining tiles using todo list from before
-            for i, future in enumerate(
-                executor.as_completed(
-                    func=_preprocess_task_wrapper,
-                    iterable=list(tasks.items()),
-                    max_submitted_tasks=dask_max_submitted_tasks,
-                    chunksize=dask_chunksize,
-                ),
-                1,
-            ):
-                task_key, result = future.result()
-                logger.debug(
-                    "preprocessing task %s/%s %s processed successfully",
-                    i,
-                    len(tasks),
-                    task_key,
+    with ExitStack() as exit_stack:
+        # create fresh Executor if it was not passed on
+        if executor is None:
+            executor = exit_stack.enter_context(
+                Executor(
+                    max_workers=workers,
+                    concurrency="dask" if dask_scheduler else "processes",
+                    start_method=multiprocessing_start_method,
+                    multiprocessing_module=multiprocessing_module,
+                    dask_scheduler=dask_scheduler,
                 )
-                process.config.set_preprocessing_task_result(task_key, result)
-                yield f"preprocessing task {task_key} finished"
-    finally:
-        if create_executor:
-            executor.close()
-
+            )
+        duration = exit_stack.enter_context(Timer())
+        logger.info(
+            "run preprocessing on %s tasks using %s workers", len(tasks), workers
+        )
+        for num_processed, future in enumerate(
+            executor.as_completed(
+                func=_preprocess_task_wrapper,
+                iterable=tasks.values(),
+                fkwargs=dict(append_data=True),
+                max_submitted_tasks=dask_max_submitted_tasks,
+                chunksize=dask_chunksize,
+            ),
+            1,
+        ):
+            result = future.result()
+            logger.debug(
+                "preprocessing task %s/%s %s processed successfully",
+                num_processed,
+                len(tasks),
+                result.task_key,
+            )
+            process.config.set_preprocessing_task_result(result.task_key, result.data)
+            yield future
     process.config.preprocessing_tasks_finished = True
 
     logger.info("%s task(s) iterated in %s", str(len(tasks)), duration)
@@ -496,18 +346,21 @@ def _run_on_single_tile(
     tile=None,
     dask_scheduler=None,
 ):
-    logger.info("run process on single tile")
-    create_executor = executor is None
-    executor = executor or Executor(
-        concurrency="dask" if dask_scheduler else None,
-        dask_scheduler=dask_scheduler,
-    )
-    try:
+    with ExitStack() as exit_stack:
+        # create fresh Executor if it was not passed on
+        if executor is None:
+            executor = exit_stack.enter_context(
+                Executor(
+                    concurrency="dask" if dask_scheduler else None,
+                    dask_scheduler=dask_scheduler,
+                )
+            )
+        logger.info("run process on single tile")
         return next(
             executor.as_completed(
                 func=_execute_and_write,
                 iterable=[
-                    TileProcess(
+                    TileTask(
                         tile=tile,
                         config=process.config,
                         skip=(
@@ -518,10 +371,7 @@ def _run_on_single_tile(
                 ],
                 fkwargs=dict(output_writer=process.config.output),
             )
-        ).result()
-    finally:
-        if create_executor:
-            executor.close()
+        )
 
 
 def _run_area(
@@ -541,7 +391,7 @@ def _run_area(
 
     # for output drivers requiring writing data in parent process
     if process.config.output.write_in_parent_process:
-        for process_info in _run_multi(
+        for future in _run_multi(
             executor=executor,
             func=_execute,
             zoom_levels=zoom_levels,
@@ -555,11 +405,11 @@ def _run_area(
             write_in_parent_process=True,
             skip_output_check=skip_output_check,
         ):
-            yield process_info
+            yield future
 
     # for output drivers which can write data in child processes
     else:
-        for process_info in _run_multi(
+        for future in _run_multi(
             executor=executor,
             func=_execute_and_write,
             fkwargs=dict(output_writer=process.config.output),
@@ -574,7 +424,7 @@ def _run_area(
             write_in_parent_process=False,
             skip_output_check=skip_output_check,
         ):
-            yield process_info
+            yield future
 
 
 def _filter_skipable(
@@ -612,118 +462,348 @@ def _run_multi(
     workers = min([workers, total_tiles])
     num_processed = 0
 
-    # here we store the parents of processed tiles so we can update overviews
-    # also in "continue" mode in case there were updates at the baselevel
-    overview_parents = set()
-
     # If an Executor is passed on, don't close after processing. If no Executor is passed on,
     # create one and properly close it afterwards.
-    create_executor = executor is None
-    executor = executor or Executor(
-        max_workers=workers,
-        concurrency="dask" if dask_scheduler else "processes",
-        start_method=multiprocessing_start_method,
-        multiprocessing_module=multiprocessing_module,
-        dask_scheduler=dask_scheduler,
-    )
-
-    try:
+    with ExitStack() as exit_stack:
+        # create fresh Executor if it was not passed on
+        if executor is None:
+            executor = exit_stack.enter_context(
+                Executor(
+                    max_workers=workers,
+                    concurrency="dask" if dask_scheduler else "processes",
+                    start_method=multiprocessing_start_method,
+                    multiprocessing_module=multiprocessing_module,
+                    dask_scheduler=dask_scheduler,
+                )
+            )
         with Timer() as duration:
             logger.info(
-                "run process on %s tiles using %s workers", total_tiles, workers
+                "run process on %s tiles using %s workers on executor %s",
+                total_tiles,
+                workers,
+                executor,
             )
+            if process.config.baselevels:
+                f = _run_multi_overviews
+            else:
+                f = _run_multi_no_overviews
+            for num_processed, future in enumerate(
+                f(
+                    zoom_levels=zoom_levels,
+                    executor=executor,
+                    func=func,
+                    process=process,
+                    skip_output_check=skip_output_check,
+                    fkwargs=fkwargs,
+                    dask_chunksize=dask_chunksize,
+                    dask_max_submitted_tasks=dask_max_submitted_tasks,
+                    write_in_parent_process=write_in_parent_process,
+                ),
+                1,
+            ):
+                logger.debug("task %s finished: %s", num_processed, future)
+                yield future
 
-            for i, zoom in enumerate(zoom_levels):
+        logger.info("%s tile(s) iterated in %s", str(num_processed), duration)
 
-                # get generator list of tiles, whether they are to be skipped and skip_info
-                # from _filter_skipable and pass on to executor
-                for future in executor.as_completed(
-                    func=func,
-                    iterable=(
-                        (
-                            TileProcess(
-                                tile=tile,
-                                config=process.config,
-                                skip=(
-                                    process.mode == "continue"
-                                    and process.config.output_reader.tiles_exist(tile)
-                                )
-                                if skip_output_check
-                                else False,
-                            ),
-                            skip,
-                            process_msg,
-                        )
-                        for tile, skip, process_msg in _filter_skipable(
-                            process=process,
-                            tiles_batches=process.get_process_tiles(
-                                zoom, batch_by="row"
-                            ),
-                            target_set=(
-                                overview_parents
-                                if process.config.baselevels and i
-                                else None
-                            ),
-                            skip_output_check=skip_output_check,
-                        )
+
+def _compute_task_graph(
+    executor=None,
+    process=None,
+    skip_output_check=False,
+    zoom_levels=None,
+    tile=None,
+    with_results=False,
+    raise_errors=False,
+    **kwargs,
+):
+    # TODO optimize memory management, e.g. delete preprocessing tasks from input
+    # once the dask graph is ready.
+    from dask.delayed import delayed
+    from distributed import as_completed
+
+    # materialize all tasks including dependencies
+    with Timer() as t:
+        coll = to_dask_collection(
+            process.task_batches(
+                zoom=zoom_levels, tile=tile, skip_output_check=skip_output_check
+            )
+        )
+    logger.debug("%s dask collection generated in %s", len(coll), t)
+
+    # send to scheduler
+    with Timer() as t:
+        futures = executor._executor.compute(coll, optimize_graph=True, traverse=True)
+    logger.debug("sent to scheduler in %s", t)
+
+    for future in as_completed(
+        futures, with_results=with_results, raise_errors=raise_errors
+    ):
+        futures.remove(future)
+        if process.config.output.write_in_parent_process:
+            yield FinishedFuture(
+                result=_write(
+                    process_info=future.result(),
+                    output_writer=process.config.output,
+                    append_output=True,
+                )
+            )
+        else:
+            yield future
+
+
+def _compute_tasks(
+    executor=None,
+    process=None,
+    workers=None,
+    zoom_levels=None,
+    tile=None,
+    skip_output_check=False,
+    **kwargs,
+):
+    num_processed = 0
+    if not process.config.preprocessing_tasks_finished:
+        tasks = process.config.preprocessing_tasks()
+        logger.info(
+            "run preprocessing on %s tasks using %s workers", len(tasks), workers
+        )
+        # process all remaining tiles using todo list from before
+        for i, future in enumerate(
+            executor.as_completed(
+                func=_preprocess_task_wrapper,
+                iterable=tasks.values(),
+                fkwargs=dict(append_data=True),
+                **kwargs,
+            ),
+            1,
+        ):
+            result = future.result()
+            logger.debug(
+                "preprocessing task %s/%s %s processed successfully",
+                i,
+                len(tasks),
+                result.task_key,
+            )
+            process.config.set_preprocessing_task_result(result.task_key, result.data)
+            yield future
+
+    # run single tile
+    if tile:
+        logger.info("run process on single tile")
+        yield next(
+            executor.as_completed(
+                func=_execute_and_write,
+                iterable=[
+                    TileTask(
+                        tile=tile,
+                        config=process.config,
+                        skip=(
+                            process.config.mode == "continue"
+                            and process.config.output_reader.tiles_exist(tile)
+                        ),
                     ),
-                    fkwargs=fkwargs,
-                    max_submitted_tasks=dask_max_submitted_tasks,
-                    chunksize=dask_chunksize,
-                    item_skip_bool=True,
-                ):
-                    # tiles which were not processed
-                    if isinstance(future, SkippedFuture):
-                        process_info = ProcessInfo(
-                            tile=future.result().tile,
-                            processed=False,
-                            process_msg=future.skip_info,
-                            written=False,
-                            write_msg="nothing written",
+                ],
+                fkwargs=dict(output_writer=process.config.output),
+            )
+        )
+
+    else:
+        # for output drivers requiring writing data in parent process
+        if process.config.output.write_in_parent_process:
+            func = _execute
+            fkwargs = dict(append_data=False)
+            write_in_parent_process = True
+
+        # for output drivers which can write data in child processes
+        else:
+            func = _execute_and_write
+            fkwargs = dict(append_data=False, output_writer=process.config.output)
+            write_in_parent_process = False
+
+        # for outputs which have overviews
+        if process.config.baselevels:
+            _process_batches = _run_multi_overviews
+        # for outputs with no overviews
+        else:
+            _process_batches = _run_multi_no_overviews
+
+        for num_processed, future in enumerate(
+            _process_batches(
+                zoom_levels=zoom_levels,
+                executor=executor,
+                func=func,
+                process=process,
+                skip_output_check=skip_output_check,
+                fkwargs=fkwargs,
+                write_in_parent_process=write_in_parent_process,
+                **kwargs,
+            ),
+            1,
+        ):
+            logger.debug(
+                "task %s finished",
+                num_processed,
+            )
+            yield future
+
+
+def _run_multi_overviews(
+    zoom_levels=None,
+    executor=None,
+    func=None,
+    process=None,
+    skip_output_check=None,
+    fkwargs=None,
+    dask_chunksize=None,
+    dask_max_submitted_tasks=None,
+    write_in_parent_process=None,
+):
+    # here we store the parents of processed tiles so we can update overviews
+    # also in "continue" mode in case there were updates at the baselevel
+    overview_parents = set()
+
+    for i, zoom in enumerate(zoom_levels):
+
+        # get generator list of tiles, whether they are to be skipped and skip_info
+        # from _filter_skipable and pass on to executor
+        for future in executor.as_completed(
+            func=func,
+            iterable=(
+                (
+                    TileTask(
+                        tile=tile,
+                        config=process.config,
+                        skip=(
+                            process.mode == "continue"
+                            and process.config.output_reader.tiles_exist(tile)
                         )
-                    # tiles which were processed
-                    else:
-                        # trigger output write for driver which require parent process for writing
-                        if write_in_parent_process:
-                            output_data, process_info = future.result()
-                            process_info = _write(
-                                process_info=process_info,
-                                output_data=output_data,
-                                output_writer=process.config.output,
-                            )
-
-                        # output already has been written, so just use task process info
-                        else:
-                            process_info = future.result()
-                            # in case of building overviews from baselevels, remember which parent
-                            # tile needs to be updated later on
-                            if (
-                                not skip_output_check
-                                and process.config.baselevels
-                                and process_info.processed
-                                and process_info.tile.zoom > 0
-                            ):
-                                overview_parents.add(process_info.tile.get_parent())
-                    num_processed += 1
-                    logger.debug(
-                        "tile %s/%s finished: %s, %s, %s",
-                        num_processed,
-                        total_tiles,
-                        process_info.tile,
-                        process_info.process_msg,
-                        process_info.write_msg,
+                        if skip_output_check
+                        else False,
+                    ),
+                    skip,
+                    process_msg,
+                )
+                for tile, skip, process_msg in _filter_skipable(
+                    process=process,
+                    tiles_batches=process.get_process_tiles(zoom, batch_by="row"),
+                    target_set=(
+                        overview_parents if process.config.baselevels and i else None
+                    ),
+                    skip_output_check=skip_output_check,
+                )
+            ),
+            fkwargs=fkwargs,
+            max_submitted_tasks=dask_max_submitted_tasks,
+            chunksize=dask_chunksize,
+            item_skip_bool=True,
+        ):
+            # tiles which were not processed
+            if isinstance(future, SkippedFuture):
+                process_info = TileProcessInfo(
+                    tile=future.result().tile,
+                    processed=False,
+                    process_msg=future.skip_info,
+                    written=False,
+                    write_msg="nothing written",
+                )
+            # tiles which were processed
+            else:
+                # trigger output write for driver which require parent process for writing
+                # the code coverage below is omitted because we don't usually calculate overviews
+                # when writing single files
+                if write_in_parent_process:  # pragma: no cover
+                    process_info = _write(
+                        process_info=future.result(),
+                        output_writer=process.config.output,
                     )
-                    try:
-                        overview_parents.remove(process_info.tile)
-                    except KeyError:
-                        pass
-                    yield process_info
 
-    finally:
-        if create_executor:
-            executor.close()
+                # output already has been written, so just use task process info
+                else:
+                    process_info = future.result()
+                    # in case of building overviews from baselevels, remember which parent
+                    # tile needs to be updated later on
+                    if (
+                        not skip_output_check
+                        and process.config.baselevels
+                        and process_info.processed
+                        and process_info.tile.zoom > 0
+                    ):
+                        overview_parents.add(process_info.tile.get_parent())
+            try:
+                overview_parents.remove(process_info.tile)
+            except KeyError:
+                pass
+            yield FinishedFuture(result=process_info)
+
+
+def _run_multi_no_overviews(
+    zoom_levels=None,
+    executor=None,
+    func=None,
+    process=None,
+    skip_output_check=None,
+    fkwargs=None,
+    dask_chunksize=None,
+    dask_max_submitted_tasks=None,
+    write_in_parent_process=None,
+):
+    # get generator list of tiles, whether they are to be skipped and skip_info
+    # from _filter_skipable and pass on to executor
+    for future in executor.as_completed(
+        func=func,
+        iterable=(
+            (
+                TileTask(
+                    tile=tile,
+                    config=process.config,
+                    skip=(
+                        process.mode == "continue"
+                        and process.config.output_reader.tiles_exist(tile)
+                    )
+                    if skip_output_check
+                    else False,
+                ),
+                skip,
+                process_msg,
+            )
+            for tile, skip, process_msg in _filter_skipable(
+                process=process,
+                tiles_batches=(
+                    batch
+                    for zoom in zoom_levels
+                    for batch in process.get_process_tiles(zoom, batch_by="row")
+                ),
+                target_set=None,
+                skip_output_check=skip_output_check,
+            )
+        ),
+        fkwargs=fkwargs,
+        max_submitted_tasks=dask_max_submitted_tasks,
+        chunksize=dask_chunksize,
+        item_skip_bool=True,
+    ):
+        # tiles which were not processed
+        if isinstance(future, SkippedFuture):
+            process_info = TileProcessInfo(
+                tile=future.result().tile,
+                processed=False,
+                process_msg=future.skip_info,
+                written=False,
+                write_msg="nothing written",
+            )
+        # tiles which were processed
+        else:
+            # trigger output write for driver which require parent process for writing
+            if write_in_parent_process:
+                process_info = _write(
+                    process_info=future.result(),
+                    output_writer=process.config.output,
+                )
 
-    logger.info("%s tile(s) iterated in %s", str(num_processed), duration)
+            # output already has been written, so just use task process info
+            else:
+                process_info = future.result()
+        yield FinishedFuture(result=process_info)
 
 
 ###############################
@@ -731,7 +811,7 @@ def _run_multi(
 ###############################
 
 
-def _execute(tile_process=None, **_):
+def _execute(tile_process=None, dependencies=None, **_):
     logger.debug(
         (tile_process.tile.id, "running on %s" % multiprocessing.current_process().name)
     )
@@ -739,65 +819,70 @@ def _execute(tile_process=None, **_):
     # skip execution if overwrite is disabled and tile exists
     if tile_process.skip:
         logger.debug((tile_process.tile.id, "tile exists, skipping"))
-        return None, ProcessInfo(
+        return TileProcessInfo(
             tile=tile_process.tile,
             processed=False,
             process_msg="output already exists",
             written=False,
             write_msg="nothing written",
+            data=None,
         )
 
     # execute on process tile
     with Timer() as duration:
         try:
-            output = tile_process.execute()
+            output = tile_process.execute(dependencies=dependencies)
         except MapcheteNodataTile:  # pragma: no cover
             output = "empty"
     processor_message = "processed in %s" % duration
     logger.debug((tile_process.tile.id, processor_message))
-    return output, ProcessInfo(
+    return TileProcessInfo(
         tile=tile_process.tile,
         processed=True,
         process_msg=processor_message,
         written=None,
         write_msg=None,
+        data=output,
     )
 
 
-def _write(process_info=None, output_data=None, output_writer=None, **_):
+def _write(process_info=None, output_writer=None, append_data=False, **_):
     if process_info.processed:
         try:
-            output_data = output_writer.streamline_output(output_data)
+            output_data = output_writer.streamline_output(process_info.data)
         except MapcheteNodataTile:
             output_data = None
         if output_data is None:
             message = "output empty, nothing written"
             logger.debug((process_info.tile.id, message))
-            return ProcessInfo(
+            return TileProcessInfo(
                 tile=process_info.tile,
                 processed=process_info.processed,
                 process_msg=process_info.process_msg,
                 written=False,
                 write_msg=message,
             )
-
         with Timer() as duration:
             output_writer.write(process_tile=process_info.tile, data=output_data)
         message = "output written in %s" % duration
         logger.debug((process_info.tile.id, message))
-        return ProcessInfo(
+        return TileProcessInfo(
             tile=process_info.tile,
             processed=process_info.processed,
             process_msg=process_info.process_msg,
             written=True,
             write_msg=message,
+            data=output_data if append_data else None,
         )
 
     return process_info
 
 
-def _execute_and_write(tile_process=None, output_writer=None, **_):
-    output_data, process_info = _execute(tile_process=tile_process)
+def _execute_and_write(
+    tile_process=None, output_writer=None, dependencies=None, append_data=False, **_
+):
     return _write(
-        process_info=process_info, output_data=output_data, output_writer=output_writer
+        process_info=_execute(tile_process=tile_process, dependencies=dependencies),
+        output_writer=output_writer,
+        append_data=append_data,
     )
diff --git a/mapchete/_tasks.py b/mapchete/_tasks.py
new file mode 100644
index 00000000..72a557c4
--- /dev/null
+++ b/mapchete/_tasks.py
@@ -0,0 +1,432 @@
+from collections import namedtuple
+from itertools import chain
+import logging
+from shapely.geometry import box, mapping
+from traceback import format_exc
+from uuid import uuid4
+
+from mapchete._timer import Timer
+from mapchete.config import get_process_func
+from mapchete.errors import (
+    MapcheteNodataTile,
+    NoTaskGeometry,
+    MapcheteProcessException,
+    MapcheteProcessOutputError,
+)
+from mapchete.io import raster
+from mapchete.io._geometry_operations import to_shape
+from mapchete.io.vector import IndexedFeatures
+from mapchete.validate import validate_bounds
+from mapchete._user_process import MapcheteProcess
+
+
+logger = logging.getLogger(__name__)
+
+
+TaskResult = namedtuple(
+    "TaskResult",
+    "task_id tile processed process_msg result",
+    defaults=(None, None, None, None, None),
+)
+
+
+class Task:
+    """Generic processing task.
+
+    Can optionally have spatial properties attached which helps building up dependencies
+    between tasks.
+    """
+
+    def __init__(
+        self,
+        id=None,
+        func=None,
+        fargs=None,
+        fkwargs=None,
+        geometry=None,
+        bounds=None,
+        dependencies=None,
+    ):
+        self.id = id or uuid4().hex
+        self.func = func
+        self.fargs = fargs or ()
+        self.fkwargs = fkwargs or {}
+        self.dependencies = dependencies or {}
+        if geometry and bounds:
+            raise ValueError("only provide one of either 'geometry' or 'bounds'")
+        elif geometry:
+            self.geometry = to_shape(geometry)
+            self.bounds = self.geometry.bounds
+        elif bounds:
+            self.bounds = validate_bounds(bounds)
+            self.geometry = box(*self.bounds)
+        else:
+            self.bounds, self.geometry = None, None
+
+    def __repr__(self):  # pragma: no cover
+        return f"Task(id={self.id}, bounds={self.bounds})"
+
+    def to_dict(self):
+        return {
+            "id": self.id,
+            "bounds": self.bounds,
+            "geometry": mapping(self.geometry) if self.geometry else None,
+            "properties": {},
+        }
+
+    def add_dependencies(self, dependencies):
+        dependencies = dependencies or {}
+        if not isinstance(dependencies, dict):
+            raise TypeError(
+                f"dependencies must be a dictionary, not {type(dependencies)}"
+            )
+        self.dependencies.update(dependencies)
+
+    def execute(self, dependencies=None):
+        return self.func(*self.fargs, **self.fkwargs)
+
+    def has_geometry(self):
+        return self.geometry is not None
+
+    @property
+    def __geo_interface__(self):
+        if self.has_geometry():
+            return mapping(self.geometry)
+        else:
+            raise NoTaskGeometry(f"{self} has no geo information assigned")
+
+
+def _execute_task_wrapper(task, **kwargs):
+    return task.execute(**kwargs)
+
+
+class TaskBatch:
+    def __init__(self, tasks=None, id=None, func=None, fkwargs=None):
+        if tasks is None:
+            raise TypeError("TaskBatch requires at least one Task")
+        self.id = id or uuid4().hex
+        self.tasks = IndexedFeatures(
+            (self._validate(t) for t in tasks), allow_non_geo_objects=True
+        )
+        self.bounds = self.tasks.bounds
+        self.func = func or _execute_task_wrapper
+        self.fkwargs = fkwargs or {}
+
+    def __repr__(self):  # pragma: no cover
+        return f"TaskBatch(id={self.id}, bounds={self.bounds})"
+
+    def __iter__(self):
+        return iter(self.tasks.values())
+
+    def __len__(self):
+        return len(self.tasks)
+
+    def items(self):
+        return self.tasks.items()
+
+    def keys(self):
+        return self.tasks.keys()
+
+    def values(self):
+        return self.tasks.values()
+
+    def intersection(self, other):
+        if isinstance(other, Task):
+            return self.tasks.filter(other.bounds)
+        elif isinstance(other, (tuple, list)):
+            return self.tasks.filter(validate_bounds(other))
+        else:
+            raise TypeError(
+                f"intersection only works with other Task instances or bounds, not {other}"
+            )
+
+    def _validate(self, item):
+        if isinstance(item, Task):
+            return item
+        else:
+            raise TypeError("TaskBatch items must be Taskss, not %s", type(item))
+
+
+class TileTask(Task):
+    """
+    Class to process on a specific process tile.
+
+    If skip is set to True, all attributes will be set to None.
+    """
+
+    def __init__(self, tile=None, id=None, config=None, skip=False, dependencies=None):
+        """Set attributes depending on baselevels or not."""
+        self.tile = (
+            config.process_pyramid.tile(*tile) if isinstance(tile, tuple) else tile
+        )
+        self.id = id or f"tile_task_{self.tile.zoom}-{self.tile.row}-{self.tile.col}"
+        self.skip = skip
+        self.config_zoom_levels = None if skip else config.zoom_levels
+        self.config_baselevels = None if skip else config.baselevels
+        self.process = None if skip else config.process
+        self.config_dir = None if skip else config.config_dir
+        if (
+            skip
+            or self.tile.zoom not in self.config_zoom_levels
+            or self.tile.zoom in self.config_baselevels
+        ):
+            self.input, self.process_func_params, self.output_params = {}, {}, {}
+        else:
+            self.input = config.get_inputs_for_tile(tile)
+            self.process_func_params = config.get_process_func_params(tile.zoom)
+            self.output_params = config.output_reader.output_params
+        self.mode = None if skip else config.mode
+        self.output_reader = (
+            None if skip or not config.baselevels else config.output_reader
+        )
+        super().__init__(
+            id=f"tile_task_{tile.zoom}-{tile.row}-{tile.col}", geometry=tile.bbox
+        )
+
+    def execute(self, dependencies=None):
+        """
+        Run the Mapchete process and return the result.
+
+        If baselevels are defined it will generate the result from the other zoom levels
+        accordingly.
+
+        Parameters
+        ----------
+        process_tile : Tile or tile index tuple
+            Member of the process tile pyramid (not necessarily the output
+            pyramid, if output has a different metatiling setting)
+
+        Returns
+        -------
+        data : NumPy array or features
+            process output
+        """
+        if self.mode not in ["memory", "continue", "overwrite"]:  # pragma: no cover
+            raise ValueError("process mode must be memory, continue or overwrite")
+
+        if self.tile.zoom not in self.config_zoom_levels:
+            raise MapcheteNodataTile
+
+        dependencies = dependencies or {}
+        process_output = self._execute(dependencies=dependencies)
+        if isinstance(process_output, str) and process_output == "empty":
+            raise MapcheteNodataTile
+        elif process_output is None:
+            raise MapcheteProcessOutputError("process output is empty")
+        return process_output
+
+    def _execute(self, dependencies=None):
+        # If baselevel is active and zoom is outside of baselevel,
+        # interpolate from other zoom levels.
+        if self.config_baselevels:
+            if self.tile.zoom < min(self.config_baselevels["zooms"]):
+                return self._interpolate_from_baselevel("lower", dependencies)
+            elif self.tile.zoom > max(self.config_baselevels["zooms"]):
+                return self._interpolate_from_baselevel("higher", dependencies)
+        # Otherwise, execute from process file.
+        process_func = get_process_func(
+            process=self.process, config_dir=self.config_dir
+        )
+        try:
+            with Timer() as duration:
+                # append dependent preprocessing task results to input objects
+                if dependencies:
+                    for task_key, task_result in dependencies.items():
+                        if not task_key.startswith("tile_task"):
+                            inp_key, task_key = task_key.split(":")[0], ":".join(
+                                task_key.split(":")[1:]
+                            )
+                            if task_key in [None, ""]:  # pragma: no cover
+                                raise KeyError(f"malformed task key: {inp_key}")
+                            for inp in self.input.values():
+                                if inp.input_key == inp_key:
+                                    inp.set_preprocessing_task_result(
+                                        task_key=task_key, result=task_result.data
+                                    )
+                # Actually run process.
+                process_data = process_func(
+                    MapcheteProcess(
+                        tile=self.tile,
+                        params=self.process_func_params,
+                        input=self.input,
+                        output_params=self.output_params,
+                    ),
+                    **self.process_func_params,
+                )
+        except MapcheteNodataTile:
+            raise
+        except Exception as e:
+            # Log process time
+            logger.exception(e)
+            logger.error((self.tile.id, "exception in user process", e, str(duration)))
+            new = MapcheteProcessException(format_exc())
+            new.old = e
+            raise new
+
+        return process_data
+
+    def _interpolate_from_baselevel(self, baselevel=None, dependencies=None):
+        # This is a special tile derived from a pyramid which has the pixelbuffer setting
+        # from the output pyramid but metatiling from the process pyramid. This is due to
+        # performance reasons as for the usual case overview tiles do not need the
+        # process pyramid pixelbuffers.
+        tile = self.config_baselevels["tile_pyramid"].tile(*self.tile.id)
+
+        # get output_tiles that intersect with process tile
+        output_tiles = (
+            list(self.output_reader.pyramid.tiles_from_bounds(tile.bounds, tile.zoom))
+            if tile.pixelbuffer > self.output_reader.pyramid.pixelbuffer
+            else self.output_reader.pyramid.intersecting(tile)
+        )
+
+        with Timer() as duration:
+            # resample from parent tile
+            if baselevel == "higher":
+                parent_tile = self.tile.get_parent()
+                process_data = raster.resample_from_array(
+                    self.output_reader.read(parent_tile),
+                    in_affine=parent_tile.affine,
+                    out_tile=self.tile,
+                    resampling=self.config_baselevels["higher"],
+                    nodata=self.output_reader.output_params["nodata"],
+                )
+            # resample from children tiles
+            elif baselevel == "lower":
+                src_tiles = {}
+                for process_info in dependencies.values():
+                    logger.debug("reading output from dependend tasks")
+                    process_tile = process_info.tile
+                    for output_tile in self.output_reader.pyramid.intersecting(
+                        process_tile
+                    ):
+                        if process_info.data is not None:
+                            src_tiles[output_tile] = raster.extract_from_array(
+                                in_raster=process_info.data,
+                                in_affine=process_tile.affine,
+                                out_tile=output_tile,
+                            )
+                if self.output_reader.pyramid.pixelbuffer:  # pragma: no cover
+                    # if there is a pixelbuffer around the output tiles, we need to read more child tiles
+                    child_tiles = [
+                        child_tile
+                        for output_tile in output_tiles
+                        for child_tile in self.output_reader.pyramid.tiles_from_bounds(
+                            output_tile.bounds, output_tile.zoom + 1
+                        )
+                    ]
+                else:
+                    child_tiles = [
+                        child_tile
+                        for output_tile in output_tiles
+                        for child_tile in output_tile.get_children()
+                    ]
+                for child_tile in child_tiles:
+                    if child_tile not in src_tiles:
+                        src_tiles[child_tile] = self.output_reader.read(child_tile)
+
+                process_data = raster.resample_from_array(
+                    in_raster=raster.create_mosaic(
+                        [(src_tile, data) for src_tile, data in src_tiles.items()],
+                        nodata=self.output_reader.output_params["nodata"],
+                    ),
+                    out_tile=self.tile,
+                    resampling=self.config_baselevels["lower"],
+                    nodata=self.output_reader.output_params["nodata"],
+                )
+        logger.debug((self.tile.id, "generated from baselevel", str(duration)))
+        return process_data
+
+
+class TileTaskBatch(TaskBatch):
+    """Combines TileTask instances of same pyramid and zoom level into one batch."""
+
+    def __init__(self, tasks=None, id=None, func=None, fkwargs=None):
+        self.id = id or uuid4().hex
+        self.bounds = None, None, None, None
+        self.tasks = {tile: item for tile, item in self._validate(tasks)}
+        self._update_bounds()
+        self.func = func or _execute_task_wrapper
+        self.fkwargs = fkwargs or {}
+
+    def _update_bounds(self):
+        for tile in self.tasks.keys():
+            left, bottom, right, top = self.bounds
+            self.bounds = (
+                tile.left if left is None else min(left, tile.left),
+                tile.bottom if bottom is None else min(bottom, tile.bottom),
+                tile.right if right is None else max(right, tile.right),
+                tile.top if top is None else max(top, tile.top),
+            )
+
+    def intersection(self, other):
+        if isinstance(other, TileTask):
+            if other.tile.zoom + 1 != self._zoom:  # pragma: no cover
+                raise ValueError("intersecting tile has to be from zoom level above")
+            return [
+                self.tasks[child]
+                for child in other.tile.get_children()
+                if child in self.tasks
+            ]
+
+        if isinstance(other, Task):
+            return self.intersection(other.bounds)
+
+        elif isinstance(other, (tuple, list)):
+            return [
+                self.tasks[tile]
+                for tile in self._tp.tiles_from_bounds(bounds=other, zoom=self._zoom)
+                if tile in self.tasks
+            ]
+
+        else:
+            raise TypeError(
+                "intersections only works with other Task instances or bounds"
+            )
+
+    def _validate(self, items):
+        self._tp = None
+        self._zoom = None
+        for item in items:
+            if not isinstance(item, TileTask):  # pragma: no cover
+                raise TypeError(
+                    "TileTaskBatch items must be TileTasks, not %s", type(item)
+                )
+            if self._tp is None:
+                self._tp = item.tile.buffered_tp
+            elif item.tile.buffered_tp != self._tp:  # pragma: no cover
+                raise TypeError("all TileTasks must derive from the same pyramid.")
+            if self._zoom is None:
+                self._zoom = item.tile.zoom
+            elif item.tile.zoom != self._zoom:  # pragma: no cover
+                raise TypeError("all TileTasks must lie on the same zoom level")
+            yield item.tile, item
+
+
+def to_dask_collection(batches):
+    from dask.delayed import delayed
+
+    tasks = {}
+    with Timer() as t:
+        previous_batch = None
+        for batch in batches:
+            if previous_batch:
+                logger.debug("previous batch had %s tasks", len(previous_batch))
+            for task in batch.values():
+                if previous_batch:
+                    dependencies = {
+                        child.id: tasks[child]
+                        for child in previous_batch.intersection(task)
+                    }
+                    logger.debug(
+                        "found %s dependencies from last batch for task %s",
+                        len(dependencies),
+                        task,
+                    )
+                else:
+                    dependencies = {}
+                tasks[task] = delayed(batch.func)(
+                    task, dependencies=dependencies, **batch.fkwargs
+                )
+            previous_batch = batch
+    logger.debug("%s tasks generated in %s", len(tasks), t)
+    return list(tasks.values())
diff --git a/mapchete/_user_process.py b/mapchete/_user_process.py
new file mode 100644
index 00000000..abb2b022
--- /dev/null
+++ b/mapchete/_user_process.py
@@ -0,0 +1,174 @@
+from mapchete.commons import clip as commons_clip
+from mapchete.commons import contours as commons_contours
+from mapchete.commons import hillshade as commons_hillshade
+from mapchete.validate import deprecated_kwargs
+
+
+class MapcheteProcess(object):
+    """
+    Process class inherited by user process script.
+
+    Its attributes and methods can be accessed via "self" from within a
+    Mapchete process Python file.
+
+    Parameters
+    ----------
+    tile : BufferedTile
+        Tile process should be run on
+    config : MapcheteConfig
+        process configuration
+    params : dictionary
+        process parameters
+
+    Attributes
+    ----------
+    identifier : string
+        process identifier
+    title : string
+        process title
+    version : string
+        process version string
+    abstract : string
+        short text describing process purpose
+    tile : BufferedTile
+        Tile process should be run on
+    tile_pyramid : BufferedTilePyramid
+        process tile pyramid
+    output_pyramid : BufferedTilePyramid
+        output tile pyramid
+    params : dictionary
+        process parameters
+    """
+
+    def __init__(
+        self, tile=None, params=None, input=None, output_params=None, config=None
+    ):
+        """Initialize Mapchete process."""
+        self.identifier = ""
+        self.title = ""
+        self.version = ""
+        self.abstract = ""
+
+        self.tile = tile
+        self.tile_pyramid = tile.tile_pyramid
+        if config is not None:
+            input = config.get_inputs_for_tile(tile)
+            params = config.params_at_zoom(tile.zoom)
+        self.params = dict(params, input=input)
+        self.input = input
+        self.output_params = output_params
+
+    def write(self, data, **kwargs):
+        """Deprecated."""
+        raise DeprecationWarning(
+            "Please return process output data instead of using self.write()."
+        )
+
+    def read(self, **kwargs):
+        """
+        Read existing output data from a previous run.
+
+        Returns
+        -------
+        process output : NumPy array (raster) or feature iterator (vector)
+        """
+        raise DeprecationWarning(
+            "Read existing output from within a process is deprecated"
+        )
+
+    @deprecated_kwargs
+    def open(self, input_id, **kwargs):
+        """
+        Open input data.
+
+        Parameters
+        ----------
+        input_id : string
+            input identifier from configuration file or file path
+        kwargs : driver specific parameters (e.g. resampling)
+
+        Returns
+        -------
+        tiled input data : InputTile
+            reprojected input data within tile
+        """
+        if input_id not in self.input:
+            raise ValueError("%s not found in config as input" % input_id)
+        return self.input[input_id]
+
+    def hillshade(self, elevation, azimuth=315.0, altitude=45.0, z=1.0, scale=1.0):
+        """
+        Calculate hillshading from elevation data.
+
+        Parameters
+        ----------
+        elevation : array
+            input elevation data
+        azimuth : float
+            horizontal angle of light source (315: North-West)
+        altitude : float
+            vertical angle of light source (90 would result in slope shading)
+        z : float
+            vertical exaggeration factor
+        scale : float
+            scale factor of pixel size units versus height units (insert 112000
+            when having elevation values in meters in a geodetic projection)
+
+        Returns
+        -------
+        hillshade : array
+        """
+        return commons_hillshade.hillshade(
+            elevation, self.tile, azimuth, altitude, z, scale
+        )
+
+    def contours(self, elevation, interval=100, field="elev", base=0):
+        """
+        Extract contour lines from elevation data.
+
+        Parameters
+        ----------
+        elevation : array
+            input elevation data
+        interval : integer
+            elevation value interval when drawing contour lines
+        field : string
+            output field name containing elevation value
+        base : integer
+            elevation base value the intervals are computed from
+
+        Returns
+        -------
+        contours : iterable
+            contours as GeoJSON-like pairs of properties and geometry
+        """
+        return commons_contours.extract_contours(
+            elevation, self.tile, interval=interval, field=field, base=base
+        )
+
+    def clip(self, array, geometries, inverted=False, clip_buffer=0):
+        """
+        Clip array by geometry.
+
+        Parameters
+        ----------
+        array : array
+            raster data to be clipped
+        geometries : iterable
+            geometries used to clip source array
+        inverted : bool
+            invert clipping (default: False)
+        clip_buffer : int
+            buffer (in pixels) geometries before applying clip
+
+        Returns
+        -------
+        clipped array : array
+        """
+        return commons_clip.clip_array_with_vector(
+            array,
+            self.tile.affine,
+            geometries,
+            inverted=inverted,
+            clip_buffer=clip_buffer * self.tile.pixel_x_size,
+        )
diff --git a/mapchete/commands/_cp.py b/mapchete/commands/_cp.py
index 7206ad50..5de037ac 100644
--- a/mapchete/commands/_cp.py
+++ b/mapchete/commands/_cp.py
@@ -11,7 +11,7 @@
 from shapely.geometry.base import BaseGeometry
 
 import mapchete
-from mapchete.io import fs_from_path, tiles_exist
+from mapchete.io import fs_from_path, tiles_exist, copy
 from mapchete.io.vector import reproject_geometry
 
 logger = logging.getLogger(__name__)
@@ -136,7 +136,7 @@ def _empty_callback(*_):
             msg = f"copy {src_metadata} to {dst_metadata}"
             logger.debug(msg)
             msg_callback(msg)
-            _copy(src_fs, src_metadata, dst_fs, dst_metadata)
+            copy(src_metadata, dst_metadata, src_fs=src_fs, dst_fs=dst_fs)
 
         with mapchete.open(
             dst_tiledir,
@@ -252,7 +252,7 @@ def _copy_tile(
 
         # copy from source to target
         dst_path = dst_mp.config.output_reader.get_path(tile)
-        _copy(src_fs, src_path, dst_fs, dst_path)
+        copy(src_path, dst_path, src_fs=src_fs, dst_fs=dst_fs)
         msg = f"{tile}: copy {src_path} to {dst_path}"
         logger.debug(msg)
         return 1, msg
@@ -260,17 +260,3 @@ def _copy_tile(
     msg = f"{tile}: source tile ({src_path}) does not exist"
     logger.debug(msg)
     return 0, msg
-
-
-def _copy(src_fs, src_path, dst_fs, dst_path):
-    # create parent directories on local filesystems
-    if dst_fs.protocol == "file":
-        dst_fs.makedirs(os.path.dirname(dst_path), exist_ok=True)
-
-    # copy either within a filesystem or between filesystems
-    if src_fs == dst_fs:
-        src_fs.copy(src_path, dst_path)
-    else:
-        with src_fs.open(src_path, "rb") as src:
-            with dst_fs.open(dst_path, "wb") as dst:
-                dst.write(src.read())
diff --git a/mapchete/commands/_execute.py b/mapchete/commands/_execute.py
index 216c527d..60c364af 100755
--- a/mapchete/commands/_execute.py
+++ b/mapchete/commands/_execute.py
@@ -2,6 +2,7 @@
 
 import logging
 from multiprocessing import cpu_count
+import traceback
 from typing import Callable, List, Tuple, Union
 import warnings
 
@@ -10,6 +11,7 @@
 
 import mapchete
 from mapchete.config import bounds_from_opts, raw_conf, raw_conf_process_pyramid
+from mapchete._processing import PreprocessingProcessInfo, TileProcessInfo
 
 logger = logging.getLogger(__name__)
 
@@ -34,6 +36,7 @@ def execute(
     dask_max_submitted_tasks=1000,
     dask_chunksize=100,
     dask_client=None,
+    dask_compute_graph=True,
     msg_callback: Callable = None,
     as_iterator: bool = False,
 ) -> mapchete.Job:
@@ -80,6 +83,8 @@ def execute(
         Number of tasks submitted to the scheduler at once. (default: 100)
     dask_client : dask.distributed.Client
         Reusable Client instance if required. Otherwise a new client will be created.
+    dask_compute_graph : bool
+        Build and compute dask graph instead of submitting tasks as preprocessing & zoom tiles batches. (default: True)
     msg_callback : Callable
         Optional callback function for process messages.
     as_iterator : bool
@@ -167,6 +172,7 @@ def _empty_callback(_):
                 zoom=None if tile else zoom,
                 dask_max_submitted_tasks=dask_max_submitted_tasks,
                 dask_chunksize=dask_chunksize,
+                dask_compute_graph=dask_compute_graph,
             ),
             executor_concurrency=concurrency,
             executor_kwargs=dict(
@@ -179,8 +185,8 @@ def _empty_callback(_):
             tiles_tasks=tiles_tasks,
         )
     # explicitly exit the mp object on failure
-    except Exception:  # pragma: no cover
-        mp.__exit__(None, None, None)
+    except Exception as exc:  # pragma: no cover
+        mp.__exit__(exc, repr(exc), traceback.format_exc())
         raise
 
 
@@ -194,25 +200,25 @@ def _process_everything(
     **kwargs,
 ):
     try:
-        for preprocessing_task_info in mp.batch_preprocessor(
-            executor=executor,
-            workers=workers,
-            dask_max_submitted_tasks=dask_max_submitted_tasks,
-            dask_chunksize=dask_chunksize,
-        ):  # pragma: no cover
-            yield preprocessing_task_info
-            msg_callback(preprocessing_task_info)
-        for process_info in mp.batch_processor(
+        for future in mp.compute(
             executor=executor,
             workers=workers,
             dask_max_submitted_tasks=dask_max_submitted_tasks,
             dask_chunksize=dask_chunksize,
             **kwargs,
         ):
-            yield process_info
-            msg_callback(
-                f"Tile {process_info.tile.id}: {process_info.process_msg}, {process_info.write_msg}"
-            )
-    # explicitly exit the mp object on success
-    finally:
+            process_info = future.result()
+            if isinstance(process_info, PreprocessingProcessInfo):
+                msg_callback(f"Task {process_info.task_key} finished")
+            elif isinstance(process_info, TileProcessInfo):
+                msg_callback(
+                    f"Task {process_info.tile.id}: {process_info.process_msg}, {process_info.write_msg}"
+                )
+            else:  # pragma: no cover
+                raise TypeError(f"unknown process info type: {type(process_info)}")
+            yield future
+        # explicitly exit the mp object on success
         mp.__exit__(None, None, None)
+    except Exception as exc:  # pragma: no cover
+        mp.__exit__(exc, repr(exc), traceback.format_exc())
+        raise
diff --git a/mapchete/config.py b/mapchete/config.py
index b69727d4..3e6d04af 100644
--- a/mapchete/config.py
+++ b/mapchete/config.py
@@ -341,25 +341,42 @@ def preprocessing_tasks_per_input(self):
         }
 
     def preprocessing_tasks(self):
+        """Get mapping of all preprocessing tasks."""
         return {
             task_key: task
-            for preprocessing_tasks in self.preprocessing_tasks_per_input().values()
-            for task_key, task in preprocessing_tasks.items()
+            for _, inp_preprocessing_tasks in self.preprocessing_tasks_per_input().items()
+            for task_key, task in inp_preprocessing_tasks.items()
         }
 
     def preprocessing_tasks_count(self):
-        """Return number of unique preprocessing tasks."""
+        """Return number of preprocessing tasks."""
         return len(self.preprocessing_tasks())
 
+    def preprocessing_task_finished(self, task_key):
+        """Return True if task of given key has already been run."""
+        inp_key, task_key = task_key.split(":")[0], ":".join(task_key.split(":")[1:])
+        try:
+            inp = self.input[inp_key]
+        except KeyError:  # pragma: no cover
+            raise KeyError(f"input {inp_key} not found")
+        return inp.preprocessing_task_finished(task_key)
+
     def set_preprocessing_task_result(self, task_key, result):
         """Append preprocessing task result to input."""
-        found = False
+        if ":" in task_key:
+            inp_key = task_key.split(":")[0]
+        else:
+            raise KeyError(
+                f"preprocessing task cannot be assigned to an input: {task_key}"
+            )
         for inp in self.input.values():
-            if task_key in inp.preprocessing_tasks:
-                found = True
-                inp.preprocessing_tasks_results[task_key] = result
-        if not found:
-            raise KeyError(f"task key {task_key} not found in any input")
+            if inp_key == inp.input_key:
+                break
+        else:  # pragma: no cover
+            raise KeyError(
+                f"task {task_key} cannot be assigned to input with key {inp_key}"
+            )
+        inp.set_preprocessing_task_result(task_key, result)
 
     @cached_property
     def zoom_levels(self):
@@ -476,7 +493,6 @@ def input(self):
                 if v is not None
             ]
         )
-
         if self._init_inputs:
             return initialize_inputs(
                 raw_inputs,
@@ -1025,6 +1041,7 @@ def initialize_inputs(
                         delimiters=delimiters,
                     ),
                     readonly=readonly,
+                    input_key=k,
                 )
             except Exception as e:
                 logger.exception(e)
@@ -1035,15 +1052,21 @@ def initialize_inputs(
         elif isinstance(v, dict):
             logger.debug("load input reader for abstract input %s", v)
             try:
+                abstract = deepcopy(v)
+                if "path" in abstract:
+                    abstract.update(
+                        path=absolute_path(path=abstract["path"], base_dir=config_dir)
+                    )
                 reader = load_input_reader(
                     dict(
-                        abstract=deepcopy(v),
+                        abstract=abstract,
                         pyramid=pyramid,
                         pixelbuffer=pyramid.pixelbuffer,
                         delimiters=delimiters,
                         conf_dir=config_dir,
                     ),
                     readonly=readonly,
+                    input_key=k,
                 )
             except Exception as e:
                 logger.exception(e)
@@ -1055,6 +1078,10 @@ def initialize_inputs(
         reader.bbox(out_crs=pyramid.crs)
         initalized_inputs[k] = reader
 
+    logger.debug(
+        "initialized inputs: %s",
+        initalized_inputs.keys(),
+    )
     return initalized_inputs
 
 
diff --git a/mapchete/errors.py b/mapchete/errors.py
index e6b307e7..03d45aa6 100644
--- a/mapchete/errors.py
+++ b/mapchete/errors.py
@@ -43,3 +43,7 @@ class MapcheteIOError(IOError):
 
 class JobCancelledError(Exception):
     """Raised when Job gets cancelled."""
+
+
+class NoTaskGeometry(TypeError):
+    """Raised when Task has no assigned geo information."""
diff --git a/mapchete/formats/__init__.py b/mapchete/formats/__init__.py
index 6711815e..aa5b0c26 100644
--- a/mapchete/formats/__init__.py
+++ b/mapchete/formats/__init__.py
@@ -98,7 +98,7 @@ def load_output_writer(output_params, readonly=False):
     raise MapcheteDriverError("no loader for driver '%s' could be found." % driver_name)
 
 
-def load_input_reader(input_params, readonly=False):
+def load_input_reader(input_params, readonly=False, input_key=None):
     """
     Return input class of driver.
 
@@ -126,7 +126,9 @@ def load_input_reader(input_params, readonly=False):
         if hasattr(driver_, "METADATA") and (
             driver_.METADATA["driver_name"] == driver_name
         ):
-            return driver_.InputData(input_params, readonly=readonly)
+            return driver_.InputData(
+                input_params, readonly=readonly, input_key=input_key
+            )
     raise MapcheteDriverError("no loader for driver '%s' could be found." % driver_name)
 
 
diff --git a/mapchete/formats/base.py b/mapchete/formats/base.py
index 6dff6257..4d341868 100644
--- a/mapchete/formats/base.py
+++ b/mapchete/formats/base.py
@@ -24,6 +24,7 @@
     prepare_array,
     read_raster_window,
 )
+from mapchete._tasks import Task
 from mapchete.io.vector import read_vector_window
 from mapchete.tile import BufferedTilePyramid
 
@@ -51,8 +52,9 @@ class InputData(object):
 
     METADATA = {"driver_name": None, "data_type": None, "mode": "r"}
 
-    def __init__(self, input_params, **kwargs):
+    def __init__(self, input_params, input_key=None, **kwargs):
         """Initialize relevant input information."""
+        self.input_key = input_key
         self.pyramid = input_params.get("pyramid")
         self.pixelbuffer = input_params.get("pixelbuffer")
         self.crs = self.pyramid.crs if self.pyramid else None
@@ -106,7 +108,9 @@ def cleanup(self):
         """Optional cleanup function called when Mapchete exits."""
         pass
 
-    def add_preprocessing_task(self, func, fargs=None, fkwargs=None, key=None):
+    def add_preprocessing_task(
+        self, func, fargs=None, fkwargs=None, key=None, geometry=None, bounds=None
+    ):
         """
         Add longer running preprocessing function to be called right before processing.
 
@@ -118,18 +122,56 @@ def add_preprocessing_task(self, func, fargs=None, fkwargs=None, key=None):
             fargs = (fargs,)
         fkwargs = fkwargs or {}
         key = f"{func}-{get_hash((func, fargs, fkwargs))}" if key is None else key
+        if self.input_key:
+            key = f"{self.input_key}:{key}"
         if key in self.preprocessing_tasks:  # pragma: no cover
             raise KeyError(f"preprocessing task with key {key} already exists")
         logger.debug(f"add preprocessing task {key, func}")
-        self.preprocessing_tasks[key] = (func, fargs, fkwargs)
+        self.preprocessing_tasks[key] = Task(
+            id=key,
+            func=func,
+            fargs=fargs,
+            fkwargs=fkwargs,
+            geometry=geometry,
+            bounds=bounds,
+        )
 
     def get_preprocessing_task_result(self, task_key):
+        """
+        Get result of preprocessing task.
+        """
+        if self.input_key and not task_key.startswith(f"{self.input_key}:"):
+            task_key = f"{self.input_key}:{task_key}"
         if task_key not in self.preprocessing_tasks:
             raise KeyError(f"task {task_key} is not a task for current input")
         if task_key not in self.preprocessing_tasks_results:
             raise ValueError(f"task {task_key} has not yet been executed")
         return self.preprocessing_tasks_results[task_key]
 
+    def set_preprocessing_task_result(self, task_key, result):
+        """
+        Set result of preprocessing task.
+        """
+        if self.input_key and not task_key.startswith(
+            f"{self.input_key}:"
+        ):  # pragma: no cover
+            task_key = f"{self.input_key}:{task_key}"
+        if task_key not in self.preprocessing_tasks:  # pragma: no cover
+            raise KeyError(f"task {task_key} is not a task for current input")
+        if task_key in self.preprocessing_tasks_results:  # pragma: no cover
+            raise KeyError(f"task {task_key} has already been set")
+        self.preprocessing_tasks_results[task_key] = result
+
+    def preprocessing_task_finished(self, task_key):
+        """
+        Return whether preprocessing task already ran.
+        """
+        if self.input_key and not task_key.startswith(f"{self.input_key}:"):
+            task_key = f"{self.input_key}:{task_key}"
+        if task_key not in self.preprocessing_tasks:  # pragma: no cover
+            raise KeyError(f"task {task_key} is not a task for current input")
+        return task_key in self.preprocessing_tasks_results
+
 
 class InputTile(object):
     """
@@ -142,6 +184,9 @@ class InputTile(object):
         driver specific parameters
     """
 
+    preprocessing_tasks_results = {}
+    input_key = None
+
     def __init__(self, tile, **kwargs):
         """Initialize."""
 
@@ -166,6 +211,12 @@ def is_empty(self):
         """
         raise NotImplementedError
 
+    def set_preprocessing_task_result(self, task_key=None, result=None):
+        """
+        Adds a preprocessing task result.
+        """
+        self.preprocessing_tasks_results[task_key] = result
+
     def __enter__(self):
         """Required for 'with' statement."""
         return self
@@ -405,7 +456,12 @@ def output_cleaned(self, process_data):
         """
         if self.METADATA["data_type"] == "raster":
             if is_numpy_or_masked_array(process_data):
-                return process_data
+                return prepare_array(
+                    process_data,
+                    masked=True,
+                    nodata=self.output_params["nodata"],
+                    dtype=self.profile()["dtype"],
+                )
             elif is_numpy_or_masked_array_with_tags(process_data):
                 data, tags = process_data
                 return self.output_cleaned(data), tags
@@ -415,7 +471,7 @@ def output_cleaned(self, process_data):
     def streamline_output(self, process_data):
         if isinstance(process_data, str) and process_data == "empty":
             raise MapcheteNodataTile
-        elif process_data is None:
+        elif process_data is None:  # pragma: no cover
             raise MapcheteProcessOutputError("process output is empty")
         elif self.output_is_valid(process_data):
             return self.output_cleaned(process_data)
diff --git a/mapchete/formats/default/gtiff.py b/mapchete/formats/default/gtiff.py
index 295598df..953b207e 100644
--- a/mapchete/formats/default/gtiff.py
+++ b/mapchete/formats/default/gtiff.py
@@ -614,7 +614,7 @@ def close(self, exc_type=None, exc_value=None, exc_traceback=None):
         """Build overviews and write file."""
         try:
             # only in case no Exception was raised
-            if not exc_type:
+            if exc_type is None:
                 # build overviews
                 if self.overviews and self.dst is not None:
                     logger.debug(
diff --git a/mapchete/formats/default/raster_file.py b/mapchete/formats/default/raster_file.py
index 2ffca02b..346d8d43 100644
--- a/mapchete/formats/default/raster_file.py
+++ b/mapchete/formats/default/raster_file.py
@@ -8,6 +8,7 @@
 from cached_property import cached_property
 from copy import deepcopy
 import logging
+import numpy.ma as ma
 import os
 import rasterio
 from rasterio.vrt import WarpedVRT
@@ -16,7 +17,12 @@
 
 from mapchete.formats import base
 from mapchete.io.vector import reproject_geometry, segmentize_geometry
-from mapchete.io.raster import read_raster_window
+from mapchete.io.raster import (
+    read_raster_window,
+    convert_raster,
+    read_raster,
+    resample_from_array,
+)
 from mapchete import io
 
 
@@ -61,11 +67,55 @@ class InputData(base.InputData):
         "mode": "r",
         "file_extensions": ["tif", "vrt", "png", "jp2"],
     }
+    _cached_path = None
+    _cache_keep = False
+    _memory_cache_active = False
 
     def __init__(self, input_params, **kwargs):
         """Initialize."""
         super().__init__(input_params, **kwargs)
-        self.path = input_params["path"]
+        self.path = (
+            input_params["abstract"]["path"]
+            if "abstract" in input_params
+            else input_params["path"]
+        )
+        self._cache_task = f"cache_{self.path}"
+        if "abstract" in input_params and "cache" in input_params["abstract"]:
+            if isinstance(input_params["abstract"]["cache"], dict):
+                if "path" in input_params["abstract"]["cache"]:
+                    self._cached_path = io.absolute_path(
+                        path=input_params["abstract"]["cache"]["path"],
+                        base_dir=input_params["conf_dir"],
+                    )
+                else:  # pragma: no cover
+                    raise ValueError("please provide a cache path")
+                # add preprocessing task to cache data
+                self.add_preprocessing_task(
+                    convert_raster,
+                    key=self._cache_task,
+                    fkwargs=dict(
+                        inp=self.path,
+                        out=self._cached_path,
+                        format=input_params["abstract"]["cache"].get("format", "COG"),
+                    ),
+                    geometry=self.bbox(),
+                )
+                self._cache_keep = input_params["abstract"]["cache"].get("keep", False)
+            elif (
+                isinstance(input_params["abstract"]["cache"], str)
+                and input_params["abstract"]["cache"] == "memory"
+            ):
+                self._memory_cache_active = True
+                self.add_preprocessing_task(
+                    read_raster,
+                    key=self._cache_task,
+                    fkwargs=dict(inp=self.path),
+                    geometry=self.bbox(),
+                )
+            else:  # pragma: no cover
+                raise ValueError(
+                    f"invalid cache configuration given: {input_params['abstract']['cache']}"
+                )
 
     @cached_property
     def profile(self):
@@ -86,7 +136,19 @@ def open(self, tile, **kwargs):
         input tile : ``InputTile``
             tile view of input data
         """
-        return InputTile(tile, self, **kwargs)
+        if self._memory_cache_active and self.preprocessing_task_finished(
+            self._cache_task
+        ):
+            in_memory_raster = self.get_preprocessing_task_result(self._cache_task)
+        else:
+            in_memory_raster = None
+        return InputTile(
+            tile,
+            self,
+            in_memory_raster=in_memory_raster,
+            cache_task_key=self._cache_task,
+            **kwargs,
+        )
 
     def bbox(self, out_crs=None):
         """
@@ -136,6 +198,15 @@ def exists(self):
         """
         return os.path.isfile(self.path)  # pragma: no cover
 
+    def cleanup(self):
+        """Cleanup when mapchete closes."""
+        if self._cached_path and not self._cache_keep:
+            logger.debug("remove cached file %s", self._cached_path)
+            try:
+                io.fs_from_path(self._cached_path).rm(self._cached_path)
+            except FileNotFoundError:
+                pass
+
 
 class InputTile(base.InputTile):
     """
@@ -150,24 +221,43 @@ class InputTile(base.InputTile):
     Attributes
     ----------
     tile : tile : ``Tile``
-    raster_file : ``InputData``
+    input_data : ``InputData``
         parent InputData object
     resampling : string
         resampling method passed on to rasterio
     """
 
-    def __init__(self, tile, raster_file, **kwargs):
+    _memory_cache_active = False
+    _in_memory_raster = None
+
+    def __init__(
+        self, tile, input_data, in_memory_raster=None, cache_task_key=None, **kwargs
+    ):
         """Initialize."""
         self.tile = tile
-        self.raster_file = raster_file
-        if io.path_is_remote(raster_file.path):
-            file_ext = os.path.splitext(raster_file.path)[1]
-            self.gdal_opts = {
-                "GDAL_DISABLE_READDIR_ON_OPEN": True,
-                "CPL_VSIL_CURL_ALLOWED_EXTENSIONS": "%s,.ovr" % file_ext,
-            }
+        self.bbox = input_data.bbox(out_crs=self.tile.crs)
+        self.profile = input_data.profile
+        self.cache_task_key = cache_task_key
+        self.input_key = input_data.input_key
+        if input_data._memory_cache_active:
+            self._memory_cache_active = True
+            self._in_memory_raster = in_memory_raster
         else:
-            self.gdal_opts = {}
+            self.path = input_data._cached_path or input_data.path
+            if io.path_is_remote(input_data.path):
+                file_ext = os.path.splitext(self.path)[1]
+                self.gdal_opts = {
+                    "GDAL_DISABLE_READDIR_ON_OPEN": True,
+                    "CPL_VSIL_CURL_ALLOWED_EXTENSIONS": "%s,.ovr" % file_ext,
+                }
+            else:
+                self.gdal_opts = {}
+
+    def __repr__(self):  # pragma: no cover
+        source = (
+            repr(self._in_memory_raster) if self._memory_cache_active else self.path
+        )
+        return f"raster_file.InputTile(tile={self.tile.id}, source={source})"
 
     def read(self, indexes=None, resampling="nearest", **kwargs):
         """
@@ -177,13 +267,33 @@ def read(self, indexes=None, resampling="nearest", **kwargs):
         -------
         data : array
         """
-        return read_raster_window(
-            self.raster_file.path,
-            self.tile,
-            indexes=self._get_band_indexes(indexes),
-            resampling=resampling,
-            gdal_opts=self.gdal_opts,
-        )
+        if self._memory_cache_active:
+            self._in_memory_raster = (
+                self._in_memory_raster
+                or self.preprocessing_tasks_results.get(self.cache_task_key)
+            )
+            if self._in_memory_raster is None:  # pragma: no cover
+                raise RuntimeError("preprocessing tasks have not yet been run")
+            return resample_from_array(
+                in_raster=ma.stack(
+                    [
+                        self._in_memory_raster.data[i - 1]
+                        for i in self._get_band_indexes(indexes)
+                    ]
+                ),
+                in_affine=self._in_memory_raster.affine,
+                in_crs=self._in_memory_raster.crs,
+                out_tile=self.tile,
+                resampling=resampling,
+            )
+        else:
+            return read_raster_window(
+                self.path,
+                self.tile,
+                indexes=self._get_band_indexes(indexes),
+                resampling=resampling,
+                gdal_opts=self.gdal_opts,
+            )
 
     def is_empty(self, indexes=None):
         """
@@ -194,9 +304,7 @@ def is_empty(self, indexes=None):
         is empty : bool
         """
         # empty if tile does not intersect with file bounding box
-        return not self.tile.bbox.intersects(
-            self.raster_file.bbox(out_crs=self.tile.crs)
-        )
+        return not self.tile.bbox.intersects(self.bbox)
 
     def _get_band_indexes(self, indexes=None):
         """Return valid band indexes."""
@@ -206,7 +314,7 @@ def _get_band_indexes(self, indexes=None):
             else:
                 return [indexes]
         else:
-            return range(1, self.raster_file.profile["count"] + 1)
+            return range(1, self.profile["count"] + 1)
 
 
 def get_segmentize_value(input_file=None, tile_pyramid=None):
diff --git a/mapchete/formats/default/vector_file.py b/mapchete/formats/default/vector_file.py
index f76c2d9c..f29ace53 100644
--- a/mapchete/formats/default/vector_file.py
+++ b/mapchete/formats/default/vector_file.py
@@ -4,12 +4,24 @@
 Currently limited by extensions .shp and .geojson but could be extended easily.
 """
 
+from cached_property import cached_property
 import fiona
+import logging
 from shapely.geometry import box
 from rasterio.crs import CRS
 
 from mapchete.formats import base
-from mapchete.io.vector import reproject_geometry, read_vector_window
+from mapchete.io.vector import (
+    reproject_geometry,
+    read_vector_window,
+    convert_vector,
+    read_vector,
+    IndexedFeatures,
+)
+from mapchete.io import fs_from_path, absolute_path
+
+
+logger = logging.getLogger(__name__)
 
 
 METADATA = {
@@ -49,11 +61,62 @@ class InputData(base.InputData):
         "mode": "r",
         "file_extensions": ["shp", "geojson"],
     }
+    _cached_path = None
+    _cache_keep = False
+    _memory_cache_active = False
 
     def __init__(self, input_params, **kwargs):
         """Initialize."""
         super().__init__(input_params, **kwargs)
-        self.path = input_params["path"]
+        self.path = (
+            input_params["abstract"]["path"]
+            if "abstract" in input_params
+            else input_params["path"]
+        )
+        self._cache_task = f"cache_{self.path}"
+        if "abstract" in input_params and "cache" in input_params["abstract"]:
+            if isinstance(input_params["abstract"]["cache"], dict):
+                if "path" in input_params["abstract"]["cache"]:
+                    self._cached_path = absolute_path(
+                        path=input_params["abstract"]["cache"]["path"],
+                        base_dir=input_params["conf_dir"],
+                    )
+                else:  # pragma: no cover
+                    raise ValueError("please provide a cache path")
+                # add preprocessing task to cache data
+                self.add_preprocessing_task(
+                    convert_vector,
+                    key=f"cache_{self.path}",
+                    fkwargs=dict(
+                        inp=self.path,
+                        out=self._cached_path,
+                        format=input_params["abstract"]["cache"].get(
+                            "format", "FlatGeobuf"
+                        ),
+                    ),
+                    geometry=self.bbox(),
+                )
+                self._cache_keep = input_params["abstract"]["cache"].get("keep", False)
+            elif (
+                isinstance(input_params["abstract"]["cache"], str)
+                and input_params["abstract"]["cache"] == "memory"
+            ):
+                self._memory_cache_active = True
+                self.add_preprocessing_task(
+                    read_vector,
+                    key=f"cache_{self.path}",
+                    fkwargs=dict(inp=self.path, index=None),
+                    geometry=self.bbox(),
+                )
+            else:  # pragma: no cover
+                raise ValueError(
+                    f"invalid cache configuration given: {input_params['abstract']['cache']}"
+                )
+
+    @cached_property
+    def in_memory_features(self):
+        """This property can be accessed once the preprocessing task is finished."""
+        return IndexedFeatures(self.get_preprocessing_task_result(f"cache_{self.path}"))
 
     def open(self, tile, **kwargs):
         """
@@ -68,7 +131,29 @@ def open(self, tile, **kwargs):
         input tile : ``InputTile``
             tile view of input data
         """
-        return InputTile(tile, self, **kwargs)
+        if self._memory_cache_active and self.preprocessing_task_finished(
+            self._cache_task
+        ):
+            tile_features = IndexedFeatures(
+                self.in_memory_features.filter(
+                    reproject_geometry(
+                        tile.bbox, src_crs=tile.crs, dst_crs=self.in_memory_features.crs
+                    ).bounds
+                ),
+                crs=self.in_memory_features.crs,
+                index=None,
+            )
+
+        else:
+            tile_features = None
+
+        return InputTile(
+            tile,
+            self,
+            in_memory_features=tile_features,
+            cache_task_key=self._cache_task,
+            **kwargs,
+        )
 
     def bbox(self, out_crs=None):
         """
@@ -93,6 +178,15 @@ def bbox(self, out_crs=None):
             bbox, src_crs=inp_crs, dst_crs=out_crs, clip_to_crs_bounds=False
         )
 
+    def cleanup(self):
+        """Cleanup when mapchete closes."""
+        if self._cached_path and not self._cache_keep:
+            logger.debug("remove cached file %s", self._cached_path)
+            try:
+                fs_from_path(self._cached_path).rm(self._cached_path)
+            except FileNotFoundError:
+                pass
+
 
 class InputTile(base.InputTile):
     """
@@ -107,15 +201,33 @@ class InputTile(base.InputTile):
     Attributes
     ----------
     tile : tile : ``Tile``
-    vector_file : string
+    input_data : string
         path to input vector file
     """
 
-    def __init__(self, tile, vector_file, **kwargs):
+    _memory_cache_active = False
+    _in_memory_features = None
+
+    def __init__(
+        self, tile, input_data, in_memory_features=None, cache_task_key=None, **kwargs
+    ):
         """Initialize."""
         self.tile = tile
-        self.vector_file = vector_file
         self._cache = {}
+        self.bbox = input_data.bbox(out_crs=self.tile.crs)
+        self.cache_task_key = cache_task_key
+        self.input_key = input_data.input_key
+        if input_data._memory_cache_active:
+            self._memory_cache_active = True
+            self._in_memory_features = in_memory_features
+        else:
+            self.path = input_data._cached_path or input_data.path
+
+    def __repr__(self):  # pragma: no cover
+        source = (
+            repr(self._in_memory_features) if self._memory_cache_active else self.path
+        )
+        return f"vector_file.InputTile(tile={self.tile.id}, source={source})"
 
     def read(self, validity_check=True, clip_to_crs_bounds=False, **kwargs):
         """
@@ -133,6 +245,13 @@ def read(self, validity_check=True, clip_to_crs_bounds=False, **kwargs):
         -------
         data : list
         """
+        if self._memory_cache_active:
+            self._in_memory_features = (
+                self._in_memory_features
+                or self.preprocessing_tasks_results.get(self.cache_task_key)
+            )
+            if self._in_memory_features is None:  # pragma: no cover
+                raise RuntimeError("preprocessing tasks have not yet been run")
         return (
             []
             if self.is_empty()
@@ -149,7 +268,7 @@ def is_empty(self):
         -------
         is empty : bool
         """
-        if not self.tile.bbox.intersects(self.vector_file.bbox()):
+        if not self.tile.bbox.intersects(self.bbox):
             return True
         return len(self._read_from_cache(True)) == 0
 
@@ -158,7 +277,9 @@ def _read_from_cache(self, validity_check=True, clip_to_crs_bounds=False):
         if checked not in self._cache:
             self._cache[checked] = list(
                 read_vector_window(
-                    self.vector_file.path,
+                    self._in_memory_features
+                    if self._memory_cache_active
+                    else self.path,
                     self.tile,
                     validity_check=validity_check,
                     clip_to_crs_bounds=clip_to_crs_bounds,
diff --git a/mapchete/io/__init__.py b/mapchete/io/__init__.py
index 8442ad35..0e31b811 100644
--- a/mapchete/io/__init__.py
+++ b/mapchete/io/__init__.py
@@ -13,6 +13,7 @@
     get_gdal_options,
 )
 from mapchete.io._path import (
+    copy,
     fs_from_path,
     path_is_remote,
     path_exists,
@@ -24,6 +25,7 @@
 
 
 __all__ = [
+    "copy",
     "fs_from_path",
     "GDAL_HTTP_OPTS",
     "get_best_zoom_level",
diff --git a/mapchete/io/_path.py b/mapchete/io/_path.py
index c0e2ee0a..0f5f8933 100644
--- a/mapchete/io/_path.py
+++ b/mapchete/io/_path.py
@@ -97,11 +97,10 @@ def makedirs(path):
     ----------
     path : path
     """
-    if not path_is_remote(path):
-        try:
-            os.makedirs(path)
-        except OSError:
-            pass
+    fs = fs_from_path(path)
+    # create parent directories on local filesystems
+    if fs.protocol == "file":
+        fs.makedirs(path, exist_ok=True)
 
 
 def tiles_exist(
@@ -110,7 +109,7 @@ def tiles_exist(
     output_tiles_batches=None,
     process_tiles=None,
     process_tiles_batches=None,
-    **kwargs
+    **kwargs,
 ):
     """
     Yield tiles and whether their output already exists or not.
@@ -294,3 +293,23 @@ def fs_from_path(path, timeout=5, session=None, username=None, password=None, **
         return fsspec.filesystem("https", auth=auth, asynchronous=False, **kwargs)
     else:
         return fsspec.filesystem("file", asynchronous=False, **kwargs)
+
+
+def copy(src_path, dst_path, src_fs=None, dst_fs=None, overwrite=False):
+    """Copy path from one place to the other."""
+    src_fs = src_fs or fs_from_path(src_path)
+    dst_fs = dst_fs or fs_from_path(dst_path)
+
+    if not overwrite and dst_fs.exists(dst_path):
+        raise IOError(f"{dst_path} already exists")
+
+    # create parent directories on local filesystems
+    makedirs(os.path.dirname(dst_path))
+
+    # copy either within a filesystem or between filesystems
+    if src_fs == dst_fs:
+        src_fs.copy(src_path, dst_path)
+    else:
+        with src_fs.open(src_path, "rb") as src:
+            with dst_fs.open(dst_path, "wb") as dst:
+                dst.write(src.read())
diff --git a/mapchete/io/raster.py b/mapchete/io/raster.py
index b28f3d89..5ab37bc0 100644
--- a/mapchete/io/raster.py
+++ b/mapchete/io/raster.py
@@ -22,7 +22,14 @@
 import warnings
 
 from mapchete.errors import MapcheteIOError
-from mapchete.io import path_is_remote, get_gdal_options, path_exists, fs_from_path
+from mapchete.io import (
+    path_is_remote,
+    get_gdal_options,
+    path_exists,
+    fs_from_path,
+    copy,
+    makedirs,
+)
 from mapchete.io._misc import MAPCHETE_IO_RETRY_SETTINGS
 from mapchete.tile import BufferedTile
 from mapchete.validate import validate_write_window_params
@@ -651,7 +658,6 @@ def resample_from_array(
     if nodataval is not None:  # pragma: no cover
         warnings.warn("'nodataval' is deprecated, please use 'nodata'")
         nodata = nodata or nodataval
-    # TODO rename function
     if isinstance(in_raster, ma.MaskedArray):
         pass
     elif isinstance(in_raster, np.ndarray):
@@ -889,7 +895,9 @@ def _get_tiles_properties(tiles):
             raise ValueError("all tiles must have the same CRS")
         if isinstance(data, np.ndarray):
             if data[0].dtype != tiles[0][1][0].dtype:
-                raise TypeError("all tile data must have the same dtype")
+                raise TypeError(
+                    f"all tile data must have the same dtype: {data[0].dtype} != {tiles[0][1][0].dtype}"
+                )
     return tile.tile_pyramid, tile.pixel_x_size, data[0].dtype
 
 
@@ -1037,3 +1045,52 @@ def _prepare_masked(data, masked=True, nodata=0, dtype=None):
         return ma.masked_values(data.astype(dtype, copy=False), nodata, copy=False)
     else:
         return ma.filled(data.astype(dtype, copy=False), nodata)
+
+
+def convert_raster(inp, out, overwrite=False, exists_ok=True, **kwargs):
+    """
+    Convert raster file to a differernt format.
+
+    When kwargs are given, the operation will be conducted by rasterio, without kwargs,
+    the file is simply copied to the destination using fsspec.
+
+    Parameters
+    ----------
+    inp : str
+        Path to input file.
+    out : str
+        Path to output file.
+    overwrite : bool
+        Overwrite output file. (default: False)
+    skip_exists : bool
+        Skip conversion if outpu already exists. (default: True)
+    kwargs : mapping
+        Creation parameters passed on to output file.
+    """
+    if path_exists(out):
+        if not exists_ok:
+            raise IOError(f"{out} already exists")
+        elif not overwrite:
+            logger.debug("output %s already exists and will not be overwritten")
+            return
+    kwargs = kwargs or {}
+    if kwargs:
+        logger.debug("convert raster file %s to %s using %s", inp, out, kwargs)
+        with rasterio.open(inp, "r") as src:
+            makedirs(os.path.dirname(out))
+            with rasterio_write(out, mode="w", **{**src.meta, **kwargs}) as dst:
+                dst.write(src.read())
+    else:
+        logger.debug("copy %s to %s", inp, out)
+        copy(inp, out, overwrite=overwrite)
+
+
+def read_raster(inp, **kwargs):
+    logger.debug("reading {inp} into memory")
+    with rasterio.open(inp, "r") as src:
+        return ReferencedRaster(
+            data=src.read(masked=True),
+            affine=src.transform,
+            bounds=src.bounds,
+            crs=src.crs,
+        )
diff --git a/mapchete/io/vector.py b/mapchete/io/vector.py
index dafeacac..f9dd5f42 100644
--- a/mapchete/io/vector.py
+++ b/mapchete/io/vector.py
@@ -1,20 +1,23 @@
 """Functions handling vector data."""
 
+from contextlib import ExitStack
 import os
 import logging
 import fiona
 from fiona.errors import DriverError, FionaError, FionaValueError
 from fiona.io import MemoryFile
+import json
 from retry import retry
 from rasterio.crs import CRS
-from shapely.geometry import box, mapping
+from shapely.geometry import box, mapping, shape
 from shapely.errors import TopologicalError
 from tilematrix import clip_geometry_to_srs_bounds
 from itertools import chain
+import warnings
 
 from mapchete.errors import GeometryTypeError, MapcheteIOError
 from mapchete.io._misc import MAPCHETE_IO_RETRY_SETTINGS
-from mapchete.io._path import fs_from_path, path_exists
+from mapchete.io._path import fs_from_path, path_exists, makedirs, copy
 from mapchete.io._geometry_operations import (
     reproject_geometry,
     segmentize_geometry,
@@ -23,6 +26,7 @@
     clean_geometry_type,
     _repair,
 )
+from mapchete.validate import validate_bounds
 
 __all__ = [
     "reproject_geometry",
@@ -35,9 +39,7 @@
 logger = logging.getLogger(__name__)
 
 
-def read_vector_window(
-    input_files, tile, validity_check=True, clip_to_crs_bounds=False
-):
+def read_vector_window(inp, tile, validity_check=True, clip_to_crs_bounds=False):
     """
     Read a window of an input vector dataset.
 
@@ -45,8 +47,8 @@ def read_vector_window(
 
     Parameters
     ----------
-    input_file : string
-        path to vector file
+    inp : string or IndexedFeatures
+        path to vector file or an IndexedFeatures instance
     tile : ``Tile``
         tile extent to read data from
     validity_check : bool
@@ -71,9 +73,7 @@ def read_vector_window(
                         validity_check=validity_check,
                         clip_to_crs_bounds=clip_to_crs_bounds,
                     )
-                    for path in (
-                        input_files if isinstance(input_files, list) else [input_files]
-                    )
+                    for path in (inp if isinstance(inp, list) else [inp])
                 ]
             )
         ]
@@ -83,13 +83,11 @@ def read_vector_window(
         raise MapcheteIOError(e)
 
 
-def _read_vector_window(
-    input_file, tile, validity_check=True, clip_to_crs_bounds=False
-):
+def _read_vector_window(inp, tile, validity_check=True, clip_to_crs_bounds=False):
     if tile.pixelbuffer and tile.is_on_edge():
         return chain.from_iterable(
             _get_reprojected_features(
-                input_file=input_file,
+                inp=inp,
                 dst_bounds=bbox.bounds,
                 dst_crs=tile.crs,
                 validity_check=validity_check,
@@ -101,7 +99,7 @@ def _read_vector_window(
         )
     else:
         features = _get_reprojected_features(
-            input_file=input_file,
+            inp=inp,
             dst_bounds=tile.bounds,
             dst_crs=tile.crs,
             validity_check=validity_check,
@@ -242,62 +240,282 @@ def __exit__(self, *args):
 @retry(
     logger=logger,
     exceptions=(DriverError, FionaError, FionaValueError),
-    **MAPCHETE_IO_RETRY_SETTINGS
+    **MAPCHETE_IO_RETRY_SETTINGS,
 )
 def _get_reprojected_features(
-    input_file=None,
+    inp=None,
     dst_bounds=None,
     dst_crs=None,
     validity_check=False,
     clip_to_crs_bounds=False,
 ):
-    logger.debug("reading %s", input_file)
-    try:
-        with fiona.open(input_file, "r") as src:
-            src_crs = CRS(src.crs)
-            # reproject tile bounding box to source file CRS for filter
-            if src_crs == dst_crs:
-                dst_bbox = box(*dst_bounds)
-            else:
-                dst_bbox = reproject_geometry(
-                    box(*dst_bounds),
-                    src_crs=dst_crs,
-                    dst_crs=src_crs,
-                    validity_check=True,
+    logger.debug("reading %s", inp)
+    with ExitStack() as exit_stack:
+        if isinstance(inp, str):
+            try:
+                src = exit_stack.enter_context(fiona.open(inp, "r"))
+                src_crs = CRS(src.crs)
+            except Exception as e:
+                if path_exists(inp):
+                    logger.error("error while reading file %s: %s", inp, e)
+                    raise e
+                else:
+                    raise FileNotFoundError("%s not found" % inp)
+        else:
+            src = inp
+            src_crs = inp.crs
+        # reproject tile bounding box to source file CRS for filter
+        if src_crs == dst_crs:
+            dst_bbox = box(*dst_bounds)
+        else:
+            dst_bbox = reproject_geometry(
+                box(*dst_bounds),
+                src_crs=dst_crs,
+                dst_crs=src_crs,
+                validity_check=True,
+            )
+        for feature in src.filter(bbox=dst_bbox.bounds):
+
+            try:
+                # check validity
+                original_geom = _repair(to_shape(feature["geometry"]))
+
+                # clip with bounds and omit if clipped geometry is empty
+                clipped_geom = original_geom.intersection(dst_bbox)
+
+                # reproject each feature to tile CRS
+                g = reproject_geometry(
+                    clean_geometry_type(
+                        clipped_geom,
+                        original_geom.geom_type,
+                        raise_exception=False,
+                    ),
+                    src_crs=src_crs,
+                    dst_crs=dst_crs,
+                    validity_check=validity_check,
+                    clip_to_crs_bounds=False,
                 )
-            for feature in src.filter(bbox=dst_bbox.bounds):
-
-                try:
-                    # check validity
-                    original_geom = _repair(to_shape(feature["geometry"]))
-
-                    # clip with bounds and omit if clipped geometry is empty
-                    clipped_geom = original_geom.intersection(dst_bbox)
-
-                    # reproject each feature to tile CRS
-                    g = reproject_geometry(
-                        clean_geometry_type(
-                            clipped_geom,
-                            original_geom.geom_type,
-                            raise_exception=False,
-                        ),
-                        src_crs=src_crs,
-                        dst_crs=dst_crs,
-                        validity_check=validity_check,
-                        clip_to_crs_bounds=False,
-                    )
-                    if not g.is_empty:
-                        yield {
-                            "properties": feature["properties"],
-                            "geometry": mapping(g),
-                        }
-                # this can be handled quietly
-                except TopologicalError as e:  # pragma: no cover
-                    logger.warning("feature omitted: %s", e)
-
-    except Exception as e:
-        if path_exists(input_file):
-            logger.error("error while reading file %s: %s", input_file, e)
-            raise e
+                if not g.is_empty:
+                    yield {
+                        "properties": feature["properties"],
+                        "geometry": mapping(g),
+                    }
+            # this can be handled quietly
+            except TopologicalError as e:  # pragma: no cover
+                logger.warning("feature omitted: %s", e)
+
+
+def bounds_intersect(bounds1, bounds2):
+    bounds1 = validate_bounds(bounds1)
+    bounds2 = validate_bounds(bounds2)
+    horizontal = (
+        # partial overlap
+        bounds1.left <= bounds2.left <= bounds1.right
+        or bounds1.left <= bounds2.right <= bounds1.right
+        # bounds 1 within bounds 2
+        or bounds2.left <= bounds1.left < bounds1.right <= bounds2.right
+        # bounds 2 within bounds 1
+        or bounds1.left <= bounds2.left < bounds2.right <= bounds1.right
+    )
+    vertical = (
+        # partial overlap
+        bounds1.bottom <= bounds2.bottom <= bounds1.top
+        or bounds1.bottom <= bounds2.top <= bounds1.top
+        # bounds 1 within bounds 2
+        or bounds2.bottom <= bounds1.bottom < bounds1.top <= bounds2.top
+        # bounds 2 within bounds 1
+        or bounds1.bottom <= bounds2.bottom < bounds2.top <= bounds1.top
+    )
+    return horizontal and vertical
+
+
+class FakeIndex:
+    """Provides a fake spatial index in case rtree is not installed."""
+
+    def __init__(self):
+        self._items = []
+
+    def insert(self, id, bounds):
+        self._items.append((id, bounds))
+
+    def intersection(self, bounds):
+        return [
+            id for id, i_bounds in self._items if bounds_intersect(i_bounds, bounds)
+        ]
+
+
+class IndexedFeatures:
+    """
+    Behaves like a mapping of GeoJSON-like objects but has a filter() method.
+
+    Parameters
+    ----------
+    features : iterable
+        Features to be indexed
+    index : string
+        Spatial index to use. Can either be "rtree" (if installed) or None.
+    """
+
+    def __init__(self, features, index="rtree", allow_non_geo_objects=False, crs=None):
+        if index == "rtree":
+            try:
+                from rtree import index
+
+                self._index = index.Index()
+            except ImportError:  # pragma: no cover
+                warnings.warn(
+                    "It is recommended to install rtree in order to significantly speed up spatial indexes."
+                )
+                self._index = FakeIndex()
+        else:
+            self._index = FakeIndex()
+
+        self.crs = features.crs if hasattr(features, "crs") else crs
+        self._items = {}
+        self._non_geo_items = set()
+        self.bounds = (None, None, None, None)
+        for feature in features:
+            if isinstance(feature, tuple):
+                id_, feature = feature
+            else:
+                id_ = self._get_feature_id(feature)
+            self._items[id_] = feature
+            bounds = self._get_feature_bounds(
+                feature, allow_non_geo_objects=allow_non_geo_objects
+            )
+            if bounds is None:
+                self._non_geo_items.add(id_)
+            else:
+                self._update_bounds(bounds)
+                self._index.insert(id_, bounds)
+
+    def __repr__(self):  # pragma: no cover
+        return f"IndexedFeatures(features={len(self)}, index={self._index.__repr__()})"
+
+    def __len__(self):
+        return len(self._items)
+
+    def __str__(self):  # pragma: no cover
+        return "IndexedFeatures([%s])" % (", ".join([str(f) for f in self]))
+
+    def __getitem__(self, key):
+        try:
+            return self._items[hash(key)]
+        except KeyError:
+            raise KeyError(f"no feature with id {key} exists")
+
+    def __iter__(self):
+        return iter(self._items.values())
+
+    def items(self):
+        return self._items.items()
+
+    def keys(self):
+        return self._items.keys()
+
+    def values(self):
+        return self._items.values()
+
+    def filter(self, bounds=None, bbox=None):
+        """
+        Return features intersecting with bounds.
+
+        Parameters
+        ----------
+        bounds : list or tuple
+            Bounding coordinates (left, bottom, right, top).
+
+        Returns
+        -------
+        features : list
+            List of features.
+        """
+        bounds = bounds or bbox
+        return [
+            self._items[id_]
+            for id_ in chain(self._index.intersection(bounds), self._non_geo_items)
+        ]
+
+    def _update_bounds(self, bounds):
+        left, bottom, right, top = self.bounds
+        self.bounds = (
+            bounds.left if left is None else min(left, bounds.left),
+            bounds.bottom if bottom is None else min(bottom, bounds.bottom),
+            bounds.right if right is None else max(right, bounds.right),
+            bounds.top if top is None else max(top, bounds.top),
+        )
+
+    def _get_feature_id(self, feature):
+        if hasattr(feature, "id"):
+            return hash(feature.id)
+        elif isinstance(feature, dict) and "id" in feature:
+            return hash(feature["id"])
         else:
-            raise FileNotFoundError("%s not found" % input_file)
+            try:
+                return hash(feature)
+            except TypeError:
+                raise TypeError("features need to have an id or have to be hashable")
+
+    def _get_feature_bounds(self, feature, allow_non_geo_objects=False):
+        try:
+            if hasattr(feature, "bounds"):
+                return validate_bounds(feature.bounds)
+            elif hasattr(feature, "__geo_interface__"):
+                return validate_bounds(shape(feature).bounds)
+            elif feature.get("bounds"):
+                return validate_bounds(feature["bounds"])
+            elif feature.get("geometry"):
+                return validate_bounds(to_shape(feature["geometry"]).bounds)
+            else:
+                raise TypeError("no bounds")
+        except Exception as exc:
+            if allow_non_geo_objects:
+                return None
+            else:
+                logger.exception(exc)
+                raise TypeError(f"cannot determine bounds from feature: {feature}")
+
+
+def convert_vector(inp, out, overwrite=False, exists_ok=True, **kwargs):
+    """
+    Convert vector file to a differernt format.
+
+    When kwargs are given, the operation will be conducted by Fiona, without kwargs,
+    the file is simply copied to the destination using fsspec.
+
+    Parameters
+    ----------
+    inp : str
+        Path to input file.
+    out : str
+        Path to output file.
+    overwrite : bool
+        Overwrite output file. (default: False)
+    skip_exists : bool
+        Skip conversion if outpu already exists. (default: True)
+    kwargs : mapping
+        Creation parameters passed on to output file.
+    """
+    if path_exists(out):
+        if not exists_ok:
+            raise IOError(f"{out} already exists")
+        elif not overwrite:
+            logger.debug("output %s already exists and will not be overwritten")
+            return
+        else:
+            fs_from_path(out).rm(out)
+    kwargs = kwargs or {}
+    if kwargs:
+        logger.debug("convert raster file %s to %s using %s", inp, out, kwargs)
+        with fiona.open(inp, "r") as src:
+            makedirs(os.path.dirname(out))
+            with fiona.open(out, mode="w", **{**src.meta, **kwargs}) as dst:
+                dst.writerecords(src)
+    else:
+        logger.debug("copy %s to %s", inp, out)
+        copy(inp, out, overwrite=overwrite)
+
+
+def read_vector(inp, index="rtree"):
+    with fiona.open(inp, "r") as src:
+        return IndexedFeatures(src, index=index)
diff --git a/mapchete/tile.py b/mapchete/tile.py
index b93e91b3..aebc711f 100644
--- a/mapchete/tile.py
+++ b/mapchete/tile.py
@@ -240,6 +240,9 @@ def __init__(self, tile, pixelbuffer=0):
         Tile.__init__(self, tile.tile_pyramid, tile.zoom, tile.row, tile.col)
         self._tile = tile
         self.pixelbuffer = pixelbuffer
+        self.buffered_tp = BufferedTilePyramid(
+            tile.tp.to_dict(), pixelbuffer=pixelbuffer
+        )
 
     @cached_property
     def left(self):
@@ -352,7 +355,7 @@ def __eq__(self, other):
         return (
             isinstance(other, self.__class__)
             and self.pixelbuffer == other.pixelbuffer
-            and self.tp == other.tp
+            and self.buffered_tp == other.buffered_tp
             and self.id == other.id
         )
 
diff --git a/requirements.txt b/requirements.txt
index aa7f6397..9d4007cb 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -21,7 +21,7 @@ pyproj
 retry>=0.9.2
 rasterio>=1.0.28
 s3fs>=0.5.2
-Shapely>=1.7.1
+Shapely>=1.7.1,<1.8.1
 tilematrix>=2021.11.0
 tqdm>=4.46.0
 werkzeug>=0.15
\ No newline at end of file
diff --git a/setup.py b/setup.py
index 43d683f4..86839d8e 100644
--- a/setup.py
+++ b/setup.py
@@ -35,7 +35,7 @@
     "pyproj",
     "retry",
     "rasterio>=1.0.28,<1.2.7",
-    "shapely",
+    "shapely<1.8.1",
     "tilematrix>=2021.11.0",
     "tqdm",
 ]
@@ -45,6 +45,7 @@
 req_http = ["fsspec[http]", "aiohttp", "requests"]
 req_s3 = ["boto3", "fsspec[s3]", "s3fs>=0.5.1"]
 req_serve = ["flask", "werkzeug>=0.15"]
+req_spatial_index = ["rtree"]
 req_stac = ["pystac"]
 req_vrt = ["lxml"]
 req_complete = (
@@ -54,6 +55,7 @@
     + req_http
     + req_s3
     + req_serve
+    + req_spatial_index
     + req_stac
     + req_vrt
 )
@@ -112,6 +114,7 @@
         "http": req_http,
         "s3": req_s3,
         "serve": req_serve,
+        "spatial_index": req_spatial_index,
         "stac": req_stac,
         "vrt": req_vrt,
     },


Test patch:
diff --git a/mapchete/testing.py b/mapchete/testing.py
index 20d413c5..40c01073 100644
--- a/mapchete/testing.py
+++ b/mapchete/testing.py
@@ -64,27 +64,30 @@ def clear_output(self):
         except FileNotFoundError:
             pass
 
-    def process_mp(self, tile=None, tile_zoom=None):
+    def process_mp(self, tile=None, tile_zoom=None, batch_preprocess=True):
         """
         Return MapcheteProcess object used to test processes.
         """
+        mp = self.mp(batch_preprocess=batch_preprocess)
         if tile:
-            tile = self.mp().config.process_pyramid.tile(*tile)
+            tile = mp.config.process_pyramid.tile(*tile)
         else:
             # just use first process tile from lowest zoom level
             tile = self.first_process_tile(zoom=tile_zoom)
         return mapchete.MapcheteProcess(
             tile=tile,
-            params=self.mp().config.params_at_zoom(tile.zoom),
-            input=self.mp().config.get_inputs_for_tile(tile),
+            params=mp.config.params_at_zoom(tile.zoom),
+            input=mp.config.get_inputs_for_tile(tile),
         )
 
-    def mp(self):
+    def mp(self, batch_preprocess=True):
         """
         Return Mapchete object from mapchete.open().
         """
         if not self._mp:
             self._mp = mapchete.open(self.dict)
+            if batch_preprocess:
+                self._mp.batch_preprocess()
         return self._mp
 
     def first_process_tile(self, zoom=None):
diff --git a/test/conftest.py b/test/conftest.py
index 93a04a41..c782f060 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -616,6 +616,33 @@ def aoi_br():
         yield example
 
 
+@pytest.fixture
+def preprocess_cache_raster_vector():
+    """Fixture for preprocess_cache_raster_vector.mapchete."""
+    with ProcessFixture(
+        os.path.join(TESTDATA_DIR, "preprocess_cache_raster_vector.mapchete"),
+    ) as example:
+        yield example
+
+
+@pytest.fixture
+def preprocess_cache_memory():
+    """Fixture for preprocess_cache_memory.mapchete."""
+    with ProcessFixture(
+        os.path.join(TESTDATA_DIR, "preprocess_cache_memory.mapchete"),
+    ) as example:
+        yield example
+
+
+@pytest.fixture
+def preprocess_cache_memory_single_file():
+    """Fixture for preprocess_cache_memory_single_file.mapchete."""
+    with ProcessFixture(
+        os.path.join(TESTDATA_DIR, "preprocess_cache_memory_single_file.mapchete"),
+    ) as example:
+        yield example
+
+
 @pytest.fixture
 def custom_grid_points():
     """Fixture for custom_grid_points.mapchete."""
diff --git a/test/test_commands.py b/test/test_commands.py
index 047cd2ab..600d0f5e 100644
--- a/test/test_commands.py
+++ b/test/test_commands.py
@@ -10,6 +10,7 @@
 
 import mapchete
 from mapchete.commands import convert, cp, execute, index, rm
+from mapchete._executor import SequentialExecutor, ConcurrentFuturesExecutor
 
 
 SCRIPTDIR = os.path.dirname(os.path.realpath(__file__))
@@ -50,6 +51,7 @@ def test_cp(mp_tmpdir, cleantopo_br, wkt_geom):
     assert len(tiles)
 
 
+@pytest.mark.remote
 def test_cp_http(mp_tmpdir, http_tiledir):
     # copy tiles and subset by bounds
     tiles = cp(
@@ -89,6 +91,37 @@ def test_execute(mp_tmpdir, cleantopo_br_metatiling_1, cleantopo_br_tif):
             assert not src.read(masked=True).mask.all()
 
 
+def test_execute_set_executor(mp_tmpdir, cleantopo_br_metatiling_1, cleantopo_br_tif):
+    zoom = 5
+    tp = TilePyramid("geodetic")
+    tiles = list(tp.tiles_from_bounds(rasterio.open(cleantopo_br_tif).bounds, zoom))
+    job = execute(cleantopo_br_metatiling_1.dict, zoom=zoom, as_iterator=True)
+
+    # invalid concurrency
+    job.set_executor_concurrency("invalid")
+    with pytest.raises(ValueError):
+        next(iter(job))
+
+    # reset to concurrent.futures
+    job.set_executor_concurrency("processes")
+    next(iter(job))
+    assert isinstance(job.executor, ConcurrentFuturesExecutor)
+
+    # with max_workers: 1 it should fall back to a sequential executor
+    job.set_executor_kwargs({"max_workers": 1})
+    next(iter(job))
+    assert isinstance(job.executor, SequentialExecutor)
+
+    # make sure everything runs smoothly from here on
+    for t in job:
+        assert t
+    assert len(tiles) == len(job)
+    mp = cleantopo_br_metatiling_1.mp()
+    for t in tiles:
+        with rasterio.open(mp.config.output.get_path(t)) as src:
+            assert not src.read(masked=True).mask.all()
+
+
 def test_execute_cancel(mp_tmpdir, cleantopo_br_metatiling_1, cleantopo_br_tif):
     zoom = 5
     job = execute(cleantopo_br_metatiling_1.dict, zoom=zoom, as_iterator=True)
@@ -121,6 +154,11 @@ def test_execute_point(mp_tmpdir, example_mapchete, dummy2_tif):
     assert len(job) == 1
 
 
+def test_execute_preprocessing_tasks_dask(preprocess_cache_raster_vector):
+    job = execute(preprocess_cache_raster_vector.dict, concurrency="dask")
+    assert len(job)
+
+
 def test_convert_geodetic(cleantopo_br_tif, mp_tmpdir):
     """Automatic geodetic tile pyramid creation of raster files."""
     job = convert(cleantopo_br_tif, mp_tmpdir, output_pyramid="geodetic")
@@ -209,14 +247,34 @@ def test_convert_single_gtiff_cog(cleantopo_br_tif, mp_tmpdir):
     """Automatic geodetic tile pyramid creation of raster files."""
     single_gtiff = os.path.join(mp_tmpdir, "single_out_cog.tif")
     job = convert(
-        cleantopo_br_tif, single_gtiff, output_pyramid="geodetic", zoom=3, cog=True
+        cleantopo_br_tif, single_gtiff, output_pyramid="geodetic", zoom=5, cog=True
     )
     assert len(job)
     with rasterio.open(single_gtiff, "r") as src:
         assert src.meta["driver"] == "GTiff"
         assert src.meta["dtype"] == "uint16"
         data = src.read(masked=True)
-        assert data.mask.any()
+        assert not data.mask.all()
+    assert cog_validate(single_gtiff, strict=True)
+
+
+def test_convert_single_gtiff_cog_dask(cleantopo_br_tif, mp_tmpdir):
+    """Automatic geodetic tile pyramid creation of raster files."""
+    single_gtiff = os.path.join(mp_tmpdir, "single_out_cog.tif")
+    job = convert(
+        cleantopo_br_tif,
+        single_gtiff,
+        output_pyramid="geodetic",
+        zoom=5,
+        cog=True,
+        concurrency="dask",
+    )
+    assert len(job)
+    with rasterio.open(single_gtiff, "r") as src:
+        assert src.meta["driver"] == "GTiff"
+        assert src.meta["dtype"] == "uint16"
+        data = src.read(masked=True)
+        assert not data.mask.all()
     assert cog_validate(single_gtiff, strict=True)
 
 
@@ -241,6 +299,7 @@ def test_convert_single_gtiff_overviews(cleantopo_br_tif, mp_tmpdir):
         assert src.overviews(1)
 
 
+@pytest.mark.remote
 def test_convert_remote_single_gtiff(http_raster, mp_tmpdir):
     """Automatic geodetic tile pyramid creation of raster files."""
     single_gtiff = os.path.join(mp_tmpdir, "single_out.tif")
@@ -279,7 +338,6 @@ def test_convert_scale_ratio(cleantopo_br_tif, mp_tmpdir):
         output_dtype="uint8",
         scale_ratio=0.003,
     )
-    print(job)
     assert len(job)
     for zoom, row, col in [(4, 15, 15), (3, 7, 7)]:
         out_file = os.path.join(*[mp_tmpdir, str(zoom), str(row), str(col) + ".tif"])
diff --git a/test/test_io.py b/test/test_io.py
index d3590085..513143c4 100644
--- a/test/test_io.py
+++ b/test/test_io.py
@@ -1,26 +1,11 @@
 """Test Mapchete io module."""
 
 import pytest
-import shutil
 import rasterio
-import tempfile
-import numpy as np
-import numpy.ma as ma
-import fiona
-from fiona.errors import DriverError
 import os
-from rasterio.crs import CRS
-from rasterio.enums import Compression
-from shapely import wkt
 from shapely.errors import TopologicalError
-from shapely.geometry import shape, box, Polygon, MultiPolygon, LineString
-from shapely.ops import unary_union
-from tilematrix import Bounds
-from itertools import product
 
 import mapchete
-from mapchete.config import MapcheteConfig
-from mapchete.errors import GeometryTypeError, MapcheteIOError
 from mapchete.io import (
     get_best_zoom_level,
     path_exists,
@@ -28,25 +13,7 @@
     read_json,
     tile_to_zoom_level,
     tiles_exist,
-)
-from mapchete.io.raster import (
-    read_raster_window,
-    write_raster_window,
-    extract_from_array,
-    resample_from_array,
-    create_mosaic,
-    ReferencedRaster,
-    prepare_array,
-    RasterWindowMemoryFile,
-    read_raster_no_crs,
-)
-from mapchete.io.vector import (
-    read_vector_window,
-    reproject_geometry,
-    clean_geometry_type,
-    segmentize_geometry,
-    write_vector_window,
-    _repair,
+    copy,
 )
 from mapchete.tile import BufferedTilePyramid
 
@@ -61,990 +28,11 @@ def test_best_zoom_level(dummy1_tif):
     assert get_best_zoom_level(dummy1_tif, "mercator")
 
 
-def test_read_raster_window(dummy1_tif, minmax_zoom):
-    """Read array with read_raster_window."""
-    zoom = 8
-    # without reproject
-    config = MapcheteConfig(minmax_zoom.dict)
-    rasterfile = config.params_at_zoom(zoom)["input"]["file1"]
-    dummy1_bbox = rasterfile.bbox()
-
-    pixelbuffer = 5
-    tile_pyramid = BufferedTilePyramid("geodetic", pixelbuffer=pixelbuffer)
-    tiles = list(tile_pyramid.tiles_from_geom(dummy1_bbox, zoom))
-    # add edge tile
-    tiles.append(tile_pyramid.tile(8, 0, 0))
-    for tile in tiles:
-        width, height = tile.shape
-        for band in read_raster_window(dummy1_tif, tile):
-            assert isinstance(band, ma.MaskedArray)
-            assert band.shape == (width, height)
-        for index in range(1, 4):
-            band = read_raster_window(dummy1_tif, tile, index)
-            assert isinstance(band, ma.MaskedArray)
-            assert band.shape == (width, height)
-        for index in [None, [1, 2, 3]]:
-            band = read_raster_window(dummy1_tif, tile, index)
-            assert isinstance(band, ma.MaskedArray)
-            assert band.ndim == 3
-            assert band.shape == (3, width, height)
-
-
-def test_read_raster_window_reproject(dummy1_3857_tif, minmax_zoom):
-    """Read array with read_raster_window."""
-    zoom = 8
-    # with reproject
-    config_raw = minmax_zoom.dict
-    config_raw["input"].update(file1=dummy1_3857_tif)
-    config = MapcheteConfig(config_raw)
-    rasterfile = config.params_at_zoom(zoom)["input"]["file1"]
-    dummy1_bbox = rasterfile.bbox()
-
-    pixelbuffer = 5
-    tile_pyramid = BufferedTilePyramid("geodetic", pixelbuffer=pixelbuffer)
-    tiles = list(tile_pyramid.tiles_from_geom(dummy1_bbox, zoom))
-    # target window out of CRS bounds
-    band = read_raster_window(dummy1_3857_tif, tile_pyramid.tile(12, 0, 0))
-    assert isinstance(band, ma.MaskedArray)
-    assert band.mask.all()
-    # not intersecting tile
-    tiles.append(tile_pyramid.tile(zoom, 1, 1))  # out of CRS bounds
-    tiles.append(tile_pyramid.tile(zoom, 16, 1))  # out of file bbox
-    for tile in tiles:
-        for band in read_raster_window(dummy1_3857_tif, tile):
-            assert isinstance(band, ma.MaskedArray)
-            assert band.shape == tile.shape
-        bands = read_raster_window(dummy1_3857_tif, tile, [1])
-        assert isinstance(bands, ma.MaskedArray)
-        assert bands.shape == tile.shape
-    # errors
-    with pytest.raises(IOError):
-        read_raster_window("nonexisting_path", tile)
-
-
-def test_read_raster_window_resampling(cleantopo_br_tif):
-    """Assert various resampling options work."""
-    tp = BufferedTilePyramid("geodetic")
-    with rasterio.open(cleantopo_br_tif, "r") as src:
-        tiles = tp.tiles_from_bounds(src.bounds, 4)
-    for tile in tiles:
-        outputs = [
-            read_raster_window(cleantopo_br_tif, tile, resampling=resampling)
-            for resampling in [
-                "nearest",
-                "bilinear",
-                "cubic",
-                "cubic_spline",
-                "lanczos",
-                "average",
-                "mode",
-            ]
-        ]
-        # resampling test:
-        assert any(
-            [not np.array_equal(w, v) for v, w in zip(outputs[:-1], outputs[1:])]
-        )
-
-
-def test_read_raster_window_partly_overlapping(cleantopo_br_tif):
-    """Read array with read_raster_window where window is bigger than file."""
-    tile = BufferedTilePyramid("geodetic").tile(4, 15, 31)
-    data = read_raster_window(cleantopo_br_tif, tile)
-    assert isinstance(data, ma.MaskedArray)
-    assert data.mask.any()
-
-
-def test_read_raster_window_mask(s2_band):
-    """No resampling artefacts on mask edges."""
-    tile = BufferedTilePyramid("geodetic").tile(zoom=13, row=1918, col=8905)
-    data = read_raster_window(
-        s2_band, tile, resampling="cubic", src_nodata=0, dst_nodata=0
-    )
-    assert data.any()
-    assert not np.where(data == 1, True, False).any()
-
-
-def test_read_raster_window_input_list(cleantopo_br):
-    process_zoom = 5
-    conf = dict(**cleantopo_br.dict)
-    conf["output"].update(metatiling=1)
-    with mapchete.open(conf) as mp:
-        mp.batch_process(process_zoom)
-        tiles = [
-            (tile, mp.config.output.get_path(tile))
-            for tile in mp.config.output_pyramid.tiles_from_bounds(
-                mp.config.bounds, process_zoom
-            )
-            if path_exists(mp.config.output.get_path(tile))
-        ]
-        upper_tile = next(mp.get_process_tiles(process_zoom - 1))
-        assert len(tiles) > 1
-        resampled = resample_from_array(
-            in_raster=create_mosaic(
-                [(tile, read_raster_window(path, tile)) for tile, path in tiles]
-            ),
-            out_tile=upper_tile,
-        )
-    resampled2 = read_raster_window(
-        [p for _, p in tiles], upper_tile, src_nodata=0, dst_nodata=0
-    )
-    assert resampled.dtype == resampled2.dtype
-    assert resampled.shape == resampled2.shape
-    assert np.array_equal(resampled.mask, resampled2.mask)
-    # TODO slight rounding errors occur
-    assert np.allclose(resampled, resampled2, rtol=0.01)
-
-
-def test_read_raster_window_retry(invalid_tif):
-    tile = BufferedTilePyramid("geodetic").tile(zoom=13, row=1918, col=8905)
-    with pytest.raises(MapcheteIOError):
-        read_raster_window(invalid_tif, tile)
-    with pytest.raises(FileNotFoundError):
-        read_raster_window("not_existing.tif", tile)
-
-
-def test_read_raster_no_crs_errors():
-    with tempfile.NamedTemporaryFile() as tmpfile:
-        with pytest.raises(MapcheteIOError):
-            read_raster_no_crs(tmpfile.name)
-
-
-def test_write_raster_window():
-    """Basic output format writing."""
-    path = tempfile.NamedTemporaryFile(delete=False).name
-    # standard tile
-    tp = BufferedTilePyramid("geodetic")
-    tile = tp.tile(5, 5, 5)
-    data = ma.masked_array(np.ones((2,) + tile.shape))
-    for out_profile in [
-        dict(
-            driver="GTiff",
-            count=2,
-            dtype="uint8",
-            compress="lzw",
-            nodata=0,
-            height=tile.height,
-            width=tile.width,
-            affine=tile.affine,
-        ),
-        dict(
-            driver="GTiff",
-            count=2,
-            dtype="uint8",
-            compress="deflate",
-            nodata=0,
-            height=tile.height,
-            width=tile.width,
-            affine=tile.affine,
-        ),
-        dict(
-            driver="PNG",
-            count=2,
-            dtype="uint8",
-            nodata=0,
-            height=tile.height,
-            width=tile.width,
-            compress=None,
-            affine=tile.affine,
-        ),
-    ]:
-        try:
-            write_raster_window(
-                in_tile=tile, in_data=data, out_profile=out_profile, out_path=path
-            )
-            with rasterio.open(path, "r") as src:
-                assert src.read().any()
-                assert src.meta["driver"] == out_profile["driver"]
-                assert src.transform == tile.affine
-                if out_profile["compress"]:
-                    assert src.compression == Compression(
-                        out_profile["compress"].upper()
-                    )
-        finally:
-            shutil.rmtree(path, ignore_errors=True)
-    # with metatiling
-    tile = BufferedTilePyramid("geodetic", metatiling=4).tile(5, 1, 1)
-    data = ma.masked_array(np.ones((2,) + tile.shape))
-    out_tile = BufferedTilePyramid("geodetic").tile(5, 5, 5)
-    out_profile = dict(
-        driver="GTiff",
-        count=2,
-        dtype="uint8",
-        compress="lzw",
-        nodata=0,
-        height=out_tile.height,
-        width=out_tile.width,
-        affine=out_tile.affine,
-    )
-    try:
-        write_raster_window(
-            in_tile=tile,
-            in_data=data,
-            out_profile=out_profile,
-            out_tile=out_tile,
-            out_path=path,
-        )
-        with rasterio.open(path, "r") as src:
-            assert src.shape == out_tile.shape
-            assert src.read().any()
-            assert src.meta["driver"] == out_profile["driver"]
-            assert src.transform == out_profile["transform"]
-    finally:
-        shutil.rmtree(path, ignore_errors=True)
-
-
-def test_write_raster_window_memory():
-    """Basic output format writing."""
-    path = "memoryfile"
-    # standard tile
-    tp = BufferedTilePyramid("geodetic")
-    tile = tp.tile(5, 5, 5)
-    data = ma.masked_array(np.ones((2,) + tile.shape))
-    for out_profile in [
-        dict(
-            driver="GTiff",
-            count=2,
-            dtype="uint8",
-            compress="lzw",
-            nodata=0,
-            height=tile.height,
-            width=tile.width,
-            affine=tile.affine,
-        ),
-        dict(
-            driver="GTiff",
-            count=2,
-            dtype="uint8",
-            compress="deflate",
-            nodata=0,
-            height=tile.height,
-            width=tile.width,
-            affine=tile.affine,
-        ),
-        dict(
-            driver="PNG",
-            count=2,
-            dtype="uint8",
-            nodata=0,
-            height=tile.height,
-            width=tile.width,
-            compress=None,
-            affine=tile.affine,
-        ),
-    ]:
-        with pytest.raises(DeprecationWarning):
-            write_raster_window(
-                in_tile=tile, in_data=data, out_profile=out_profile, out_path=path
-            )
-
-
-def test_raster_window_memoryfile():
-    """Use context manager for rasterio MemoryFile."""
-    tp = BufferedTilePyramid("geodetic")
-    tile = tp.tile(5, 5, 5)
-    data = ma.masked_array(np.ones((2,) + tile.shape))
-    for out_profile in [
-        dict(
-            driver="GTiff",
-            count=2,
-            dtype="uint8",
-            compress="lzw",
-            nodata=0,
-            height=tile.height,
-            width=tile.width,
-            affine=tile.affine,
-        ),
-        dict(
-            driver="GTiff",
-            count=2,
-            dtype="uint8",
-            compress="deflate",
-            nodata=0,
-            height=tile.height,
-            width=tile.width,
-            affine=tile.affine,
-        ),
-        dict(
-            driver="PNG",
-            count=2,
-            dtype="uint8",
-            nodata=0,
-            height=tile.height,
-            width=tile.width,
-            compress=None,
-            affine=tile.affine,
-        ),
-    ]:
-        with RasterWindowMemoryFile(
-            in_tile=tile, in_data=data, out_profile=out_profile
-        ) as memfile:
-            with memfile.open() as src:
-                assert src.read().any()
-                assert src.meta["driver"] == out_profile["driver"]
-                assert src.transform == tile.affine
-                if out_profile["compress"]:
-                    assert src.compression == Compression(
-                        out_profile["compress"].upper()
-                    )
-
-
-def test_write_raster_window_errors():
-    """Basic output format writing."""
-    tile = BufferedTilePyramid("geodetic").tile(5, 5, 5)
-    data = ma.masked_array(np.ndarray((1, 1)))
-    profile = {}
-    path = ""
-    # in_tile
-    with pytest.raises(TypeError):
-        write_raster_window(
-            in_tile="invalid tile",
-            in_data=data,
-            out_profile=profile,
-            out_tile=tile,
-            out_path=path,
-        )
-    # out_tile
-    with pytest.raises(TypeError):
-        write_raster_window(
-            in_tile=tile,
-            in_data=data,
-            out_profile=profile,
-            out_tile="invalid tile",
-            out_path=path,
-        )
-    # in_data
-    with pytest.raises(TypeError):
-        write_raster_window(
-            in_tile=tile,
-            in_data="invalid data",
-            out_profile=profile,
-            out_tile=tile,
-            out_path=path,
-        )
-    # out_profile
-    with pytest.raises(TypeError):
-        write_raster_window(
-            in_tile=tile,
-            in_data=data,
-            out_profile="invalid profile",
-            out_tile=tile,
-            out_path=path,
-        )
-    # out_path
-    with pytest.raises(TypeError):
-        write_raster_window(
-            in_tile=tile, in_data=data, out_profile=profile, out_tile=tile, out_path=999
-        )
-    # cannot write
-    with pytest.raises(ValueError):
-        write_raster_window(
-            in_tile=tile,
-            in_data=data,
-            out_profile=profile,
-            out_tile=tile,
-            out_path="/invalid_path",
-        )
-
-
-def test_extract_from_array():
-    """Extract subdata from array."""
-    in_tile = BufferedTilePyramid("geodetic", metatiling=4).tile(5, 5, 5)
-    shape = (in_tile.shape[0] // 2, in_tile.shape[1])
-    data = ma.masked_array(np.concatenate([np.ones(shape), np.ones(shape) * 2]))
-    # intersecting at top
-    out_tile = BufferedTilePyramid("geodetic").tile(5, 20, 20)
-    out_array = extract_from_array(
-        in_raster=data, in_affine=in_tile.affine, out_tile=out_tile
-    )
-    assert isinstance(out_array, np.ndarray)
-    assert np.all(np.where(out_array == 1, True, False))
-    # intersecting at bottom
-    out_tile = BufferedTilePyramid("geodetic").tile(5, 22, 20)
-    out_array = extract_from_array(
-        in_raster=data, in_affine=in_tile.affine, out_tile=out_tile
-    )
-    assert isinstance(out_array, np.ndarray)
-    assert np.all(np.where(out_array == 2, True, False))
-    # not intersecting
-    out_tile = BufferedTilePyramid("geodetic").tile(5, 15, 20)
-    with pytest.raises(ValueError):
-        out_array = extract_from_array(
-            in_raster=data, in_affine=in_tile.affine, out_tile=out_tile
-        )
-
-
-def test_resample_from_array():
-    """Resample array using rasterio reproject()."""
-    in_tile = BufferedTilePyramid("geodetic").tile(5, 5, 5)
-    in_data = np.ones(in_tile.shape)
-    # tile from next toom level
-    out_tile = BufferedTilePyramid("geodetic").tile(6, 10, 10)
-    out_array = resample_from_array(in_data, in_tile.affine, out_tile)
-    assert isinstance(out_array, ma.masked_array)
-    assert np.all(np.where(out_array == 1, True, False))
-    # not intersecting tile
-    out_tile = BufferedTilePyramid("geodetic").tile(7, 0, 0)
-    out_array = resample_from_array(in_data, in_tile.affine, out_tile)
-    assert isinstance(out_array, ma.masked_array)
-    assert out_array.mask.all()
-    # data as tuple
-    in_data = (np.ones(in_tile.shape[1:]),)
-    out_tile = BufferedTilePyramid("geodetic").tile(6, 10, 10)
-    out_array = resample_from_array(in_data, in_tile.affine, out_tile)
-    # deprecated
-    resample_from_array(in_data, in_tile.affine, out_tile, nodata=-9999)
-    # errors
-    with pytest.raises(TypeError):
-        in_data = "invalid_type"
-        resample_from_array(in_data, in_tile.affine, out_tile)
-    with pytest.raises(TypeError):
-        in_data = np.ones(in_tile.shape[0])
-        resample_from_array(in_data, in_tile.affine, out_tile)
-
-
-def test_create_mosaic_errors():
-    """Check error handling of create_mosaic()."""
-    tp_geo = BufferedTilePyramid("geodetic")
-    tp_mer = BufferedTilePyramid("mercator")
-    geo_tile = tp_geo.tile(1, 0, 0)
-    geo_tile_data = np.ndarray(geo_tile.shape)
-    mer_tile = tp_mer.tile(1, 1, 0)
-    mer_tile_data = np.ndarray(mer_tile.shape)
-    # tiles error
-    with pytest.raises(TypeError):
-        create_mosaic("invalid tiles")
-    with pytest.raises(TypeError):
-        create_mosaic(["invalid tiles"])
-    # CRS error
-    with pytest.raises(ValueError):
-        create_mosaic([(geo_tile, geo_tile_data), (mer_tile, mer_tile_data)])
-    # zoom error
-    with pytest.raises(ValueError):
-        diff_zoom = tp_geo.tile(2, 1, 0)
-        diff_zoom_data = np.ndarray(diff_zoom.shape)
-        create_mosaic([(geo_tile, geo_tile_data), (diff_zoom, diff_zoom_data)])
-    # tile data error
-    with pytest.raises(TypeError):
-        # for one tile
-        create_mosaic([(geo_tile, None)])
-    with pytest.raises(TypeError):
-        # for multiple tiles
-        create_mosaic([(geo_tile, None), (geo_tile, None)])
-    # tile data type error
-    with pytest.raises(TypeError):
-        diff_type = tp_geo.tile(1, 1, 0)
-        diff_type_data = np.ndarray(diff_zoom.shape).astype("int")
-        create_mosaic([(geo_tile, geo_tile_data), (diff_type, diff_type_data)])
-    # no tiles
-    with pytest.raises(ValueError):
-        create_mosaic(tiles=[])
-
-
-def test_create_mosaic():
-    """Create mosaic from tiles."""
-    tp = BufferedTilePyramid("geodetic")
-    # quick return mosaic if there is just one tile
-    tile = tp.tile(3, 3, 3)
-    data = np.ones(tile.shape)
-    mosaic = create_mosaic([(tile, data)])
-    assert isinstance(mosaic, ReferencedRaster)
-    assert np.array_equal(data, mosaic.data)
-    assert tile.affine == mosaic.affine
-    zoom = 5
-    # multiple tiles on top left corner of tile matrix
-    for pixelbuffer in [0, 10]:
-        tp = BufferedTilePyramid("geodetic", pixelbuffer=pixelbuffer)
-        tiles = [
-            (tp.tile(zoom, row, col), np.ones(tp.tile(zoom, row, col).shape))
-            for row, col in product(range(4), range(4))
-        ]
-        # 4x4 top left tiles from zoom 5 equal top left tile from zoom 3
-        # also use tile generator
-        mosaic = create_mosaic((t for t in tiles))
-        assert isinstance(mosaic, ReferencedRaster)
-        assert np.all(np.where(mosaic.data == 1, True, False))
-        mosaic_bbox = box(
-            mosaic.affine[2],
-            mosaic.affine[5] + mosaic.data.shape[1] * mosaic.affine[4],
-            mosaic.affine[2] + mosaic.data.shape[2] * mosaic.affine[0],
-            mosaic.affine[5],
-        )
-        control_bbox = box(*unary_union([t.bbox for t, _ in tiles]).bounds)
-        assert mosaic_bbox.equals(control_bbox)
-    # multiple tiles on bottom right corner of tile matrix
-    for pixelbuffer in [0, 10]:
-        tp = BufferedTilePyramid("geodetic", pixelbuffer=pixelbuffer)
-        tiles = [
-            (tp.tile(zoom, row, col), np.ones(tp.tile(zoom, row, col).shape))
-            for row, col in product(
-                range(tp.matrix_height(zoom) - 4, tp.matrix_height(zoom)),
-                range(tp.matrix_width(zoom) - 4, tp.matrix_width(zoom)),
-            )
-        ]
-        # 4x4 top left tiles from zoom 5 equal top left tile from zoom 3
-        # also use tile generator
-        mosaic = create_mosaic((t for t in tiles))
-        assert isinstance(mosaic, ReferencedRaster)
-        assert np.all(np.where(mosaic.data == 1, True, False))
-        mosaic_bbox = box(
-            mosaic.affine[2],
-            mosaic.affine[5] + mosaic.data.shape[1] * mosaic.affine[4],
-            mosaic.affine[2] + mosaic.data.shape[2] * mosaic.affine[0],
-            mosaic.affine[5],
-        )
-        control_bbox = box(*unary_union([t.bbox for t, _ in tiles]).bounds)
-        assert mosaic_bbox.equals(control_bbox)
-
-
-def test_create_mosaic_antimeridian():
-    """Create mosaic using tiles on opposing antimeridian sides."""
-    zoom = 5
-    row = 0
-    pixelbuffer = 5
-    tp = BufferedTilePyramid("geodetic", pixelbuffer=pixelbuffer)
-    west = tp.tile(zoom, row, 0)
-    east = tp.tile(zoom, row, tp.matrix_width(zoom) - 1)
-    mosaic = create_mosaic(
-        [
-            (west, np.ones(west.shape).astype("uint8")),
-            (east, np.ones(east.shape).astype("uint8") * 2),
-        ]
-    )
-    assert isinstance(mosaic, ReferencedRaster)
-
-    # Huge array gets initialized because the two tiles are on opposing sides of the
-    # projection area. The below test should pass if the tiles are stitched together next
-    # to each other.
-    assert mosaic.data.shape == (1, west.height, west.width * 2 - 2 * pixelbuffer)
-    assert mosaic.data[0][0][0] == 2
-    assert mosaic.data[0][0][-1] == 1
-
-    # If tiles from opposing sides from Antimeridian are mosaicked it will happen that the
-    # output mosaic exceeds the CRS bounds (obviously). In such a case the mosaicking
-    # function shall make sure that the larger part of the output mosaic shall be inside
-    # the CRS bounds.
-
-    # (1) mosaic crosses Antimeridian in the West, larger part is on Western hemisphere:
-    tiles_ids = [
-        # Western hemisphere tiles
-        (zoom, row, 0),
-        (zoom, row, 1),
-        # Eastern hemisphere tile
-        (zoom, row, tp.matrix_width(zoom) - 1),
-    ]
-    tiles = [
-        (tp.tile(*tile_id), np.ones(tp.tile(*tile_id).shape)) for tile_id in tiles_ids
-    ]
-    mosaic = create_mosaic(tiles)
-    control_bounds = Bounds(
-        # Eastern tile has to be shifted
-        -(360 - tp.tile(*tiles_ids[2]).left),
-        tp.tile(*tiles_ids[2]).bottom,
-        tp.tile(*tiles_ids[1]).right,
-        tp.tile(*tiles_ids[1]).top,
-    )
-    assert mosaic.bounds == control_bounds
-
-    # (2) mosaic crosses Antimeridian in the West, larger part is on Eastern hemisphere:
-    tiles_ids = [
-        # Western hemisphere tile
-        (zoom, row, 0),
-        # Eastern hemisphere tiles
-        (zoom, row, tp.matrix_width(zoom) - 1),
-        (zoom, row, tp.matrix_width(zoom) - 2),
-    ]
-    tiles = [
-        (tp.tile(*tile_id), np.ones(tp.tile(*tile_id).shape)) for tile_id in tiles_ids
-    ]
-    mosaic = create_mosaic(tiles)
-    control_bounds = Bounds(
-        tp.tile(*tiles_ids[2]).left,
-        tp.tile(*tiles_ids[2]).bottom,
-        # Western tile has to be shifted
-        360 + tp.tile(*tiles_ids[0]).right,
-        tp.tile(*tiles_ids[0]).top,
-    )
-    assert mosaic.bounds == control_bounds
-
-
-def test_prepare_array_iterables():
-    """Convert iterable data into a proper array."""
-    # input is iterable
-    # iterable contains arrays
-    data = [np.zeros((1, 1))]
-    # output ndarray
-    output = prepare_array(data, masked=False)
-    assert isinstance(output, np.ndarray)
-    assert not isinstance(output, ma.masked_array)
-    assert output.shape == (1, 1, 1)
-    # output masked array
-    output = prepare_array(data)
-    assert isinstance(output, ma.masked_array)
-    assert output.shape == (1, 1, 1)
-    # iterable contains masked arrays
-    data = [ma.empty((1, 1))]
-    output = prepare_array(data, masked=False)
-    assert isinstance(output, np.ndarray)
-    assert not isinstance(output, ma.masked_array)
-    assert output.shape == (1, 1, 1)
-    # output masked array
-    output = prepare_array(data)
-    assert isinstance(output, ma.masked_array)
-    assert output.shape == (1, 1, 1)
-    # iterable contains masked arrays with full mask
-    data = [ma.masked_array(data=np.ones((1, 1)), mask=np.ones((1, 1)))]
-    output = prepare_array(data, masked=False)
-    assert isinstance(output, np.ndarray)
-    assert not isinstance(output, ma.masked_array)
-    assert output.shape == (1, 1, 1)
-    # output masked array
-    output = prepare_array(data)
-    assert isinstance(output, ma.masked_array)
-    assert output.shape == (1, 1, 1)
-
-
-def test_prepare_array_maskedarrays():
-    """Convert masked array data into a proper array."""
-    # input is ma.masked_array
-    data = ma.empty((1, 1, 1))
-    # output ndarray
-    output = prepare_array(data, masked=False)
-    assert isinstance(output, np.ndarray)
-    assert not isinstance(output, ma.masked_array)
-    assert output.shape == (1, 1, 1)
-    # output masked array
-    output = prepare_array(data)
-    assert isinstance(output, ma.masked_array)
-    assert output.shape == (1, 1, 1)
-    # input is ma.masked_array with full mask
-    data = ma.masked_array(data=np.ones((1, 1, 1)), mask=np.ones((1, 1, 1)))
-    # output ndarray
-    output = prepare_array(data, masked=False)
-    assert isinstance(output, np.ndarray)
-    assert not isinstance(output, ma.masked_array)
-    assert output.shape == (1, 1, 1)
-    # output masked array
-    output = prepare_array(data)
-    assert isinstance(output, ma.masked_array)
-    assert output.shape == (1, 1, 1)
-
-
-def test_prepare_array_ndarrays():
-    """Convert ndarray data into a proper array."""
-    # input is np.ndarray
-    data = np.zeros((1, 1, 1))
-    # output ndarray
-    output = prepare_array(data, masked=False)
-    assert isinstance(output, np.ndarray)
-    assert not isinstance(output, ma.masked_array)
-    assert output.shape == (1, 1, 1)
-    # output masked array
-    output = prepare_array(data)
-    assert isinstance(output, ma.masked_array)
-    assert output.shape == (1, 1, 1)
-    # input is 2D np.ndarray
-    data = np.zeros((1, 1))
-    # output ndarray
-    output = prepare_array(data, masked=False)
-    assert isinstance(output, np.ndarray)
-    assert not isinstance(output, ma.masked_array)
-    assert output.shape == (1, 1, 1)
-    # output masked array
-    output = prepare_array(data)
-    assert isinstance(output, ma.masked_array)
-    assert output.shape == (1, 1, 1)
-
-
-def test_prepare_array_errors():
-    """Convert ndarray data into a proper array."""
-    # input is iterable
-    data = [None]
-    try:
-        prepare_array(data)
-        raise Exception()
-    except ValueError:
-        pass
-    # input is not array
-    data = 5
-    try:
-        prepare_array(data)
-        raise Exception()
-    except ValueError:
-        pass
-
-
-def test_read_vector_window(geojson, landpoly_3857):
-    """Read vector data from read_vector_window."""
-    zoom = 4
-    config = MapcheteConfig(geojson.dict)
-    vectorfile = config.params_at_zoom(zoom)["input"]["file1"]
-    pixelbuffer = 5
-    tile_pyramid = BufferedTilePyramid("geodetic", pixelbuffer=pixelbuffer)
-    tiles = tile_pyramid.tiles_from_geom(vectorfile.bbox(), zoom)
-    feature_count = 0
-    for tile in tiles:
-        for feature in read_vector_window(vectorfile.path, tile):
-            assert "properties" in feature
-            assert shape(feature["geometry"]).is_valid
-            feature_count += 1
-    assert feature_count
-    # into different CRS
-    raw_config = geojson.dict
-    raw_config["input"].update(file1=landpoly_3857)
-    config = MapcheteConfig(raw_config)
-    vectorfile = config.params_at_zoom(zoom)["input"]["file1"]
-    pixelbuffer = 5
-    tile_pyramid = BufferedTilePyramid("geodetic", pixelbuffer=pixelbuffer)
-    tiles = tile_pyramid.tiles_from_geom(vectorfile.bbox(), zoom)
-    feature_count = 0
-    for tile in tiles:
-        for feature in read_vector_window(vectorfile.path, tile):
-            assert "properties" in feature
-            assert shape(feature["geometry"]).is_valid
-            feature_count += 1
-    assert feature_count
-
-
-def test_read_vector_window_errors(invalid_geojson):
-    with pytest.raises(FileNotFoundError):
-        read_vector_window(
-            "invalid_path", BufferedTilePyramid("geodetic").tile(0, 0, 0)
-        )
-    with pytest.raises(MapcheteIOError):
-        read_vector_window(
-            invalid_geojson, BufferedTilePyramid("geodetic").tile(0, 0, 0)
-        )
-
-
-def test_reproject_geometry(landpoly):
-    """Reproject geometry."""
-    with fiona.open(landpoly, "r") as src:
-        for feature in src:
-
-            # WGS84 to Spherical Mercator
-            out_geom = reproject_geometry(
-                shape(feature["geometry"]), CRS(src.crs), CRS().from_epsg(3857)
-            )
-            assert out_geom.is_valid
-
-            # WGS84 to LAEA
-            out_geom = reproject_geometry(
-                shape(feature["geometry"]), CRS(src.crs), CRS().from_epsg(3035)
-            )
-            assert out_geom.is_valid
-
-            # WGS84 to WGS84
-            out_geom = reproject_geometry(
-                shape(feature["geometry"]), CRS(src.crs), CRS().from_epsg(4326)
-            )
-            assert out_geom.is_valid
-
-
-def test_reproject_geometry_latlon2mercator():
-    # WGS84 bounds to Spherical Mercator
-    big_box = box(-180, -90, 180, 90)
-    reproject_geometry(big_box, CRS().from_epsg(4326), CRS().from_epsg(3857))
-
-    # WGS84 bounds to Spherical Mercator raising clip error
-    with pytest.raises(RuntimeError):
-        reproject_geometry(
-            big_box, CRS().from_epsg(4326), CRS().from_epsg(3857), error_on_clip=True
-        )
-    outside_box = box(-180, 87, 180, 90)
-    assert reproject_geometry(
-        outside_box,
-        CRS().from_epsg(4326),
-        CRS().from_epsg(3857),
-    ).is_valid
-
-
-def test_reproject_geometry_empty_geom():
-    # empty geometry
-    assert reproject_geometry(
-        Polygon(), CRS().from_epsg(4326), CRS().from_epsg(3857)
-    ).is_empty
-    assert reproject_geometry(
-        Polygon(), CRS().from_epsg(4326), CRS().from_epsg(4326)
-    ).is_empty
-
-
-def test_reproject_geometry_latlon2mercator_epsg():
-    # CRS parameter
-    big_box = box(-180, -90, 180, 90)
-    assert reproject_geometry(big_box, 4326, 3857) == reproject_geometry(
-        big_box, "4326", "3857"
-    )
-    with pytest.raises(TypeError):
-        reproject_geometry(big_box, 1.0, 1.0)
-
-
-def test_reproject_geometry_clip_crs_bounds_epsg():
-    bbox = wkt.loads(
-        "Polygon ((6.05416952699480682 49.79497943046440867, 6.04100166381764581 50.01055350300158864, 5.70657677854139056 50.00153486687963778, 5.72122668311700089 49.78602894452072292, 6.05416952699480682 49.79497943046440867))"
-    )
-    dst_crs = "EPSG:32632"
-
-    # reproject to UTM
-    bbox_utm = reproject_geometry(bbox, 4326, dst_crs, clip_to_crs_bounds=True)
-    assert bbox_utm.is_valid
-    assert bbox_utm.area
-    # revert to WGS84 and test geometry clipping
-    bbox_wgs84 = reproject_geometry(bbox_utm, dst_crs, 4326)
-    assert bbox_wgs84.area < bbox.area
-
-    # reproject to UTM but don't clip
-    bbox_utm = reproject_geometry(bbox, 4326, dst_crs, clip_to_crs_bounds=False)
-    assert bbox_utm.is_valid
-    assert bbox_utm.area
-    # make sure geometry was not clipped
-    bbox_wgs84 = reproject_geometry(bbox_utm, dst_crs, 4326)
-    assert bbox_wgs84.intersects(bbox)
-    assert (
-        bbox_wgs84.intersection(bbox).area
-        == pytest.approx(bbox_wgs84.area)
-        == pytest.approx(bbox.area)
-    )
-
-    # reproject to UTM don't clip but segmentize
-    bbox_utm = reproject_geometry(
-        bbox, 4326, dst_crs, clip_to_crs_bounds=False, segmentize=True
-    )
-    assert bbox_utm.is_valid
-    assert bbox_utm.area
-    # make sure geometry was not clipped
-    bbox_wgs84 = reproject_geometry(bbox_utm, dst_crs, 4326)
-    assert bbox_wgs84.intersects(bbox)
-    assert (
-        bbox_wgs84.intersection(bbox).area
-        == pytest.approx(bbox_wgs84.area)
-        == pytest.approx(bbox.area)
-    )
-
-
-def test_reproject_geometry_clip_crs_bounds_proj():
-    bbox = wkt.loads(
-        "Polygon ((6.05416952699480682 49.79497943046440867, 6.04100166381764581 50.01055350300158864, 5.70657677854139056 50.00153486687963778, 5.72122668311700089 49.78602894452072292, 6.05416952699480682 49.79497943046440867))"
-    )
-    dst_crs = "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs"
-
-    # reproject to UTM
-    bbox_utm = reproject_geometry(bbox, 4326, dst_crs, clip_to_crs_bounds=True)
-    assert bbox_utm.is_valid
-    assert bbox_utm.area
-    # revert to WGS84 and test geometry clipping
-    bbox_wgs84 = reproject_geometry(bbox_utm, dst_crs, 4326)
-    # NOTE: on some proj versions (TBD), pyproj cannot detect the CRS bounds of a CRS passed on by a proj string
-    # assert bbox_wgs84.area < bbox.area
-
-    # reproject to UTM but don't clip
-    bbox_utm = reproject_geometry(bbox, 4326, dst_crs, clip_to_crs_bounds=False)
-    assert bbox_utm.is_valid
-    assert bbox_utm.area
-    # make sure geometry was not clipped
-    bbox_wgs84 = reproject_geometry(bbox_utm, dst_crs, 4326)
-    assert bbox_wgs84.intersects(bbox)
-    assert (
-        bbox_wgs84.intersection(bbox).area
-        == pytest.approx(bbox_wgs84.area)
-        == pytest.approx(bbox.area)
-    )
-
-    # reproject to UTM don't clip but segmentize
-    bbox_utm = reproject_geometry(
-        bbox, 4326, dst_crs, clip_to_crs_bounds=False, segmentize=True
-    )
-    assert bbox_utm.is_valid
-    assert bbox_utm.area
-    # make sure geometry was not clipped
-    bbox_wgs84 = reproject_geometry(bbox_utm, dst_crs, 4326)
-    assert bbox_wgs84.intersects(bbox)
-    assert (
-        bbox_wgs84.intersection(bbox).area
-        == pytest.approx(bbox_wgs84.area)
-        == pytest.approx(bbox.area)
-    )
-
-
-def test_repair_geometry():
-    # invalid LineString
-    l = LineString([(0, 0), (0, 0), (0, 0)])
-    with pytest.raises(TopologicalError):
-        _repair(l)
-
-
-def test_write_vector_window_errors(landpoly):
-    with fiona.open(landpoly) as src:
-        feature = next(iter(src))
-    with pytest.raises((DriverError, ValueError, TypeError)):
-        write_vector_window(
-            in_data=["invalid", feature],
-            out_tile=BufferedTilePyramid("geodetic").tile(0, 0, 0),
-            out_path="/invalid_path",
-            out_schema=dict(geometry="Polygon", properties=dict()),
-        )
-
-
-def test_segmentize_geometry():
-    """Segmentize function."""
-    # Polygon
-    polygon = box(-18, -9, 18, 9)
-    out = segmentize_geometry(polygon, 1)
-    assert out.is_valid
-    # wrong type
-    with pytest.raises(TypeError):
-        segmentize_geometry(polygon.centroid, 1)
-
-
-def test_clean_geometry_type(geometrycollection):
-    """Filter and break up geometries."""
-    polygon = box(-18, -9, 18, 9)
-    # invalid type
-    with pytest.raises(TypeError):
-        clean_geometry_type(polygon, "invalid_type")
-
-    # don't return geometry
-    with pytest.raises(GeometryTypeError):
-        clean_geometry_type(polygon, "LineString", raise_exception=True)
-
-    # return geometry as is
-    assert clean_geometry_type(polygon, "Polygon").geom_type == "Polygon"
-    assert clean_geometry_type(polygon, "MultiPolygon").geom_type == "Polygon"
-
-    # don't allow multipart geometries
-    with pytest.raises(GeometryTypeError):
-        clean_geometry_type(
-            MultiPolygon([polygon]),
-            "Polygon",
-            allow_multipart=False,
-            raise_exception=True,
-        )
-
-    # multipolygons from geometrycollection
-    result = clean_geometry_type(
-        geometrycollection, "Polygon", allow_multipart=True, raise_exception=False
-    )
-    assert result.geom_type == "MultiPolygon"
-    assert not result.is_empty
-
-    # polygons from geometrycollection
-    result = clean_geometry_type(
-        geometrycollection, "Polygon", allow_multipart=False, raise_exception=False
-    )
-    assert result.geom_type == "GeometryCollection"
-    assert result.is_empty
-
-
 @pytest.mark.remote
 def test_s3_path_exists(s2_band_remote):
     assert path_exists(s2_band_remote)
 
 
-@pytest.mark.remote
-def test_s3_read_raster_window(s2_band_remote):
-    tile = BufferedTilePyramid("geodetic").tile(10, 276, 1071)
-    assert read_raster_window(s2_band_remote, tile).any()
-
-
 @pytest.mark.remote
 def test_remote_path_exists(http_raster):
     assert path_exists(http_raster)
@@ -1286,6 +274,24 @@ def test_tiles_exist_s3(gtiff_s3):
         assert set(output_tiles) == existing.union(not_existing)
 
 
+def test_copy(cleantopo_br_tif, tmpdir):
+    out = os.path.join(tmpdir, "copied.tif")
+
+    # copy and verify file is valid
+    copy(cleantopo_br_tif, out)
+    with rasterio.open(out) as src:
+        assert not src.read(masked=True).mask.all()
+
+    # try to copy again, catching the IOError
+    with pytest.raises(IOError):
+        copy(cleantopo_br_tif, out)
+
+    # copy again but overwrite
+    copy(cleantopo_br_tif, out, overwrite=True)
+    with rasterio.open(out) as src:
+        assert not src.read(masked=True).mask.all()
+
+
 def test_custom_grid_points(custom_grid_points):
     mp = custom_grid_points.process_mp(tile=(3, 1245, 37))
     with mp.open("inp") as points:
diff --git a/test/test_io_raster.py b/test/test_io_raster.py
new file mode 100644
index 00000000..69e0210e
--- /dev/null
+++ b/test/test_io_raster.py
@@ -0,0 +1,805 @@
+import pytest
+import shutil
+import rasterio
+import tempfile
+import numpy as np
+import numpy.ma as ma
+import os
+from rasterio.enums import Compression
+from shapely.geometry import box
+from shapely.ops import unary_union
+from tilematrix import Bounds
+from itertools import product
+
+
+import mapchete
+from mapchete.config import MapcheteConfig
+from mapchete.errors import MapcheteIOError
+from mapchete.io import path_exists
+from mapchete.io.raster import (
+    read_raster_window,
+    write_raster_window,
+    extract_from_array,
+    resample_from_array,
+    create_mosaic,
+    ReferencedRaster,
+    prepare_array,
+    RasterWindowMemoryFile,
+    read_raster_no_crs,
+    convert_raster,
+)
+from mapchete.tile import BufferedTilePyramid
+
+
+def test_read_raster_window(dummy1_tif, minmax_zoom):
+    """Read array with read_raster_window."""
+    zoom = 8
+    # without reproject
+    config = MapcheteConfig(minmax_zoom.dict)
+    rasterfile = config.params_at_zoom(zoom)["input"]["file1"]
+    dummy1_bbox = rasterfile.bbox()
+
+    pixelbuffer = 5
+    tile_pyramid = BufferedTilePyramid("geodetic", pixelbuffer=pixelbuffer)
+    tiles = list(tile_pyramid.tiles_from_geom(dummy1_bbox, zoom))
+    # add edge tile
+    tiles.append(tile_pyramid.tile(8, 0, 0))
+    for tile in tiles:
+        width, height = tile.shape
+        for band in read_raster_window(dummy1_tif, tile):
+            assert isinstance(band, ma.MaskedArray)
+            assert band.shape == (width, height)
+        for index in range(1, 4):
+            band = read_raster_window(dummy1_tif, tile, index)
+            assert isinstance(band, ma.MaskedArray)
+            assert band.shape == (width, height)
+        for index in [None, [1, 2, 3]]:
+            band = read_raster_window(dummy1_tif, tile, index)
+            assert isinstance(band, ma.MaskedArray)
+            assert band.ndim == 3
+            assert band.shape == (3, width, height)
+
+
+def test_read_raster_window_reproject(dummy1_3857_tif, minmax_zoom):
+    """Read array with read_raster_window."""
+    zoom = 8
+    # with reproject
+    config_raw = minmax_zoom.dict
+    config_raw["input"].update(file1=dummy1_3857_tif)
+    config = MapcheteConfig(config_raw)
+    rasterfile = config.params_at_zoom(zoom)["input"]["file1"]
+    dummy1_bbox = rasterfile.bbox()
+
+    pixelbuffer = 5
+    tile_pyramid = BufferedTilePyramid("geodetic", pixelbuffer=pixelbuffer)
+    tiles = list(tile_pyramid.tiles_from_geom(dummy1_bbox, zoom))
+    # target window out of CRS bounds
+    band = read_raster_window(dummy1_3857_tif, tile_pyramid.tile(12, 0, 0))
+    assert isinstance(band, ma.MaskedArray)
+    assert band.mask.all()
+    # not intersecting tile
+    tiles.append(tile_pyramid.tile(zoom, 1, 1))  # out of CRS bounds
+    tiles.append(tile_pyramid.tile(zoom, 16, 1))  # out of file bbox
+    for tile in tiles:
+        for band in read_raster_window(dummy1_3857_tif, tile):
+            assert isinstance(band, ma.MaskedArray)
+            assert band.shape == tile.shape
+        bands = read_raster_window(dummy1_3857_tif, tile, [1])
+        assert isinstance(bands, ma.MaskedArray)
+        assert bands.shape == tile.shape
+    # errors
+    with pytest.raises(IOError):
+        read_raster_window("nonexisting_path", tile)
+
+
+def test_read_raster_window_resampling(cleantopo_br_tif):
+    """Assert various resampling options work."""
+    tp = BufferedTilePyramid("geodetic")
+    with rasterio.open(cleantopo_br_tif, "r") as src:
+        tiles = tp.tiles_from_bounds(src.bounds, 4)
+    for tile in tiles:
+        outputs = [
+            read_raster_window(cleantopo_br_tif, tile, resampling=resampling)
+            for resampling in [
+                "nearest",
+                "bilinear",
+                "cubic",
+                "cubic_spline",
+                "lanczos",
+                "average",
+                "mode",
+            ]
+        ]
+        # resampling test:
+        assert any(
+            [not np.array_equal(w, v) for v, w in zip(outputs[:-1], outputs[1:])]
+        )
+
+
+def test_read_raster_window_partly_overlapping(cleantopo_br_tif):
+    """Read array with read_raster_window where window is bigger than file."""
+    tile = BufferedTilePyramid("geodetic").tile(4, 15, 31)
+    data = read_raster_window(cleantopo_br_tif, tile)
+    assert isinstance(data, ma.MaskedArray)
+    assert data.mask.any()
+
+
+def test_read_raster_window_mask(s2_band):
+    """No resampling artefacts on mask edges."""
+    tile = BufferedTilePyramid("geodetic").tile(zoom=13, row=1918, col=8905)
+    data = read_raster_window(
+        s2_band, tile, resampling="cubic", src_nodata=0, dst_nodata=0
+    )
+    assert data.any()
+    assert not np.where(data == 1, True, False).any()
+
+
+def test_read_raster_window_input_list(cleantopo_br):
+    process_zoom = 5
+    conf = dict(**cleantopo_br.dict)
+    conf["output"].update(metatiling=1)
+    with mapchete.open(conf) as mp:
+        mp.batch_process(process_zoom)
+        tiles = [
+            (tile, mp.config.output.get_path(tile))
+            for tile in mp.config.output_pyramid.tiles_from_bounds(
+                mp.config.bounds, process_zoom
+            )
+            if path_exists(mp.config.output.get_path(tile))
+        ]
+        upper_tile = next(mp.get_process_tiles(process_zoom - 1))
+        assert len(tiles) > 1
+        resampled = resample_from_array(
+            in_raster=create_mosaic(
+                [(tile, read_raster_window(path, tile)) for tile, path in tiles]
+            ),
+            out_tile=upper_tile,
+        )
+    resampled2 = read_raster_window(
+        [p for _, p in tiles], upper_tile, src_nodata=0, dst_nodata=0
+    )
+    assert resampled.dtype == resampled2.dtype
+    assert resampled.shape == resampled2.shape
+    assert np.array_equal(resampled.mask, resampled2.mask)
+    # TODO slight rounding errors occur
+    assert np.allclose(resampled, resampled2, rtol=0.01)
+
+
+def test_read_raster_window_retry(invalid_tif):
+    tile = BufferedTilePyramid("geodetic").tile(zoom=13, row=1918, col=8905)
+    with pytest.raises(MapcheteIOError):
+        read_raster_window(invalid_tif, tile)
+    with pytest.raises(FileNotFoundError):
+        read_raster_window("not_existing.tif", tile)
+
+
+def test_read_raster_no_crs_errors():
+    with tempfile.NamedTemporaryFile() as tmpfile:
+        with pytest.raises(MapcheteIOError):
+            read_raster_no_crs(tmpfile.name)
+
+
+def test_write_raster_window():
+    """Basic output format writing."""
+    path = tempfile.NamedTemporaryFile(delete=False).name
+    # standard tile
+    tp = BufferedTilePyramid("geodetic")
+    tile = tp.tile(5, 5, 5)
+    data = ma.masked_array(np.ones((2,) + tile.shape))
+    for out_profile in [
+        dict(
+            driver="GTiff",
+            count=2,
+            dtype="uint8",
+            compress="lzw",
+            nodata=0,
+            height=tile.height,
+            width=tile.width,
+            affine=tile.affine,
+        ),
+        dict(
+            driver="GTiff",
+            count=2,
+            dtype="uint8",
+            compress="deflate",
+            nodata=0,
+            height=tile.height,
+            width=tile.width,
+            affine=tile.affine,
+        ),
+        dict(
+            driver="PNG",
+            count=2,
+            dtype="uint8",
+            nodata=0,
+            height=tile.height,
+            width=tile.width,
+            compress=None,
+            affine=tile.affine,
+        ),
+    ]:
+        try:
+            write_raster_window(
+                in_tile=tile, in_data=data, out_profile=out_profile, out_path=path
+            )
+            with rasterio.open(path, "r") as src:
+                assert src.read().any()
+                assert src.meta["driver"] == out_profile["driver"]
+                assert src.transform == tile.affine
+                if out_profile["compress"]:
+                    assert src.compression == Compression(
+                        out_profile["compress"].upper()
+                    )
+        finally:
+            shutil.rmtree(path, ignore_errors=True)
+    # with metatiling
+    tile = BufferedTilePyramid("geodetic", metatiling=4).tile(5, 1, 1)
+    data = ma.masked_array(np.ones((2,) + tile.shape))
+    out_tile = BufferedTilePyramid("geodetic").tile(5, 5, 5)
+    out_profile = dict(
+        driver="GTiff",
+        count=2,
+        dtype="uint8",
+        compress="lzw",
+        nodata=0,
+        height=out_tile.height,
+        width=out_tile.width,
+        affine=out_tile.affine,
+    )
+    try:
+        write_raster_window(
+            in_tile=tile,
+            in_data=data,
+            out_profile=out_profile,
+            out_tile=out_tile,
+            out_path=path,
+        )
+        with rasterio.open(path, "r") as src:
+            assert src.shape == out_tile.shape
+            assert src.read().any()
+            assert src.meta["driver"] == out_profile["driver"]
+            assert src.transform == out_profile["transform"]
+    finally:
+        shutil.rmtree(path, ignore_errors=True)
+
+
+def test_write_raster_window_memory():
+    """Basic output format writing."""
+    path = "memoryfile"
+    # standard tile
+    tp = BufferedTilePyramid("geodetic")
+    tile = tp.tile(5, 5, 5)
+    data = ma.masked_array(np.ones((2,) + tile.shape))
+    for out_profile in [
+        dict(
+            driver="GTiff",
+            count=2,
+            dtype="uint8",
+            compress="lzw",
+            nodata=0,
+            height=tile.height,
+            width=tile.width,
+            affine=tile.affine,
+        ),
+        dict(
+            driver="GTiff",
+            count=2,
+            dtype="uint8",
+            compress="deflate",
+            nodata=0,
+            height=tile.height,
+            width=tile.width,
+            affine=tile.affine,
+        ),
+        dict(
+            driver="PNG",
+            count=2,
+            dtype="uint8",
+            nodata=0,
+            height=tile.height,
+            width=tile.width,
+            compress=None,
+            affine=tile.affine,
+        ),
+    ]:
+        with pytest.raises(DeprecationWarning):
+            write_raster_window(
+                in_tile=tile, in_data=data, out_profile=out_profile, out_path=path
+            )
+
+
+def test_raster_window_memoryfile():
+    """Use context manager for rasterio MemoryFile."""
+    tp = BufferedTilePyramid("geodetic")
+    tile = tp.tile(5, 5, 5)
+    data = ma.masked_array(np.ones((2,) + tile.shape))
+    for out_profile in [
+        dict(
+            driver="GTiff",
+            count=2,
+            dtype="uint8",
+            compress="lzw",
+            nodata=0,
+            height=tile.height,
+            width=tile.width,
+            affine=tile.affine,
+        ),
+        dict(
+            driver="GTiff",
+            count=2,
+            dtype="uint8",
+            compress="deflate",
+            nodata=0,
+            height=tile.height,
+            width=tile.width,
+            affine=tile.affine,
+        ),
+        dict(
+            driver="PNG",
+            count=2,
+            dtype="uint8",
+            nodata=0,
+            height=tile.height,
+            width=tile.width,
+            compress=None,
+            affine=tile.affine,
+        ),
+    ]:
+        with RasterWindowMemoryFile(
+            in_tile=tile, in_data=data, out_profile=out_profile
+        ) as memfile:
+            with memfile.open() as src:
+                assert src.read().any()
+                assert src.meta["driver"] == out_profile["driver"]
+                assert src.transform == tile.affine
+                if out_profile["compress"]:
+                    assert src.compression == Compression(
+                        out_profile["compress"].upper()
+                    )
+
+
+def test_write_raster_window_errors():
+    """Basic output format writing."""
+    tile = BufferedTilePyramid("geodetic").tile(5, 5, 5)
+    data = ma.masked_array(np.ndarray((1, 1)))
+    profile = {}
+    path = ""
+    # in_tile
+    with pytest.raises(TypeError):
+        write_raster_window(
+            in_tile="invalid tile",
+            in_data=data,
+            out_profile=profile,
+            out_tile=tile,
+            out_path=path,
+        )
+    # out_tile
+    with pytest.raises(TypeError):
+        write_raster_window(
+            in_tile=tile,
+            in_data=data,
+            out_profile=profile,
+            out_tile="invalid tile",
+            out_path=path,
+        )
+    # in_data
+    with pytest.raises(TypeError):
+        write_raster_window(
+            in_tile=tile,
+            in_data="invalid data",
+            out_profile=profile,
+            out_tile=tile,
+            out_path=path,
+        )
+    # out_profile
+    with pytest.raises(TypeError):
+        write_raster_window(
+            in_tile=tile,
+            in_data=data,
+            out_profile="invalid profile",
+            out_tile=tile,
+            out_path=path,
+        )
+    # out_path
+    with pytest.raises(TypeError):
+        write_raster_window(
+            in_tile=tile, in_data=data, out_profile=profile, out_tile=tile, out_path=999
+        )
+    # cannot write
+    with pytest.raises(ValueError):
+        write_raster_window(
+            in_tile=tile,
+            in_data=data,
+            out_profile=profile,
+            out_tile=tile,
+            out_path="/invalid_path",
+        )
+
+
+def test_extract_from_array():
+    """Extract subdata from array."""
+    in_tile = BufferedTilePyramid("geodetic", metatiling=4).tile(5, 5, 5)
+    shape = (in_tile.shape[0] // 2, in_tile.shape[1])
+    data = ma.masked_array(np.concatenate([np.ones(shape), np.ones(shape) * 2]))
+    # intersecting at top
+    out_tile = BufferedTilePyramid("geodetic").tile(5, 20, 20)
+    out_array = extract_from_array(
+        in_raster=data, in_affine=in_tile.affine, out_tile=out_tile
+    )
+    assert isinstance(out_array, np.ndarray)
+    assert np.all(np.where(out_array == 1, True, False))
+    # intersecting at bottom
+    out_tile = BufferedTilePyramid("geodetic").tile(5, 22, 20)
+    out_array = extract_from_array(
+        in_raster=data, in_affine=in_tile.affine, out_tile=out_tile
+    )
+    assert isinstance(out_array, np.ndarray)
+    assert np.all(np.where(out_array == 2, True, False))
+    # not intersecting
+    out_tile = BufferedTilePyramid("geodetic").tile(5, 15, 20)
+    with pytest.raises(ValueError):
+        out_array = extract_from_array(
+            in_raster=data, in_affine=in_tile.affine, out_tile=out_tile
+        )
+
+
+def test_resample_from_array():
+    """Resample array using rasterio reproject()."""
+    in_tile = BufferedTilePyramid("geodetic").tile(5, 5, 5)
+    in_data = np.ones(in_tile.shape)
+    # tile from next toom level
+    out_tile = BufferedTilePyramid("geodetic").tile(6, 10, 10)
+    out_array = resample_from_array(in_data, in_tile.affine, out_tile)
+    assert isinstance(out_array, ma.masked_array)
+    assert np.all(np.where(out_array == 1, True, False))
+    # not intersecting tile
+    out_tile = BufferedTilePyramid("geodetic").tile(7, 0, 0)
+    out_array = resample_from_array(in_data, in_tile.affine, out_tile)
+    assert isinstance(out_array, ma.masked_array)
+    assert out_array.mask.all()
+    # data as tuple
+    in_data = (np.ones(in_tile.shape[1:]),)
+    out_tile = BufferedTilePyramid("geodetic").tile(6, 10, 10)
+    out_array = resample_from_array(in_data, in_tile.affine, out_tile)
+    # deprecated
+    resample_from_array(in_data, in_tile.affine, out_tile, nodata=-9999)
+    # errors
+    with pytest.raises(TypeError):
+        in_data = "invalid_type"
+        resample_from_array(in_data, in_tile.affine, out_tile)
+    with pytest.raises(TypeError):
+        in_data = np.ones(in_tile.shape[0])
+        resample_from_array(in_data, in_tile.affine, out_tile)
+
+
+def test_create_mosaic_errors():
+    """Check error handling of create_mosaic()."""
+    tp_geo = BufferedTilePyramid("geodetic")
+    tp_mer = BufferedTilePyramid("mercator")
+    geo_tile = tp_geo.tile(1, 0, 0)
+    geo_tile_data = np.ndarray(geo_tile.shape)
+    mer_tile = tp_mer.tile(1, 1, 0)
+    mer_tile_data = np.ndarray(mer_tile.shape)
+    # tiles error
+    with pytest.raises(TypeError):
+        create_mosaic("invalid tiles")
+    with pytest.raises(TypeError):
+        create_mosaic(["invalid tiles"])
+    # CRS error
+    with pytest.raises(ValueError):
+        create_mosaic([(geo_tile, geo_tile_data), (mer_tile, mer_tile_data)])
+    # zoom error
+    with pytest.raises(ValueError):
+        diff_zoom = tp_geo.tile(2, 1, 0)
+        diff_zoom_data = np.ndarray(diff_zoom.shape)
+        create_mosaic([(geo_tile, geo_tile_data), (diff_zoom, diff_zoom_data)])
+    # tile data error
+    with pytest.raises(TypeError):
+        # for one tile
+        create_mosaic([(geo_tile, None)])
+    with pytest.raises(TypeError):
+        # for multiple tiles
+        create_mosaic([(geo_tile, None), (geo_tile, None)])
+    # tile data type error
+    with pytest.raises(TypeError):
+        diff_type = tp_geo.tile(1, 1, 0)
+        diff_type_data = np.ndarray(diff_zoom.shape).astype("int")
+        create_mosaic([(geo_tile, geo_tile_data), (diff_type, diff_type_data)])
+    # no tiles
+    with pytest.raises(ValueError):
+        create_mosaic(tiles=[])
+
+
+def test_create_mosaic():
+    """Create mosaic from tiles."""
+    tp = BufferedTilePyramid("geodetic")
+    # quick return mosaic if there is just one tile
+    tile = tp.tile(3, 3, 3)
+    data = np.ones(tile.shape)
+    mosaic = create_mosaic([(tile, data)])
+    assert isinstance(mosaic, ReferencedRaster)
+    assert np.array_equal(data, mosaic.data)
+    assert tile.affine == mosaic.affine
+    zoom = 5
+    # multiple tiles on top left corner of tile matrix
+    for pixelbuffer in [0, 10]:
+        tp = BufferedTilePyramid("geodetic", pixelbuffer=pixelbuffer)
+        tiles = [
+            (tp.tile(zoom, row, col), np.ones(tp.tile(zoom, row, col).shape))
+            for row, col in product(range(4), range(4))
+        ]
+        # 4x4 top left tiles from zoom 5 equal top left tile from zoom 3
+        # also use tile generator
+        mosaic = create_mosaic((t for t in tiles))
+        assert isinstance(mosaic, ReferencedRaster)
+        assert np.all(np.where(mosaic.data == 1, True, False))
+        mosaic_bbox = box(
+            mosaic.affine[2],
+            mosaic.affine[5] + mosaic.data.shape[1] * mosaic.affine[4],
+            mosaic.affine[2] + mosaic.data.shape[2] * mosaic.affine[0],
+            mosaic.affine[5],
+        )
+        control_bbox = box(*unary_union([t.bbox for t, _ in tiles]).bounds)
+        assert mosaic_bbox.equals(control_bbox)
+    # multiple tiles on bottom right corner of tile matrix
+    for pixelbuffer in [0, 10]:
+        tp = BufferedTilePyramid("geodetic", pixelbuffer=pixelbuffer)
+        tiles = [
+            (tp.tile(zoom, row, col), np.ones(tp.tile(zoom, row, col).shape))
+            for row, col in product(
+                range(tp.matrix_height(zoom) - 4, tp.matrix_height(zoom)),
+                range(tp.matrix_width(zoom) - 4, tp.matrix_width(zoom)),
+            )
+        ]
+        # 4x4 top left tiles from zoom 5 equal top left tile from zoom 3
+        # also use tile generator
+        mosaic = create_mosaic((t for t in tiles))
+        assert isinstance(mosaic, ReferencedRaster)
+        assert np.all(np.where(mosaic.data == 1, True, False))
+        mosaic_bbox = box(
+            mosaic.affine[2],
+            mosaic.affine[5] + mosaic.data.shape[1] * mosaic.affine[4],
+            mosaic.affine[2] + mosaic.data.shape[2] * mosaic.affine[0],
+            mosaic.affine[5],
+        )
+        control_bbox = box(*unary_union([t.bbox for t, _ in tiles]).bounds)
+        assert mosaic_bbox.equals(control_bbox)
+
+
+def test_create_mosaic_antimeridian():
+    """Create mosaic using tiles on opposing antimeridian sides."""
+    zoom = 5
+    row = 0
+    pixelbuffer = 5
+    tp = BufferedTilePyramid("geodetic", pixelbuffer=pixelbuffer)
+    west = tp.tile(zoom, row, 0)
+    east = tp.tile(zoom, row, tp.matrix_width(zoom) - 1)
+    mosaic = create_mosaic(
+        [
+            (west, np.ones(west.shape).astype("uint8")),
+            (east, np.ones(east.shape).astype("uint8") * 2),
+        ]
+    )
+    assert isinstance(mosaic, ReferencedRaster)
+
+    # Huge array gets initialized because the two tiles are on opposing sides of the
+    # projection area. The below test should pass if the tiles are stitched together next
+    # to each other.
+    assert mosaic.data.shape == (1, west.height, west.width * 2 - 2 * pixelbuffer)
+    assert mosaic.data[0][0][0] == 2
+    assert mosaic.data[0][0][-1] == 1
+
+    # If tiles from opposing sides from Antimeridian are mosaicked it will happen that the
+    # output mosaic exceeds the CRS bounds (obviously). In such a case the mosaicking
+    # function shall make sure that the larger part of the output mosaic shall be inside
+    # the CRS bounds.
+
+    # (1) mosaic crosses Antimeridian in the West, larger part is on Western hemisphere:
+    tiles_ids = [
+        # Western hemisphere tiles
+        (zoom, row, 0),
+        (zoom, row, 1),
+        # Eastern hemisphere tile
+        (zoom, row, tp.matrix_width(zoom) - 1),
+    ]
+    tiles = [
+        (tp.tile(*tile_id), np.ones(tp.tile(*tile_id).shape)) for tile_id in tiles_ids
+    ]
+    mosaic = create_mosaic(tiles)
+    control_bounds = Bounds(
+        # Eastern tile has to be shifted
+        -(360 - tp.tile(*tiles_ids[2]).left),
+        tp.tile(*tiles_ids[2]).bottom,
+        tp.tile(*tiles_ids[1]).right,
+        tp.tile(*tiles_ids[1]).top,
+    )
+    assert mosaic.bounds == control_bounds
+
+    # (2) mosaic crosses Antimeridian in the West, larger part is on Eastern hemisphere:
+    tiles_ids = [
+        # Western hemisphere tile
+        (zoom, row, 0),
+        # Eastern hemisphere tiles
+        (zoom, row, tp.matrix_width(zoom) - 1),
+        (zoom, row, tp.matrix_width(zoom) - 2),
+    ]
+    tiles = [
+        (tp.tile(*tile_id), np.ones(tp.tile(*tile_id).shape)) for tile_id in tiles_ids
+    ]
+    mosaic = create_mosaic(tiles)
+    control_bounds = Bounds(
+        tp.tile(*tiles_ids[2]).left,
+        tp.tile(*tiles_ids[2]).bottom,
+        # Western tile has to be shifted
+        360 + tp.tile(*tiles_ids[0]).right,
+        tp.tile(*tiles_ids[0]).top,
+    )
+    assert mosaic.bounds == control_bounds
+
+
+def test_prepare_array_iterables():
+    """Convert iterable data into a proper array."""
+    # input is iterable
+    # iterable contains arrays
+    data = [np.zeros((1, 1))]
+    # output ndarray
+    output = prepare_array(data, masked=False)
+    assert isinstance(output, np.ndarray)
+    assert not isinstance(output, ma.masked_array)
+    assert output.shape == (1, 1, 1)
+    # output masked array
+    output = prepare_array(data)
+    assert isinstance(output, ma.masked_array)
+    assert output.shape == (1, 1, 1)
+    # iterable contains masked arrays
+    data = [ma.empty((1, 1))]
+    output = prepare_array(data, masked=False)
+    assert isinstance(output, np.ndarray)
+    assert not isinstance(output, ma.masked_array)
+    assert output.shape == (1, 1, 1)
+    # output masked array
+    output = prepare_array(data)
+    assert isinstance(output, ma.masked_array)
+    assert output.shape == (1, 1, 1)
+    # iterable contains masked arrays with full mask
+    data = [ma.masked_array(data=np.ones((1, 1)), mask=np.ones((1, 1)))]
+    output = prepare_array(data, masked=False)
+    assert isinstance(output, np.ndarray)
+    assert not isinstance(output, ma.masked_array)
+    assert output.shape == (1, 1, 1)
+    # output masked array
+    output = prepare_array(data)
+    assert isinstance(output, ma.masked_array)
+    assert output.shape == (1, 1, 1)
+
+
+def test_prepare_array_maskedarrays():
+    """Convert masked array data into a proper array."""
+    # input is ma.masked_array
+    data = ma.empty((1, 1, 1))
+    # output ndarray
+    output = prepare_array(data, masked=False)
+    assert isinstance(output, np.ndarray)
+    assert not isinstance(output, ma.masked_array)
+    assert output.shape == (1, 1, 1)
+    # output masked array
+    output = prepare_array(data)
+    assert isinstance(output, ma.masked_array)
+    assert output.shape == (1, 1, 1)
+    # input is ma.masked_array with full mask
+    data = ma.masked_array(data=np.ones((1, 1, 1)), mask=np.ones((1, 1, 1)))
+    # output ndarray
+    output = prepare_array(data, masked=False)
+    assert isinstance(output, np.ndarray)
+    assert not isinstance(output, ma.masked_array)
+    assert output.shape == (1, 1, 1)
+    # output masked array
+    output = prepare_array(data)
+    assert isinstance(output, ma.masked_array)
+    assert output.shape == (1, 1, 1)
+
+
+def test_prepare_array_ndarrays():
+    """Convert ndarray data into a proper array."""
+    # input is np.ndarray
+    data = np.zeros((1, 1, 1))
+    # output ndarray
+    output = prepare_array(data, masked=False)
+    assert isinstance(output, np.ndarray)
+    assert not isinstance(output, ma.masked_array)
+    assert output.shape == (1, 1, 1)
+    # output masked array
+    output = prepare_array(data)
+    assert isinstance(output, ma.masked_array)
+    assert output.shape == (1, 1, 1)
+    # input is 2D np.ndarray
+    data = np.zeros((1, 1))
+    # output ndarray
+    output = prepare_array(data, masked=False)
+    assert isinstance(output, np.ndarray)
+    assert not isinstance(output, ma.masked_array)
+    assert output.shape == (1, 1, 1)
+    # output masked array
+    output = prepare_array(data)
+    assert isinstance(output, ma.masked_array)
+    assert output.shape == (1, 1, 1)
+
+
+def test_prepare_array_errors():
+    """Convert ndarray data into a proper array."""
+    # input is iterable
+    data = [None]
+    try:
+        prepare_array(data)
+        raise Exception()
+    except ValueError:
+        pass
+    # input is not array
+    data = 5
+    try:
+        prepare_array(data)
+        raise Exception()
+    except ValueError:
+        pass
+
+
+@pytest.mark.remote
+def test_s3_read_raster_window(s2_band_remote):
+    tile = BufferedTilePyramid("geodetic").tile(10, 276, 1071)
+    assert read_raster_window(s2_band_remote, tile).any()
+
+
+def test_convert_raster_copy(cleantopo_br_tif, tmpdir):
+    out = os.path.join(tmpdir, "copied.tif")
+
+    # copy
+    convert_raster(cleantopo_br_tif, out)
+    with rasterio.open(out) as src:
+        assert not src.read(masked=True).mask.all()
+
+    # raise error if output exists
+    with pytest.raises(IOError):
+        convert_raster(cleantopo_br_tif, out, exists_ok=False)
+
+    # do nothing if output exists
+    convert_raster(cleantopo_br_tif, out)
+    with rasterio.open(out) as src:
+        assert not src.read(masked=True).mask.all()
+
+
+def test_convert_raster_overwrite(cleantopo_br_tif, tmpdir):
+    out = os.path.join(tmpdir, "copied.tif")
+
+    # write an invalid file
+    with open(out, "w") as dst:
+        dst.write("invalid")
+
+    # overwrite
+    convert_raster(cleantopo_br_tif, out, overwrite=True)
+    with rasterio.open(out) as src:
+        assert not src.read(masked=True).mask.all()
+
+
+def test_convert_raster_other_format_copy(cleantopo_br_tif, tmpdir):
+    out = os.path.join(tmpdir, "copied.jp2")
+
+    convert_raster(cleantopo_br_tif, out, driver="JP2OpenJPEG")
+    with rasterio.open(out) as src:
+        assert not src.read(masked=True).mask.all()
+
+    # raise error if output exists
+    with pytest.raises(IOError):
+        convert_raster(cleantopo_br_tif, out, exists_ok=False)
+
+
+def test_convert_raster_other_format_overwrite(cleantopo_br_tif, tmpdir):
+    out = os.path.join(tmpdir, "copied.jp2")
+
+    # write an invalid file
+    with open(out, "w") as dst:
+        dst.write("invalid")
+
+    # overwrite
+    convert_raster(cleantopo_br_tif, out, driver="JP2OpenJPEG", overwrite=True)
+    with rasterio.open(out) as src:
+        assert not src.read(masked=True).mask.all()
diff --git a/test/test_io_vector.py b/test/test_io_vector.py
new file mode 100644
index 00000000..88696b4b
--- /dev/null
+++ b/test/test_io_vector.py
@@ -0,0 +1,494 @@
+import pytest
+import fiona
+from fiona.errors import DriverError
+import os
+from rasterio.crs import CRS
+from shapely import wkt
+from shapely.errors import TopologicalError
+from shapely.geometry import shape, box, Polygon, MultiPolygon, LineString, mapping
+
+from mapchete.config import MapcheteConfig
+from mapchete.errors import GeometryTypeError, MapcheteIOError
+from mapchete.io.vector import (
+    read_vector_window,
+    reproject_geometry,
+    clean_geometry_type,
+    segmentize_geometry,
+    write_vector_window,
+    _repair,
+    convert_vector,
+    IndexedFeatures,
+    bounds_intersect,
+)
+from mapchete.tile import BufferedTilePyramid
+
+
+def test_read_vector_window(geojson, landpoly_3857):
+    """Read vector data from read_vector_window."""
+    zoom = 4
+    config = MapcheteConfig(geojson.dict)
+    vectorfile = config.params_at_zoom(zoom)["input"]["file1"]
+    pixelbuffer = 5
+    tile_pyramid = BufferedTilePyramid("geodetic", pixelbuffer=pixelbuffer)
+    tiles = tile_pyramid.tiles_from_geom(vectorfile.bbox(), zoom)
+    feature_count = 0
+    for tile in tiles:
+        for feature in read_vector_window(vectorfile.path, tile):
+            assert "properties" in feature
+            assert shape(feature["geometry"]).is_valid
+            feature_count += 1
+    assert feature_count
+    # into different CRS
+    raw_config = geojson.dict
+    raw_config["input"].update(file1=landpoly_3857)
+    config = MapcheteConfig(raw_config)
+    vectorfile = config.params_at_zoom(zoom)["input"]["file1"]
+    pixelbuffer = 5
+    tile_pyramid = BufferedTilePyramid("geodetic", pixelbuffer=pixelbuffer)
+    tiles = tile_pyramid.tiles_from_geom(vectorfile.bbox(), zoom)
+    feature_count = 0
+    for tile in tiles:
+        for feature in read_vector_window(vectorfile.path, tile):
+            assert "properties" in feature
+            assert shape(feature["geometry"]).is_valid
+            feature_count += 1
+    assert feature_count
+
+
+def test_read_vector_window_errors(invalid_geojson):
+    with pytest.raises(FileNotFoundError):
+        read_vector_window(
+            "invalid_path", BufferedTilePyramid("geodetic").tile(0, 0, 0)
+        )
+    with pytest.raises(MapcheteIOError):
+        read_vector_window(
+            invalid_geojson, BufferedTilePyramid("geodetic").tile(0, 0, 0)
+        )
+
+
+def test_reproject_geometry(landpoly):
+    """Reproject geometry."""
+    with fiona.open(landpoly, "r") as src:
+        for feature in src:
+
+            # WGS84 to Spherical Mercator
+            out_geom = reproject_geometry(
+                shape(feature["geometry"]), CRS(src.crs), CRS().from_epsg(3857)
+            )
+            assert out_geom.is_valid
+
+            # WGS84 to LAEA
+            out_geom = reproject_geometry(
+                shape(feature["geometry"]), CRS(src.crs), CRS().from_epsg(3035)
+            )
+            assert out_geom.is_valid
+
+            # WGS84 to WGS84
+            out_geom = reproject_geometry(
+                shape(feature["geometry"]), CRS(src.crs), CRS().from_epsg(4326)
+            )
+            assert out_geom.is_valid
+
+
+def test_reproject_geometry_latlon2mercator():
+    # WGS84 bounds to Spherical Mercator
+    big_box = box(-180, -90, 180, 90)
+    reproject_geometry(big_box, CRS().from_epsg(4326), CRS().from_epsg(3857))
+
+    # WGS84 bounds to Spherical Mercator raising clip error
+    with pytest.raises(RuntimeError):
+        reproject_geometry(
+            big_box, CRS().from_epsg(4326), CRS().from_epsg(3857), error_on_clip=True
+        )
+    outside_box = box(-180, 87, 180, 90)
+    assert reproject_geometry(
+        outside_box,
+        CRS().from_epsg(4326),
+        CRS().from_epsg(3857),
+    ).is_valid
+
+
+def test_reproject_geometry_empty_geom():
+    # empty geometry
+    assert reproject_geometry(
+        Polygon(), CRS().from_epsg(4326), CRS().from_epsg(3857)
+    ).is_empty
+    assert reproject_geometry(
+        Polygon(), CRS().from_epsg(4326), CRS().from_epsg(4326)
+    ).is_empty
+
+
+def test_reproject_geometry_latlon2mercator_epsg():
+    # CRS parameter
+    big_box = box(-180, -90, 180, 90)
+    assert reproject_geometry(big_box, 4326, 3857) == reproject_geometry(
+        big_box, "4326", "3857"
+    )
+    with pytest.raises(TypeError):
+        reproject_geometry(big_box, 1.0, 1.0)
+
+
+def test_reproject_geometry_clip_crs_bounds_epsg():
+    bbox = wkt.loads(
+        "Polygon ((6.05416952699480682 49.79497943046440867, 6.04100166381764581 50.01055350300158864, 5.70657677854139056 50.00153486687963778, 5.72122668311700089 49.78602894452072292, 6.05416952699480682 49.79497943046440867))"
+    )
+    dst_crs = "EPSG:32632"
+
+    # reproject to UTM
+    bbox_utm = reproject_geometry(bbox, 4326, dst_crs, clip_to_crs_bounds=True)
+    assert bbox_utm.is_valid
+    assert bbox_utm.area
+    # revert to WGS84 and test geometry clipping
+    bbox_wgs84 = reproject_geometry(bbox_utm, dst_crs, 4326)
+    assert bbox_wgs84.area < bbox.area
+
+    # reproject to UTM but don't clip
+    bbox_utm = reproject_geometry(bbox, 4326, dst_crs, clip_to_crs_bounds=False)
+    assert bbox_utm.is_valid
+    assert bbox_utm.area
+    # make sure geometry was not clipped
+    bbox_wgs84 = reproject_geometry(bbox_utm, dst_crs, 4326)
+    assert bbox_wgs84.intersects(bbox)
+    assert (
+        bbox_wgs84.intersection(bbox).area
+        == pytest.approx(bbox_wgs84.area)
+        == pytest.approx(bbox.area)
+    )
+
+    # reproject to UTM don't clip but segmentize
+    bbox_utm = reproject_geometry(
+        bbox, 4326, dst_crs, clip_to_crs_bounds=False, segmentize=True
+    )
+    assert bbox_utm.is_valid
+    assert bbox_utm.area
+    # make sure geometry was not clipped
+    bbox_wgs84 = reproject_geometry(bbox_utm, dst_crs, 4326)
+    assert bbox_wgs84.intersects(bbox)
+    assert (
+        bbox_wgs84.intersection(bbox).area
+        == pytest.approx(bbox_wgs84.area)
+        == pytest.approx(bbox.area)
+    )
+
+
+def test_reproject_geometry_clip_crs_bounds_proj():
+    bbox = wkt.loads(
+        "Polygon ((6.05416952699480682 49.79497943046440867, 6.04100166381764581 50.01055350300158864, 5.70657677854139056 50.00153486687963778, 5.72122668311700089 49.78602894452072292, 6.05416952699480682 49.79497943046440867))"
+    )
+    dst_crs = "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs"
+
+    # reproject to UTM
+    bbox_utm = reproject_geometry(bbox, 4326, dst_crs, clip_to_crs_bounds=True)
+    assert bbox_utm.is_valid
+    assert bbox_utm.area
+    # revert to WGS84 and test geometry clipping
+    bbox_wgs84 = reproject_geometry(bbox_utm, dst_crs, 4326)
+    # NOTE: on some proj versions (TBD), pyproj cannot detect the CRS bounds of a CRS passed on by a proj string
+    # assert bbox_wgs84.area < bbox.area
+
+    # reproject to UTM but don't clip
+    bbox_utm = reproject_geometry(bbox, 4326, dst_crs, clip_to_crs_bounds=False)
+    assert bbox_utm.is_valid
+    assert bbox_utm.area
+    # make sure geometry was not clipped
+    bbox_wgs84 = reproject_geometry(bbox_utm, dst_crs, 4326)
+    assert bbox_wgs84.intersects(bbox)
+    assert (
+        bbox_wgs84.intersection(bbox).area
+        == pytest.approx(bbox_wgs84.area)
+        == pytest.approx(bbox.area)
+    )
+
+    # reproject to UTM don't clip but segmentize
+    bbox_utm = reproject_geometry(
+        bbox, 4326, dst_crs, clip_to_crs_bounds=False, segmentize=True
+    )
+    assert bbox_utm.is_valid
+    assert bbox_utm.area
+    # make sure geometry was not clipped
+    bbox_wgs84 = reproject_geometry(bbox_utm, dst_crs, 4326)
+    assert bbox_wgs84.intersects(bbox)
+    assert (
+        bbox_wgs84.intersection(bbox).area
+        == pytest.approx(bbox_wgs84.area)
+        == pytest.approx(bbox.area)
+    )
+
+
+def test_repair_geometry():
+    # invalid LineString
+    l = LineString([(0, 0), (0, 0), (0, 0)])
+    with pytest.raises(TopologicalError):
+        _repair(l)
+
+
+def test_write_vector_window_errors(landpoly):
+    with fiona.open(landpoly) as src:
+        feature = next(iter(src))
+    with pytest.raises((DriverError, ValueError, TypeError)):
+        write_vector_window(
+            in_data=["invalid", feature],
+            out_tile=BufferedTilePyramid("geodetic").tile(0, 0, 0),
+            out_path="/invalid_path",
+            out_schema=dict(geometry="Polygon", properties=dict()),
+        )
+
+
+def test_segmentize_geometry():
+    """Segmentize function."""
+    # Polygon
+    polygon = box(-18, -9, 18, 9)
+    out = segmentize_geometry(polygon, 1)
+    assert out.is_valid
+    # wrong type
+    with pytest.raises(TypeError):
+        segmentize_geometry(polygon.centroid, 1)
+
+
+def test_clean_geometry_type(geometrycollection):
+    """Filter and break up geometries."""
+    polygon = box(-18, -9, 18, 9)
+    # invalid type
+    with pytest.raises(TypeError):
+        clean_geometry_type(polygon, "invalid_type")
+
+    # don't return geometry
+    with pytest.raises(GeometryTypeError):
+        clean_geometry_type(polygon, "LineString", raise_exception=True)
+
+    # return geometry as is
+    assert clean_geometry_type(polygon, "Polygon").geom_type == "Polygon"
+    assert clean_geometry_type(polygon, "MultiPolygon").geom_type == "Polygon"
+
+    # don't allow multipart geometries
+    with pytest.raises(GeometryTypeError):
+        clean_geometry_type(
+            MultiPolygon([polygon]),
+            "Polygon",
+            allow_multipart=False,
+            raise_exception=True,
+        )
+
+    # multipolygons from geometrycollection
+    result = clean_geometry_type(
+        geometrycollection, "Polygon", allow_multipart=True, raise_exception=False
+    )
+    assert result.geom_type == "MultiPolygon"
+    assert not result.is_empty
+
+    # polygons from geometrycollection
+    result = clean_geometry_type(
+        geometrycollection, "Polygon", allow_multipart=False, raise_exception=False
+    )
+    assert result.geom_type == "GeometryCollection"
+    assert result.is_empty
+
+
+def test_convert_vector_copy(aoi_br_geojson, tmpdir):
+    out = os.path.join(tmpdir, "copied.geojson")
+
+    # copy
+    convert_vector(aoi_br_geojson, out)
+    with fiona.open(out) as src:
+        assert list(iter(src))
+
+    # raise error if output exists
+    with pytest.raises(IOError):
+        convert_vector(aoi_br_geojson, out, exists_ok=False)
+
+    # do nothing if output exists
+    convert_vector(aoi_br_geojson, out)
+    with fiona.open(out) as src:
+        assert list(iter(src))
+
+
+def test_convert_vector_overwrite(aoi_br_geojson, tmpdir):
+    out = os.path.join(tmpdir, "copied.geojson")
+
+    # write an invalid file
+    with open(out, "w") as dst:
+        dst.write("invalid")
+
+    # overwrite
+    convert_vector(aoi_br_geojson, out, overwrite=True)
+    with fiona.open(out) as src:
+        assert list(iter(src))
+
+
+def test_convert_vector_other_format_copy(aoi_br_geojson, tmpdir):
+    out = os.path.join(tmpdir, "copied.gpkg")
+
+    convert_vector(aoi_br_geojson, out, driver="GPKG")
+    with fiona.open(out) as src:
+        assert list(iter(src))
+
+    # raise error if output exists
+    with pytest.raises(IOError):
+        convert_vector(aoi_br_geojson, out, exists_ok=False)
+
+
+def test_convert_vector_other_format_overwrite(aoi_br_geojson, tmpdir):
+    out = os.path.join(tmpdir, "copied.gkpk")
+
+    # write an invalid file
+    with open(out, "w") as dst:
+        dst.write("invalid")
+
+    # overwrite
+    convert_vector(aoi_br_geojson, out, driver="GPKG", overwrite=True)
+    with fiona.open(out) as src:
+        assert list(iter(src))
+
+
+def test_indexed_features(landpoly):
+    with fiona.open(landpoly) as src:
+        some_id = next(iter(src))["id"]
+        features = IndexedFeatures(src)
+
+    assert features[some_id]
+    with pytest.raises(KeyError):
+        features["invalid_key"]
+
+    assert features.items()
+    assert features.keys()
+
+    for f in features:
+        assert isinstance(f, dict)
+        assert "properties" in f
+        assert "geometry" in f
+
+
+def test_indexed_features_bounds():
+    feature = {"bounds": (0, 1, 2, 3)}
+    assert IndexedFeatures([(0, feature)])
+
+    feature = {"geometry": mapping(box(0, 1, 2, 3))}
+    assert IndexedFeatures([(0, feature)])
+
+    feature = {"geometry": box(0, 1, 2, 3)}
+    assert IndexedFeatures([(0, feature)])
+
+    # no featuere id
+    with pytest.raises(TypeError):
+        IndexedFeatures([feature])
+
+    class Foo:
+        __geo_interface__ = mapping(box(0, 1, 2, 3))
+        id = 0
+
+    feature = Foo()
+    assert IndexedFeatures([(0, feature)])
+
+    # no bounds
+    feature = {"properties": {"foo": "bar"}, "id": 0}
+    with pytest.raises(TypeError):
+        IndexedFeatures([feature])
+
+
+def test_bounds_intersect():
+    b1 = (0, 0, 2, 2)
+    b2 = (1, 1, 3, 3)
+    b3 = (0, 1, 2, 3)
+    b4 = (1, 0, 3, 2)
+    b5 = (0, -1, 2, 1)
+    b6 = (0, -2, 2, 0)
+    b7 = (0, 0, 2, 2)
+    b8 = (0.5, 0.5, 1.5, 1.5)
+
+    assert bounds_intersect(b1, b1)
+    assert bounds_intersect(b1, b2)
+    assert bounds_intersect(b1, b3)
+    assert bounds_intersect(b1, b4)
+    assert bounds_intersect(b1, b5)
+    assert bounds_intersect(b1, b6)
+    assert bounds_intersect(b1, b7)
+    assert bounds_intersect(b1, b8)
+
+    intersecting = [
+        (-121.2451171875, 75.7177734375, -120.849609375, 76.025390625),
+        (-119.3994140625, 75.498046875, -117.4658203125, 76.11328125),
+        (-122.958984375, 75.8056640625, -115.400390625, 77.5634765625),
+        (-125.859375, 71.0595703125, -115.3125, 74.5751953125),
+        (-115.048828125, 77.6953125, -113.5986328125, 78.0908203125),
+        (-114.8291015625, 76.6845703125, -113.4228515625, 76.904296875),
+        (-113.3349609375, 77.34375, -109.5556640625, 78.134765625),
+        (-109.423828125, 67.8955078125, -109.2041015625, 68.0712890625),
+        (-113.3349609375, 78.2666015625, -109.2041015625, 78.75),
+        (-108.017578125, 73.5205078125, -107.578125, 73.6083984375),
+        (-117.685546875, 74.3994140625, -105.380859375, 76.8603515625),
+        (-105.1171875, 68.37890625, -104.4140625, 68.5986328125),
+        (-106.083984375, 77.080078125, -103.974609375, 77.783203125),
+        (-104.8974609375, 75.0146484375, -103.5791015625, 75.4541015625),
+        (-104.677734375, 76.3330078125, -103.0078125, 76.6845703125),
+        (-102.2607421875, 68.5546875, -101.689453125, 68.818359375),
+        (-102.5244140625, 77.6953125, -100.9423828125, 77.9150390625),
+        (-119.1357421875, 68.466796875, -100.8544921875, 73.7841796875),
+        (-101.6015625, 76.552734375, -100.4150390625, 76.7724609375),
+        (-100.634765625, 70.4443359375, -100.1953125, 70.6640625),
+        (-100.5908203125, 68.6865234375, -99.9755859375, 69.169921875),
+        (-105.64453125, 77.783203125, -98.9208984375, 79.4091796875),
+        (-100.1513671875, 79.6728515625, -98.6572265625, 80.15625),
+        (-99.404296875, 73.828125, -97.6025390625, 74.1357421875),
+        (-104.4580078125, 74.970703125, -97.3388671875, 76.6845703125),
+        (-97.7783203125, 74.443359375, -97.2509765625, 74.619140625),
+        (-97.119140625, 72.9052734375, -96.5478515625, 73.212890625),
+        (-96.9873046875, 75.3662109375, -96.416015625, 75.5859375),
+        (-102.744140625, 71.279296875, -96.2841796875, 73.9599609375),
+        (-96.591796875, 69.345703125, -96.064453125, 69.5654296875),
+        (-96.3720703125, 75.4541015625, -95.9326171875, 75.673828125),
+        (-95.9765625, 69.3017578125, -95.361328125, 69.6533203125),
+        (-99.5361328125, 68.466796875, -95.2294921875, 69.9169921875),
+        (-95.7568359375, 74.4873046875, -95.2294921875, 74.6630859375),
+        (-98.3935546875, 77.783203125, -94.833984375, 78.837890625),
+        (-94.658203125, 78.1787109375, -94.306640625, 78.2666015625),
+        (-94.8779296875, 75.76171875, -94.2626953125, 75.9814453125),
+        (-96.6357421875, 74.619140625, -93.3837890625, 75.6298828125),
+        (-96.2841796875, 77.431640625, -93.1201171875, 77.8271484375),
+        (-91.23046875, 77.1240234375, -90.703125, 77.255859375),
+        (-90.0, 68.5986328125, -89.9560546875, 68.642578125),
+        (-90.0439453125, 68.6865234375, -89.9560546875, 68.818359375),
+        (-90.087890625, 71.8505859375, -89.9560546875, 72.0703125),
+        (-96.8994140625, 74.53125, -89.9560546875, 77.255859375),
+        (-90.5712890625, 76.46484375, -89.9560546875, 76.8603515625),
+        (-91.1865234375, 77.2119140625, -89.9560546875, 77.6513671875),
+        (-96.7236328125, 78.134765625, -89.9560546875, 81.38671875),
+        (-168.134765625, 13.6669921875, -89.9560546875, 74.1796875),
+    ]
+    for b in intersecting:
+        assert bounds_intersect(b, (-135.0, 60, -90, 80))
+
+
+def test_bounds_not_intersect():
+    b1 = (-3, -3, -2, -2)
+    b2 = (1, 1, 3, 3)
+    b3 = (0, 1, 2, 3)
+    b4 = (1, 0, 3, 2)
+    b5 = (0, -1, 2, 1)
+    b6 = (0, -2, 2, 0)
+    b7 = (0, 0, 2, 2)
+    b8 = (0.5, 0.5, 1.5, 1.5)
+
+    assert not bounds_intersect(b1, b2)
+    assert not bounds_intersect(b1, b3)
+    assert not bounds_intersect(b1, b4)
+    assert not bounds_intersect(b1, b5)
+    assert not bounds_intersect(b1, b6)
+    assert not bounds_intersect(b1, b8)
+    assert not bounds_intersect(b1, b7)
+
+
+def test_indexed_features_fakeindex(landpoly):
+    with fiona.open(landpoly) as src:
+        features = list(src)
+        idx = IndexedFeatures(features)
+        fake_idx = IndexedFeatures(features, index=None)
+    bounds = (-135.0, 60, -90, 80)
+    assert len(idx.filter(bounds)) == len(fake_idx.filter(bounds))
+
+
+def test_indexed_features_polygon(aoi_br_geojson):
+    with fiona.open(aoi_br_geojson) as src:
+        index = IndexedFeatures(src)
+    tp = BufferedTilePyramid("geodetic")
+    for tile in tp.tiles_from_bounds(bounds=index.bounds, zoom=5):
+        assert len(index.filter(tile.bounds)) == 1
diff --git a/test/test_mapchete.py b/test/test_mapchete.py
index ccf3aa99..45d21805 100644
--- a/test/test_mapchete.py
+++ b/test/test_mapchete.py
@@ -22,6 +22,11 @@
 from mapchete.io.raster import create_mosaic, _shift_required
 from mapchete.errors import MapcheteProcessOutputError
 from mapchete.tile import BufferedTilePyramid, count_tiles
+from mapchete._executor import DaskExecutor
+from mapchete._processing import (
+    PreprocessingProcessInfo,
+    TileProcessInfo,
+)
 
 
 def test_empty_execute(mp_tmpdir, cleantopo_br):
@@ -186,6 +191,34 @@ def test_baselevels(mp_tmpdir, baselevels):
         )
 
 
+def test_baselevels_dask(mp_tmpdir, baselevels):
+    """Baselevel interpolation."""
+    with mapchete.open(baselevels.dict, mode="continue") as mp:
+        # process data before getting baselevels
+        with DaskExecutor() as executor:
+            mp.batch_process(executor=executor)
+
+        # get tile from lower zoom level
+        for tile in mp.get_process_tiles(4):
+            data = mp.get_raw_output(tile)
+            assert not data.mask.all()
+            # write for next zoom level
+            mp.write(tile, data)
+            assert not mp.get_raw_output(tile.get_parent()).mask.all()
+
+        # get tile from higher zoom level
+        tile = next(mp.get_process_tiles(6))
+        # process and save
+        mp.write(tile, mp.get_raw_output(tile))
+        # read from baselevel
+        assert any(
+            [
+                not mp.get_raw_output(upper_tile).mask.all()
+                for upper_tile in tile.get_children()
+            ]
+        )
+
+
 def test_baselevels_custom_nodata(mp_tmpdir, baselevels_custom_nodata):
     """Baselevel interpolation."""
     fill_value = -32768.0
@@ -377,7 +410,7 @@ def test_processing(mp_tmpdir, cleantopo_br, cleantopo_tl):
                             testfile.meta["transform"], mosaic.affine
                         ):
                             assert file_item == mosaic_item
-                        band = testfile.read(1, masked=True)
+                        band = testfile.read(masked=True)
                         assert band.shape == mosaic.data.shape
                         assert ma.allclose(band, mosaic.data)
                         assert ma.allclose(band.mask, mosaic.data.mask)
@@ -594,3 +627,140 @@ def test_bufferedtiles():
     assert a != tp_buffered.tile(5, 5, 5)
 
     assert a.get_neighbors() != a.get_neighbors(connectedness=4)
+
+
+def test_compute_dask_graph(preprocess_cache_memory):
+    with preprocess_cache_memory.mp(batch_preprocess=False) as mp:
+        preprocessing_tasks = 0
+        tile_tasks = 0
+        for future in mp.compute(concurrency="dask", dask_compute_graph=True):
+            result = future.result()
+            if isinstance(result, PreprocessingProcessInfo):
+                assert result.data is not None
+                preprocessing_tasks += 1
+            else:
+                assert isinstance(result, TileProcessInfo)
+                tile_tasks += 1
+                if result.written:
+                    assert result.data is not None
+                else:
+                    assert result.data is None
+    assert tile_tasks == 20
+    assert preprocessing_tasks == 2
+
+
+def test_compute_dask(preprocess_cache_memory):
+    with preprocess_cache_memory.mp(batch_preprocess=False) as mp:
+        preprocessing_tasks = 0
+        tile_tasks = 0
+        for future in mp.compute(concurrency="dask", dask_compute_graph=False):
+            result = future.result()
+            if isinstance(result, PreprocessingProcessInfo):
+                preprocessing_tasks += 1
+            else:
+                assert isinstance(result, TileProcessInfo)
+                tile_tasks += 1
+                assert result.data is None
+    assert tile_tasks == 20
+    assert preprocessing_tasks == 2
+
+
+def test_compute_threads(preprocess_cache_memory):
+    with preprocess_cache_memory.mp(batch_preprocess=False) as mp:
+        preprocessing_tasks = 0
+        tile_tasks = 0
+        for future in mp.compute(concurrency="threads", dask_compute_graph=False):
+            result = future.result()
+            if isinstance(result, PreprocessingProcessInfo):
+                preprocessing_tasks += 1
+            else:
+                assert isinstance(result, TileProcessInfo)
+                tile_tasks += 1
+                assert result.data is None
+    assert tile_tasks == 20
+    assert preprocessing_tasks == 2
+
+
+def test_compute_processes(preprocess_cache_memory):
+    with preprocess_cache_memory.mp(batch_preprocess=False) as mp:
+        preprocessing_tasks = 0
+        tile_tasks = 0
+        for future in mp.compute(concurrency="processes", dask_compute_graph=False):
+            result = future.result()
+            if isinstance(result, PreprocessingProcessInfo):
+                preprocessing_tasks += 1
+            else:
+                assert isinstance(result, TileProcessInfo)
+                tile_tasks += 1
+                assert result.data is None
+    assert tile_tasks == 20
+    assert preprocessing_tasks == 2
+
+
+def test_compute_dask_graph_single_file(preprocess_cache_memory_single_file):
+    with preprocess_cache_memory_single_file.mp(batch_preprocess=False) as mp:
+        preprocessing_tasks = 0
+        tile_tasks = 0
+        for future in mp.compute(concurrency="dask", dask_compute_graph=True):
+            result = future.result()
+            if isinstance(result, PreprocessingProcessInfo):
+                assert result.data is not None
+                preprocessing_tasks += 1
+            else:
+                assert isinstance(result, TileProcessInfo)
+                tile_tasks += 1
+                assert result.data is None
+    assert tile_tasks == 12
+    assert preprocessing_tasks == 2
+    with rasterio.open(mp.config.output.path) as src:
+        assert not src.read(masked=True).mask.all()
+
+
+def test_compute_dask_single_file(preprocess_cache_memory_single_file):
+    with preprocess_cache_memory_single_file.mp(batch_preprocess=False) as mp:
+        preprocessing_tasks = 0
+        tile_tasks = 0
+        for future in mp.compute(concurrency="dask", dask_compute_graph=False):
+            result = future.result()
+            if isinstance(result, PreprocessingProcessInfo):
+                preprocessing_tasks += 1
+            else:
+                assert isinstance(result, TileProcessInfo)
+                tile_tasks += 1
+                assert result.data is None
+    assert tile_tasks == 12
+    assert preprocessing_tasks == 2
+    with rasterio.open(mp.config.output.path) as src:
+        assert not src.read(masked=True).mask.all()
+
+
+def test_compute_threads_single_file(preprocess_cache_memory_single_file):
+    with preprocess_cache_memory_single_file.mp(batch_preprocess=False) as mp:
+        preprocessing_tasks = 0
+        tile_tasks = 0
+        for future in mp.compute(concurrency="threads", dask_compute_graph=False):
+            result = future.result()
+            if isinstance(result, PreprocessingProcessInfo):
+                preprocessing_tasks += 1
+            else:
+                assert isinstance(result, TileProcessInfo)
+                tile_tasks += 1
+                assert result.data is None
+    assert tile_tasks == 12
+    assert preprocessing_tasks == 2
+
+
+def test_compute_processes_single_file(preprocess_cache_memory_single_file):
+    with preprocess_cache_memory_single_file.mp(batch_preprocess=False) as mp:
+        preprocessing_tasks = 0
+        tile_tasks = 0
+        for future in mp.compute(concurrency="processes", dask_compute_graph=False):
+            result = future.result()
+            if isinstance(result, PreprocessingProcessInfo):
+                preprocessing_tasks += 1
+            else:
+                assert isinstance(result, TileProcessInfo)
+                tile_tasks += 1
+                assert result.data is None
+    assert tile_tasks == 12
+    assert preprocessing_tasks == 2
diff --git a/test/test_preprocessing.py b/test/test_preprocessing.py
index 3f3fd184..622b7b1e 100644
--- a/test/test_preprocessing.py
+++ b/test/test_preprocessing.py
@@ -1,5 +1,8 @@
 import mapchete
+from mapchete.io import fs_from_path
+from mapchete._tasks import Task
 import pytest
+import rasterio
 
 
 def _trivial_func(arg, kwarg=None):
@@ -80,10 +83,138 @@ def test_run_preprocessing_tasks(example_mapchete):
             == 3
         )
 
-        # number of unique tasks should be 2
-        assert mp.config.preprocessing_tasks_count() == 2
+        assert mp.config.preprocessing_tasks_count() == 3
 
         mp.batch_preprocess()
+        assert mp.config.preprocessing_task_finished(f"{inp1.input_key}:test_task")
         assert inp1.get_preprocessing_task_result("test_task") == "foofoobar"
+        assert mp.config.preprocessing_task_finished(
+            f"{inp1.input_key}:test_other_task"
+        )
         assert inp1.get_preprocessing_task_result("test_other_task") == "barfoobar"
+        assert mp.config.preprocessing_task_finished(f"{inp2.input_key}:test_task")
         assert inp2.get_preprocessing_task_result("test_task") == "foofoobar"
+
+
+def test_preprocess_cache_raster_vector(preprocess_cache_raster_vector):
+    # tile containing only features: 5 29 62
+    mp = preprocess_cache_raster_vector.process_mp((5, 29, 62))
+    with mp.open("clip") as clip:
+        assert clip.path.endswith("cache/aoi_br.fgb")
+        assert clip.read()
+    with mp.open("inp") as raster:
+        assert raster.path.endswith("cache/cleantopo_br.tif")
+        assert raster.read().mask.all()
+
+    # tile containing features and raster: 5 30 62
+    mp = preprocess_cache_raster_vector.process_mp((5, 30, 62))
+    with mp.open("clip") as clip:
+        assert clip.path.endswith("cache/aoi_br.fgb")
+        assert clip.read()
+    with mp.open("inp") as raster:
+        assert raster.path.endswith("cache/cleantopo_br.tif")
+        assert not raster.read().mask.all()
+
+
+def test_preprocess_cache_raster_vector_tasks(preprocess_cache_raster_vector):
+    with preprocess_cache_raster_vector.mp() as mp:
+        for i in ["clip", "inp"]:
+            input_data = mp.config.input_at_zoom(key=i, zoom=5)
+            for task in input_data.preprocessing_tasks.values():
+                assert isinstance(task, Task)
+                assert task.has_geometry()
+
+
+def test_preprocessing_tasks_dependencies(preprocess_cache_memory):
+    with preprocess_cache_memory.mp() as mp:
+        for i in ["clip", "inp"]:
+            input_data = mp.config.input_at_zoom(key=i, zoom=5)
+            for task in input_data.preprocessing_tasks.values():
+                assert isinstance(task, Task)
+                assert task.has_geometry()
+        list(mp.compute())
+
+        out_path = mp.config.output_reader.path
+        total_tifs = len(
+            [
+                file
+                for directory in fs_from_path(out_path).walk(out_path)
+                for file in directory[2]
+                if file.endswith(".tif")
+            ]
+        )
+        assert total_tifs == 9
+
+
+def test_preprocessing_tasks_dependencies_dask(preprocess_cache_memory):
+    with preprocess_cache_memory.mp(batch_preprocess=False) as mp:
+        for i in ["clip", "inp"]:
+            input_data = mp.config.input_at_zoom(key=i, zoom=5)
+            for task in input_data.preprocessing_tasks.values():
+                assert isinstance(task, Task)
+                assert task.has_geometry()
+
+        list(mp.compute(concurrency="dask"))
+
+        out_path = mp.config.output_reader.path
+        total_tifs = [
+            f"{directory[0]}/{file}"
+            for directory in fs_from_path(out_path).walk(out_path)
+            for file in directory[2]
+            if file.endswith(".tif")
+        ]
+        assert len(total_tifs) == 9
+
+
+def test_preprocessing_tasks_dependencies_single_tile(preprocess_cache_memory):
+    with preprocess_cache_memory.mp(batch_preprocess=False) as mp:
+        for i in ["clip", "inp"]:
+            input_data = mp.config.input_at_zoom(key=i, zoom=5)
+            for task in input_data.preprocessing_tasks.values():
+                assert isinstance(task, Task)
+                assert task.has_geometry()
+
+        tile = (5, 31, 63)
+        list(mp.compute(tile=tile))
+
+        out_path = mp.config.output_reader.get_path(
+            mp.config.process_pyramid.tile(*tile)
+        )
+        with rasterio.open(out_path) as src:
+            assert not src.read(masked=True).mask.all()
+
+        out_path = mp.config.output_reader.path
+        total_tifs = [
+            f"{directory[0]}/{file}"
+            for directory in fs_from_path(out_path).walk(out_path)
+            for file in directory[2]
+            if file.endswith(".tif")
+        ]
+        assert len(total_tifs) == 1
+
+
+def test_preprocessing_tasks_dependencies_single_tile_dask(preprocess_cache_memory):
+    with preprocess_cache_memory.mp(batch_preprocess=False) as mp:
+        for i in ["clip", "inp"]:
+            input_data = mp.config.input_at_zoom(key=i, zoom=5)
+            for task in input_data.preprocessing_tasks.values():
+                assert isinstance(task, Task)
+                assert task.has_geometry()
+
+        tile = (5, 31, 63)
+        list(mp.compute(concurrency="dask", tile=tile))
+
+        out_path = mp.config.output_reader.get_path(
+            mp.config.process_pyramid.tile(*tile)
+        )
+        with rasterio.open(out_path) as src:
+            assert not src.read(masked=True).mask.all()
+
+        out_path = mp.config.output_reader.path
+        total_tifs = [
+            f"{directory[0]}/{file}"
+            for directory in fs_from_path(out_path).walk(out_path)
+            for file in directory[2]
+            if file.endswith(".tif")
+        ]
+        assert len(total_tifs) == 1
diff --git a/test/test_tasks.py b/test/test_tasks.py
new file mode 100644
index 00000000..186fce9a
--- /dev/null
+++ b/test/test_tasks.py
@@ -0,0 +1,149 @@
+from itertools import chain
+from typing import Type
+import pytest
+from shapely.geometry import shape
+
+
+from mapchete.errors import NoTaskGeometry
+from mapchete._tasks import Task, TileTask, TaskBatch, TileTaskBatch, to_dask_collection
+
+
+def test_task_geo_interface():
+    task = Task(bounds=(0, 1, 2, 3))
+    assert not shape(task).is_empty
+
+
+def test_task_errors():
+    # provide geometry and bounds at the same time
+    with pytest.raises(ValueError):
+        Task(geometry="foo", bounds="bar")
+
+    # task has no geo information
+    with pytest.raises(NoTaskGeometry):
+        Task().__geo_interface__
+
+    # invalid dependencies
+    with pytest.raises(TypeError):
+        Task().add_dependencies("invalid")
+
+
+def test_task_dict():
+    task_dict = Task(bounds=(0, 1, 2, 3)).to_dict()
+    assert "geometry" in task_dict
+    assert "properties" in task_dict
+    assert "id" in task_dict
+    assert "bounds" in task_dict
+
+
+def test_task_dependencies():
+    task = Task()
+    task.add_dependencies({"foo": "bar"})
+    assert "foo" in task.dependencies
+
+
+def test_tile_task_geo_interface(example_mapchete):
+    task = TileTask(
+        tile=example_mapchete.first_process_tile(), config=example_mapchete.mp().config
+    )
+    assert not shape(task).is_empty
+
+
+def test_task_batches():
+    batch = TaskBatch(
+        (Task(func=str, fargs=(i,), bounds=(0, 1, 2, 3)) for i in range(10))
+    )
+    assert batch.items()
+    assert batch.keys()
+    assert batch.values()
+
+    for task in batch:
+        assert isinstance(task, Task)
+
+    other_task = Task(bounds=(0, 1, 2, 3))
+    assert len(batch.intersection(other_task)) == 10
+    assert len(batch.intersection((0, 1, 2, 3))) == 10
+
+
+def test_task_batches_errors():
+    with pytest.raises(TypeError):
+        TaskBatch()
+
+    with pytest.raises(TypeError):
+        TaskBatch(["invalid"])
+
+    batch = TaskBatch(
+        (Task(func=str, fargs=(i,), bounds=(0, 1, 2, 3)) for i in range(10))
+    )
+    with pytest.raises(TypeError):
+        batch.intersection("invalid")
+
+
+def test_tile_task_batch(dem_to_hillshade):
+    tile_task_batch = TileTaskBatch(
+        (
+            TileTask(tile=process_tile, config=dem_to_hillshade.mp().config)
+            for process_tile in dem_to_hillshade.mp().get_process_tiles(zoom=5)
+        )
+    )
+
+    for task in tile_task_batch:
+        assert isinstance(task, TileTask)
+
+    other_tile = dem_to_hillshade.first_process_tile().get_parent()
+    tile_task = TileTask(tile=other_tile, config=dem_to_hillshade.mp().config)
+    assert len(tile_task_batch.intersection(tile_task)) == 4
+
+    task = Task(bounds=other_tile.bounds)
+    assert len(tile_task_batch.intersection(task)) == 4
+
+    assert len(tile_task_batch.intersection(task.bounds)) == 4
+
+    with pytest.raises(TypeError):
+        tile_task_batch.intersection("invalid")
+
+
+def test_task_batches_to_dask_graph(dem_to_hillshade):
+    preprocessing_batch = TaskBatch(
+        (Task(func=str, fargs=(i,), bounds=(0, 1, 2, 3)) for i in range(10))
+    )
+    assert preprocessing_batch.items()
+    assert preprocessing_batch.keys()
+    assert preprocessing_batch.values()
+    zoom_batches = (
+        TileTaskBatch(
+            (
+                TileTask(tile=process_tile, config=dem_to_hillshade.mp().config)
+                for process_tile in dem_to_hillshade.mp().get_process_tiles(zoom=zoom)
+            )
+        )
+        for zoom in reversed(dem_to_hillshade.mp().config.zoom_levels)
+    )
+    collection = to_dask_collection((preprocessing_batch, *zoom_batches))
+    import dask
+
+    dask.compute(collection)
+
+
+def test_task_batches_mixed_geometries():
+    batch = TaskBatch(
+        chain(
+            (Task(func=str, fargs=(i,), bounds=(0, 1, 2, 3)) for i in range(10)),
+            (Task(func=str, fargs=(i,)) for i in range(10)),
+        )
+    )
+    assert len(batch.items()) == 20
+    assert len(batch.keys()) == 20
+    assert len(batch.values()) == 20
+
+    for task in batch:
+        assert isinstance(task, Task)
+
+    # other task intersecting with all tasks
+    other_task = Task(bounds=(0, 1, 2, 3))
+    assert len(batch.intersection(other_task)) == 20
+    assert len(batch.intersection((0, 1, 2, 3))) == 20
+
+    # other task not intersecting with spatial tasks
+    other_task = Task(bounds=(3, 4, 5, 6))
+    assert len(batch.intersection(other_task)) == 10
+    assert len(batch.intersection((3, 4, 5, 6))) == 10
diff --git a/test/testdata/preprocess_cache_memory.mapchete b/test/testdata/preprocess_cache_memory.mapchete
new file mode 100644
index 00000000..f25b2bc7
--- /dev/null
+++ b/test/testdata/preprocess_cache_memory.mapchete
@@ -0,0 +1,22 @@
+process: mapchete.processes.convert
+zoom_levels: 
+    min: 0
+    max: 5
+baselevels:
+    min: 5
+pyramid:
+    grid: geodetic
+input:
+    inp: 
+        format: raster_file
+        path: cleantopo_br.tif
+        cache: memory
+    clip:
+        format: vector_file
+        path: aoi_br.geojson
+        cache: memory
+output:
+    dtype: uint16
+    bands: 1
+    format: GTiff
+    path: tmp/preprocess_cache_raster_vector
diff --git a/test/testdata/preprocess_cache_memory_single_file.mapchete b/test/testdata/preprocess_cache_memory_single_file.mapchete
new file mode 100644
index 00000000..d886a76b
--- /dev/null
+++ b/test/testdata/preprocess_cache_memory_single_file.mapchete
@@ -0,0 +1,18 @@
+process: mapchete.processes.convert
+zoom_levels: 5
+pyramid:
+    grid: geodetic
+input:
+    inp: 
+        format: raster_file
+        path: cleantopo_br.tif
+        cache: memory
+    clip:
+        format: vector_file
+        path: aoi_br.geojson
+        cache: memory
+output:
+    dtype: uint16
+    bands: 1
+    format: GTiff
+    path: tmp/preprocess_cache_raster_vector_single_file.tif
diff --git a/test/testdata/preprocess_cache_raster_vector.mapchete b/test/testdata/preprocess_cache_raster_vector.mapchete
new file mode 100644
index 00000000..714a88ea
--- /dev/null
+++ b/test/testdata/preprocess_cache_raster_vector.mapchete
@@ -0,0 +1,24 @@
+process: mapchete.processes.convert
+zoom_levels: 5
+pyramid:
+    grid: geodetic
+input:
+    inp: 
+        format: raster_file
+        path: cleantopo_br.tif
+        cache:
+            path: tmp/preprocess_cache_raster_vector/cache/cleantopo_br.tif
+            format: COG
+            keep: false
+    clip:
+        format: vector_file
+        path: aoi_br.geojson
+        cache:
+            path: tmp/preprocess_cache_raster_vector/cache/aoi_br.fgb
+            format: FlatGeobuf
+            keep: false
+output:
+    dtype: uint16
+    bands: 1
+    format: GTiff
+    path: tmp/preprocess_cache_raster_vector