Repository: ungarj/mapchete (language: python)

Issue:
provide STAC metadata for all generated products
see https://github.com/radiantearth/stac-spec

details like how to represent the generated pyramids in best way to be discussed
add concurrent processes warnings on Python<3.7

use "spawn" as default multiprocessing method
it seems `fsspec` has issues when used with multiprocessing fork: https://github.com/dask/gcsfs/issues/379


Discussion / hints:
Article to introduce STAC
* https://medium.com/radiant-earth-insights/spatiotemporal-asset-catalog-stac-community-advances-and-0-6-1-release-bddcd967d86f

STAC Validators:
* https://staclint.com/
* https://github.com/sparkgeo/stac-validator
keep #153 in mind by adding custom metadata items for STAC in mapchete file
Main issue is how to handle tile directories, i.e. potentially millions of files with predictable paths. Creating one asset per file does not scale and is probably also not intended by the STAC spec but it can be implemented quickly.

Short-term solution: one asset per file
Long-term solution: one asset per tile directory

@8snit @Scartography any comments?
Use structure required by GDAL 3.3: https://gdal.org/drivers/raster/stacta.html

To make this work with mapchete we have to
* [ ] provide an option to automatically add STAC metadata to new outputs, and
* [ ] provide a tool which creates the STAC metadata file for existing outputs

For the first point we can add the required metadata in the `output` section of a mapchete file:
```
output:
  path: testdata/tmp/example
  format: GTiff
  dtype: float32
  bands: 3
  stac:
    id: test_output  # optional
    properties:  # optional
      whatever: property
    eo:bands:  # optional for STAC extensions
      - name: B04
        common_name: red
      - name: B03
        common_name: green
      - name: B02
        common_name: blue
    links:  # optional
      - rel: parent
        href: /path/to/parent/catalog.json
      - rel: root
        href: /path/to/root/catalog.json
```
This would then add an <id>.json file to the tile directory or next to the single file.

For the second point we need a CLI tool and additional options for existing tools to add and/or update STAC metadata:
```
# to add or update STAC json
$ mapchete stac <mapchete_file or tile_directory> --stac-item-id <some-id>

# to add STAC json while using convert; per default an existing stac json should be used
$ mapchete convert <input> <output> --stac-item-id <some-id>
```

The STAC item id should default to the tile directory path when using a tile directory and the file basename (without extension) when using a single file output.



Fix patch:
diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
index d8c19316..45cf02f4 100644
--- a/.github/workflows/python-package.yml
+++ b/.github/workflows/python-package.yml
@@ -15,7 +15,6 @@ jobs:
       fail-fast: false
       matrix:
         python-version: [ 3.7, 3.8 ]
-        rasterio-version: [ 1.1.4, "latest" ]
         experimental: [ false ]
     runs-on: ubuntu-latest
     continue-on-error: ${{ matrix.experimental }}
@@ -30,15 +29,12 @@ jobs:
 
     - name: Install dependencies
       env:
-        RASTERIO_VERSION: ${{ matrix.rasterio-version }}
         CURL_CA_BUNDLE: /etc/ssl/certs/ca-certificates.crt
       run: |
-        sudo apt-add-repository -y ppa:ubuntugis/ubuntugis-unstable && sudo apt-get -y update && sudo apt-get install -y gdal-bin python-tk libgdal-dev
+        sudo apt-add-repository -y ppa:ubuntugis/ubuntugis-unstable && sudo apt-get -y update && sudo apt-get install -y gdal-bin python-tk libgdal-dev libproj-dev
         python -m pip install --upgrade pip
-        pip install fiona shapely --no-binary :all:
-        if [ $RASTERIO_VERSION == "latest" ]; then pip install rasterio --no-binary :all:; else pip install rasterio==$RASTERIO_VERSION --no-binary :all:; fi;
-        pip install -r requirements.txt -r test/requirements.txt
         pip install -e .[complete]
+        pip install -r test/requirements.txt
 
     - name: Lint with flake8
       run: |
diff --git a/.gitignore b/.gitignore
index 8aee8931..06c9e218 100644
--- a/.gitignore
+++ b/.gitignore
@@ -10,4 +10,5 @@ venv/
 .eggs/
 htmlcov/
 .pytest*
-.env
\ No newline at end of file
+.env
+dask-worker-space/
\ No newline at end of file
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index c585986d..e4e1cd70 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -2,6 +2,64 @@
 Changelog
 #########
 
+
+-----------------
+0.43 - 2021-09-17
+-----------------
+
+* core
+
+  * `mapchete.io.vector.reproject_geometry()`:
+
+    * use `pyproj` to determine CRS bounds to clip geometries when reprojecting
+    * enable geometry segmentation before geometry is clipped (`segmentize_on_clip=False` and `segmentize_fraction=100` args)
+  * suppress `rasterio` warnings when reading rasters (too many `rasterio.errors.NodataShadowWarning`s)
+
+* packaging
+
+  * add `pyproj` to dependencies
+
+
+-----------------
+0.42 - 2021-08-27
+-----------------
+
+* core
+
+  * add option for input drivers to let processing-heavy tasks be handled by ``mapchete.Executor`` by implementing ``InputData.add_preprocessing_task()`` and ``InputData.get_preprocessing_task_result()`` methods
+  * check futures submitted to ``mapchete.Executor`` for exceptions before yielding
+  * remove deprecated ``multi`` (now called ``workers``), ``distributed`` and ``max_chunksize`` arguments
+
+* CLI
+
+  * remove deprecated ``--max-chunksize`` option
+  * replace "tiles" with "tasks" in progress
+
+
+-----------------
+0.41 - 2021-08-17
+-----------------
+
+* core
+
+  * add ``mapchete.commands`` package
+  * add ``dask`` as optional ``Executor``
+  * expose futures in ``Executor`` class to facilitate job cancellation
+  * use ``spawn`` as default multiprocessing start method (#351)
+  * fix missing worker logs (#185)
+  * rename ``mapchete.cli.utils`` to ``mapchete.cli.options``
+  * enable providing process code from within process config
+
+* packaging
+
+  * updated API doc
+
+* CLI
+
+  * CLI: change ``--multi`` option to ``--worker``
+  * enable optional concurrency for ``mapchete cp``
+
+
 -----------------
 0.40 - 2021-06-24
 -----------------
diff --git a/LICENSE b/LICENSE
index ec16c22a..106e0c8b 100644
--- a/LICENSE
+++ b/LICENSE
@@ -1,6 +1,6 @@
 The MIT License (MIT)
 
-Copyright (c) 2015 - 2020 EOX IT Services
+Copyright (c) 2015 - 2021 EOX IT Services
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal
diff --git a/README.rst b/README.rst
index decbed8c..e162fea0 100644
--- a/README.rst
+++ b/README.rst
@@ -193,6 +193,6 @@ License
 
 MIT License
 
-Copyright (c) 2015 - 2020 `EOX IT Services`_
+Copyright (c) 2015 - 2021 `EOX IT Services`_
 
 .. _`EOX IT Services`: https://eox.at/
diff --git a/doc/source/apidoc/mapchete.cli.default.convert.rst b/doc/source/apidoc/mapchete.cli.default.convert.rst
index 84c68f1d..99eb8eca 100644
--- a/doc/source/apidoc/mapchete.cli.default.convert.rst
+++ b/doc/source/apidoc/mapchete.cli.default.convert.rst
@@ -2,6 +2,6 @@ mapchete.cli.default.convert module
 ===================================
 
 .. automodule:: mapchete.cli.default.convert
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.cli.default.cp.rst b/doc/source/apidoc/mapchete.cli.default.cp.rst
new file mode 100644
index 00000000..29fd3e37
--- /dev/null
+++ b/doc/source/apidoc/mapchete.cli.default.cp.rst
@@ -0,0 +1,7 @@
+mapchete.cli.default.cp module
+==============================
+
+.. automodule:: mapchete.cli.default.cp
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.cli.default.create.rst b/doc/source/apidoc/mapchete.cli.default.create.rst
index 7a9ff23d..856efe13 100644
--- a/doc/source/apidoc/mapchete.cli.default.create.rst
+++ b/doc/source/apidoc/mapchete.cli.default.create.rst
@@ -2,6 +2,6 @@ mapchete.cli.default.create module
 ==================================
 
 .. automodule:: mapchete.cli.default.create
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.cli.default.execute.rst b/doc/source/apidoc/mapchete.cli.default.execute.rst
index 2777b362..434167c7 100644
--- a/doc/source/apidoc/mapchete.cli.default.execute.rst
+++ b/doc/source/apidoc/mapchete.cli.default.execute.rst
@@ -2,6 +2,6 @@ mapchete.cli.default.execute module
 ===================================
 
 .. automodule:: mapchete.cli.default.execute
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.cli.default.formats.rst b/doc/source/apidoc/mapchete.cli.default.formats.rst
index 6c23f340..df63b832 100644
--- a/doc/source/apidoc/mapchete.cli.default.formats.rst
+++ b/doc/source/apidoc/mapchete.cli.default.formats.rst
@@ -2,6 +2,6 @@ mapchete.cli.default.formats module
 ===================================
 
 .. automodule:: mapchete.cli.default.formats
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.cli.default.index.rst b/doc/source/apidoc/mapchete.cli.default.index.rst
index cdd1d74d..e20cf940 100644
--- a/doc/source/apidoc/mapchete.cli.default.index.rst
+++ b/doc/source/apidoc/mapchete.cli.default.index.rst
@@ -2,6 +2,6 @@ mapchete.cli.default.index module
 =================================
 
 .. automodule:: mapchete.cli.default.index
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.cli.default.processes.rst b/doc/source/apidoc/mapchete.cli.default.processes.rst
index fd9809ff..681847d0 100644
--- a/doc/source/apidoc/mapchete.cli.default.processes.rst
+++ b/doc/source/apidoc/mapchete.cli.default.processes.rst
@@ -2,6 +2,6 @@ mapchete.cli.default.processes module
 =====================================
 
 .. automodule:: mapchete.cli.default.processes
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.cli.default.rm.rst b/doc/source/apidoc/mapchete.cli.default.rm.rst
new file mode 100644
index 00000000..73820fdc
--- /dev/null
+++ b/doc/source/apidoc/mapchete.cli.default.rm.rst
@@ -0,0 +1,7 @@
+mapchete.cli.default.rm module
+==============================
+
+.. automodule:: mapchete.cli.default.rm
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.cli.default.rst b/doc/source/apidoc/mapchete.cli.default.rst
index 3363bb6f..31ee79db 100644
--- a/doc/source/apidoc/mapchete.cli.default.rst
+++ b/doc/source/apidoc/mapchete.cli.default.rst
@@ -5,19 +5,22 @@ Submodules
 ----------
 
 .. toctree::
+   :maxdepth: 4
 
    mapchete.cli.default.convert
+   mapchete.cli.default.cp
    mapchete.cli.default.create
    mapchete.cli.default.execute
    mapchete.cli.default.formats
    mapchete.cli.default.index
    mapchete.cli.default.processes
+   mapchete.cli.default.rm
    mapchete.cli.default.serve
 
 Module contents
 ---------------
 
 .. automodule:: mapchete.cli.default
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.cli.default.serve.rst b/doc/source/apidoc/mapchete.cli.default.serve.rst
index 75c55859..f456ca98 100644
--- a/doc/source/apidoc/mapchete.cli.default.serve.rst
+++ b/doc/source/apidoc/mapchete.cli.default.serve.rst
@@ -2,6 +2,6 @@ mapchete.cli.default.serve module
 =================================
 
 .. automodule:: mapchete.cli.default.serve
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.cli.main.rst b/doc/source/apidoc/mapchete.cli.main.rst
index ffceb2bb..53f9b791 100644
--- a/doc/source/apidoc/mapchete.cli.main.rst
+++ b/doc/source/apidoc/mapchete.cli.main.rst
@@ -2,6 +2,6 @@ mapchete.cli.main module
 ========================
 
 .. automodule:: mapchete.cli.main
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.cli.options.rst b/doc/source/apidoc/mapchete.cli.options.rst
new file mode 100644
index 00000000..bbbcf996
--- /dev/null
+++ b/doc/source/apidoc/mapchete.cli.options.rst
@@ -0,0 +1,7 @@
+mapchete.cli.options module
+===========================
+
+.. automodule:: mapchete.cli.options
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.cli.rst b/doc/source/apidoc/mapchete.cli.rst
index 0e976aba..1685ca23 100644
--- a/doc/source/apidoc/mapchete.cli.rst
+++ b/doc/source/apidoc/mapchete.cli.rst
@@ -5,21 +5,23 @@ Subpackages
 -----------
 
 .. toctree::
+   :maxdepth: 4
 
-    mapchete.cli.default
+   mapchete.cli.default
 
 Submodules
 ----------
 
 .. toctree::
+   :maxdepth: 4
 
    mapchete.cli.main
-   mapchete.cli.utils
+   mapchete.cli.options
 
 Module contents
 ---------------
 
 .. automodule:: mapchete.cli
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.cli.utils.rst b/doc/source/apidoc/mapchete.cli.utils.rst
index 075f5fe3..a1c7c7f7 100644
--- a/doc/source/apidoc/mapchete.cli.utils.rst
+++ b/doc/source/apidoc/mapchete.cli.utils.rst
@@ -2,6 +2,6 @@ mapchete.cli.utils module
 =========================
 
 .. automodule:: mapchete.cli.utils
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.commands.rst b/doc/source/apidoc/mapchete.commands.rst
new file mode 100644
index 00000000..e0bbc757
--- /dev/null
+++ b/doc/source/apidoc/mapchete.commands.rst
@@ -0,0 +1,10 @@
+mapchete.commands package
+=========================
+
+Module contents
+---------------
+
+.. automodule:: mapchete.commands
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.commons.clip.rst b/doc/source/apidoc/mapchete.commons.clip.rst
index adfed0c2..a6d0ab1a 100644
--- a/doc/source/apidoc/mapchete.commons.clip.rst
+++ b/doc/source/apidoc/mapchete.commons.clip.rst
@@ -2,6 +2,6 @@ mapchete.commons.clip module
 ============================
 
 .. automodule:: mapchete.commons.clip
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.commons.contours.rst b/doc/source/apidoc/mapchete.commons.contours.rst
index 56fb50dd..5092f78a 100644
--- a/doc/source/apidoc/mapchete.commons.contours.rst
+++ b/doc/source/apidoc/mapchete.commons.contours.rst
@@ -2,6 +2,6 @@ mapchete.commons.contours module
 ================================
 
 .. automodule:: mapchete.commons.contours
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.commons.hillshade.rst b/doc/source/apidoc/mapchete.commons.hillshade.rst
index 38f9b3c0..002b06ce 100644
--- a/doc/source/apidoc/mapchete.commons.hillshade.rst
+++ b/doc/source/apidoc/mapchete.commons.hillshade.rst
@@ -2,6 +2,6 @@ mapchete.commons.hillshade module
 =================================
 
 .. automodule:: mapchete.commons.hillshade
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.commons.rst b/doc/source/apidoc/mapchete.commons.rst
index 081bd4b2..891673ef 100644
--- a/doc/source/apidoc/mapchete.commons.rst
+++ b/doc/source/apidoc/mapchete.commons.rst
@@ -5,6 +5,7 @@ Submodules
 ----------
 
 .. toctree::
+   :maxdepth: 4
 
    mapchete.commons.clip
    mapchete.commons.contours
@@ -14,6 +15,6 @@ Module contents
 ---------------
 
 .. automodule:: mapchete.commons
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.config.rst b/doc/source/apidoc/mapchete.config.rst
index 8ef00a47..1e0562d2 100644
--- a/doc/source/apidoc/mapchete.config.rst
+++ b/doc/source/apidoc/mapchete.config.rst
@@ -2,6 +2,6 @@ mapchete.config module
 ======================
 
 .. automodule:: mapchete.config
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.errors.rst b/doc/source/apidoc/mapchete.errors.rst
index 1947953e..e0d0e2bd 100644
--- a/doc/source/apidoc/mapchete.errors.rst
+++ b/doc/source/apidoc/mapchete.errors.rst
@@ -2,6 +2,6 @@ mapchete.errors module
 ======================
 
 .. automodule:: mapchete.errors
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.formats.base.rst b/doc/source/apidoc/mapchete.formats.base.rst
index 0baf06f4..0047faf3 100644
--- a/doc/source/apidoc/mapchete.formats.base.rst
+++ b/doc/source/apidoc/mapchete.formats.base.rst
@@ -2,6 +2,6 @@ mapchete.formats.base module
 ============================
 
 .. automodule:: mapchete.formats.base
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.formats.default.flatgeobuf.rst b/doc/source/apidoc/mapchete.formats.default.flatgeobuf.rst
new file mode 100644
index 00000000..ac42b5ae
--- /dev/null
+++ b/doc/source/apidoc/mapchete.formats.default.flatgeobuf.rst
@@ -0,0 +1,7 @@
+mapchete.formats.default.flatgeobuf module
+==========================================
+
+.. automodule:: mapchete.formats.default.flatgeobuf
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.formats.default.geobuf.rst b/doc/source/apidoc/mapchete.formats.default.geobuf.rst
new file mode 100644
index 00000000..70588a24
--- /dev/null
+++ b/doc/source/apidoc/mapchete.formats.default.geobuf.rst
@@ -0,0 +1,7 @@
+mapchete.formats.default.geobuf module
+======================================
+
+.. automodule:: mapchete.formats.default.geobuf
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.formats.default.geojson.rst b/doc/source/apidoc/mapchete.formats.default.geojson.rst
index b1b20d42..9c0e6111 100644
--- a/doc/source/apidoc/mapchete.formats.default.geojson.rst
+++ b/doc/source/apidoc/mapchete.formats.default.geojson.rst
@@ -2,6 +2,6 @@ mapchete.formats.default.geojson module
 =======================================
 
 .. automodule:: mapchete.formats.default.geojson
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.formats.default.gtiff.rst b/doc/source/apidoc/mapchete.formats.default.gtiff.rst
index 49b66a18..d6907f9d 100644
--- a/doc/source/apidoc/mapchete.formats.default.gtiff.rst
+++ b/doc/source/apidoc/mapchete.formats.default.gtiff.rst
@@ -2,6 +2,6 @@ mapchete.formats.default.gtiff module
 =====================================
 
 .. automodule:: mapchete.formats.default.gtiff
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.formats.default.mapchete_input.rst b/doc/source/apidoc/mapchete.formats.default.mapchete_input.rst
index 5c572a76..11d84fab 100644
--- a/doc/source/apidoc/mapchete.formats.default.mapchete_input.rst
+++ b/doc/source/apidoc/mapchete.formats.default.mapchete_input.rst
@@ -2,6 +2,6 @@ mapchete.formats.default.mapchete\_input module
 ===============================================
 
 .. automodule:: mapchete.formats.default.mapchete_input
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.formats.default.png.rst b/doc/source/apidoc/mapchete.formats.default.png.rst
index dc4358e7..9bf633c6 100644
--- a/doc/source/apidoc/mapchete.formats.default.png.rst
+++ b/doc/source/apidoc/mapchete.formats.default.png.rst
@@ -2,6 +2,6 @@ mapchete.formats.default.png module
 ===================================
 
 .. automodule:: mapchete.formats.default.png
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.formats.default.png_hillshade.rst b/doc/source/apidoc/mapchete.formats.default.png_hillshade.rst
index e5f2aca9..ad0bc7b1 100644
--- a/doc/source/apidoc/mapchete.formats.default.png_hillshade.rst
+++ b/doc/source/apidoc/mapchete.formats.default.png_hillshade.rst
@@ -2,6 +2,6 @@ mapchete.formats.default.png\_hillshade module
 ==============================================
 
 .. automodule:: mapchete.formats.default.png_hillshade
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.formats.default.raster_file.rst b/doc/source/apidoc/mapchete.formats.default.raster_file.rst
index fa66bb8d..4cea8760 100644
--- a/doc/source/apidoc/mapchete.formats.default.raster_file.rst
+++ b/doc/source/apidoc/mapchete.formats.default.raster_file.rst
@@ -2,6 +2,6 @@ mapchete.formats.default.raster\_file module
 ============================================
 
 .. automodule:: mapchete.formats.default.raster_file
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.formats.default.rst b/doc/source/apidoc/mapchete.formats.default.rst
index 413052b5..acac96d6 100644
--- a/doc/source/apidoc/mapchete.formats.default.rst
+++ b/doc/source/apidoc/mapchete.formats.default.rst
@@ -5,7 +5,10 @@ Submodules
 ----------
 
 .. toctree::
+   :maxdepth: 4
 
+   mapchete.formats.default.flatgeobuf
+   mapchete.formats.default.geobuf
    mapchete.formats.default.geojson
    mapchete.formats.default.gtiff
    mapchete.formats.default.mapchete_input
@@ -19,6 +22,6 @@ Module contents
 ---------------
 
 .. automodule:: mapchete.formats.default
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.formats.default.tile_directory.rst b/doc/source/apidoc/mapchete.formats.default.tile_directory.rst
index 3867c687..4b83da4c 100644
--- a/doc/source/apidoc/mapchete.formats.default.tile_directory.rst
+++ b/doc/source/apidoc/mapchete.formats.default.tile_directory.rst
@@ -2,6 +2,6 @@ mapchete.formats.default.tile\_directory module
 ===============================================
 
 .. automodule:: mapchete.formats.default.tile_directory
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.formats.default.vector_file.rst b/doc/source/apidoc/mapchete.formats.default.vector_file.rst
index d18b36f1..d25edb3e 100644
--- a/doc/source/apidoc/mapchete.formats.default.vector_file.rst
+++ b/doc/source/apidoc/mapchete.formats.default.vector_file.rst
@@ -2,6 +2,6 @@ mapchete.formats.default.vector\_file module
 ============================================
 
 .. automodule:: mapchete.formats.default.vector_file
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.formats.drivers.rst b/doc/source/apidoc/mapchete.formats.drivers.rst
index 97a4989c..d3da5074 100644
--- a/doc/source/apidoc/mapchete.formats.drivers.rst
+++ b/doc/source/apidoc/mapchete.formats.drivers.rst
@@ -2,6 +2,6 @@ mapchete.formats.drivers module
 ===============================
 
 .. automodule:: mapchete.formats.drivers
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.formats.rst b/doc/source/apidoc/mapchete.formats.rst
index 393d555c..0bb40123 100644
--- a/doc/source/apidoc/mapchete.formats.rst
+++ b/doc/source/apidoc/mapchete.formats.rst
@@ -5,13 +5,15 @@ Subpackages
 -----------
 
 .. toctree::
+   :maxdepth: 4
 
-    mapchete.formats.default
+   mapchete.formats.default
 
 Submodules
 ----------
 
 .. toctree::
+   :maxdepth: 4
 
    mapchete.formats.base
    mapchete.formats.drivers
@@ -20,6 +22,6 @@ Module contents
 ---------------
 
 .. automodule:: mapchete.formats
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.index.rst b/doc/source/apidoc/mapchete.index.rst
index f9c179ce..4f966b89 100644
--- a/doc/source/apidoc/mapchete.index.rst
+++ b/doc/source/apidoc/mapchete.index.rst
@@ -2,6 +2,6 @@ mapchete.index module
 =====================
 
 .. automodule:: mapchete.index
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.io.raster.rst b/doc/source/apidoc/mapchete.io.raster.rst
index 20c491d6..8071613f 100644
--- a/doc/source/apidoc/mapchete.io.raster.rst
+++ b/doc/source/apidoc/mapchete.io.raster.rst
@@ -2,6 +2,6 @@ mapchete.io.raster module
 =========================
 
 .. automodule:: mapchete.io.raster
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.io.rst b/doc/source/apidoc/mapchete.io.rst
index c6b4a279..02ba1764 100644
--- a/doc/source/apidoc/mapchete.io.rst
+++ b/doc/source/apidoc/mapchete.io.rst
@@ -5,6 +5,7 @@ Submodules
 ----------
 
 .. toctree::
+   :maxdepth: 4
 
    mapchete.io.raster
    mapchete.io.vector
@@ -13,6 +14,6 @@ Module contents
 ---------------
 
 .. automodule:: mapchete.io
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.io.vector.rst b/doc/source/apidoc/mapchete.io.vector.rst
index 7ae78eaa..a2113ee6 100644
--- a/doc/source/apidoc/mapchete.io.vector.rst
+++ b/doc/source/apidoc/mapchete.io.vector.rst
@@ -2,6 +2,6 @@ mapchete.io.vector module
 =========================
 
 .. automodule:: mapchete.io.vector
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.log.rst b/doc/source/apidoc/mapchete.log.rst
index 36bf698f..4cd4f933 100644
--- a/doc/source/apidoc/mapchete.log.rst
+++ b/doc/source/apidoc/mapchete.log.rst
@@ -2,6 +2,6 @@ mapchete.log module
 ===================
 
 .. automodule:: mapchete.log
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.processes.contours.rst b/doc/source/apidoc/mapchete.processes.contours.rst
index 9d689fb6..cf2df7d8 100644
--- a/doc/source/apidoc/mapchete.processes.contours.rst
+++ b/doc/source/apidoc/mapchete.processes.contours.rst
@@ -2,6 +2,6 @@ mapchete.processes.contours module
 ==================================
 
 .. automodule:: mapchete.processes.contours
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.processes.convert.rst b/doc/source/apidoc/mapchete.processes.convert.rst
index a8a5b70c..e30a25d4 100644
--- a/doc/source/apidoc/mapchete.processes.convert.rst
+++ b/doc/source/apidoc/mapchete.processes.convert.rst
@@ -2,6 +2,6 @@ mapchete.processes.convert module
 =================================
 
 .. automodule:: mapchete.processes.convert
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.processes.examples.example_process.rst b/doc/source/apidoc/mapchete.processes.examples.example_process.rst
index d0de5253..cc2ef372 100644
--- a/doc/source/apidoc/mapchete.processes.examples.example_process.rst
+++ b/doc/source/apidoc/mapchete.processes.examples.example_process.rst
@@ -2,6 +2,6 @@ mapchete.processes.examples.example\_process module
 ===================================================
 
 .. automodule:: mapchete.processes.examples.example_process
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.processes.examples.rst b/doc/source/apidoc/mapchete.processes.examples.rst
index 607c9815..39860691 100644
--- a/doc/source/apidoc/mapchete.processes.examples.rst
+++ b/doc/source/apidoc/mapchete.processes.examples.rst
@@ -5,6 +5,7 @@ Submodules
 ----------
 
 .. toctree::
+   :maxdepth: 4
 
    mapchete.processes.examples.example_process
 
@@ -12,6 +13,6 @@ Module contents
 ---------------
 
 .. automodule:: mapchete.processes.examples
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.processes.hillshade.rst b/doc/source/apidoc/mapchete.processes.hillshade.rst
index e330936e..8049c7a9 100644
--- a/doc/source/apidoc/mapchete.processes.hillshade.rst
+++ b/doc/source/apidoc/mapchete.processes.hillshade.rst
@@ -2,6 +2,6 @@ mapchete.processes.hillshade module
 ===================================
 
 .. automodule:: mapchete.processes.hillshade
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.processes.rst b/doc/source/apidoc/mapchete.processes.rst
index 0b3135d6..887d4ce6 100644
--- a/doc/source/apidoc/mapchete.processes.rst
+++ b/doc/source/apidoc/mapchete.processes.rst
@@ -5,13 +5,15 @@ Subpackages
 -----------
 
 .. toctree::
+   :maxdepth: 4
 
-    mapchete.processes.examples
+   mapchete.processes.examples
 
 Submodules
 ----------
 
 .. toctree::
+   :maxdepth: 4
 
    mapchete.processes.contours
    mapchete.processes.convert
@@ -21,6 +23,6 @@ Module contents
 ---------------
 
 .. automodule:: mapchete.processes
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.rst b/doc/source/apidoc/mapchete.rst
index 062cfe6d..588a578d 100644
--- a/doc/source/apidoc/mapchete.rst
+++ b/doc/source/apidoc/mapchete.rst
@@ -5,17 +5,20 @@ Subpackages
 -----------
 
 .. toctree::
+   :maxdepth: 4
 
-    mapchete.cli
-    mapchete.commons
-    mapchete.formats
-    mapchete.io
-    mapchete.processes
+   mapchete.cli
+   mapchete.commands
+   mapchete.commons
+   mapchete.formats
+   mapchete.io
+   mapchete.processes
 
 Submodules
 ----------
 
 .. toctree::
+   :maxdepth: 4
 
    mapchete.config
    mapchete.errors
@@ -28,6 +31,6 @@ Module contents
 ---------------
 
 .. automodule:: mapchete
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.tile.rst b/doc/source/apidoc/mapchete.tile.rst
index 80b55b94..f58e04bd 100644
--- a/doc/source/apidoc/mapchete.tile.rst
+++ b/doc/source/apidoc/mapchete.tile.rst
@@ -2,6 +2,6 @@ mapchete.tile module
 ====================
 
 .. automodule:: mapchete.tile
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/doc/source/apidoc/mapchete.validate.rst b/doc/source/apidoc/mapchete.validate.rst
index 1bff8c2a..02673e34 100644
--- a/doc/source/apidoc/mapchete.validate.rst
+++ b/doc/source/apidoc/mapchete.validate.rst
@@ -2,6 +2,6 @@ mapchete.validate module
 ========================
 
 .. automodule:: mapchete.validate
-    :members:
-    :undoc-members:
-    :show-inheritance:
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/mapchete/__init__.py b/mapchete/__init__.py
index eaca058b..65c606aa 100644
--- a/mapchete/__init__.py
+++ b/mapchete/__init__.py
@@ -1,13 +1,14 @@
 import logging
 
 from mapchete._core import open, Mapchete
-from mapchete._processing import MapcheteProcess
+from mapchete._executor import Executor
+from mapchete._processing import Job, MapcheteProcess, ProcessInfo
 from mapchete.tile import count_tiles
 from mapchete._timer import Timer
 
 
-__all__ = ["open", "count_tiles", "Mapchete", "MapcheteProcess", "Timer"]
-__version__ = "0.40"
+__all__ = ["open", "count_tiles", "Mapchete", "MapcheteProcess", "ProcessInfo", "Timer"]
+__version__ = "0.43"
 
 
 logger = logging.getLogger(__name__)
diff --git a/mapchete/_core.py b/mapchete/_core.py
index 7d4719fc..f5685e78 100644
--- a/mapchete/_core.py
+++ b/mapchete/_core.py
@@ -5,14 +5,16 @@
 import multiprocessing
 import os
 import threading
+import warnings
 
-from mapchete.config import MapcheteConfig
+from mapchete.config import MapcheteConfig, MULTIPROCESSING_DEFAULT_START_METHOD
 from mapchete.errors import MapcheteNodataTile
 from mapchete.formats import read_output_metadata
 from mapchete.io import fs_from_path, tiles_exist
 from mapchete._processing import (
     _run_on_single_tile,
     _run_area,
+    _preprocess,
     ProcessInfo,
     TileProcess,
 )
@@ -76,15 +78,15 @@ def open(some_input, with_cache=False, fs=None, fs_kwargs=None, **kwargs):
                 path=some_input,
                 fs=fs,
                 fs_kwargs=fs_kwargs,
-                **kwargs
+                **kwargs,
             ),
             config_dir=os.getcwd(),
             zoom_levels=kwargs.get("zoom"),
         )
         kwargs.update(mode="readonly")
         return Mapchete(MapcheteConfig(config, **kwargs))
-    else:
-        return Mapchete(MapcheteConfig(some_input, **kwargs), with_cache=with_cache)
+
+    return Mapchete(MapcheteConfig(some_input, **kwargs), with_cache=with_cache)
 
 
 class Mapchete(object):
@@ -155,9 +157,9 @@ def get_process_tiles(self, zoom=None):
             ):
                 yield tile
         else:
-            for zoom in reversed(self.config.zoom_levels):
+            for i in reversed(self.config.zoom_levels):
                 for tile in self.config.process_pyramid.tiles_from_geom(
-                    self.config.area_at_zoom(zoom), zoom
+                    self.config.area_at_zoom(i), i
                 ):
                     yield tile
 
@@ -184,15 +186,75 @@ def skip_tiles(self, tiles=None):
             for tile in tiles:
                 yield (tile, False)
 
+    def batch_preprocessor(
+        self,
+        dask_scheduler=None,
+        workers=None,
+        executor=None,
+    ):
+        """
+        Run all required preprocessing steps and yield over results.
+
+        The task count can be determined by self.config.preprocessing_tasks_count().
+
+        Parameters
+        ----------
+
+        dask_schedulter : str
+            URL to a dask scheduler if distributed execution is desired
+        workers : int
+            number of workers to be used for local processing
+        executor : mapchete.Executor
+            optional executor class to be used for processing
+
+        """
+        # process everything using executor and yield from results
+        yield from _preprocess(
+            self.config.preprocessing_tasks(),
+            process=self,
+            dask_scheduler=dask_scheduler,
+            workers=workers,
+            executor=executor,
+        )
+
+    def batch_preprocess(
+        self,
+        dask_scheduler=None,
+        workers=None,
+        executor=None,
+    ):
+        """
+        Run all required preprocessing steps.
+
+        Parameters
+        ----------
+
+        dask_schedulter : str
+            URL to a dask scheduler if distributed execution is desired
+        workers : int
+            number of workers to be used for local processing
+        executor : mapchete.Executor
+            optional executor class to be used for processing
+        """
+        list(
+            self.batch_preprocessor(
+                dask_scheduler=dask_scheduler,
+                workers=workers,
+                executor=executor,
+            )
+        )
+
     def batch_process(
         self,
         zoom=None,
         tile=None,
+        dask_scheduler=None,
         multi=None,
-        max_chunksize=1,
+        workers=None,
         multiprocessing_module=None,
-        multiprocessing_start_method="fork",
+        multiprocessing_start_method=MULTIPROCESSING_DEFAULT_START_METHOD,
         skip_output_check=False,
+        executor=None,
     ):
         """
         Process a large batch of tiles.
@@ -207,30 +269,31 @@ def batch_process(
         tile : tuple
             zoom, row and column of tile to be processed (cannot be used with
             zoom)
-        multi : int
+        workers : int
             number of workers (default: number of CPU cores)
-        max_chunksize : int
-            maximum number of process tiles to be queued for each worker;
-            (default: 1)
         multiprocessing_module : module
             either Python's standard 'multiprocessing' or Celery's 'billiard' module
             (default: multiprocessing)
         multiprocessing_start_method : str
             "fork", "forkserver" or "spawn"
-            (default: "fork")
+            (default: "spawn")
         skip_output_check : bool
             skip checking whether process tiles already have existing output before
             starting to process;
+        executor : mapchete.Executor
+            optional executor class to be used for processing
         """
         list(
             self.batch_processor(
                 zoom=zoom,
                 tile=tile,
-                multi=multi or multiprocessing.cpu_count(),
-                max_chunksize=max_chunksize,
-                multiprocessing_module=multiprocessing_module or multiprocessing,
+                dask_scheduler=dask_scheduler,
+                workers=workers,
+                multi=multi,
+                multiprocessing_module=multiprocessing_module,
                 multiprocessing_start_method=multiprocessing_start_method,
                 skip_output_check=skip_output_check,
+                executor=executor,
             )
         )
 
@@ -238,11 +301,13 @@ def batch_processor(
         self,
         zoom=None,
         tile=None,
+        dask_scheduler=None,
         multi=None,
-        max_chunksize=1,
+        workers=None,
         multiprocessing_module=None,
-        multiprocessing_start_method="fork",
+        multiprocessing_start_method=MULTIPROCESSING_DEFAULT_START_METHOD,
         skip_output_check=False,
+        executor=None,
     ):
         """
         Process a large batch of tiles and yield report messages per tile.
@@ -257,63 +322,97 @@ def batch_processor(
             zoom)
         multi : int
             number of workers (default: number of CPU cores)
-        max_chunksize : int
-            maximum number of process tiles to be queued for each worker;
-            (default: 1)
         multiprocessing_module : module
             either Python's standard 'multiprocessing' or Celery's 'billiard' module
             (default: multiprocessing)
         multiprocessing_start_method : str
             "fork", "forkserver" or "spawn"
-            (default: "fork")
+            (default: "spawn")
         skip_output_check : bool
             skip checking whether process tiles already have existing output before
             starting to process;
+        executor : mapchete.Executor
+            optional executor class to be used for processing
         """
+        if multi is not None:  # pragma: no cover
+            warnings.warn(
+                DeprecationWarning(
+                    "the 'multi' keyword is deprecated and should be called 'workers'"
+                )
+            )
+            workers = workers or multi
+
         if zoom and tile:
             raise ValueError("use either zoom or tile")
 
         # run single tile
         if tile:
             yield _run_on_single_tile(
-                process=self, tile=self.config.process_pyramid.tile(*tuple(tile))
+                executor=executor,
+                process=self,
+                tile=self.config.process_pyramid.tile(*tuple(tile)),
             )
         # run area
         else:
             for process_info in _run_area(
                 process=self,
                 zoom_levels=list(_get_zoom_level(zoom, self)),
-                multi=multi or multiprocessing.cpu_count(),
-                max_chunksize=max_chunksize,
+                dask_scheduler=dask_scheduler,
+                workers=workers or multiprocessing.cpu_count(),
                 multiprocessing_module=multiprocessing_module or multiprocessing,
                 multiprocessing_start_method=multiprocessing_start_method,
                 skip_output_check=skip_output_check,
+                executor=executor,
             ):
                 yield process_info
 
-    def count_tiles(self, minzoom, maxzoom, init_zoom=0):
+    def count_tasks(self, minzoom=None, maxzoom=None, init_zoom=0):
+        """
+        Count all preprocessing tasks and tiles at given zoom levels.
+
+        Parameters
+        ----------
+        minzoom : int
+            limits minimum process zoom
+        maxzoom : int
+            limits maximum process zoom
+        init_zoom : int
+            initial zoom level used for tile count algorithm
+
+        Returns
+        -------
+        number of tasks
+        """
+        return self.config.preprocessing_tasks_count() + self.count_tiles(
+            minzoom=minzoom, maxzoom=maxzoom, init_zoom=0
+        )
+
+    def count_tiles(self, minzoom=None, maxzoom=None, init_zoom=0):
         """
         Count number of tiles intersecting with process area at given zoom levels.
 
         Parameters
         ----------
-        geometry : shapely geometry
-        pyramid : TilePyramid
         minzoom : int
+            limits minimum process zoom
         maxzoom : int
+            limits maximum process zoom
         init_zoom : int
+            initial zoom level used for tile count algorithm
 
         Returns
         -------
         number of tiles
         """
+        minzoom = min(self.config.init_zoom_levels) if minzoom is None else minzoom
+        maxzoom = max(self.config.init_zoom_levels) if maxzoom is None else maxzoom
         if (minzoom, maxzoom) not in self._count_tiles_cache:
             self._count_tiles_cache[(minzoom, maxzoom)] = count_tiles(
                 self.config.area_at_zoom(),
                 self.config.process_pyramid,
                 minzoom,
                 maxzoom,
-                init_zoom=0,
+                init_zoom=init_zoom,
             )
         return self._count_tiles_cache[(minzoom, maxzoom)]
 
@@ -335,6 +434,7 @@ def execute(self, process_tile, raise_nodata=False):
             process output
         """
         process_tile = validate_tile(process_tile, self.config.process_pyramid)
+        self.batch_preprocess()
         try:
             return self.config.output.streamline_output(
                 TileProcess(tile=process_tile, config=self.config).execute()
@@ -342,8 +442,7 @@ def execute(self, process_tile, raise_nodata=False):
         except MapcheteNodataTile:
             if raise_nodata:  # pragma: no cover
                 raise
-            else:
-                return self.config.output.empty(process_tile)
+            return self.config.output.empty(process_tile)
 
     def read(self, output_tile):
         """
@@ -543,8 +642,10 @@ def __exit__(self, exc_type, exc_value, exc_traceback):
         # run input drivers cleanup
         for ip in self.config.input.values():
             if ip is not None:
+                logger.debug(f"running cleanup on {ip}...")
                 ip.cleanup()
         # run output driver cleanup
+        logger.debug(f"closing output driver {self.config.output}...")
         self.config.output.close(
             exc_type=exc_type, exc_value=exc_value, exc_traceback=exc_traceback
         )
diff --git a/mapchete/_executor.py b/mapchete/_executor.py
new file mode 100644
index 00000000..2c4ccae9
--- /dev/null
+++ b/mapchete/_executor.py
@@ -0,0 +1,305 @@
+from cached_property import cached_property
+import concurrent.futures
+from functools import partial
+from itertools import chain
+import logging
+import multiprocessing
+import os
+import sys
+import warnings
+
+from mapchete.config import MULTIPROCESSING_DEFAULT_START_METHOD
+from mapchete.log import set_log_level
+
+
+logger = logging.getLogger(__name__)
+
+
+class Executor:
+    """
+    Executor factory for dask and concurrent.futures executor
+
+    Will move into the mapchete core package.
+    """
+
+    def __new__(self, *args, concurrency=None, **kwargs):
+        if concurrency == "dask":
+            try:
+                return DaskExecutor(*args, **kwargs)
+            except ImportError as e:  # pragma: no cover
+                raise ImportError(
+                    f"this feature requires the mapchete[dask] extra: {e}"
+                )
+        elif concurrency is None or kwargs.get("max_workers") == 1:
+            return SequentialExecutor(*args, **kwargs)
+        elif concurrency in ["processes", "threads"]:
+            return ConcurrentFuturesExecutor(*args, concurrency=concurrency, **kwargs)
+        else:  # pragma: no cover
+            raise ValueError(
+                f"concurrency must be one of None, 'processes', 'threads' or 'dask', not {concurrency}"
+            )
+
+
+class _ExecutorBase:
+    """Define base methods and properties of executors."""
+
+    futures = None
+    _as_completed = None
+    _executor = None
+    _executor_cls = None
+    _executor_args = ()
+    _executor_kwargs = {}
+
+    def as_completed(self, func, iterable, fargs=None, fkwargs=None):
+        """Yield finished tasks."""
+        fargs = fargs or ()
+        fkwargs = fkwargs or {}
+        logger.debug("submitting tasks to executor")
+        futures = [
+            self._executor.submit(func, *chain([item], fargs), **fkwargs)
+            for item in iterable
+        ]
+        self.futures.extend(futures)
+        logger.debug(f"added {len(futures)} tasks")
+        for future in self._as_completed(futures):
+            if future.exception():  # pragma: no cover
+                raise future.exception()
+            yield future
+
+    def map(self, func, iterable, fargs=None, fkwargs=None):  # pragma: no cover
+        raise NotImplementedError
+
+    def cancel(self):
+        logger.debug(f"cancel {len(self.futures)} futures...")
+        for future in self.futures:
+            future.cancel()
+        logger.debug(f"{len(self.futures)} futures cancelled")
+
+    def close(self):  # pragma: no cover
+        self.__exit__(None, None, None)
+
+    def _as_completed(self, *args, **kwargs):  # pragma: no cover
+        raise NotImplementedError()
+
+    @cached_property
+    def _executor(self):
+        return self._executor_cls(*self._executor_args, **self._executor_kwargs)
+
+    def __enter__(self):
+        """Enter context manager."""
+        return self
+
+    def __exit__(self, *args):
+        """Exit context manager."""
+        logger.debug(f"closing executor {self._executor}...")
+        try:
+            self._executor.close()
+        except Exception:
+            self._executor.__exit__(*args)
+        logger.debug(f"closed executor {self._executor}")
+
+    def __repr__(self):  # pragma: no cover
+        return f"<Executor ({self._executor_cls})>"
+
+
+class DaskExecutor(_ExecutorBase):
+    """Execute tasks using dask cluster."""
+
+    def __init__(
+        self,
+        *args,
+        address=None,
+        dask_scheduler=None,
+        dask_client=None,
+        max_workers=None,
+        **kwargs,
+    ):
+        from dask.distributed import Client, LocalCluster
+
+        self.futures = []
+        self._executor_client = dask_client
+        if self._executor_client:  # pragma: no cover
+            logger.debug(f"using existing dask client: {dask_client}")
+        else:
+            local_cluster_kwargs = dict(
+                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),
+            )
+            logger.debug(
+                f"starting dask.distributed.Client with kwargs {self._executor_kwargs}"
+            )
+
+    def cancel(self):
+        logger.debug(f"cancel {len(self.futures)} futures...")
+        for future in self.futures:
+            future.cancel()
+        logger.debug(f"{len(self.futures)} futures cancelled")
+
+    def map(self, func, iterable, fargs=None, fkwargs=None):
+        fargs = fargs or []
+        fkwargs = fkwargs or {}
+        return [
+            f.result()
+            for f in self._executor.map(partial(func, *fargs, **fkwargs), iterable)
+        ]
+
+    def _as_completed(self, futures):
+        from dask.distributed import as_completed
+
+        if futures:
+            for future in as_completed(futures):
+                yield future
+
+    @cached_property
+    def _executor(self):
+        return self._executor_client or self._executor_cls(
+            *self._executor_args, **self._executor_kwargs
+        )
+
+    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")
+        else:
+            logger.debug(f"closing executor {self._executor}...")
+            try:
+                self._executor.close()
+            except Exception:  # pragma: no cover
+                self._executor.__exit__(*args)
+            logger.debug(f"closed executor {self._executor}")
+
+
+class ConcurrentFuturesExecutor(_ExecutorBase):
+    """Execute tasks using concurrent.futures."""
+
+    def __init__(
+        self,
+        *args,
+        max_workers=None,
+        concurrency="processes",
+        **kwargs,
+    ):
+        """Set attributes."""
+
+        self.max_workers = max_workers or os.cpu_count()
+        self.futures = []
+        self._executor_kwargs = dict(
+            max_workers=self.max_workers,
+        )
+        if sys.version_info >= (3, 7):
+            self._executor_kwargs.update(
+                initializer=set_log_level,
+                initargs=(logger.getEffectiveLevel(),),
+            )
+        else:  # pragma: no cover
+            warnings.warn(UserWarning("worker logs are not available on python<3.7"))
+        if concurrency == "processes":
+            self._executor_cls = concurrent.futures.ProcessPoolExecutor
+            if sys.version_info >= (3, 7):
+                start_method = (
+                    kwargs.get("multiprocessing_start_method")
+                    or MULTIPROCESSING_DEFAULT_START_METHOD
+                )
+                self._executor_kwargs.update(
+                    mp_context=multiprocessing.get_context(method=start_method)
+                )
+                if start_method != "spawn":  # pragma: no cover
+                    warnings.warn(
+                        UserWarning(
+                            "using a start method other than 'spawn' can cause mapchete to hang"
+                        )
+                    )
+            else:  # pragma: no cover
+                warnings.warn(
+                    UserWarning(
+                        "multiprocessing start method cannot be set on python<3.7"
+                    )
+                )
+        elif concurrency == "threads":
+            self._executor_cls = concurrent.futures.ThreadPoolExecutor
+        else:  # pragma: no cover
+            raise ValueError("concurrency must either be 'processes' or 'threads'")
+        logger.debug(
+            f"init ConcurrentFuturesExecutor using {concurrency} with {self.max_workers} workers"
+        )
+
+    def _as_completed(self, futures):
+        """Yield finished tasks."""
+        for future in concurrent.futures.as_completed(futures):
+            yield future
+
+    def map(self, func, iterable, fargs=None, fkwargs=None):
+        fargs = fargs or []
+        fkwargs = fkwargs or {}
+        return list(self._executor.map(partial(func, *fargs, **fkwargs), iterable))
+
+
+class SequentialExecutor(_ExecutorBase):
+    """Execute tasks sequentially in single process."""
+
+    def __init__(self, *args, **kwargs):
+        """Set attributes."""
+        logger.debug("init SequentialExecutor")
+        self.futures = []
+        self._cancel = False
+
+    def as_completed(self, func, iterable, fargs=None, fkwargs=None):
+        """Yield finished tasks."""
+        fargs = fargs or []
+        fkwargs = fkwargs or {}
+        for i in iterable:
+            if self._cancel:
+                return
+            yield FakeFuture(func, fargs=[i, *fargs], fkwargs=fkwargs)
+
+    def map(self, func, iterable, fargs=None, fkwargs=None):
+        fargs = fargs or []
+        fkwargs = fkwargs or {}
+        return list(map(partial(func, *fargs, **fkwargs), iterable))
+
+    def cancel(self):
+        self._cancel = True
+
+    def __exit__(self, *args):
+        """Exit context manager."""
+        logger.debug("SequentialExecutor closed")
+
+    def __repr__(self):  # pragma: no cover
+        """Return string representation."""
+        return "SequentialExecutor"
+
+
+class FakeFuture:
+    """Wrapper class to mimick future interface."""
+
+    def __init__(self, func, fargs=None, fkwargs=None):
+        """Set attributes."""
+        fargs = fargs or []
+        fkwargs = fkwargs or {}
+        try:
+            self._result, self._exception = func(*fargs, **fkwargs), None
+        except Exception as e:  # pragma: no cover
+            self._result, self._exception = None, e
+
+    def result(self):
+        """Return task result."""
+        if self._exception:
+            logger.exception(self._exception)
+            raise self._exception
+        else:
+            return self._result
+
+    def exception(self):
+        """Raise task exception if any."""
+        return self._exception
+
+    def cancelled(self):  # pragma: no cover
+        """Sequential futures cannot be cancelled."""
+        return False
+
+    def __repr__(self):  # pragma: no cover
+        """Return string representation."""
+        return f"FakeFuture(result={self._result}, exception={self._exception})"
diff --git a/mapchete/_processing.py b/mapchete/_processing.py
index dd51d967..1c214f94 100644
--- a/mapchete/_processing.py
+++ b/mapchete/_processing.py
@@ -1,15 +1,17 @@
 from collections import namedtuple
-from functools import partial
+from concurrent.futures._base import CancelledError
 from itertools import chain
 import logging
 import multiprocessing
 import os
 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
 from mapchete.errors import MapcheteNodataTile, MapcheteProcessException
 from mapchete.io import raster
 from mapchete._timer import Timer
@@ -21,6 +23,69 @@
 ProcessInfo = namedtuple("ProcessInfo", "tile processed process_msg written write_msg")
 
 
+class Job:
+    """
+    Wraps the output of a processing function into a generator with known length.
+
+    This class also exposes the internal Executor.cancel() function in order to cancel all remaining
+    tasks/futures.
+
+    Will move into the mapchete core package.
+    """
+
+    def __init__(
+        self,
+        func: Generator,
+        fargs: tuple = None,
+        fkwargs: dict = None,
+        as_iterator: bool = False,
+        total: int = None,
+        executor_concurrency: str = "processes",
+        executor_kwargs: dict = None,
+    ):
+        self.func = func
+        self.fargs = fargs or ()
+        self.fkwargs = fkwargs or {}
+        self.status = "pending"
+        self.executor = None
+        self.executor_concurrency = executor_concurrency
+        self.executor_kwargs = executor_kwargs or {}
+        self._total = total
+        self._as_iterator = as_iterator
+        if not as_iterator:
+            self._results = list(self._run())
+
+    def _run(self):
+        if self._total == 0:
+            return
+        with Executor(
+            concurrency=self.executor_concurrency, **self.executor_kwargs
+        ) as self.executor:
+            self.status = "running"
+            yield from self.func(*self.fargs, executor=self.executor, **self.fkwargs)
+            self.status = "finished"
+
+    def cancel(self):
+        if self._as_iterator:
+            # requires client and futures
+            if self.executor is None:  # pragma: no cover
+                raise ValueError("nothing to cancel because no executor is running")
+            self.executor.cancel()
+            self.status = "cancelled"
+
+    def __len__(self):
+        return self._total
+
+    def __iter__(self):
+        if self._as_iterator:
+            yield from self._run()
+        else:
+            return self._results
+
+    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.
@@ -36,7 +101,7 @@ def __init__(self, tile=None, config=None, skip=False):
         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_path = None if skip else config.process_path
+        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:
             self.input, self.process_func_params, self.output_params = {}, {}, {}
@@ -85,7 +150,7 @@ def _execute(self):
                 return self._interpolate_from_baselevel("higher")
         # Otherwise, execute from process file.
         process_func = get_process_func(
-            process_path=self.process_path, config_dir=self.config_dir
+            process=self.process, config_dir=self.config_dir
         )
         try:
             with Timer() as t:
@@ -97,7 +162,7 @@ def _execute(self):
                         input=self.input,
                         output_params=self.output_params,
                     ),
-                    **self.process_func_params
+                    **self.process_func_params,
                 )
         except MapcheteNodataTile:
             raise
@@ -343,101 +408,66 @@ def clip(self, array, geometries, inverted=False, clip_buffer=0):
         )
 
 
-#############################################################
-# wrappers helping to abstract multiprocessing and billiard #
-#############################################################
-class Executor:
-    """Wrapper class to be used with multiprocessing or billiard."""
+#######################
+# batch preprocessing #
+#######################
 
-    def __init__(
-        self,
-        start_method="fork",
-        max_workers=None,
-        multiprocessing_module=multiprocessing,
-    ):
-        """Set attributes."""
-        self.start_method = start_method
-        self.max_workers = max_workers or os.cpu_count()
-        self.multiprocessing_module = multiprocessing_module
-        if self.max_workers != 1:
-            logger.debug(
-                "init %s Executor with start_method %s and %s workers",
-                self.multiprocessing_module,
-                self.start_method,
-                self.max_workers,
-            )
-            self._pool = self.multiprocessing_module.get_context(
-                self.start_method
-            ).Pool(self.max_workers)
 
-    def as_completed(
-        self, func=None, iterable=None, fargs=None, fkwargs=None, chunksize=1
-    ):
-        """Yield finished tasks."""
-        fargs = fargs or []
-        fkwargs = fkwargs or {}
-        if self.max_workers == 1:
-            for i in iterable:
-                yield _exception_wrapper(func, fargs, fkwargs, i)
-        else:
-            logger.debug(
-                "open multiprocessing.Pool and %s %s workers",
-                self.start_method,
-                self.max_workers,
-            )
-            for finished_task in self._pool.imap_unordered(
-                partial(_exception_wrapper, func, fargs, fkwargs),
-                iterable,
-                chunksize=chunksize or 1,
-            ):
-                yield finished_task
-
-    def __enter__(self):
-        """Enter context manager."""
-        if self.max_workers != 1:
-            self._pool.__enter__()
-        return self
+def _preprocess_task_wrapper(task_tuple):
+    task_key, (func, fargs, fkwargs) = task_tuple
+    return task_key, func(*fargs, **fkwargs)
 
-    def __exit__(self, *args):
-        """Exit context manager."""
-        if self.max_workers != 1:
-            logger.debug("closing %s and workers", self._pool)
-            self._pool.__exit__(*args)
-            logger.debug("%s closed", self._pool)
 
+def _preprocess(
+    tasks,
+    process=None,
+    dask_scheduler=None,
+    workers=None,
+    multiprocessing_module=None,
+    multiprocessing_start_method=None,
+    executor=None,
+):
+    # If preprocessing tasks already finished, don't run them again.
+    if process.config.preprocessing_tasks_finished:  # pragma: no cover
+        return
+
+    # 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 t:
+            logger.debug(
+                "run preprocessing on %s tasks using %s workers", len(tasks), workers
+            )
 
-class FinishedTask:
-    """Wrapper class to encapsulate exceptions."""
-
-    def __init__(self, func, fargs=None, fkwargs=None):
-        """Set attributes."""
-        fargs = fargs or []
-        fkwargs = fkwargs or {}
-        try:
-            self._result, self._exception = func(*fargs, **fkwargs), None
-        except Exception as e:  # pragma: no cover
-            self._result, self._exception = None, e
-
-    def result(self):
-        """Return task result."""
-        if self._exception:
-            logger.exception(self._exception)
-            raise self._exception
-        else:
-            return self._result
-
-    def exception(self):
-        """Raise task exception if any."""
-        return self._exception
-
-    def __repr__(self):
-        """Return string representation."""
-        return "FinishedTask(result=%s, exception=%s)" % (self._result, self._exception)
+            # process all remaining tiles using todo list from before
+            for i, future in enumerate(
+                executor.as_completed(
+                    func=_preprocess_task_wrapper,
+                    iterable=[(k, v) for k, v in tasks.items()],
+                ),
+                1,
+            ):
+                task_key, result = future.result()
+                logger.debug(
+                    f"preprocessing task {i}/{len(tasks)} {task_key} processed successfully"
+                )
+                process.config.set_preprocessing_task_result(task_key, result)
+                yield f"preprocessing task {task_key} finished"
+    finally:
+        if create_executor:
+            executor.close()
 
+    process.config.preprocessing_tasks_finished = True
 
-def _exception_wrapper(func, fargs, fkwargs, i):
-    """Wrap function around FinishedTask object."""
-    return FinishedTask(func, list(chain([i], fargs)), fkwargs)
+    logger.info("%s task(s) iterated in %s", str(len(tasks)), t)
 
 
 ###########################
@@ -445,7 +475,7 @@ def _exception_wrapper(func, fargs, fkwargs, i):
 ###########################
 
 
-def _run_on_single_tile(process=None, tile=None):
+def _run_on_single_tile(executor=None, process=None, tile=None):
     logger.debug("run process on single tile")
     return _execute_and_write(
         tile_process=TileProcess(
@@ -461,10 +491,11 @@ def _run_on_single_tile(process=None, tile=None):
 
 
 def _run_area(
+    executor=None,
     process=None,
     zoom_levels=None,
-    multi=None,
-    max_chunksize=None,
+    dask_scheduler=None,
+    workers=None,
     multiprocessing_module=None,
     multiprocessing_start_method=None,
     skip_output_check=False,
@@ -475,13 +506,14 @@ 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(
+            executor=executor,
             func=_execute,
             zoom_levels=zoom_levels,
             process=process,
-            multi=multi,
+            dask_scheduler=dask_scheduler,
+            workers=workers,
             multiprocessing_start_method=multiprocessing_start_method,
             multiprocessing_module=multiprocessing_module,
-            max_chunksize=max_chunksize,
             write_in_parent_process=True,
             skip_output_check=skip_output_check,
         ):
@@ -490,14 +522,15 @@ def _run_area(
     # for output drivers which can write data in child processes
     else:
         for process_info in _run_multi(
+            executor=executor,
             func=_execute_and_write,
             fkwargs=dict(output_writer=process.config.output),
             zoom_levels=zoom_levels,
             process=process,
-            multi=multi,
+            dask_scheduler=dask_scheduler,
+            workers=workers,
             multiprocessing_start_method=multiprocessing_start_method,
             multiprocessing_module=multiprocessing_module,
-            max_chunksize=max_chunksize,
             write_in_parent_process=False,
             skip_output_check=skip_output_check,
         ):
@@ -520,52 +553,117 @@ def _filter_skipable(process=None, tiles=None, todo=None, target_set=None):
 
 
 def _run_multi(
+    executor=None,
     func=None,
     zoom_levels=None,
     process=None,
-    multi=None,
+    dask_scheduler=None,
+    workers=None,
     multiprocessing_start_method=None,
     multiprocessing_module=None,
-    max_chunksize=None,
     write_in_parent_process=False,
     fkwargs=None,
     skip_output_check=False,
 ):
     total_tiles = process.count_tiles(min(zoom_levels), max(zoom_levels))
-    workers = min([multi, total_tiles])
+    workers = min([workers, total_tiles])
     num_processed = 0
-    logger.debug("run process on %s tiles using %s workers", total_tiles, workers)
 
     # 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()
 
-    with Timer() as t, Executor(
+    # 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,
-    ) as executor:
-
-        for i, zoom in enumerate(zoom_levels):
-
-            if skip_output_check:  # pragma: no cover
-                # don't check outputs and simply proceed
-                todo = process.get_process_tiles(zoom)
-            else:
-                # check which process output already exists and which process tiles need
-                # to be added to todo list
-                todo = set()
-                logger.debug("check skippable tiles")
-                for process_info in _filter_skipable(
-                    process=process,
-                    tiles=process.get_process_tiles(zoom),
-                    todo=todo,
-                    target_set=(
-                        overview_parents if process.config.baselevels and i else None
+        dask_scheduler=dask_scheduler,
+    )
+
+    try:
+        with Timer() as t:
+            logger.debug(
+                "run process on %s tiles using %s workers", total_tiles, workers
+            )
+
+            for i, zoom in enumerate(zoom_levels):
+
+                if skip_output_check:  # pragma: no cover
+                    # don't check outputs and simply proceed
+                    todo = process.get_process_tiles(zoom)
+                else:
+                    # check which process output already exists and which process tiles need
+                    # to be added to todo list
+                    todo = set()
+                    logger.debug("check skippable tiles")
+                    for process_info in _filter_skipable(
+                        process=process,
+                        tiles=process.get_process_tiles(zoom),
+                        todo=todo,
+                        target_set=(
+                            overview_parents
+                            if process.config.baselevels and i
+                            else None
+                        ),
+                    ):
+                        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,
+                        )
+                        yield process_info
+
+                # process all remaining tiles using todo list from before
+                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,
+                        )
+                        for tile in todo
                     ),
+                    fkwargs=fkwargs,
                 ):
+                    # 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.info(
+                    logger.debug(
                         "tile %s/%s finished: %s, %s, %s",
                         num_processed,
                         total_tiles,
@@ -575,60 +673,11 @@ def _run_multi(
                     )
                     yield process_info
 
-            # process all remaining tiles using todo list from before
-            for task 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,
-                    )
-                    for tile in todo
-                ),
-                fkwargs=fkwargs,
-                chunksize=max_chunksize,
-            ):
-                # trigger output write for driver which require parent process for writing
-                if write_in_parent_process:
-                    output_data, process_info = task.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 = task.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.info(
-                    "tile %s/%s finished: %s, %s, %s",
-                    num_processed,
-                    total_tiles,
-                    process_info.tile,
-                    process_info.process_msg,
-                    process_info.write_msg,
-                )
-                yield process_info
+    finally:
+        if create_executor:
+            executor.close()
 
-    logger.debug("%s tile(s) iterated in %s", str(num_processed), t)
+    logger.info("%s tile(s) iterated in %s", str(num_processed), t)
 
 
 ###############################
@@ -636,7 +685,7 @@ def _run_multi(
 ###############################
 
 
-def _execute(tile_process=None):
+def _execute(tile_process=None, **kwargs):
     logger.debug(
         (tile_process.tile.id, "running on %s" % multiprocessing.current_process().name)
     )
@@ -670,7 +719,7 @@ def _execute(tile_process=None):
         )
 
 
-def _write(process_info=None, output_data=None, output_writer=None):
+def _write(process_info=None, output_data=None, output_writer=None, **kwargs):
     if process_info.processed:
         try:
             output_data = output_writer.streamline_output(output_data)
@@ -702,7 +751,7 @@ def _write(process_info=None, output_data=None, output_writer=None):
         return process_info
 
 
-def _execute_and_write(tile_process=None, output_writer=None):
+def _execute_and_write(tile_process=None, output_writer=None, **kwargs):
     output_data, process_info = _execute(tile_process=tile_process)
     return _write(
         process_info=process_info, output_data=output_data, output_writer=output_writer
diff --git a/mapchete/cli/default/convert.py b/mapchete/cli/default/convert.py
index a550884b..acfc80de 100644
--- a/mapchete/cli/default/convert.py
+++ b/mapchete/cli/default/convert.py
@@ -1,30 +1,15 @@
 import click
-import fiona
-import logging
-from multiprocessing import cpu_count
-import os
-from pprint import pformat
-import rasterio
-from rasterio.enums import Resampling
 from rasterio.dtypes import dtype_ranges
+from rasterio.enums import Resampling
 from rasterio.rio.options import creation_options
-from shapely.geometry import box
-import sys
 import tilematrix
+import tqdm
 
-from mapchete.cli import utils
-from mapchete.config import raw_conf, raw_conf_output_pyramid
-from mapchete.formats import (
-    driver_from_file,
-    available_output_formats,
-    available_input_formats,
-)
-from mapchete.io import read_json, get_best_zoom_level
-from mapchete.io.vector import reproject_geometry
-from mapchete.tile import BufferedTilePyramid
-from mapchete.validate import validate_zooms
+import mapchete
+from mapchete import commands
+from mapchete.cli import options
+from mapchete.formats import available_output_formats
 
-logger = logging.getLogger(__name__)
 OUTPUT_FORMATS = available_output_formats()
 
 
@@ -37,21 +22,20 @@ def _validate_bidx(ctx, param, bidx):
 
 
 @click.command(help="Convert outputs or other geodata.")
-@utils.arg_input
-@utils.arg_output
-@utils.opt_zoom
-@utils.opt_bounds
-@utils.opt_bounds_crs
-@utils.opt_area
-@utils.opt_area_crs
-@utils.opt_point
-@utils.opt_point_crs
-@utils.opt_wkt_geometry
+@options.arg_tiledir
+@options.arg_output
+@options.opt_zoom
+@options.opt_bounds
+@options.opt_bounds_crs
+@options.opt_area
+@options.opt_area_crs
+@options.opt_point
+@options.opt_point_crs
 @click.option(
     "--clip-geometry",
     "-c",
     type=click.Path(exists=True),
-    help="Clip output by geometry",
+    help="Clip output by geometry.",
 )
 @click.option("--bidx", callback=_validate_bidx, help="Band indexes to copy.")
 @click.option(
@@ -66,7 +50,7 @@ def _validate_bidx(ctx, param, bidx):
 )
 @click.option(
     "--output-format",
-    type=click.Choice(available_output_formats()),
+    type=click.Choice(OUTPUT_FORMATS),
     help="Output format.",
 )
 @click.option(
@@ -92,7 +76,7 @@ def _validate_bidx(ctx, param, bidx):
     default=0.0,
     help="Scaling offset (for raster output only).",
 )
-@utils.opt_resampling_method
+@options.opt_resampling_method
 @click.option(
     "--overviews", is_flag=True, help="Generate overviews (single GTiff output only)."
 )
@@ -107,286 +91,63 @@ def _validate_bidx(ctx, param, bidx):
     is_flag=True,
     help="Write a valid COG. This will automatically generate verviews. (GTiff only)",
 )
-@utils.opt_overwrite
-@utils.opt_verbose
-@utils.opt_no_pbar
-@utils.opt_debug
-@utils.opt_multi
-@utils.opt_logfile
-@utils.opt_vrt
-@utils.opt_idx_out_dir
+@options.opt_overwrite
+@options.opt_verbose
+@options.opt_no_pbar
+@options.opt_debug
+@options.opt_workers
+@options.opt_multi
+@options.opt_concurrency
+@options.opt_logfile
+@options.opt_vrt
+@options.opt_idx_out_dir
 def convert(
-    input_,
+    tiledir,
     output,
-    zoom=None,
-    bounds=None,
-    bounds_crs=None,
-    area=None,
-    area_crs=None,
-    point=None,
-    point_crs=None,
-    wkt_geometry=None,
-    clip_geometry=None,
-    bidx=None,
-    output_pyramid=None,
-    output_metatiling=None,
-    output_format=None,
-    output_dtype=None,
-    output_geometry_type=None,
-    creation_options=None,
-    scale_ratio=None,
-    scale_offset=None,
-    resampling_method=None,
-    overviews=False,
-    overviews_resampling_method=None,
-    cog=False,
-    overwrite=False,
-    logfile=None,
-    verbose=False,
-    no_pbar=False,
-    debug=False,
-    multi=None,
+    *args,
     vrt=False,
     idx_out_dir=None,
+    debug=False,
+    no_pbar=False,
+    verbose=False,
+    logfile=None,
+    **kwargs,
 ):
-    try:
-        input_info = _get_input_info(input_)
-        output_info = _get_output_info(output)
-    except Exception as e:
-        raise click.BadArgumentUsage(e)
-
-    # collect mapchete configuration
-    mapchete_config = dict(
-        process="mapchete.processes.convert",
-        input=dict(inp=input_, clip=clip_geometry),
-        pyramid=(
-            dict(
-                grid=output_pyramid,
-                metatiling=(
-                    output_metatiling
-                    or (
-                        input_info["pyramid"].get("metatiling", 1)
-                        if input_info["pyramid"]
-                        else 1
-                    )
-                ),
-                pixelbuffer=(
-                    input_info["pyramid"].get("pixelbuffer", 0)
-                    if input_info["pyramid"]
-                    else 0
-                ),
-            )
-            if output_pyramid
-            else input_info["pyramid"]
-        ),
-        output=dict(
-            {
-                k: v
-                for k, v in input_info["output_params"].items()
-                if k not in ["delimiters", "bounds", "mode"]
-            },
-            path=output,
-            format=(
-                output_format
-                or output_info["driver"]
-                or input_info["output_params"]["format"]
-            ),
-            dtype=output_dtype or input_info["output_params"].get("dtype"),
-            **creation_options,
-            **dict(overviews=True, overviews_resampling=overviews_resampling_method)
-            if overviews
-            else dict(),
-        ),
-        config_dir=os.getcwd(),
-        zoom_levels=zoom or input_info["zoom_levels"],
-        scale_ratio=scale_ratio,
-        scale_offset=scale_offset,
-        resampling=resampling_method,
-        band_indexes=bidx,
-    )
-
-    # assert all required information is there
-    if mapchete_config["output"]["format"] is None:
-        # this happens if input file is e.g. JPEG2000 and output is a tile directory
-        raise click.BadOptionUsage("output-format", "Output format required.")
-    if mapchete_config["output"]["format"] == "GTiff":
-        mapchete_config["output"].update(cog=cog)
-    output_type = OUTPUT_FORMATS[mapchete_config["output"]["format"]]["data_type"]
-    if bidx is not None:
-        mapchete_config["output"].update(bands=len(bidx))
-    if mapchete_config["pyramid"] is None:
-        raise click.BadOptionUsage("output-pyramid", "Output pyramid required.")
-    elif mapchete_config["zoom_levels"] is None:
-        try:
-            mapchete_config.update(
-                zoom_levels=dict(
-                    min=0,
-                    max=get_best_zoom_level(input_, mapchete_config["pyramid"]["grid"]),
-                )
-            )
-        except Exception:
-            raise click.BadOptionUsage("zoom", "Zoom levels required.")
-    elif input_info["input_type"] != output_type:
-        raise click.BadArgumentUsage(
-            "Output format type (%s) is incompatible with input format (%s)."
-            % (output_type, input_info["input_type"])
+    with mapchete.Timer() as t:
+        job = commands.convert(
+            tiledir,
+            output,
+            *args,
+            as_iterator=True,
+            msg_callback=tqdm.tqdm.write if verbose else None,
+            **kwargs,
         )
-    if output_metatiling:
-        mapchete_config["output"].update(metatiling=output_metatiling)
-    if input_info["output_params"].get("schema") and output_geometry_type:
-        mapchete_config["output"]["schema"].update(geometry=output_geometry_type)
-
-    # determine process bounds
-    out_pyramid = BufferedTilePyramid.from_dict(mapchete_config["pyramid"])
-    inp_bounds = (
-        bounds
-        or reproject_geometry(
-            box(*input_info["bounds"]),
-            src_crs=input_info["crs"],
-            dst_crs=out_pyramid.crs,
-        ).bounds
-        if input_info["bounds"]
-        else out_pyramid.bounds
-    )
-    # if clip-geometry is available, intersect determined bounds with clip bounds
-    if clip_geometry:
-        clip_intersection = _clip_bbox(
-            clip_geometry, dst_crs=out_pyramid.crs
-        ).intersection(box(*inp_bounds))
-        if clip_intersection.is_empty:
-            click.echo(
-                "Process area is empty: clip bounds don't intersect with input bounds."
-            )
+        if not len(job):
             return
-    # add process bounds and output type
-    mapchete_config.update(
-        bounds=(clip_intersection.bounds if clip_geometry else inp_bounds),
-        bounds_crs=bounds_crs,
-        clip_to_output_dtype=mapchete_config["output"].get("dtype", None),
-    )
-    logger.debug("temporary config generated: %s", pformat(mapchete_config))
-
-    utils._process_area(
-        debug=debug,
-        mapchete_config=mapchete_config,
-        mode="overwrite" if overwrite else "continue",
-        zoom=zoom,
-        wkt_geometry=wkt_geometry,
-        point=point,
-        point_crs=point_crs,
-        bounds=bounds,
-        bounds_crs=bounds_crs,
-        area=area,
-        area_crs=area_crs,
-        multi=multi or cpu_count(),
-        verbose_dst=open(os.devnull, "w") if debug or not verbose else sys.stdout,
-        no_pbar=no_pbar,
-        vrt=vrt,
-        idx_out_dir=idx_out_dir,
-    )
-
-
-def _clip_bbox(clip_geometry, dst_crs=None):
-    with fiona.open(clip_geometry) as src:
-        return reproject_geometry(box(*src.bounds), src_crs=src.crs, dst_crs=dst_crs)
-
-
-def _get_input_info(input_):
-
-    # assuming single file if path has a file extension
-    if os.path.splitext(input_)[1]:
-        driver = driver_from_file(input_)
-
-        # single file input can be a mapchete file or a rasterio/fiona file
-        if driver == "Mapchete":
-            logger.debug("input is mapchete file")
-            input_info = _input_mapchete_info(input_)
-
-        elif driver == "raster_file":
-            # this should be readable by rasterio
-            logger.debug("input is raster_file")
-            input_info = _input_rasterio_info(input_)
-
-        elif driver == "vector_file":
-            # this should be readable by Fiona
-            input_info = _input_fiona_info(input_)
-        else:  # pragma: no cover
-            raise NotImplementedError("driver %s is not supported" % driver)
-
-    # assuming tile directory
-    else:
-        logger.debug("input is tile directory")
-        input_info = _input_tile_directory_info(input_)
-
-    return input_info
-
-
-def _input_mapchete_info(input_):
-    conf = raw_conf(input_)
-    output_params = conf["output"]
-    pyramid = raw_conf_output_pyramid(conf)
-    return dict(
-        output_params=output_params,
-        pyramid=pyramid.to_dict(),
-        crs=pyramid.crs,
-        zoom_levels=validate_zooms(conf["zoom_levels"], expand=False),
-        pixel_size=None,
-        input_type=OUTPUT_FORMATS[output_params["format"]]["data_type"],
-        bounds=conf.get("bounds"),
-    )
-
-
-def _input_rasterio_info(input_):
-    with rasterio.open(input_) as src:
-        return dict(
-            output_params=dict(
-                bands=src.meta["count"],
-                dtype=src.meta["dtype"],
-                format=src.driver if src.driver in available_input_formats() else None,
-            ),
-            pyramid=None,
-            crs=src.crs,
-            zoom_levels=None,
-            pixel_size=src.transform[0],
-            input_type="raster",
-            bounds=src.bounds,
+        list(
+            tqdm.tqdm(
+                job,
+                unit="task",
+                disable=debug or no_pbar,
+            )
         )
-
-
-def _input_fiona_info(input_):
-    with fiona.open(input_) as src:
-        return dict(
-            output_params=dict(
-                schema=src.schema,
-                format=src.driver if src.driver in available_input_formats() else None,
-            ),
-            pyramid=None,
-            crs=src.crs,
-            zoom_levels=None,
-            input_type="vector",
-            bounds=src.bounds,
+        tqdm.tqdm.write(f"processing {tiledir} finished in {t}")
+
+    if vrt:
+        tqdm.tqdm.write("creating VRT(s)")
+        list(
+            tqdm.tqdm(
+                commands.index(
+                    output,
+                    *args,
+                    vrt=vrt,
+                    idx_out_dir=idx_out_dir,
+                    as_iterator=True,
+                    msg_callback=tqdm.tqdm.write if verbose else None,
+                    **kwargs,
+                ),
+                unit="tile",
+                disable=debug or no_pbar,
+            )
         )
-
-
-def _input_tile_directory_info(input_):
-    conf = read_json(os.path.join(input_, "metadata.json"))
-    pyramid = BufferedTilePyramid.from_dict(conf["pyramid"])
-    return dict(
-        output_params=conf["driver"],
-        pyramid=pyramid.to_dict(),
-        crs=pyramid.crs,
-        zoom_levels=None,
-        pixel_size=None,
-        input_type=OUTPUT_FORMATS[conf["driver"]["format"]]["data_type"],
-        bounds=None,
-    )
-
-
-def _get_output_info(output):
-    _, file_ext = os.path.splitext(output)
-    if not file_ext:
-        return dict(type="TileDirectory", driver=None)
-    elif file_ext == ".tif":
-        return dict(type="SingleFile", driver="GTiff")
-    else:
-        raise TypeError("Could not determine output from extension: %s", file_ext)
+        tqdm.tqdm.write(f"index(es) creation for {tiledir} finished")
diff --git a/mapchete/cli/default/cp.py b/mapchete/cli/default/cp.py
index 98539923..e532ea61 100644
--- a/mapchete/cli/default/cp.py
+++ b/mapchete/cli/default/cp.py
@@ -1,161 +1,50 @@
 import click
-import fiona
-import json
-import logging
-import os
-from shapely.geometry import box, shape
-from shapely.ops import unary_union
-from tilematrix import TilePyramid
 import tqdm
 
-import mapchete
-from mapchete.cli import utils
-from mapchete.config import _guess_geometry
-from mapchete.formats import read_output_metadata
-from mapchete.io import fs_from_path, tiles_exist
-from mapchete.io.vector import reproject_geometry
-
-logger = logging.getLogger(__name__)
+from mapchete import commands
+from mapchete.cli import options
 
 
 @click.command(help="Copy TileDirectory from one source to another.")
-@utils.arg_input
-@utils.arg_output
-@utils.opt_zoom
-@utils.opt_area
-@utils.opt_area_crs
-@utils.opt_bounds
-@utils.opt_bounds_crs
-@utils.opt_wkt_geometry
-@utils.opt_overwrite
-@utils.opt_verbose
-@utils.opt_no_pbar
-@utils.opt_debug
-@utils.opt_logfile
-@utils.opt_multi
-@utils.opt_http_username
-@utils.opt_http_password
-@utils.opt_src_fs_opts
-@utils.opt_dst_fs_opts
-def cp(
-    input_,
-    output,
-    zoom=None,
-    area=None,
-    area_crs=None,
-    bounds=None,
-    bounds_crs=None,
-    wkt_geometry=None,
-    overwrite=False,
-    verbose=False,
-    no_pbar=False,
-    debug=False,
-    logfile=None,
-    multi=None,
-    username=None,
-    password=None,
-    src_fs_opts=None,
-    dst_fs_opts=None,
-):
+@options.arg_src_tiledir
+@options.arg_dst_tiledir
+@options.opt_zoom
+@options.opt_area
+@options.opt_area_crs
+@options.opt_bounds
+@options.opt_bounds_crs
+@options.opt_overwrite
+@options.opt_verbose
+@options.opt_no_pbar
+@options.opt_debug
+@options.opt_logfile
+@options.opt_workers
+@options.opt_dask_scheduler
+@options.opt_concurrency
+@options.opt_http_username
+@options.opt_http_password
+@options.opt_src_fs_opts
+@options.opt_dst_fs_opts
+def cp(*args, debug=False, no_pbar=False, verbose=False, logfile=None, **kwargs):
     """Copy TileDirectory."""
-    if zoom is None:  # pragma: no cover
-        raise click.UsageError("zoom level(s) required")
-
-    src_fs = fs_from_path(input_, username=username, password=password, **src_fs_opts)
-    dst_fs = fs_from_path(output, username=username, password=password, **dst_fs_opts)
-
-    # open source tile directory
-    with mapchete.open(
-        input_,
-        zoom=zoom,
-        area=area,
-        area_crs=area_crs,
-        bounds=bounds,
-        bounds_crs=bounds_crs,
-        wkt_geometry=wkt_geometry,
-        username=username,
-        password=password,
-        fs=src_fs,
-        fs_kwargs=src_fs_opts,
-    ) as src_mp:
-        tp = src_mp.config.output_pyramid
-
-        # copy metadata to destination if necessary
-        src_metadata = os.path.join(input_, "metadata.json")
-        dst_metadata = os.path.join(output, "metadata.json")
-        if not dst_fs.exists(dst_metadata):
-            logger.debug(f"copy {src_metadata} to {dst_metadata}")
-            _copy(src_fs, src_metadata, dst_fs, dst_metadata)
-
-        with mapchete.open(
-            output,
-            zoom=zoom,
-            area=area,
-            area_crs=area_crs,
-            bounds=bounds,
-            bounds_crs=bounds_crs,
-            wkt_geometry=wkt_geometry,
-            username=username,
-            password=password,
-            fs=dst_fs,
-            fs_kwargs=dst_fs_opts,
-        ) as dst_mp:
-            for z in range(min(zoom), max(zoom) + 1):
-                click.echo(f"copy zoom {z}...")
-                # materialize all tiles
-                aoi_geom = src_mp.config.area_at_zoom(z)
-                tiles = [
-                    t
-                    for t in tp.tiles_from_geom(aoi_geom, z)
-                    # this is required to omit tiles touching the config area
-                    if aoi_geom.intersection(t.bbox).area
-                ]
-
-                # check which source tiles exist
-                logger.debug("looking for existing source tiles...")
-                src_tiles_exist = {
-                    tile: exists
-                    for tile, exists in tiles_exist(
-                        config=src_mp.config, output_tiles=tiles, multi=multi
-                    )
-                }
-
-                logger.debug("looking for existing destination tiles...")
-                # chech which destination tiles exist
-                dst_tiles_exist = {
-                    tile: exists
-                    for tile, exists in tiles_exist(
-                        config=dst_mp.config, output_tiles=tiles, multi=multi
-                    )
-                }
-
-                # copy
-                for tile in tqdm.tqdm(tiles, unit="tile", disable=debug or no_pbar):
-                    src_path = src_mp.config.output_reader.get_path(tile)
-                    # only copy if source tile exists
-                    if src_tiles_exist[tile]:
-                        # skip if destination tile exists and overwrite is deactivated
-                        if dst_tiles_exist[tile] and not overwrite:
-                            logger.debug(f"{tile}: destination tile exists")
-                            continue
-                        # copy from source to target
-                        else:
-                            dst_path = dst_mp.config.output_reader.get_path(tile)
-                            logger.debug(f"{tile}: copy {src_path} to {dst_path}")
-                            _copy(src_fs, src_path, dst_fs, dst_path)
-                    else:
-                        logger.debug(f"{tile}: source tile ({src_path}) does not exist")
-
-
-def _copy(src_fs, src_path, dst_fs, dst_path):
-    # create parent directories on local filesystems
-    if dst_fs.protocol == "file":
-        dst_fs.mkdir(os.path.dirname(dst_path), create_parents=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())
+    # handle deprecated options
+    for x in ["password", "username"]:
+        if kwargs.get(x):  # pragma: no cover
+            raise click.BadOptionUsage(
+                x,
+                f"'--{x} foo' is deprecated. You should use '--src-fs-opts {x}=foo' instead.",
+            )
+        kwargs.pop(x)
+    # copy
+    list(
+        tqdm.tqdm(
+            commands.cp(
+                *args,
+                as_iterator=True,
+                msg_callback=tqdm.tqdm.write if verbose else None,
+                **kwargs,
+            ),
+            unit="tile",
+            disable=debug or no_pbar,
+        )
+    )
diff --git a/mapchete/cli/default/create.py b/mapchete/cli/default/create.py
index 1486dc47..dbc7e099 100644
--- a/mapchete/cli/default/create.py
+++ b/mapchete/cli/default/create.py
@@ -7,7 +7,7 @@
 from shutil import copyfile
 from oyaml import dump
 
-from mapchete.cli import utils
+from mapchete.cli import options
 
 FORMAT_MANDATORY = {
     "GTiff": {"bands": None, "dtype": None},
@@ -29,12 +29,12 @@
 
 
 @click.command(help="Create a new process.")
-@utils.arg_create_mapchete_file
-@utils.arg_process_file
-@utils.arg_out_format
-@utils.opt_out_path
-@utils.opt_pyramid_type
-@utils.opt_force
+@options.arg_create_mapchete_file
+@options.arg_process_file
+@options.arg_out_format
+@options.opt_out_path
+@options.opt_pyramid_type
+@options.opt_force
 def create(
     mapchete_file,
     process_file,
diff --git a/mapchete/cli/default/execute.py b/mapchete/cli/default/execute.py
old mode 100755
new mode 100644
index a863c664..069d71fb
--- a/mapchete/cli/default/execute.py
+++ b/mapchete/cli/default/execute.py
@@ -1,108 +1,83 @@
-"""Command line utility to execute a Mapchete process."""
-
 import click
-import logging
-from multiprocessing import cpu_count
-import os
-import sys
 import tqdm
 
-from mapchete.cli import utils
-from mapchete.config import raw_conf_process_pyramid
-
-
-# workaround for https://github.com/tqdm/tqdm/issues/481
-tqdm.monitor_interval = 0
-
-logger = logging.getLogger(__name__)
+import mapchete
+from mapchete import commands
+from mapchete.cli import options
 
 
 @click.command(help="Execute a process.")
-@utils.arg_mapchete_files
-@utils.opt_zoom
-@utils.opt_bounds
-@utils.opt_bounds_crs
-@utils.opt_area
-@utils.opt_area_crs
-@utils.opt_point
-@utils.opt_point_crs
-@utils.opt_wkt_geometry
-@utils.opt_tile
-@utils.opt_overwrite
-@utils.opt_multi
-@utils.opt_input_file
-@utils.opt_logfile
-@utils.opt_verbose
-@utils.opt_no_pbar
-@utils.opt_debug
-@utils.opt_max_chunksize
-@utils.opt_multiprocessing_start_method
-@utils.opt_vrt
-@utils.opt_idx_out_dir
+@options.arg_mapchete_files
+@options.opt_zoom
+@options.opt_bounds
+@options.opt_bounds_crs
+@options.opt_area
+@options.opt_area_crs
+@options.opt_point
+@options.opt_point_crs
+@options.opt_tile
+@options.opt_overwrite
+@options.opt_concurrency
+@options.opt_workers
+@options.opt_multi
+@options.opt_dask_scheduler
+@options.opt_input_file
+@options.opt_logfile
+@options.opt_verbose
+@options.opt_no_pbar
+@options.opt_debug
+@options.opt_vrt
+@options.opt_idx_out_dir
 def execute(
     mapchete_files,
-    zoom=None,
-    bounds=None,
-    bounds_crs=None,
-    area=None,
-    area_crs=None,
-    point=None,
-    point_crs=None,
-    wkt_geometry=None,
-    tile=None,
-    overwrite=False,
-    multi=None,
-    input_file=None,
-    logfile=None,
-    verbose=False,
-    no_pbar=False,
-    debug=False,
-    max_chunksize=None,
-    multiprocessing_start_method=None,
+    *args,
     vrt=False,
     idx_out_dir=None,
+    debug=False,
+    no_pbar=False,
+    verbose=False,
+    logfile=None,
+    input_file=None,
+    **kwargs,
 ):
-    """Execute a Mapchete process."""
-    mode = "overwrite" if overwrite else "continue"
-    # send verbose messages to /dev/null if not activated
-    verbose_dst = open(os.devnull, "w") if debug or not verbose else sys.stdout
-
+    if input_file is not None:  # pragma: no cover
+        raise click.BadOptionUsage(
+            "input-file",
+            "'--input-file' is deprecated.",
+        )
     for mapchete_file in mapchete_files:
-        tqdm.tqdm.write("preparing to process %s" % mapchete_file, file=verbose_dst)
-        # process single tile
-        if tile:
-            utils._process_single_tile(
-                raw_conf_process_pyramid=raw_conf_process_pyramid,
-                mapchete_config=mapchete_file,
-                tile=tile,
-                mode=mode,
-                input_file=input_file,
-                debug=debug,
-                verbose_dst=verbose_dst,
-                vrt=vrt,
-                idx_out_dir=idx_out_dir,
-                no_pbar=no_pbar,
+        tqdm.tqdm.write(f"preparing to process {mapchete_file}")
+        with mapchete.Timer() as t:
+            list(
+                tqdm.tqdm(
+                    commands.execute(
+                        mapchete_file,
+                        *args,
+                        as_iterator=True,
+                        msg_callback=tqdm.tqdm.write if verbose else None,
+                        **kwargs,
+                    ),
+                    unit="task",
+                    disable=debug or no_pbar,
+                )
             )
-        # process area
-        else:
-            utils._process_area(
-                debug=debug,
-                mapchete_config=mapchete_file,
-                mode=mode,
-                zoom=zoom,
-                wkt_geometry=wkt_geometry,
-                point=point,
-                point_crs=point_crs,
-                bounds=bounds,
-                bounds_crs=bounds_crs,
-                area=area,
-                area_crs=area_crs,
-                input_file=input_file,
-                multi=multi or cpu_count(),
-                verbose_dst=verbose_dst,
-                max_chunksize=max_chunksize,
-                multiprocessing_start_method=multiprocessing_start_method,
-                no_pbar=no_pbar,
-                vrt=vrt,
-                idx_out_dir=idx_out_dir,
+            tqdm.tqdm.write(f"processing {mapchete_file} finished in {t}")
+
+        if vrt:
+            tqdm.tqdm.write("creating VRT(s)")
+            list(
+                tqdm.tqdm(
+                    commands.index(
+                        mapchete_file,
+                        *args,
+                        vrt=vrt,
+                        idx_out_dir=idx_out_dir,
+                        as_iterator=True,
+                        msg_callback=tqdm.tqdm.write if verbose else None,
+                        **kwargs,
+                    ),
+                    unit="tile",
+                    disable=debug or no_pbar,
+                )
             )
+            tqdm.tqdm.write(f"index(es) creation for {mapchete_file} finished")
diff --git a/mapchete/cli/default/formats.py b/mapchete/cli/default/formats.py
index 85682eb1..059e377a 100644
--- a/mapchete/cli/default/formats.py
+++ b/mapchete/cli/default/formats.py
@@ -2,14 +2,14 @@
 
 import click
 
-from mapchete.cli import utils
+from mapchete.cli import options
 from mapchete.formats import available_input_formats, available_output_formats
 
 
 @click.command(help="List available input and/or output formats.")
-@utils.opt_input_formats
-@utils.opt_output_formats
-@utils.opt_debug
+@options.opt_input_formats
+@options.opt_output_formats
+@options.opt_debug
 def formats(input_formats, output_formats, debug=False):
     """List input and/or output formats."""
     if input_formats == output_formats:
diff --git a/mapchete/cli/default/index.py b/mapchete/cli/default/index.py
index 1babab07..e9ca6e3a 100644
--- a/mapchete/cli/default/index.py
+++ b/mapchete/cli/default/index.py
@@ -8,8 +8,8 @@
 import tqdm
 
 import mapchete
-from mapchete.cli import utils
-from mapchete.index import zoom_index_gen
+from mapchete.cli import options
+from mapchete import commands
 
 
 # workaround for https://github.com/tqdm/tqdm/issues/481
@@ -19,146 +19,51 @@
 
 
 @click.command(help="Create index of output tiles.")
-@utils.arg_inputs
-@utils.opt_idx_out_dir
-@utils.opt_geojson
-@utils.opt_gpkg
-@utils.opt_shp
-@utils.opt_vrt
-@utils.opt_txt
-@utils.opt_fieldname
-@utils.opt_basepath
-@utils.opt_for_gdal
-@utils.opt_zoom
-@utils.opt_bounds
-@utils.opt_bounds_crs
-@utils.opt_area
-@utils.opt_area_crs
-@utils.opt_point
-@utils.opt_point_crs
-@utils.opt_wkt_geometry
-@utils.opt_tile
-@utils.opt_http_username
-@utils.opt_http_password
-@utils.opt_verbose
-@utils.opt_no_pbar
-@utils.opt_debug
-@utils.opt_logfile
-def index(
-    inputs,
-    idx_out_dir=None,
-    geojson=False,
-    gpkg=False,
-    shp=False,
-    vrt=False,
-    txt=False,
-    fieldname=None,
-    basepath=None,
-    for_gdal=False,
-    zoom=None,
-    bounds=None,
-    bounds_crs=None,
-    area=None,
-    area_crs=None,
-    point=None,
-    point_crs=None,
-    wkt_geometry=None,
-    tile=None,
-    username=None,
-    password=None,
-    verbose=False,
-    no_pbar=False,
-    debug=False,
-    logfile=None,
-):
-    if not any([geojson, gpkg, shp, txt, vrt]):
-        raise click.MissingParameter(
-            """At least one of '--geojson', '--gpkg', '--shp', '--vrt' or '--txt'"""
-            """must be provided.""",
-            param_type="option",
+@options.arg_tiledir
+@options.opt_idx_out_dir
+@options.opt_geojson
+@options.opt_gpkg
+@options.opt_shp
+@options.opt_vrt
+@options.opt_txt
+@options.opt_fieldname
+@options.opt_basepath
+@options.opt_for_gdal
+@options.opt_zoom
+@options.opt_bounds
+@options.opt_bounds_crs
+@options.opt_area
+@options.opt_area_crs
+@options.opt_point
+@options.opt_point_crs
+@options.opt_tile
+@options.opt_http_username
+@options.opt_http_password
+@options.opt_fs_opts
+@options.opt_verbose
+@options.opt_no_pbar
+@options.opt_debug
+@options.opt_logfile
+def index(*args, debug=False, no_pbar=False, verbose=False, logfile=None, **kwargs):
+    """Create various index files from process output."""
+    # handle deprecated options
+    for x in ["password", "username"]:
+        if kwargs.get(x):  # pragma: no cover
+            raise click.BadOptionUsage(
+                x,
+                f"'--{x} foo' is deprecated. You should use '--fs-opts {x}=foo' instead.",
+            )
+        kwargs.pop(x)
+    list(
+        tqdm.tqdm(
+            commands.index(
+                *args,
+                as_iterator=True,
+                msg_callback=tqdm.tqdm.write if verbose else None,
+                **kwargs,
+            ),
+            unit="tile",
+            disable=debug or no_pbar,
         )
-    # send verbose messages to /dev/null if not activated
-    verbose_dst = open(os.devnull, "w") if debug or not verbose else sys.stdout
-
-    for input_ in inputs:
-
-        tqdm.tqdm.write("create index(es) for %s" % input_, file=verbose_dst)
-
-        with click_spinner.spinner(disable=debug) as spinner:
-
-            # process single tile
-            if tile:
-                with mapchete.open(
-                    input_, mode="readonly", username=username, password=password
-                ) as mp:
-                    spinner.stop()
-                    tile = mp.config.process_pyramid.tile(*tile)
-                    for tile in tqdm.tqdm(
-                        zoom_index_gen(
-                            mp=mp,
-                            zoom=tile.zoom,
-                            out_dir=idx_out_dir
-                            if idx_out_dir
-                            else mp.config.output.path,
-                            geojson=geojson,
-                            gpkg=gpkg,
-                            shapefile=shp,
-                            vrt=vrt,
-                            txt=txt,
-                            fieldname=fieldname,
-                            basepath=basepath,
-                            for_gdal=for_gdal,
-                        ),
-                        total=mp.count_tiles(tile.zoom, tile.zoom),
-                        unit="tile",
-                        disable=debug or no_pbar,
-                    ):
-                        logger.debug("%s indexed", tile)
-
-            # process area
-            else:
-                with mapchete.open(
-                    input_,
-                    mode="readonly",
-                    zoom=zoom,
-                    wkt_geometry=wkt_geometry,
-                    point=point,
-                    point_crs=point_crs,
-                    bounds=bounds,
-                    bounds_crs=bounds_crs,
-                    area=area,
-                    area_crs=area_crs,
-                    username=username,
-                    password=password,
-                ) as mp:
-                    spinner.stop()
-                    logger.debug("process bounds: %s", mp.config.init_bounds)
-                    logger.debug("process area: %s", mp.config.init_area)
-                    logger.debug("process zooms: %s", mp.config.init_zoom_levels)
-                    logger.debug("fieldname: %s", fieldname)
-                    for tile in tqdm.tqdm(
-                        zoom_index_gen(
-                            mp=mp,
-                            zoom=mp.config.init_zoom_levels,
-                            out_dir=(
-                                idx_out_dir if idx_out_dir else mp.config.output.path
-                            ),
-                            geojson=geojson,
-                            gpkg=gpkg,
-                            shapefile=shp,
-                            vrt=vrt,
-                            txt=txt,
-                            fieldname=fieldname,
-                            basepath=basepath,
-                            for_gdal=for_gdal,
-                        ),
-                        total=mp.count_tiles(
-                            min(mp.config.init_zoom_levels),
-                            max(mp.config.init_zoom_levels),
-                        ),
-                        unit="tile",
-                        disable=debug or no_pbar,
-                    ):
-                        logger.debug("%s indexed", tile)
-
-        tqdm.tqdm.write("index(es) creation for %s finished" % input_, file=verbose_dst)
+    )
+    tqdm.tqdm.write(f"index(es) creation for {kwargs.get('tiledir')} finished")
diff --git a/mapchete/cli/default/processes.py b/mapchete/cli/default/processes.py
index c84bf6b2..17d95427 100644
--- a/mapchete/cli/default/processes.py
+++ b/mapchete/cli/default/processes.py
@@ -3,7 +3,7 @@
 import click
 import logging
 
-from mapchete.cli import utils
+from mapchete.cli import options
 from mapchete.processes import process_names_docstrings
 
 logger = logging.getLogger(__name__)
@@ -13,7 +13,7 @@
 @click.option(
     "--process_name", "-n", type=click.STRING, help="Print docstring of process."
 )
-@utils.opt_debug
+@options.opt_debug
 def processes(process_name=None, docstrings=False, debug=False):
     """List available processes."""
     processes = process_names_docstrings(process_name=process_name)
diff --git a/mapchete/cli/default/rm.py b/mapchete/cli/default/rm.py
index 2a42d94b..7571c84d 100644
--- a/mapchete/cli/default/rm.py
+++ b/mapchete/cli/default/rm.py
@@ -1,98 +1,50 @@
 import click
-import logging
 import tqdm
 
-import mapchete
-from mapchete.cli import utils
-from mapchete.io import fs_from_path, rm, tiles_exist
-
-logger = logging.getLogger(__name__)
-
-
-@click.command("rm", help="Remove tiles from TileDirectory.")
-@utils.arg_input
-@utils.opt_zoom
-@utils.opt_area
-@utils.opt_area_crs
-@utils.opt_bounds
-@utils.opt_bounds_crs
-@utils.opt_wkt_geometry
-@utils.opt_multi
-@utils.opt_verbose
-@utils.opt_no_pbar
-@utils.opt_debug
-@utils.opt_logfile
-@utils.opt_force
-@utils.opt_fs_opts
-def rm_(
-    input_,
-    zoom=None,
-    area=None,
-    area_crs=None,
-    bounds=None,
-    bounds_crs=None,
-    wkt_geometry=None,
-    multi=None,
-    verbose=False,
-    no_pbar=False,
+from mapchete import commands
+from mapchete.cli import options
+
+
+@click.command(help="Remove tiles from TileDirectory.")
+@options.arg_tiledir
+@options.opt_zoom
+@options.opt_area
+@options.opt_area_crs
+@options.opt_bounds
+@options.opt_bounds_crs
+@options.opt_multi
+@options.opt_verbose
+@options.opt_no_pbar
+@options.opt_debug
+@options.opt_logfile
+@options.opt_force
+@options.opt_fs_opts
+def rm(
+    *args,
+    force=False,
     debug=False,
+    no_pbar=False,
+    verbose=False,
     logfile=None,
-    force=None,
-    fs_opts=None,
+    **kwargs,
 ):
-    """Copy TileDirectory."""
-    if zoom is None:  # pragma: no cover
-        raise click.UsageError("zoom level(s) required")
-
-    src_fs = fs_from_path(input_, **fs_opts)
-
-    # open source tile directory
-    with mapchete.open(
-        input_,
-        zoom=zoom,
-        area=area,
-        area_crs=area_crs,
-        bounds=bounds,
-        bounds_crs=bounds_crs,
-        wkt_geometry=wkt_geometry,
-        fs=src_fs,
-        fs_kwargs=fs_opts,
-        mode="readonly",
-    ) as src_mp:
-        tp = src_mp.config.output_pyramid
-
-        tiles = {}
-        for z in range(min(zoom), max(zoom) + 1):
-            tiles[z] = []
-            # check which source tiles exist
-            logger.debug(f"looking for existing source tiles in zoom {z}...")
-            for tile, exists in tiles_exist(
-                config=src_mp.config,
-                output_tiles=[
-                    t
-                    for t in tp.tiles_from_geom(src_mp.config.area_at_zoom(z), z)
-                    # this is required to omit tiles touching the config area
-                    if src_mp.config.area_at_zoom(z).intersection(t.bbox).area
-                ],
-                multi=multi,
-            ):
-                if exists:
-                    tiles[z].append(tile)
-
-        total_tiles = sum([len(v) for v in tiles.values()])
-
-        if total_tiles:
-            if force or click.confirm(
-                f"Do you want to delete {total_tiles} tiles?", abort=True
-            ):
-                # remove
-                rm(
-                    [
-                        src_mp.config.output_reader.get_path(tile)
-                        for zoom_tiles in tiles.values()
-                        for tile in zoom_tiles
-                    ],
-                    fs=src_fs,
+    """Remove tiles from TileDirectory."""
+    tiles_to_delete = commands.rm(
+        *args,
+        as_iterator=True,
+        msg_callback=tqdm.tqdm.write if verbose else None,
+        **kwargs,
+    )
+    if len(tiles_to_delete):
+        if force or click.confirm(
+            f"Do you want to delete {len(tiles_to_delete)} tiles?", abort=True
+        ):
+            list(
+                tqdm.tqdm(
+                    tiles_to_delete,
+                    unit="tile",
+                    disable=debug or no_pbar,
                 )
-        else:  # pragma: no cover
-            click.echo("No tiles found to delete.")
+            )
+    else:  # pragma: no cover
+        tqdm.tqdm.write("No tiles found to delete.")
diff --git a/mapchete/cli/default/serve.py b/mapchete/cli/default/serve.py
index 1ac9e438..c0369b75 100755
--- a/mapchete/cli/default/serve.py
+++ b/mapchete/cli/default/serve.py
@@ -9,24 +9,24 @@
 from rasterio.io import MemoryFile
 
 import mapchete
-from mapchete.cli import utils
+from mapchete.cli import options
 from mapchete.tile import BufferedTilePyramid
 
 logger = logging.getLogger(__name__)
 
 
 @click.command(help="Serve a process on localhost.")
-@utils.arg_mapchete_files
-@utils.opt_port
-@utils.opt_internal_cache
-@utils.opt_zoom
-@utils.opt_bounds
-@utils.opt_overwrite
-@utils.opt_readonly
-@utils.opt_memory
-@utils.opt_input_file
-@utils.opt_debug
-@utils.opt_logfile
+@options.arg_mapchete_files
+@options.opt_port
+@options.opt_internal_cache
+@options.opt_zoom
+@options.opt_bounds
+@options.opt_overwrite
+@options.opt_readonly
+@options.opt_memory
+@options.opt_input_file
+@options.opt_debug
+@options.opt_logfile
 def serve(
     mapchete_files,
     port=None,
diff --git a/mapchete/cli/default/stac.py b/mapchete/cli/default/stac.py
new file mode 100644
index 00000000..4dd3d159
--- /dev/null
+++ b/mapchete/cli/default/stac.py
@@ -0,0 +1,129 @@
+import click
+import fsspec
+import logging
+import json
+import os
+import oyaml as yaml
+
+from mapchete.cli import options
+from mapchete.config import raw_conf, raw_conf_output_pyramid
+from mapchete.formats import read_output_metadata
+from mapchete.stac import tile_directory_stac_item
+from mapchete.validate import validate_zooms
+
+
+logger = logging.getLogger(__name__)
+
+
+@click.group()
+def stac():
+    pass
+
+
+@stac.command(help="Create STAC item metadata.")
+@options.arg_input
+@click.option("--item-id", "-i", type=click.STRING, help="Unique item ID.")
+@click.option(
+    "--item-metadata",
+    "-m",
+    type=click.Path(),
+    help="Optional additional item metadata to be appended. Must be a YAML file.",
+)
+@options.opt_zoom
+@click.option("--item-path", "-p", type=click.Path(), help="Path of output STAC item.")
+@click.option("--asset-basepath", type=click.Path(), help="Alternative asset basepath.")
+@click.option(
+    "--indent",
+    type=click.INT,
+    default=4,
+    help="Indentation for output JSON. (default: 4)",
+)
+@options.opt_bounds
+@options.opt_bounds_crs
+@options.opt_force
+@options.opt_debug
+def create_item(
+    input_,
+    item_id=None,
+    item_metadata=None,
+    asset_basepath=None,
+    zoom=None,
+    bounds=None,
+    bounds_crs=None,
+    item_path=None,
+    indent=None,
+    force=None,
+    **kwargs,
+):
+    (
+        tile_pyramid,
+        default_basepath,
+        default_id,
+        default_bounds,
+        default_bounds_crs,
+        default_zoom,
+        default_item_metadata,
+    ) = output_info(input_)
+
+    if default_zoom:
+        zoom = zoom or validate_zooms(default_zoom)
+
+    if zoom is None:  # pragma: no cover
+        raise ValueError("zoom must be set")
+    min_zoom, max_zoom = min(zoom), max(zoom)
+
+    if item_metadata:  # pragma: no cover
+        with fsspec.open(item_metadata) as src:
+            metadata = yaml.safe_load(src.read())
+    else:
+        metadata = default_item_metadata or {}
+
+    item_id = item_id or metadata.get("id", default_id)
+    logger.debug("use item ID %s", item_id)
+    item_path = item_path or os.path.join(default_basepath, f"{item_id}.json")
+    item = tile_directory_stac_item(
+        item_id=item_id,
+        item_metadata=metadata,
+        tile_pyramid=tile_pyramid,
+        min_zoom=min_zoom,
+        max_zoom=max_zoom,
+        bounds=bounds or default_bounds,
+        bounds_crs=bounds_crs or default_bounds_crs,
+        item_path=item_path,
+        asset_basepath=asset_basepath,
+        relative_paths=None,
+        bands_type=None,
+        crs_unit_to_meter=1,
+    )
+    logger.debug("item_path: %s", item_path)
+    item_json = json.dumps(item.to_dict(), indent=indent)
+    click.echo(item_json)
+    if force or click.confirm(f"Write output to {item_path}?", abort=True):
+        with fsspec.open(item_path, "w") as dst:
+            dst.write(item_json)
+
+
+def output_info(inp):
+    if inp.endswith(".mapchete"):
+        conf = raw_conf(inp)
+        default_basepath = os.path.dirname(conf["output"]["path"].strip("/") + "/")
+        return (
+            raw_conf_output_pyramid(conf),
+            default_basepath,
+            os.path.basename(default_basepath),
+            conf.get("bounds"),
+            conf.get("bounds_crs"),
+            conf.get("zoom_levels"),
+            conf["output"].get("stac"),
+        )
+
+    default_basepath = inp.strip("/")
+    return (
+        read_output_metadata(os.path.join(inp, "metadata.json"))["pyramid"],
+        default_basepath,
+        os.path.basename(default_basepath),
+        None,
+        None,
+        None,
+        None,
+    )
diff --git a/mapchete/cli/utils.py b/mapchete/cli/options.py
similarity index 57%
rename from mapchete/cli/utils.py
rename to mapchete/cli/options.py
index 7e5ebd14..372f807c 100644
--- a/mapchete/cli/utils.py
+++ b/mapchete/cli/options.py
@@ -8,7 +8,11 @@
 import tqdm
 
 import mapchete
-from mapchete.config import raw_conf, bounds_from_opts
+from mapchete.config import (
+    raw_conf,
+    bounds_from_opts,
+    MULTIPROCESSING_DEFAULT_START_METHOD,
+)
 from mapchete.formats import available_output_formats
 from mapchete.index import zoom_index_gen
 from mapchete.log import set_log_level, setup_logfile
@@ -19,23 +23,6 @@
 
 
 MULTIPROCESSING_START_METHODS = get_all_start_methods()
-MULTIPROCESSING_START_METHOD_DEFAULT = (
-    "fork" if "fork" in MULTIPROCESSING_START_METHODS else "spawn"
-)
-
-
-# verbose stdout writer #
-#########################
-def write_verbose_msg(process_info, dst):
-    tqdm.tqdm.write(
-        "Tile %s: %s, %s"
-        % (
-            tuple(process_info.tile.id),
-            process_info.process_msg,
-            process_info.write_msg,
-        ),
-        file=dst,
-    )
 
 
 # click callbacks #
@@ -68,14 +55,6 @@ def _validate_mapchete_files(ctx, param, mapchete_files):
     return mapchete_files
 
 
-def _validate_inputs(ctx, param, inputs):
-    if len(inputs) == 0:
-        raise click.MissingParameter(
-            "at least one mapchete file or path to Tile Directory required"
-        )
-    return inputs
-
-
 def _set_debug_log_level(ctx, param, debug):
     if debug:
         set_log_level(logging.DEBUG)
@@ -132,6 +111,10 @@ def _cb_key_val(ctx, param, value):
         return out
 
 
+def _cb_none_concurrency(ctx, param, value):
+    return None if value == "none" else value
+
+
 # click arguments #
 ###################
 arg_mapchete_file = click.argument("mapchete_file", type=click.Path(exists=True))
@@ -149,23 +132,22 @@ def _cb_key_val(ctx, param, value):
 arg_input_raster = click.argument("input_raster", type=click.Path(exists=True))
 arg_out_dir = click.argument("output_dir", type=click.Path())
 arg_input = click.argument("input_", metavar="INPUT", type=click.STRING)
-arg_inputs = click.argument(
-    "inputs", metavar="INPUTS", nargs=-1, callback=_validate_inputs
-)
 arg_output = click.argument("output", type=click.STRING)
+arg_src_tiledir = click.argument("src_tiledir", type=click.STRING)
+arg_dst_tiledir = click.argument("dst_tiledir", type=click.STRING)
+arg_tiledir = click.argument("tiledir", type=click.STRING)
 
 
 # click options #
 #################
 opt_out_path = click.option(
     "--out-path",
-    "-op",
     type=click.Path(),
     default=os.path.join(os.getcwd(), "output"),
-    help="Process output path.",
+    help="Output path.",
 )
 opt_idx_out_dir = click.option(
-    "--idx-out-dir", "-od", type=click.Path(), help="Index output directory."
+    "--idx-out-dir", type=click.Path(), help="Index output directory."
 )
 opt_input_file = click.option(
     "--input-file",
@@ -232,7 +214,13 @@ def _cb_key_val(ctx, param, value):
     "--multi",
     "-m",
     type=click.INT,
-    help="Number of concurrent processes.",
+    help="Number of workers when processing concurrently.",
+)
+opt_workers = click.option(
+    "--workers",
+    "-w",
+    type=click.INT,
+    help="Number of workers when processing concurrently.",
 )
 opt_force = click.option(
     "--force", "-f", is_flag=True, help="Overwrite if files already exist."
@@ -255,20 +243,13 @@ def _cb_key_val(ctx, param, value):
     callback=_set_debug_log_level,
     help="Deactivate progress bar and print debug log output.",
 )
-opt_max_chunksize = click.option(
-    "--max-chunksize",
-    "-c",
-    type=click.INT,
-    default=1,
-    help="Maximum number of process tiles to be queued for each  worker. (default: 1)",
-)
 opt_multiprocessing_start_method = click.option(
     "--multiprocessing-start-method",
     type=click.Choice(MULTIPROCESSING_START_METHODS),
-    default=MULTIPROCESSING_START_METHOD_DEFAULT,
+    default=MULTIPROCESSING_DEFAULT_START_METHOD,
     help=(
         "Method used by multiprocessing module to start child workers. Availability of "
-        f"methods depends on OS (default: {MULTIPROCESSING_START_METHOD_DEFAULT})"
+        f"methods depends on OS (default: {MULTIPROCESSING_DEFAULT_START_METHOD})"
     ),
 )
 opt_input_formats = click.option(
@@ -356,189 +337,29 @@ def _cb_key_val(ctx, param, value):
     metavar="NAME=VALUE",
     multiple=True,
     callback=_cb_key_val,
-    help="Configuration options for source fsspec filesystem. ",
+    help="Configuration options for source fsspec filesystem.",
 )
 opt_dst_fs_opts = click.option(
     "--dst-fs-opts",
     metavar="NAME=VALUE",
     multiple=True,
     callback=_cb_key_val,
-    help="Configuration options for destination fsspec filesystem. ",
+    help="Configuration options for destination fsspec filesystem.",
 )
 opt_fs_opts = click.option(
     "--fs-opts",
     metavar="NAME=VALUE",
     multiple=True,
     callback=_cb_key_val,
-    help="Configuration options for destination fsspec filesystem. ",
+    help="Configuration options for destination fsspec filesystem.",
+)
+opt_dask_scheduler = click.option(
+    "--dask-scheduler", type=click.STRING, help="Address for dask scheduler."
+)
+opt_concurrency = click.option(
+    "--concurrency",
+    type=click.Choice(["processes", "threads", "dask", "none"]),
+    default="processes",
+    callback=_cb_none_concurrency,
+    help="Decide which Executor to use for concurrent processing.",
 )
-
-
-# convenience processing functions #
-####################################
-def _process_single_tile(
-    debug=None,
-    raw_conf_process_pyramid=None,
-    mapchete_config=None,
-    tile=None,
-    mode=None,
-    input_file=None,
-    verbose_dst=None,
-    vrt=None,
-    idx_out_dir=None,
-    no_pbar=None,
-):
-    with click_spinner.spinner(disable=debug) as spinner:
-        with mapchete.Timer() as t:
-            tile = raw_conf_process_pyramid(raw_conf(mapchete_config)).tile(*tile)
-            with mapchete.open(
-                mapchete_config,
-                mode=mode,
-                bounds=tile.bounds,
-                zoom=tile.zoom,
-                single_input_file=input_file,
-            ) as mp:
-                spinner.stop()
-                tqdm.tqdm.write("processing 1 tile", file=verbose_dst)
-
-                # run process on tile
-                for result in mp.batch_processor(tile=tile):
-                    write_verbose_msg(result, dst=verbose_dst)
-
-            tqdm.tqdm.write(
-                (
-                    "processing %s finished in %s" % (mapchete_config, t)
-                    if isinstance(mapchete_config, str)
-                    else "processing finished in %s" % t
-                ),
-                file=verbose_dst,
-            )
-
-            # write VRT index
-            if vrt:
-                with mapchete.Timer() as t_vrt:
-                    tqdm.tqdm.write("creating VRT", file=verbose_dst)
-                    for tile in tqdm.tqdm(
-                        zoom_index_gen(
-                            mp=mp,
-                            zoom=tile.zoom,
-                            out_dir=idx_out_dir or mp.config.output.path,
-                            vrt=vrt,
-                        ),
-                        total=mp.count_tiles(tile.zoom, tile.zoom),
-                        unit="tile",
-                        disable=debug or no_pbar,
-                    ):
-                        logger.debug("%s indexed", tile)
-
-                    tqdm.tqdm.write(
-                        (
-                            "VRT(s) for %s created in %s" % (mapchete_config, t_vrt)
-                            if isinstance(mapchete_config, str)
-                            else "VRT(s) created in %s" % t_vrt
-                        ),
-                        file=verbose_dst,
-                    )
-
-
-def _process_area(
-    debug=None,
-    mapchete_config=None,
-    mode=None,
-    zoom=None,
-    wkt_geometry=None,
-    point=None,
-    point_crs=None,
-    bounds=None,
-    bounds_crs=None,
-    area=None,
-    area_crs=None,
-    input_file=None,
-    multi=None,
-    verbose_dst=None,
-    max_chunksize=None,
-    multiprocessing_start_method=None,
-    no_pbar=None,
-    vrt=None,
-    idx_out_dir=None,
-):
-    multi = multi or cpu_count()
-    with click_spinner.spinner(disable=debug) as spinner:
-        with mapchete.Timer() as t:
-            with mapchete.open(
-                mapchete_config,
-                mode=mode,
-                zoom=zoom,
-                bounds=bounds_from_opts(
-                    wkt_geometry=wkt_geometry,
-                    point=point,
-                    point_crs=point_crs,
-                    bounds=bounds,
-                    bounds_crs=bounds_crs,
-                    raw_conf=raw_conf(mapchete_config),
-                ),
-                area=area,
-                area_crs=area_crs,
-                single_input_file=input_file,
-            ) as mp:
-                spinner.stop()
-                tiles_count = mp.count_tiles(
-                    min(mp.config.init_zoom_levels), max(mp.config.init_zoom_levels)
-                )
-
-                tqdm.tqdm.write(
-                    "processing %s tile(s) on %s worker(s)" % (tiles_count, multi),
-                    file=verbose_dst,
-                )
-
-                # run process on tiles
-                for process_info in tqdm.tqdm(
-                    mp.batch_processor(
-                        multi=multi,
-                        zoom=zoom,
-                        max_chunksize=max_chunksize,
-                        multiprocessing_start_method=multiprocessing_start_method,
-                    ),
-                    total=tiles_count,
-                    unit="tile",
-                    disable=debug or no_pbar,
-                ):
-                    write_verbose_msg(process_info, dst=verbose_dst)
-
-            tqdm.tqdm.write(
-                (
-                    "processing %s finished in %s" % (mapchete_config, t)
-                    if isinstance(mapchete_config, str)
-                    else "processing finished in %s" % t
-                ),
-                file=verbose_dst,
-            )
-
-            # write VRT index
-            if vrt:
-                with mapchete.Timer() as t_vrt:
-                    tqdm.tqdm.write("creating VRT(s)", file=verbose_dst)
-                    for tile in tqdm.tqdm(
-                        zoom_index_gen(
-                            mp=mp,
-                            zoom=mp.config.init_zoom_levels,
-                            out_dir=idx_out_dir or mp.config.output.path,
-                            vrt=vrt,
-                        ),
-                        total=mp.count_tiles(
-                            min(mp.config.init_zoom_levels),
-                            max(mp.config.init_zoom_levels),
-                        ),
-                        unit="tile",
-                        disable=debug or no_pbar,
-                    ):
-                        logger.debug("%s indexed", tile)
-
-                    tqdm.tqdm.write(
-                        (
-                            "VRT(s) for %s created in %s" % (mapchete_config, t_vrt)
-                            if isinstance(mapchete_config, str)
-                            else "VRT(s) created in %s" % t_vrt
-                        ),
-                        file=verbose_dst,
-                    )
diff --git a/mapchete/commands/__init__.py b/mapchete/commands/__init__.py
new file mode 100644
index 00000000..6a3d5012
--- /dev/null
+++ b/mapchete/commands/__init__.py
@@ -0,0 +1,12 @@
+"""
+This package contains easy to access functions which otherwise would have to be called via the CLI.
+This should make the use from within other scripts, notebooks, etc. easier.
+"""
+from mapchete.commands._convert import convert
+from mapchete.commands._cp import cp
+from mapchete.commands._execute import execute
+from mapchete.commands._index import index
+from mapchete.commands._rm import rm
+
+
+__all__ = ["convert", "cp", "execute", "index", "rm"]
diff --git a/mapchete/commands/_convert.py b/mapchete/commands/_convert.py
new file mode 100644
index 00000000..fa08494b
--- /dev/null
+++ b/mapchete/commands/_convert.py
@@ -0,0 +1,411 @@
+import fiona
+import logging
+from multiprocessing import cpu_count
+import os
+from pprint import pformat
+import rasterio
+from rasterio.crs import CRS
+from shapely.geometry import box
+from shapely.geometry.base import BaseGeometry
+import tilematrix
+from typing import Callable, List, Tuple, Union
+import warnings
+
+import mapchete
+from mapchete.commands._execute import execute
+from mapchete.config import raw_conf, raw_conf_output_pyramid
+from mapchete.formats import (
+    driver_from_file,
+    available_output_formats,
+    available_input_formats,
+)
+from mapchete.io import read_json, get_best_zoom_level
+from mapchete.io.vector import reproject_geometry
+from mapchete.tile import BufferedTilePyramid
+from mapchete.validate import validate_zooms
+
+logger = logging.getLogger(__name__)
+OUTPUT_FORMATS = available_output_formats()
+
+
+def convert(
+    tiledir: str,
+    output: str,
+    zoom: Union[int, List[int]] = None,
+    area: Union[BaseGeometry, str, dict] = None,
+    area_crs: Union[CRS, str] = None,
+    bounds: Tuple[float] = None,
+    bounds_crs: Union[CRS, str] = None,
+    point: Tuple[float, float] = None,
+    point_crs: Tuple[float, float] = None,
+    tile: Tuple[int, int, int] = None,
+    overwrite: bool = False,
+    concurrency: str = "processes",
+    dask_scheduler: str = None,
+    dask_client=None,
+    workers: int = None,
+    multi: int = None,
+    clip_geometry: str = None,
+    bidx: List[int] = None,
+    output_pyramid: str = None,
+    output_metatiling: int = None,
+    output_format: str = None,
+    output_dtype: str = None,
+    output_geometry_type: str = None,
+    creation_options: dict = None,
+    scale_ratio: float = None,
+    scale_offset: float = None,
+    resampling_method: str = "nearest",
+    overviews: bool = False,
+    overviews_resampling_method: str = "cubic_spline",
+    cog: bool = False,
+    msg_callback: Callable = None,
+    as_iterator: bool = False,
+) -> mapchete.Job:
+    """
+    Convert mapchete outputs or other geodata.
+
+    This is a wrapper around the mapchete.processes.convert process which helps generating tiled
+    outputs for raster and vector data or single COGs from TileDirectory raster inputs.
+
+    It also supports clipping of the input by a vector dataset.
+
+    If only a subset of a TileDirectory is desired, please see the mapchete.commands.cp command.
+
+    Parameters
+    ----------
+    tiledir : str
+        Path to TileDirectory or mapchete config.
+    output : str
+        Path to output.
+    zoom : integer or list of integers
+        Single zoom, minimum and maximum zoom or a list of zoom levels.
+    area : str, dict, BaseGeometry
+        Geometry to override bounds or area provided in process configuration. Can be either a
+        WKT string, a GeoJSON mapping, a shapely geometry or a path to a Fiona-readable file.
+    area_crs : CRS or str
+        CRS of area (default: process CRS).
+    bounds : tuple
+        Override bounds or area provided in process configuration.
+    bounds_crs : CRS or str
+        CRS of area (default: process CRS).
+    point : iterable
+        X and y coordinates of point whose corresponding process tile bounds will be used.
+    point_crs : str or CRS
+        CRS of point (defaults to process pyramid CRS).
+    tile : tuple
+        Zoom, row and column of tile to be processed (cannot be used with zoom)
+    overwrite : bool
+        Overwrite existing output.
+    workers : int
+        Number of execution workers when processing concurrently.
+    concurrency : str
+        Concurrency to be used. Could either be "processes", "threads" or "dask".
+    dask_scheduler : str
+        URL to dask scheduler if required.
+    dask_client : dask.distributed.Client
+        Reusable Client instance if required. Otherwise a new client will be created.
+    clip_geometry : str
+        Path to Fiona-readable file by which output will be clipped.
+    bidx : list of integers
+        Band indexes to read from source.
+    output_pyramid : str
+        Output pyramid to write to.
+    output_metatiling : int
+        Output metatiling.
+    output_format : str
+        Output format. Can be any raster or vector format available by mapchete.
+    output_dtype : str
+        Output data type (for raster output only).
+    output_geometry_type :
+        Output geometry type (for vector output only).
+    creation_options : dict
+        Output driver specific creation options.
+    scale_ratio : float
+        Scaling factor (for raster output only).
+    scale_offset : float
+        Scaling offset (for raster output only).
+    resampling_method : str
+        Resampling method used. (default: nearest).
+    overviews : bool
+        Generate overviews (single GTiff output only).
+    overviews_resampling_method : str
+        Resampling method used for overviews. (default: cubic_spline)
+    cog : bool
+        Write a valid COG. This will automatically generate verviews. (GTiff only)
+    msg_callback : Callable
+        Optional callback function for process messages.
+    as_iterator : bool
+        Returns as generator but with a __len__() property.
+
+    Returns
+    -------
+    mapchete.Job instance either with already processed items or a generator with known length.
+
+    Examples
+    --------
+    >>> convert("foo", "bar")
+
+    This will run the whole conversion process.
+
+    >>> for i in convert("foo", "bar", as_iterator=True):
+    >>>     print(i)
+
+    This will return a generator where through iteration, tiles are copied.
+
+    >>> list(tqdm.tqdm(convert("foo", "bar", as_iterator=True)))
+
+    Usage within a process bar.
+    """
+
+    def _empty_callback(*args, **kwargs):
+        pass
+
+    msg_callback = msg_callback or _empty_callback
+    if multi is not None:  # pragma: no cover
+        warnings.warn("The 'multi' parameter is deprecated and is now named 'workers'")
+    workers = workers or multi or cpu_count()
+    creation_options = creation_options or {}
+    bidx = [bidx] if isinstance(bidx, int) else bidx
+    try:
+        input_info = _get_input_info(tiledir)
+        output_info = _get_output_info(output)
+    except Exception as e:
+        raise ValueError(e)
+
+    # collect mapchete configuration
+    mapchete_config = dict(
+        process="mapchete.processes.convert",
+        input=dict(inp=tiledir, clip=clip_geometry),
+        pyramid=(
+            dict(
+                grid=output_pyramid,
+                metatiling=(
+                    output_metatiling
+                    or (
+                        input_info["pyramid"].get("metatiling", 1)
+                        if input_info["pyramid"]
+                        else 1
+                    )
+                ),
+                pixelbuffer=(
+                    input_info["pyramid"].get("pixelbuffer", 0)
+                    if input_info["pyramid"]
+                    else 0
+                ),
+            )
+            if output_pyramid
+            else input_info["pyramid"]
+        ),
+        output=dict(
+            {
+                k: v
+                for k, v in input_info["output_params"].items()
+                if k not in ["delimiters", "bounds", "mode"]
+            },
+            path=output,
+            format=(
+                output_format
+                or output_info["driver"]
+                or input_info["output_params"]["format"]
+            ),
+            dtype=output_dtype or input_info["output_params"].get("dtype"),
+            **creation_options,
+            **dict(overviews=True, overviews_resampling=overviews_resampling_method)
+            if overviews
+            else dict(),
+        ),
+        config_dir=os.getcwd(),
+        zoom_levels=zoom or input_info["zoom_levels"],
+        scale_ratio=scale_ratio,
+        scale_offset=scale_offset,
+        resampling=resampling_method,
+        band_indexes=bidx,
+    )
+
+    # assert all required information is there
+    if mapchete_config["output"]["format"] is None:
+        # this happens if input file is e.g. JPEG2000 and output is a tile directory
+        raise ValueError("Output format required.")
+    if mapchete_config["output"]["format"] == "GTiff":
+        mapchete_config["output"].update(cog=cog)
+    output_type = OUTPUT_FORMATS[mapchete_config["output"]["format"]]["data_type"]
+    if bidx is not None:
+        mapchete_config["output"].update(bands=len(bidx))
+    if mapchete_config["pyramid"] is None:
+        raise ValueError("Output pyramid required.")
+    elif mapchete_config["zoom_levels"] is None:
+        try:
+            mapchete_config.update(
+                zoom_levels=dict(
+                    min=0,
+                    max=get_best_zoom_level(
+                        tiledir, mapchete_config["pyramid"]["grid"]
+                    ),
+                )
+            )
+        except Exception:
+            raise ValueError("Zoom levels required.")
+    elif input_info["input_type"] != output_type:
+        raise ValueError(
+            f"Output format type ({output_type}) is incompatible with input format ({input_info['input_type']})."
+        )
+    if output_metatiling:
+        mapchete_config["output"].update(metatiling=output_metatiling)
+    if input_info["output_params"].get("schema") and output_geometry_type:
+        mapchete_config["output"]["schema"].update(geometry=output_geometry_type)
+
+    # determine process bounds
+    out_pyramid = BufferedTilePyramid.from_dict(mapchete_config["pyramid"])
+    inp_bounds = (
+        bounds
+        or reproject_geometry(
+            box(*input_info["bounds"]),
+            src_crs=input_info["crs"],
+            dst_crs=out_pyramid.crs,
+        ).bounds
+        if input_info["bounds"]
+        else out_pyramid.bounds
+    )
+    # if clip-geometry is available, intersect determined bounds with clip bounds
+    if clip_geometry:
+        clip_intersection = _clip_bbox(
+            clip_geometry, dst_crs=out_pyramid.crs
+        ).intersection(box(*inp_bounds))
+        if clip_intersection.is_empty:
+            msg_callback(
+                "Process area is empty: clip bounds don't intersect with input bounds."
+            )
+            # this returns a Job with an empty iterator
+            return mapchete.Job(None, (), as_iterator=as_iterator, total=0)
+    # add process bounds and output type
+    mapchete_config.update(
+        bounds=(clip_intersection.bounds if clip_geometry else inp_bounds),
+        bounds_crs=bounds_crs,
+        clip_to_output_dtype=mapchete_config["output"].get("dtype", None),
+    )
+    logger.debug(f"temporary config generated: {pformat(mapchete_config)}")
+
+    return execute(
+        mapchete_config=mapchete_config,
+        mode="overwrite" if overwrite else "continue",
+        zoom=zoom,
+        point=point,
+        point_crs=point_crs,
+        bounds=bounds,
+        bounds_crs=bounds_crs,
+        area=area,
+        area_crs=area_crs,
+        concurrency=concurrency,
+        dask_scheduler=dask_scheduler,
+        dask_client=dask_client,
+        workers=workers,
+        as_iterator=as_iterator,
+        msg_callback=msg_callback,
+    )
+
+
+def _clip_bbox(clip_geometry, dst_crs=None):
+    with fiona.open(clip_geometry) as src:
+        return reproject_geometry(box(*src.bounds), src_crs=src.crs, dst_crs=dst_crs)
+
+
+def _get_input_info(tiledir):
+
+    # assuming single file if path has a file extension
+    if os.path.splitext(tiledir)[1]:
+        driver = driver_from_file(tiledir)
+
+        # single file input can be a mapchete file or a rasterio/fiona file
+        if driver == "Mapchete":
+            logger.debug("input is mapchete file")
+            input_info = _input_mapchete_info(tiledir)
+
+        elif driver == "raster_file":
+            # this should be readable by rasterio
+            logger.debug("input is raster_file")
+            input_info = _input_rasterio_info(tiledir)
+
+        elif driver == "vector_file":
+            # this should be readable by Fiona
+            input_info = _input_fiona_info(tiledir)
+        else:  # pragma: no cover
+            raise NotImplementedError(f"driver {driver} is not supported")
+
+    # assuming tile directory
+    else:
+        logger.debug("input is tile directory")
+        input_info = _input_tile_directory_info(tiledir)
+
+    return input_info
+
+
+def _input_mapchete_info(tiledir):
+    conf = raw_conf(tiledir)
+    output_params = conf["output"]
+    pyramid = raw_conf_output_pyramid(conf)
+    return dict(
+        output_params=output_params,
+        pyramid=pyramid.to_dict(),
+        crs=pyramid.crs,
+        zoom_levels=validate_zooms(conf["zoom_levels"], expand=False),
+        pixel_size=None,
+        input_type=OUTPUT_FORMATS[output_params["format"]]["data_type"],
+        bounds=conf.get("bounds"),
+    )
+
+
+def _input_rasterio_info(tiledir):
+    with rasterio.open(tiledir) as src:
+        return dict(
+            output_params=dict(
+                bands=src.meta["count"],
+                dtype=src.meta["dtype"],
+                format=src.driver if src.driver in available_input_formats() else None,
+            ),
+            pyramid=None,
+            crs=src.crs,
+            zoom_levels=None,
+            pixel_size=src.transform[0],
+            input_type="raster",
+            bounds=src.bounds,
+        )
+
+
+def _input_fiona_info(tiledir):
+    with fiona.open(tiledir) as src:
+        return dict(
+            output_params=dict(
+                schema=src.schema,
+                format=src.driver if src.driver in available_input_formats() else None,
+            ),
+            pyramid=None,
+            crs=src.crs,
+            zoom_levels=None,
+            input_type="vector",
+            bounds=src.bounds,
+        )
+
+
+def _input_tile_directory_info(tiledir):
+    conf = read_json(os.path.join(tiledir, "metadata.json"))
+    pyramid = BufferedTilePyramid.from_dict(conf["pyramid"])
+    return dict(
+        output_params=conf["driver"],
+        pyramid=pyramid.to_dict(),
+        crs=pyramid.crs,
+        zoom_levels=None,
+        pixel_size=None,
+        input_type=OUTPUT_FORMATS[conf["driver"]["format"]]["data_type"],
+        bounds=None,
+    )
+
+
+def _get_output_info(output):
+    _, file_ext = os.path.splitext(output)
+    if not file_ext:
+        return dict(type="TileDirectory", driver=None)
+    elif file_ext == ".tif":
+        return dict(type="SingleFile", driver="GTiff")
+    else:
+        raise TypeError(f"Could not determine output from extension: {file_ext}")
diff --git a/mapchete/commands/_cp.py b/mapchete/commands/_cp.py
new file mode 100644
index 00000000..7ca253fc
--- /dev/null
+++ b/mapchete/commands/_cp.py
@@ -0,0 +1,268 @@
+import fiona
+import json
+import logging
+from multiprocessing import cpu_count
+import os
+from rasterio.crs import CRS
+from shapely.geometry import box, shape
+from shapely.geometry.base import BaseGeometry
+from shapely.ops import unary_union
+from tilematrix import TilePyramid
+from typing import Callable, List, Tuple, Union
+import warnings
+
+import mapchete
+from mapchete.config import _guess_geometry
+from mapchete.formats import read_output_metadata
+from mapchete.io import fs_from_path, tiles_exist
+from mapchete.io.vector import reproject_geometry
+
+logger = logging.getLogger(__name__)
+
+
+def cp(
+    src_tiledir: str,
+    dst_tiledir: str,
+    zoom: Union[int, List[int]] = None,
+    area: Union[BaseGeometry, str, dict] = None,
+    area_crs: Union[CRS, str] = None,
+    bounds: Tuple[float] = None,
+    bounds_crs: Union[CRS, str] = None,
+    overwrite: bool = False,
+    workers: int = None,
+    multi: int = None,
+    concurrency: str = None,
+    dask_scheduler: str = None,
+    dask_client=None,
+    src_fs_opts: dict = None,
+    dst_fs_opts: dict = None,
+    msg_callback: Callable = None,
+    as_iterator: bool = False,
+) -> mapchete.Job:
+    """
+    Copy TileDirectory from source to destination.
+
+    Parameters
+    ----------
+    src_tiledir : str
+        Source TileDirectory or mapchete file.
+    dst_tiledir : str
+        Destination TileDirectory.
+    zoom : integer or list of integers
+        Single zoom, minimum and maximum zoom or a list of zoom levels.
+    area : str, dict, BaseGeometry
+        Geometry to override bounds or area provided in process configuration. Can be either a
+        WKT string, a GeoJSON mapping, a shapely geometry or a path to a Fiona-readable file.
+    area_crs : CRS or str
+        CRS of area (default: process CRS).
+    bounds : tuple
+        Override bounds or area provided in process configuration.
+    bounds_crs : CRS or str
+        CRS of area (default: process CRS).
+    overwrite : bool
+        Overwrite existing output.
+    workers : int
+        Number of threads used to check whether tiles exist.
+    concurrency : str
+        Concurrency to be used. Could either be "processes", "threads" or "dask".
+    dask_scheduler : str
+        URL to dask scheduler if required.
+    dask_client : dask.distributed.Client
+        Reusable Client instance if required. Otherwise a new client will be created.
+    src_fs_opts : dict
+        Configuration options for source fsspec filesystem.
+    dst_fs_opts : dict
+        Configuration options for destination fsspec filesystem.
+    msg_callback : Callable
+        Optional callback function for process messages.
+    as_iterator : bool
+        Returns as generator but with a __len__() property.
+
+    Returns
+    -------
+    mapchete.Job instance either with already processed items or a generator with known length.
+
+    Examples
+    --------
+    >>> cp("foo", "bar", zoom=5)
+
+    This will run the whole copy process.
+
+    >>> for i in cp("foo", "bar", zoom=5, as_iterator=True):
+    >>>     print(i)
+
+    This will return a generator where through iteration, tiles are copied.
+
+    >>> list(tqdm.tqdm(cp("foo", "bar", zoom=5, as_iterator=True)))
+
+    Usage within a process bar.
+    """
+
+    def _empty_callback(*args):
+        pass
+
+    msg_callback = msg_callback or _empty_callback
+    if multi is not None:  # pragma: no cover
+        warnings.warn("The 'multi' parameter is deprecated and is now named 'workers'")
+    workers = workers or multi or cpu_count()
+    src_fs_opts = src_fs_opts or {}
+    dst_fs_opts = dst_fs_opts or {}
+    if zoom is None:  # pragma: no cover
+        raise ValueError("zoom level(s) required")
+
+    src_fs = fs_from_path(src_tiledir, **src_fs_opts)
+    dst_fs = fs_from_path(dst_tiledir, **dst_fs_opts)
+
+    # open source tile directory
+    with mapchete.open(
+        src_tiledir,
+        zoom=zoom,
+        area=area,
+        area_crs=area_crs,
+        bounds=bounds,
+        bounds_crs=bounds_crs,
+        fs=src_fs,
+        fs_kwargs=src_fs_opts,
+    ) as src_mp:
+        tp = src_mp.config.output_pyramid
+
+        # copy metadata to destination if necessary
+        src_metadata = os.path.join(src_tiledir, "metadata.json")
+        dst_metadata = os.path.join(dst_tiledir, "metadata.json")
+        if not dst_fs.exists(dst_metadata):
+            msg = f"copy {src_metadata} to {dst_metadata}"
+            logger.debug(msg)
+            msg_callback(msg)
+            _copy(src_fs, src_metadata, dst_fs, dst_metadata)
+
+        with mapchete.open(
+            dst_tiledir,
+            zoom=zoom,
+            area=area,
+            area_crs=area_crs,
+            bounds=bounds,
+            bounds_crs=bounds_crs,
+            fs=dst_fs,
+            fs_kwargs=dst_fs_opts,
+        ) as dst_mp:
+            return mapchete.Job(
+                _copy_tiles,
+                fargs=(
+                    msg_callback,
+                    src_mp,
+                    dst_mp,
+                    tp,
+                    workers,
+                    src_fs,
+                    dst_fs,
+                    overwrite,
+                ),
+                executor_concurrency=concurrency,
+                executor_kwargs=dict(
+                    max_workers=workers,
+                    dask_scheduler=dask_scheduler,
+                    dask_client=dask_client,
+                ),
+                as_iterator=as_iterator,
+                total=src_mp.count_tiles(),
+            )
+
+
+def _copy_tiles(
+    msg_callback,
+    src_mp,
+    dst_mp,
+    tp,
+    workers,
+    src_fs,
+    dst_fs,
+    overwrite,
+    executor=None,
+):
+    for z in src_mp.config.init_zoom_levels:
+        msg_callback(f"copy tiles for zoom {z}...")
+        # materialize all tiles
+        aoi_geom = src_mp.config.area_at_zoom(z)
+        tiles = [
+            t
+            for t in tp.tiles_from_geom(aoi_geom, z)
+            # this is required to omit tiles touching the config area
+            if aoi_geom.intersection(t.bbox).area
+        ]
+
+        # check which source tiles exist
+        logger.debug("looking for existing source tiles...")
+        src_tiles_exist = {
+            tile: exists
+            for tile, exists in tiles_exist(
+                config=src_mp.config, output_tiles=tiles, multi=workers
+            )
+        }
+
+        logger.debug("looking for existing destination tiles...")
+        # chech which destination tiles exist
+        dst_tiles_exist = {
+            tile: exists
+            for tile, exists in tiles_exist(
+                config=dst_mp.config, output_tiles=tiles, multi=workers
+            )
+        }
+
+        # copy
+        copied = 0
+        for future in executor.as_completed(
+            _copy_tile,
+            tiles,
+            fargs=(
+                src_mp,
+                dst_mp,
+                src_tiles_exist,
+                dst_tiles_exist,
+                src_fs,
+                dst_fs,
+                overwrite,
+            ),
+        ):
+            c, msg = future.result()
+            copied += c
+            yield msg
+
+        msg_callback(f"{copied} tiles copied")
+
+
+def _copy_tile(
+    tile, src_mp, dst_mp, src_tiles_exist, dst_tiles_exist, src_fs, dst_fs, overwrite
+):
+    src_path = src_mp.config.output_reader.get_path(tile)
+    # only copy if source tile exists
+    if src_tiles_exist[tile]:
+        # skip if destination tile exists and overwrite is deactivated
+        if dst_tiles_exist[tile] and not overwrite:
+            msg = f"{tile}: destination tile exists"
+            logger.debug(msg)
+            return 0, msg
+        # copy from source to target
+        else:
+            dst_path = dst_mp.config.output_reader.get_path(tile)
+            _copy(src_fs, src_path, dst_fs, dst_path)
+            msg = f"{tile}: copy {src_path} to {dst_path}"
+            logger.debug(msg)
+            return 1, msg
+    else:
+        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
new file mode 100755
index 00000000..96bb5829
--- /dev/null
+++ b/mapchete/commands/_execute.py
@@ -0,0 +1,187 @@
+import logging
+from multiprocessing import cpu_count
+from rasterio.crs import CRS
+from shapely.geometry.base import BaseGeometry
+from typing import Callable, List, Tuple, Union
+import warnings
+
+import mapchete
+from mapchete.config import bounds_from_opts, raw_conf, raw_conf_process_pyramid
+
+logger = logging.getLogger(__name__)
+
+
+def execute(
+    mapchete_config: Union[str, dict],
+    zoom: Union[int, List[int]] = None,
+    area: Union[BaseGeometry, str, dict] = None,
+    area_crs: Union[CRS, str] = None,
+    bounds: Tuple[float] = None,
+    bounds_crs: Union[CRS, str] = None,
+    point: Tuple[float, float] = None,
+    point_crs: Tuple[float, float] = None,
+    tile: Tuple[int, int, int] = None,
+    overwrite: bool = False,
+    mode: str = "continue",
+    concurrency: str = "processes",
+    workers: int = None,
+    multi: int = None,
+    multiprocessing_start_method: str = None,
+    dask_scheduler: str = None,
+    dask_client=None,
+    msg_callback: Callable = None,
+    as_iterator: bool = False,
+) -> mapchete.Job:
+    """
+    Execute a Mapchete process.
+
+    Parameters
+    ----------
+    mapchete_config : str or dict
+        Mapchete configuration as file path or dictionary.
+    zoom : integer or list of integers
+        Single zoom, minimum and maximum zoom or a list of zoom levels.
+    area : str, dict, BaseGeometry
+        Geometry to override bounds or area provided in process configuration. Can be either a
+        WKT string, a GeoJSON mapping, a shapely geometry or a path to a Fiona-readable file.
+    area_crs : CRS or str
+        CRS of area (default: process CRS).
+    bounds : tuple
+        Override bounds or area provided in process configuration.
+    bounds_crs : CRS or str
+        CRS of area (default: process CRS).
+    point : iterable
+        X and y coordinates of point whose corresponding process tile bounds will be used.
+    point_crs : str or CRS
+        CRS of point (defaults to process pyramid CRS).
+    tile : tuple
+        Zoom, row and column of tile to be processed (cannot be used with zoom)
+    overwrite : bool
+        Overwrite existing output.
+    mode : str
+        Set process mode. One of "readonly", "continue" or "overwrite".
+    workers : int
+        Number of execution workers when processing concurrently.
+    multiprocessing_start_method : str
+        Method used by multiprocessing module to start child workers. Availability of methods
+        depends on OS.
+    concurrency : str
+        Concurrency to be used. Could either be "processes", "threads" or "dask".
+    dask_scheduler : str
+        URL to dask scheduler if required.
+    dask_client : dask.distributed.Client
+        Reusable Client instance if required. Otherwise a new client will be created.
+    msg_callback : Callable
+        Optional callback function for process messages.
+    as_iterator : bool
+        Returns as generator but with a __len__() property.
+
+    Returns
+    -------
+    mapchete.Job instance either with already processed items or a generator with known length.
+
+    Examples
+    --------
+    >>> execute("foo")
+
+    This will run the whole execute process.
+
+    >>> for i in execute("foo", as_iterator=True):
+    >>>     print(i)
+
+    This will return a generator where through iteration, tiles are copied.
+
+    >>> list(tqdm.tqdm(execute("foo", as_iterator=True)))
+
+    Usage within a process bar.
+    """
+    mode = "overwrite" if overwrite else mode
+
+    def _empty_callback(*args):
+        pass
+
+    msg_callback = msg_callback or _empty_callback
+    if multi is not None:  # pragma: no cover
+        warnings.warn("The 'multi' parameter is deprecated and is now named 'workers'")
+    workers = workers or multi or cpu_count()
+
+    if tile:
+        tile = raw_conf_process_pyramid(raw_conf(mapchete_config)).tile(*tile)
+        bounds = tile.bounds
+        zoom = tile.zoom
+    else:
+        bounds = bounds_from_opts(
+            point=point,
+            point_crs=point_crs,
+            bounds=bounds,
+            bounds_crs=bounds_crs,
+            raw_conf=raw_conf(mapchete_config),
+        )
+
+    # be careful opening mapchete not as context manager
+    mp = mapchete.open(
+        mapchete_config,
+        mode=mode,
+        bounds=bounds,
+        zoom=zoom,
+        area=area,
+        area_crs=area_crs,
+    )
+    try:
+        if tile:
+            tasks_count = 1 + mp.config.preprocessing_tasks_count()
+        else:
+            tasks_count = mp.count_tasks()
+        msg_callback(f"processing {tasks_count} task(s) on {workers} worker(s)")
+        # automatically use dask Executor if dask scheduler is defined
+        if dask_scheduler or dask_client:  # pragma: no cover
+            concurrency = "dask"
+        # use sequential Executor if only one tile or only one worker is defined
+        elif tasks_count == 1 or workers == 1:
+            logger.debug(
+                f"using sequential Executor because there is only one {'task' if tasks_count == 1 else 'worker'}"
+            )
+            concurrency = None
+        return mapchete.Job(
+            _process_everything,
+            fargs=(
+                msg_callback,
+                mp,
+            ),
+            fkwargs=dict(
+                tile=tile,
+                workers=workers,
+                zoom=None if tile else zoom,
+            ),
+            executor_concurrency=concurrency,
+            executor_kwargs=dict(
+                dask_scheduler=dask_scheduler,
+                dask_client=dask_client,
+                multiprocessing_start_method=multiprocessing_start_method,
+            ),
+            as_iterator=as_iterator,
+            total=tasks_count,
+        )
+    # explicitly exit the mp object on failure
+    except Exception:  # pragma: no cover
+        mp.__exit__(None, None, None)
+        raise
+
+
+def _process_everything(msg_callback, mp, executor=None, workers=None, **kwargs):
+    try:
+        for preprocessing_task_info in mp.batch_preprocessor(
+            executor=executor, workers=workers
+        ):  # pragma: no cover
+            yield preprocessing_task_info
+            msg_callback(preprocessing_task_info)
+        for process_info in mp.batch_processor(
+            executor=executor, workers=workers, **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:
+        mp.__exit__(None, None, None)
diff --git a/mapchete/commands/_index.py b/mapchete/commands/_index.py
new file mode 100644
index 00000000..49130dde
--- /dev/null
+++ b/mapchete/commands/_index.py
@@ -0,0 +1,153 @@
+import logging
+import os
+from rasterio.crs import CRS
+from shapely.geometry.base import BaseGeometry
+import sys
+import tqdm
+from typing import Callable, List, Tuple, Union
+
+import mapchete
+from mapchete.index import zoom_index_gen
+
+logger = logging.getLogger(__name__)
+
+
+def index(
+    tiledir: str,
+    idx_out_dir: str = None,
+    geojson: bool = False,
+    gpkg: bool = False,
+    shp: bool = False,
+    vrt: bool = False,
+    txt: bool = False,
+    fieldname: str = "location",
+    basepath: str = None,
+    for_gdal: bool = False,
+    zoom: Union[int, List[int]] = None,
+    area: Union[BaseGeometry, str, dict] = None,
+    area_crs: Union[CRS, str] = None,
+    bounds: Tuple[float] = None,
+    bounds_crs: Union[CRS, str] = None,
+    point: Tuple[float, float] = None,
+    point_crs: Tuple[float, float] = None,
+    tile: Tuple[int, int, int] = None,
+    fs_opts: dict = None,
+    msg_callback: Callable = None,
+    as_iterator: bool = False,
+    **kwargs,
+) -> mapchete.Job:
+    """
+    Create one or more indexes from a TileDirectory.
+
+    Parameters
+    ----------
+    tiledir : str
+        Source TileDirectory or mapchete file.
+    idx_out_dir : str
+        Alternative output dir for index. Defaults to TileDirectory path.
+    geojson : bool
+        Activate GeoJSON output.
+    gpkg : bool
+        Activate GeoPackage output.
+    shp : bool
+        Activate Shapefile output.
+    vrt : bool
+        Activate VRT output.
+    txt : bool
+        Activate TXT output.
+    fieldname : str
+        Field name which contains paths of tiles (default: "location").
+    basepath : str
+        Use custom base path for absolute paths instead of output path.
+    for_gdal : bool
+        Use GDAL compatible remote paths, i.e. add "/vsicurl/" before path.
+    zoom : integer or list of integers
+        Single zoom, minimum and maximum zoom or a list of zoom levels.
+    area : str, dict, BaseGeometry
+        Geometry to override bounds or area provided in process configuration. Can be either a
+        WKT string, a GeoJSON mapping, a shapely geometry or a path to a Fiona-readable file.
+    area_crs : CRS or str
+        CRS of area (default: process CRS).
+    bounds : tuple
+        Override bounds or area provided in process configuration.
+    bounds_crs : CRS or str
+        CRS of area (default: process CRS).
+    point : iterable
+        X and y coordinates of point whose corresponding process tile bounds will be used.
+    point_crs : str or CRS
+        CRS of point (defaults to process pyramid CRS).
+    tile : tuple
+        Zoom, row and column of tile to be processed (cannot be used with zoom)
+    fs_opts : dict
+        Configuration options for fsspec filesystem.
+    msg_callback : Callable
+        Optional callback function for process messages.
+    as_iterator : bool
+        Returns as generator but with a __len__() property.
+
+    Returns
+    -------
+    mapchete.Job instance either with already processed items or a generator with known length.
+
+    Examples
+    --------
+    >>> index("foo", vrt=True, zoom=5)
+
+    This will run the whole index process.
+
+    >>> for i in index("foo", vrt=True, zoom=5, as_iterator=True):
+    >>>     print(i)
+
+    This will return a generator where through iteration, tiles are copied.
+
+    >>> list(tqdm.tqdm(index("foo", vrt=True, zoom=5, as_iterator=True)))
+
+    Usage within a process bar.
+    """
+
+    if not any([geojson, gpkg, shp, txt, vrt]):
+        raise ValueError(
+            """At least one of '--geojson', '--gpkg', '--shp', '--vrt' or '--txt'"""
+            """must be provided."""
+        )
+
+    def _empty_callback(*args):
+        pass
+
+    msg_callback = msg_callback or _empty_callback
+    fs_opts = fs_opts or {}
+
+    msg_callback(f"create index(es) for {tiledir}")
+    # process single tile
+    with mapchete.open(
+        tiledir,
+        mode="readonly",
+        fs_kwargs=fs_opts,
+        zoom=tile[0] if tile else zoom,
+        point=point,
+        point_crs=point_crs,
+        bounds=bounds,
+        bounds_crs=bounds_crs,
+        area=area,
+        area_crs=area_crs,
+    ) as mp:
+
+        return mapchete.Job(
+            zoom_index_gen,
+            fkwargs=dict(
+                mp=mp,
+                zoom=None if tile else mp.config.init_zoom_levels,
+                tile=tile,
+                out_dir=idx_out_dir if idx_out_dir else mp.config.output.path,
+                geojson=geojson,
+                gpkg=gpkg,
+                shapefile=shp,
+                vrt=vrt,
+                txt=txt,
+                fieldname=fieldname,
+                basepath=basepath,
+                for_gdal=for_gdal,
+            ),
+            as_iterator=as_iterator,
+            total=1 if tile else mp.count_tiles(),
+        )
diff --git a/mapchete/commands/_rm.py b/mapchete/commands/_rm.py
new file mode 100644
index 00000000..415f3af0
--- /dev/null
+++ b/mapchete/commands/_rm.py
@@ -0,0 +1,160 @@
+import logging
+from rasterio.crs import CRS
+from shapely.geometry.base import BaseGeometry
+from typing import Callable, List, Tuple, Union
+
+import mapchete
+from mapchete.io import fs_from_path, tiles_exist
+
+logger = logging.getLogger(__name__)
+
+
+def rm(
+    tiledir: str,
+    zoom: Union[int, List[int]] = None,
+    area: Union[BaseGeometry, str, dict] = None,
+    area_crs: Union[CRS, str] = None,
+    bounds: Tuple[float] = None,
+    bounds_crs: Union[CRS, str] = None,
+    overwrite: bool = False,
+    multi: int = None,
+    fs_opts: dict = None,
+    msg_callback: Callable = None,
+    as_iterator: bool = False,
+) -> mapchete.Job:
+    """
+    Remove tiles from TileDirectory.
+
+    Parameters
+    ----------
+    tiledir : str
+        TileDirectory or mapchete file.
+    zoom : integer or list of integers
+        Single zoom, minimum and maximum zoom or a list of zoom levels.
+    area : str, dict, BaseGeometry
+        Geometry to override bounds or area provided in process configuration. Can be either a
+        WKT string, a GeoJSON mapping, a shapely geometry or a path to a Fiona-readable file.
+    area_crs : CRS or str
+        CRS of area (default: process CRS).
+    bounds : tuple
+        Override bounds or area provided in process configuration.
+    bounds_crs : CRS or str
+        CRS of area (default: process CRS).
+    multi : int
+        Number of threads used to check whether tiles exist.
+    fs_opts : dict
+        Configuration options for fsspec filesystem.
+    msg_callback : Callable
+        Optional callback function for process messages.
+    as_iterator : bool
+        Returns as generator but with a __len__() property.
+
+    Returns
+    -------
+    mapchete.Job instance either with already processed items or a generator with known length.
+
+    Examples
+    --------
+    >>> rm("foo", zoom=5)
+
+    This will run the whole rm process.
+
+    >>> for i in rm("foo", zoom=5, as_iterator=True):
+    >>>     print(i)
+
+    This will return a generator where through iteration, tiles are removed.
+
+    >>> list(tqdm.tqdm(rm("foo", zoom=5, as_iterator=True)))
+
+    Usage within a process bar.
+
+    """
+
+    def _empty_callback(*args):
+        pass
+
+    msg_callback = msg_callback or _empty_callback
+    fs_opts = fs_opts or {}
+    if zoom is None:  # pragma: no cover
+        raise ValueError("zoom level(s) required")
+
+    fs = fs_from_path(tiledir, **fs_opts)
+
+    with mapchete.open(
+        tiledir,
+        zoom=zoom,
+        area=area,
+        area_crs=area_crs,
+        bounds=bounds,
+        bounds_crs=bounds_crs,
+        fs=fs,
+        fs_kwargs=fs_opts,
+        mode="readonly",
+    ) as mp:
+        tp = mp.config.output_pyramid
+
+        tiles = {}
+        for z in mp.config.init_zoom_levels:
+            tiles[z] = []
+            # check which source tiles exist
+            logger.debug(f"looking for existing source tiles in zoom {z}...")
+            for tile, exists in tiles_exist(
+                config=mp.config,
+                output_tiles=[
+                    t
+                    for t in tp.tiles_from_geom(mp.config.area_at_zoom(z), z)
+                    # this is required to omit tiles touching the config area
+                    if mp.config.area_at_zoom(z).intersection(t.bbox).area
+                ],
+                multi=multi,
+            ):
+                if exists:
+                    tiles[z].append(tile)
+
+        paths = [
+            mp.config.output_reader.get_path(tile)
+            for zoom_tiles in tiles.values()
+            for tile in zoom_tiles
+        ]
+        return mapchete.Job(
+            _rm,
+            fargs=(
+                paths,
+                fs,
+                msg_callback,
+            ),
+            as_iterator=as_iterator,
+            total=len(paths),
+        )
+
+
+def _rm(paths, fs, msg_callback, recursive=False, **kwargs):
+    """
+    Remove one or multiple paths from file system.
+
+    Note: all paths have to be from the same file system!
+
+    Parameters
+    ----------
+    paths : str or list
+    fs : fsspec.FileSystem
+    """
+    logger.debug(f"got {len(paths)} path(s) on {fs}")
+
+    # s3fs enables multiple paths as input, so let's use this:
+    if "s3" in fs.protocol:  # pragma: no cover
+        fs.rm(paths, recursive=recursive)
+        for path in paths:
+            msg = f"deleted {path}"
+            logger.debug(msg)
+            yield msg
+
+    # otherwise, just iterate through the paths
+    else:
+        for path in paths:
+            fs.rm(path, recursive=recursive)
+            msg = f"deleted {path}"
+            logger.debug(msg)
+            yield msg
+
+    msg_callback(f"{len(paths)} tiles deleted")
diff --git a/mapchete/config.py b/mapchete/config.py
index 4b360670..c33957f9 100644
--- a/mapchete/config.py
+++ b/mapchete/config.py
@@ -14,6 +14,7 @@
 from collections import OrderedDict
 from copy import deepcopy
 import fiona
+import hashlib
 import importlib
 import inspect
 import logging
@@ -26,6 +27,7 @@
 from shapely.geometry.base import BaseGeometry
 from shapely.ops import cascaded_union
 import sys
+from tempfile import NamedTemporaryFile
 from tilematrix._funcs import Bounds
 import warnings
 
@@ -58,7 +60,7 @@
 
 # parameters which have to be provided in the configuration and their types
 _MANDATORY_PARAMETERS = [
-    ("process", (str, type(None))),  # path to .py file or module path
+    ("process", (str, list, type(None))),  # path to .py file or module path
     ("pyramid", dict),  # process pyramid definition
     ("input", (dict, type(None))),  # files & other types
     ("output", dict),  # process output parameters
@@ -91,8 +93,11 @@
     "metatiling",
     "delimiters",
     "mode",
+    "stac",
 ]
 
+MULTIPROCESSING_DEFAULT_START_METHOD = "spawn"
+
 
 class MapcheteConfig(object):
     """
@@ -197,6 +202,7 @@ def __init__(
         if mode not in ["memory", "continue", "readonly", "overwrite"]:
             raise MapcheteConfigError("unknown mode %s" % mode)
         self.mode = mode
+        self.preprocessing_tasks_finished = False
 
         # (1) assert mandatory params are available
         try:
@@ -206,7 +212,7 @@ def __init__(
 
         # (2) check user process
         self.config_dir = self._raw["config_dir"]
-        self.process_name = self.process_path = self._raw["process"]
+        self.process_name = self.process_path = self.process = self._raw["process"]
         if self.mode != "readonly":
             logger.debug("validating process code")
             self.process_func
@@ -322,6 +328,40 @@ def __init__(
         logger.debug("prepare output")
         self.output.prepare(process_area=self.area_at_zoom())
 
+    def input_at_zoom(self, key=None, zoom=None):
+        if zoom is None:  # pragma: no cover
+            raise ValueError("zoom not provided")
+        return self.input[get_hash(self._params_at_zoom[zoom]["input"][key])]
+
+    def preprocessing_tasks_per_input(self):
+        """Get all preprocessing tasks defined by the input drivers."""
+        return {
+            k: inp.preprocessing_tasks
+            for k, inp in self.input.items()
+            if inp is not None
+        }
+
+    def preprocessing_tasks(self):
+        return {
+            task_key: task
+            for preprocessing_tasks in self.preprocessing_tasks_per_input().values()
+            for task_key, task in preprocessing_tasks.items()
+        }
+
+    def preprocessing_tasks_count(self):
+        """Return number of unique preprocessing tasks."""
+        return len(self.preprocessing_tasks())
+
+    def set_preprocessing_task_result(self, task_key, result):
+        """Append preprocessing task result to input."""
+        found = False
+        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")
+
     @cached_property
     def zoom_levels(self):
         """Process zoom levels as defined in the configuration."""
@@ -548,7 +588,7 @@ def process_func(self):
             )
         else:
             return get_process_func(
-                process_path=self.process_path,
+                process=self.process,
                 config_dir=self.config_dir,
                 run_compile=True,
             )
@@ -590,7 +630,9 @@ def params_at_zoom(self, zoom):
         zoom level dependent process configuration
         """
         if zoom not in self.init_zoom_levels:
-            raise ValueError("zoom level not available with current configuration")
+            raise ValueError(
+                f"zoom level {zoom} not available with current configuration: {self.init_zoom_levels}"
+            )
         out = OrderedDict(
             self._params_at_zoom[zoom], input=OrderedDict(), output=self.output
         )
@@ -630,7 +672,9 @@ def area_at_zoom(self, zoom=None):
             return self._cache_full_process_area
         else:
             if zoom not in self.init_zoom_levels:
-                raise ValueError("zoom level not available with current configuration")
+                raise ValueError(
+                    f"zoom level {zoom} not available with current configuration: {self.init_zoom_levels}"
+                )
             return self._area_at_zoom(zoom)
 
     def _area_at_zoom(self, zoom):
@@ -799,12 +843,9 @@ def process_bounds(self, zoom=None):
         return self.bounds_at_zoom(zoom)
 
 
-def get_hash(x):
+def get_hash(x, length=16):
     """Return hash of x."""
-    if isinstance(x, str):
-        return hash(x)
-    elif isinstance(x, dict):
-        return hash(yaml.dump(x))
+    return hashlib.sha224(yaml.dump(dict(key=x)).encode()).hexdigest()[:length]
 
 
 def get_zoom_levels(process_zoom_levels=None, init_zoom_levels=None):
@@ -991,11 +1032,11 @@ def bounds_from_opts(
         return
 
 
-def get_process_func(process_path=None, config_dir=None, run_compile=False):
+def get_process_func(process=None, config_dir=None, run_compile=False):
     """Import and return process function."""
-    logger.debug("get process function from %s", process_path)
+    logger.debug(f"get process function from {process}")
     process_module = _load_process_module(
-        process_path=process_path, config_dir=config_dir, run_compile=run_compile
+        process=process, config_dir=config_dir, run_compile=run_compile
     )
     try:
         if hasattr(process_module, "Process"):
@@ -1006,37 +1047,51 @@ def get_process_func(process_path=None, config_dir=None, run_compile=False):
         if hasattr(process_module, "execute"):
             return process_module.execute
         else:
-            raise ImportError("No execute() function found in %s" % process_path)
+            raise ImportError("No execute() function found in %s" % process)
     except ImportError as e:
         raise MapcheteProcessImportError(e)
 
 
-def _load_process_module(process_path=None, config_dir=None, run_compile=False):
-    if process_path.endswith(".py"):
-        module_path = os.path.join(config_dir, process_path)
-        if not os.path.isfile(module_path):
-            raise MapcheteConfigError(f"{module_path} is not available")
-        try:
-            if run_compile:
-                py_compile.compile(module_path, doraise=True)
-            module_name = os.path.splitext(os.path.basename(module_path))[0]
-            # load module
-            spec = importlib.util.spec_from_file_location(module_name, module_path)
-            module = importlib.util.module_from_spec(spec)
-            spec.loader.exec_module(module)
-            # required to make imported module available using multiprocessing
-            sys.modules[module_name] = module
-            # configure process file logger
-            add_module_logger(module.__name__)
-        except py_compile.PyCompileError as e:
-            raise MapcheteProcessSyntaxError(e)
-        except ImportError as e:
-            raise MapcheteProcessImportError(e)
-    else:
-        try:
-            module = importlib.import_module(process_path)
-        except ImportError as e:
-            raise MapcheteProcessImportError(e)
+def _load_process_module(process=None, config_dir=None, run_compile=False):
+    tmpfile = None
+    try:
+        if isinstance(process, list):
+            tmpfile = NamedTemporaryFile(suffix=".py")
+            logger.debug(f"writing process code to temporary file {tmpfile.name}")
+            with open(tmpfile.name, "w") as dst:
+                for line in process:
+                    dst.write(line + "\n")
+            process = tmpfile.name
+        if process.endswith(".py"):
+            module_path = absolute_path(path=process, base_dir=config_dir)
+            if not os.path.isfile(module_path):
+                raise MapcheteConfigError(f"{module_path} is not available")
+            try:
+                if run_compile:
+                    py_compile.compile(module_path, doraise=True)
+                module_name = os.path.splitext(os.path.basename(module_path))[0]
+                # load module
+                spec = importlib.util.spec_from_file_location(module_name, module_path)
+                module = importlib.util.module_from_spec(spec)
+                spec.loader.exec_module(module)
+                # required to make imported module available using multiprocessing
+                sys.modules[module_name] = module
+                # configure process file logger
+                add_module_logger(module.__name__)
+            except py_compile.PyCompileError as e:
+                raise MapcheteProcessSyntaxError(e)
+            except ImportError as e:
+                raise MapcheteProcessImportError(e)
+        else:
+            try:
+                module = importlib.import_module(process)
+            except ImportError as e:
+                raise MapcheteProcessImportError(e)
+        logger.debug(f"return process func: {module}")
+    finally:
+        if tmpfile:
+            logger.debug(f"removing {tmpfile.name}")
+            tmpfile.close()
     return module
 
 
diff --git a/mapchete/formats/base.py b/mapchete/formats/base.py
index f0ae20dc..dbface2e 100644
--- a/mapchete/formats/base.py
+++ b/mapchete/formats/base.py
@@ -14,6 +14,7 @@
 import types
 import warnings
 
+from mapchete.config import get_hash
 from mapchete.errors import MapcheteProcessOutputError, MapcheteNodataTile
 from mapchete.formats import write_output_metadata
 from mapchete.io import makedirs, path_exists
@@ -55,6 +56,10 @@ def __init__(self, input_params, **kwargs):
         self.pyramid = input_params.get("pyramid")
         self.pixelbuffer = input_params.get("pixelbuffer")
         self.crs = self.pyramid.crs if self.pyramid else None
+        # collect preprocessing tasks to be run by the Executor
+        self.preprocessing_tasks = {}
+        # storage location of all preprocessing tasks
+        self.preprocessing_tasks_results = {}
 
     def open(self, tile, **kwargs):
         """
@@ -101,6 +106,30 @@ def cleanup(self):
         """Optional cleanup function called when Mapchete exits."""
         pass
 
+    def add_preprocessing_task(self, func, fargs=None, fkwargs=None, key=None):
+        """
+        Add longer running preprocessing function to be called right before processing.
+
+        Applied correctly this will speed up process initialization and if multiple tasks
+        are required they run in parallel as they are being passed on the Executor.
+        """
+        fargs = fargs or ()
+        if not isinstance(fargs, (tuple, list)):
+            fargs = (fargs,)
+        fkwargs = fkwargs or {}
+        key = f"{func}-{get_hash((func, fargs, fkwargs))}" if key is None else 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)
+
+    def get_preprocessing_task_result(self, 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]
+
 
 class InputTile(object):
     """
diff --git a/mapchete/formats/default/tile_directory.py b/mapchete/formats/default/tile_directory.py
index 8b046295..cb86e8b2 100644
--- a/mapchete/formats/default/tile_directory.py
+++ b/mapchete/formats/default/tile_directory.py
@@ -14,7 +14,7 @@
     read_output_metadata,
 )
 from mapchete.io import path_exists, absolute_path, tile_to_zoom_level
-from mapchete.io.vector import reproject_geometry
+from mapchete.io.vector import reproject_geometry, segmentize_geometry
 from mapchete.tile import BufferedTilePyramid
 
 
@@ -187,6 +187,7 @@ def bbox(self, out_crs=None):
             box(*self._bounds),
             src_crs=self.td_pyramid.crs,
             dst_crs=self.pyramid.crs if out_crs is None else out_crs,
+            segmentize_on_clip=True,
         )
 
 
diff --git a/mapchete/index.py b/mapchete/index.py
index 4fdf4c10..74ae3455 100644
--- a/mapchete/index.py
+++ b/mapchete/index.py
@@ -52,6 +52,7 @@ def zoom_index_gen(
     mp=None,
     out_dir=None,
     zoom=None,
+    tile=None,
     geojson=False,
     gpkg=False,
     shapefile=False,
@@ -61,6 +62,7 @@ def zoom_index_gen(
     basepath=None,
     for_gdal=True,
     threading=False,
+    **kwargs,
 ):
     """
     Generate indexes for given zoom level.
@@ -91,6 +93,9 @@ def zoom_index_gen(
         use GDAL compatible remote paths, i.e. add "/vsicurl/" before path
         (default: True)
     """
+    if tile and zoom:  # pragma: no cover
+        raise ValueError("tile and zoom cannot be used at the same time")
+    zoom = tile[0] if tile else zoom
     for zoom in get_zoom_levels(process_zoom_levels=zoom):
         with ExitStack() as es:
             # get index writers for all enabled formats
@@ -149,28 +154,36 @@ def zoom_index_gen(
 
             # all output tiles for given process area
             logger.debug("determine affected output tiles")
-            output_tiles = set(
-                [
-                    tile
-                    for tile in mp.config.output_pyramid.tiles_from_geom(
-                        mp.config.area_at_zoom(zoom), zoom
+            logger.debug(f"HANSE {zoom} {tile}")
+            if tile:
+                output_tiles = set(
+                    mp.config.output_pyramid.intersecting(
+                        mp.config.process_pyramid.tile(*tile)
                     )
-                    # this is required to omit tiles touching the config area
-                    if tile.bbox.intersection(mp.config.area_at_zoom(zoom)).area
-                ]
-            )
-
+                )
+            else:
+                output_tiles = set(
+                    [
+                        t
+                        for t in mp.config.output_pyramid.tiles_from_geom(
+                            mp.config.area_at_zoom(zoom), zoom
+                        )
+                        # this is required to omit tiles touching the config area
+                        if t.bbox.intersection(mp.config.area_at_zoom(zoom)).area
+                    ]
+                )
+            logger.debug(f"HERBERT {zoom} {output_tiles}")
             # check which tiles exist in any index
             logger.debug("check which tiles exist in index(es)")
             existing_in_any_index = set(
-                tile
-                for tile in output_tiles
+                t
+                for t in output_tiles
                 if any(
                     [
                         i.entry_exists(
-                            tile=tile,
+                            tile=t,
                             path=_tile_path(
-                                orig_path=mp.config.output.get_path(tile),
+                                orig_path=mp.config.output.get_path(t),
                                 basepath=basepath,
                                 for_gdal=for_gdal,
                             ),
@@ -186,46 +199,46 @@ def zoom_index_gen(
                 )
             )
             # tiles which do not exist in any index
-            for tile, output_exists in tiles_exist(
+            for t, output_exists in tiles_exist(
                 mp.config, output_tiles=output_tiles.difference(existing_in_any_index)
             ):
                 tile_path = _tile_path(
-                    orig_path=mp.config.output.get_path(tile),
+                    orig_path=mp.config.output.get_path(t),
                     basepath=basepath,
                     for_gdal=for_gdal,
                 )
                 indexes = [
                     i
                     for i in index_writers
-                    if not i.entry_exists(tile=tile, path=tile_path)
+                    if not i.entry_exists(tile=t, path=tile_path)
                 ]
                 if indexes and output_exists:
                     logger.debug("%s exists", tile_path)
                     logger.debug("write to %s indexes" % len(indexes))
                     for index in indexes:
-                        index.write(tile, tile_path)
+                        index.write(t, tile_path)
                 # yield tile for progress information
-                yield tile
+                yield t
 
             # tiles which exist in at least one index
-            for tile in existing_in_any_index:
+            for t in existing_in_any_index:
                 tile_path = _tile_path(
-                    orig_path=mp.config.output.get_path(tile),
+                    orig_path=mp.config.output.get_path(t),
                     basepath=basepath,
                     for_gdal=for_gdal,
                 )
                 indexes = [
                     i
                     for i in index_writers
-                    if not i.entry_exists(tile=tile, path=tile_path)
+                    if not i.entry_exists(tile=t, path=tile_path)
                 ]
                 if indexes:
                     logger.debug("%s exists", tile_path)
                     logger.debug("write to %s indexes" % len(indexes))
                     for index in indexes:
-                        index.write(tile, tile_path)
+                        index.write(t, tile_path)
                 # yield tile for progress information
-                yield tile
+                yield t
 
 
 def _index_file_path(out_dir, zoom, ext):
diff --git a/mapchete/io/__init__.py b/mapchete/io/__init__.py
index 00b66cdf..8442ad35 100644
--- a/mapchete/io/__init__.py
+++ b/mapchete/io/__init__.py
@@ -16,7 +16,6 @@
     fs_from_path,
     path_is_remote,
     path_exists,
-    rm,
     tiles_exist,
     absolute_path,
     relative_path,
@@ -35,7 +34,6 @@
     "tiles_exist",
     "absolute_path",
     "relative_path",
-    "rm",
     "makedirs",
     "write_json",
     "read_json",
diff --git a/mapchete/io/_geometry_operations.py b/mapchete/io/_geometry_operations.py
index d915d4e0..bcda9c76 100644
--- a/mapchete/io/_geometry_operations.py
+++ b/mapchete/io/_geometry_operations.py
@@ -1,4 +1,5 @@
 from fiona.transform import transform_geom
+import pyproj
 from rasterio.crs import CRS
 from shapely.errors import TopologicalError
 from shapely.geometry import (
@@ -35,6 +36,8 @@ def reproject_geometry(
     src_crs=None,
     dst_crs=None,
     error_on_clip=False,
+    segmentize_on_clip=False,
+    segmentize_fraction=100,
     validity_check=True,
     antimeridian_cutting=False,
 ):
@@ -89,20 +92,38 @@ def _reproject_geom(geometry, src_crs, dst_crs):
     # geometry needs to be clipped to its CRS bounds
     elif (
         dst_crs.is_epsg_code
-        and dst_crs.get("init") in CRS_BOUNDS  # just in case for an CRS with EPSG code
-        and dst_crs.get("init")  # if CRS has defined bounds
+        and dst_crs.get("init")
         != "epsg:4326"  # and is not WGS84 (does not need clipping)
+        and (
+            dst_crs.get("init") in CRS_BOUNDS
+            or pyproj.CRS(dst_crs.to_epsg()).area_of_use.bounds
+        )
     ):
         wgs84_crs = CRS().from_epsg(4326)
         # get dst_crs boundaries
-        crs_bbox = box(*CRS_BOUNDS[dst_crs.get("init")])
+        crs_bbox = box(
+            *CRS_BOUNDS.get(
+                dst_crs.get("init"), pyproj.CRS(dst_crs.to_epsg()).area_of_use.bounds
+            )
+        )
         # reproject geometry to WGS84
         geometry_4326 = _reproject_geom(geometry, src_crs, wgs84_crs)
         # raise error if geometry has to be clipped
         if error_on_clip and not geometry_4326.within(crs_bbox):
             raise RuntimeError("geometry outside target CRS bounds")
+
+        clipped = crs_bbox.intersection(geometry_4326)
+
+        # segmentize clipped geometry using one 100th of with or height depending on
+        # which is shorter
+        if segmentize_on_clip:
+            height = clipped.bounds[3] - clipped.bounds[1]
+            width = clipped.bounds[2] - clipped.bounds[0]
+            segmentize_value = min([height, width]) / segmentize_fraction
+            clipped = segmentize_geometry(clipped, segmentize_value)
+
         # clip geometry dst_crs boundaries and return
-        return _reproject_geom(crs_bbox.intersection(geometry_4326), wgs84_crs, dst_crs)
+        return _reproject_geom(clipped, wgs84_crs, dst_crs)
 
     # return without clipping if destination CRS does not have defined bounds
     else:
diff --git a/mapchete/io/_json.py b/mapchete/io/_json.py
index edb1d058..2907a363 100644
--- a/mapchete/io/_json.py
+++ b/mapchete/io/_json.py
@@ -3,7 +3,7 @@
 import logging
 import os
 
-from mapchete.io._path import fs_from_path, makedirs, path_exists
+from mapchete.io._path import fs_from_path, path_exists
 from mapchete.io._misc import get_boto3_bucket
 
 
@@ -15,7 +15,7 @@ def write_json(path, params, fs=None, **kwargs):
     logger.debug(f"write {params} to {path}")
     fs = fs or fs_from_path(path, **kwargs)
     # using python 3.7 or higher we can use the fs.mkdir() call
-    makedirs(os.path.dirname(path))
+    fs.mkdirs(os.path.dirname(path), exist_ok=True)
     with fs.open(path, "w") as dst:
         json.dump(params, dst, sort_keys=True, indent=4)
 
diff --git a/mapchete/io/_path.py b/mapchete/io/_path.py
index 1730e803..48bbf06a 100644
--- a/mapchete/io/_path.py
+++ b/mapchete/io/_path.py
@@ -100,31 +100,6 @@ def makedirs(path):
             pass
 
 
-def rm(paths, fs=None, recursive=False):
-    """
-    Remove one or multiple paths from file system.
-
-    Note: all paths have to be from the same file system!
-
-    Parameters
-    ----------
-    paths : str or list
-    fs : fsspec.FileSystem
-    """
-    paths = [paths] if isinstance(paths, str) else paths
-    fs = fs or fs_from_path(paths[0])
-    logger.debug(f"got {len(paths)} path(s) on {fs}")
-
-    # s3fs enables multiple paths as input, so let's use this:
-    if "s3" in fs.protocol:  # pragma: no cover
-        fs.rm(paths, recursive=recursive)
-
-    # otherwise, just iterate through the paths
-    else:
-        for path in paths:
-            fs.rm(path, recursive=recursive)
-
-
 def tiles_exist(config=None, output_tiles=None, process_tiles=None, multi=None):
     """
     Yield tiles and whether their output already exists or not.
@@ -270,7 +245,6 @@ def _s3_tiledirectories(
 
 def fs_from_path(path, timeout=5, session=None, username=None, password=None, **kwargs):
     """Guess fsspec FileSystem from path and initialize using the desired options."""
-    logger.debug(f"create new FileSystem for path {path}")
     if path.startswith("s3://"):
         return fsspec.filesystem(
             "s3",
diff --git a/mapchete/io/raster.py b/mapchete/io/raster.py
index 41220ec8..eb355239 100644
--- a/mapchete/io/raster.py
+++ b/mapchete/io/raster.py
@@ -298,12 +298,14 @@ def _read(
             ),
             resampling=Resampling[resampling],
         ) as vrt:
-            return vrt.read(
-                window=vrt.window(*dst_bounds),
-                out_shape=dst_shape,
-                indexes=indexes,
-                masked=True,
-            )
+            with warnings.catch_warnings():
+                warnings.simplefilter("ignore")
+                return vrt.read(
+                    window=vrt.window(*dst_bounds),
+                    out_shape=dst_shape,
+                    indexes=indexes,
+                    masked=True,
+                )
 
     if isinstance(input_file, str):
         logger.debug("got file path %s", input_file)
diff --git a/mapchete/processes/__init__.py b/mapchete/processes/__init__.py
index 490ba50a..bb58e52f 100644
--- a/mapchete/processes/__init__.py
+++ b/mapchete/processes/__init__.py
@@ -37,7 +37,7 @@ def _import():
     return imported
 
 
-def process_names_docstrings(process_name):
+def process_names_docstrings(process_name=None):
     """
     Return registered process module names and docstrings.
 
diff --git a/mapchete/stac.py b/mapchete/stac.py
new file mode 100644
index 00000000..0b738bfd
--- /dev/null
+++ b/mapchete/stac.py
@@ -0,0 +1,250 @@
+from collections import OrderedDict
+import datetime
+import os
+from shapely.geometry import box, mapping
+
+from mapchete.io.vector import reproject_geometry
+
+
+OUT_PIXEL_SIZE = 0.28e-3
+UNIT_TO_METER = {"mercator": 1, "geodetic": 111319.4907932732}
+KNOWN_MATRIX_PROPERTIES = {
+    "geodetic": {
+        "name": "WorldCRS84Quad",
+        "url": "http://schemas.opengis.net/tms/1.0/json/examples/WorldCRS84Quad.json",
+        "title": "CRS84 for the World",
+        "crs": "http://www.opengis.net/def/crs/OGC/1.3/CRS84",
+        "supportedCRS": "http://www.opengis.net/def/crs/OGC/1.3/CRS84",
+        "wellKnownScaleSet": "http://www.opengis.net/def/wkss/OGC/1.0/GoogleCRS84Quad",
+    },
+    "mercator": {
+        "name": "WebMercatorQuad",
+        "url": "http://schemas.opengis.net/tms/1.0/json/examples/WebMercatorQuad.json",
+        "title": "Google Maps Compatible for the World",
+        "crs": "http://www.opengis.net/def/crs/EPSG/0/3857",
+        "supportedCRS": "http://www.opengis.net/def/crs/EPSG/0/3857",
+        "wellKnownScaleSet": "http://www.opengis.net/def/wkss/OGC/1.0/GoogleMapsCompatible",
+    },
+}
+
+
+def tile_directory_stac_item(
+    item_id=None,
+    tile_pyramid=None,
+    max_zoom=None,
+    item_path=None,
+    asset_basepath=None,
+    relative_paths=False,
+    min_zoom=0,
+    item_metadata=None,
+    bounds=None,
+    bounds_crs=None,
+    bands_type="image/tiff; application=geotiff",
+    crs_unit_to_meter=1,
+):
+    """
+    Create STAC item from a Tile Directory using the tiled-assets.
+
+    Currently only works for one asset.
+
+    Parameters
+    ----------
+    item_id : str
+        Unique item ID. (required)
+    tile_pyramid : mapchete.tile.BufferedTilePyramid
+        Tile pyramid of tiled asset. (required)
+    max_zoom : int
+        Maximum zoom of tile directory. (required)
+    item_path : str
+        Path of stac item.
+    asset_basepath : str
+        Base path of asset.
+    relative_paths : bool
+        Interpret asset path schema as being relative to item. If not set, one of item_path
+        or asset_basepath has to be set. (default: False)
+    min_zoom : int
+        Minimum zoom of tile directory. (default: 0)
+    item_metadata : dict
+        Optional additional item metadata to be appended.
+    bounds : list or tuple
+        Tile directory bounds. (default: tile_pyramid.bounds)
+    bounds_crs : CRS
+        CRS of bounds. (default: tile_pyramid.crs)
+    bands_type : str
+        Media typ of tile directory tiles. (default: "image/tiff; application=geotiff")
+    crs_unit_to_meter : int or float
+        Factor to convert CRS units to meter if tile pyramid grid is not "geodetic" or "mercator". (default: 1)
+
+    Returns
+    -------
+    pystac.Item
+    """
+    try:
+        import pystac
+        from pystac.version import get_stac_version
+    except ImportError:  # pragma: no cover
+        raise ImportError(
+            "dependencies for extra mapchete[stac] is required for this feature"
+        )
+
+    if item_id is None:
+        raise ValueError("item_id must be set")
+    if max_zoom is None:
+        raise ValueError("max_zoom must be set")
+    if tile_pyramid is None:
+        raise ValueError("tile_pyramid must be set")
+
+    item_metadata = _cleanup_datetime(item_metadata or {})
+    timestamp = (
+        item_metadata.get("properties", {}).get("start_datetime")
+        or item_metadata.get("properties", {}).get("end_datetime")
+        or str(datetime.datetime.utcnow())
+    )
+    tp_grid = tile_pyramid.grid.type
+    bands_schema = "{TileMatrix}/{TileRow}/{TileCol}.tif"
+    # thumbnail_href = thumbnail_href or "0/0/0.tif"
+    # thumbnail_type = thumbnail_type or "image/tiff; application=geotiff"
+    if asset_basepath:
+        bands_schema = os.path.join(asset_basepath, bands_schema)
+        # thumbnail_href = os.path.join(asset_basepath, thumbnail_href)
+    elif not relative_paths:
+        if item_path is None:
+            raise ValueError("either alternative_basepath or item_path must be set")
+        bands_schema = os.path.join(os.path.dirname(item_path), bands_schema)
+        # thumbnail_href = os.path.join(os.path.dirname(item_path), thumbnail_href)
+
+    # use bounds provided or fall back to tile pyramid bounds
+    bounds = bounds or tile_pyramid.bounds
+    bounds_crs = bounds_crs or tile_pyramid.crs
+
+    # bounds in tilepyramid CRS
+    left, bottom, right, top = reproject_geometry(
+        box(*bounds), src_crs=bounds_crs, dst_crs=tile_pyramid.crs
+    ).bounds
+
+    # bounds in lat/lon
+    geometry_4326 = reproject_geometry(
+        box(*bounds), src_crs=bounds_crs, dst_crs="EPSG:4326"
+    )
+    bounds_4326 = [*geometry_4326.bounds]
+
+    # tiles:tile_matrix_set object:
+    # http://schemas.opengis.net/tms/1.0/json/tms-schema.json
+    if tp_grid in KNOWN_MATRIX_PROPERTIES:
+        tile_matrix_set_identifier = KNOWN_MATRIX_PROPERTIES[tp_grid]["name"]
+        crs = KNOWN_MATRIX_PROPERTIES[tp_grid]["crs"]
+        tile_matrix_set = {
+            "type": "TileMatrixSetType",
+            "title": KNOWN_MATRIX_PROPERTIES[tp_grid]["title"],
+            "identifier": tile_matrix_set_identifier,
+            "supportedCRS": KNOWN_MATRIX_PROPERTIES[tp_grid]["supportedCRS"],
+            "wellKnownScaleSet": KNOWN_MATRIX_PROPERTIES[tp_grid]["wellKnownScaleSet"],
+        }
+    else:
+        tile_matrix_set_identifier = "custom"  # ???
+        authority, code = tile_pyramid.crs.to_authority()
+        crs = f"urn:ogc:def:crs:{authority}::{code}"
+        # from http://schemas.opengis.net/tms/1.0/json/tms-schema.json
+        # "required": ["type", "identifier", "supportedCRS", "tileMatrix"],
+        tile_matrix_set = {
+            "type": "TileMatrixSetType",
+            "identifier": tile_matrix_set_identifier,
+            "supportedCRS": crs,
+        }
+    tile_matrix_set.update(
+        tileMatrix=[
+            {
+                "type": "TileMatrixType",
+                "identifier": str(zoom),
+                "scaleDenominator": _scale(
+                    tp_grid, tile_pyramid.pixel_x_size(zoom), crs_unit_to_meter
+                ),
+                "topLeftCorner": [tile_pyramid.bounds.left, tile_pyramid.bounds.top],
+                "tileWidth": tile_pyramid.tile_width(zoom),
+                "tileHeight": tile_pyramid.tile_height(zoom),
+                "matrixWidth": tile_pyramid.matrix_width(zoom),
+                "matrixHeight": tile_pyramid.matrix_height(zoom),
+            }
+            for zoom in range(min_zoom, max_zoom + 1)
+        ],
+        boundingBox={
+            "type": "BoundingBoxType",
+            "crs": crs,
+            "lowerCorner": [tile_pyramid.bounds.left, tile_pyramid.bounds.bottom],
+            "upperCorner": [tile_pyramid.bounds.right, tile_pyramid.bounds.top],
+        },
+    )
+
+    # tiles:tile_matrix_links object:
+    tile_matrix_links = {
+        "url": f"#{tile_matrix_set_identifier}",
+        "limits": {
+            str(zoom): {
+                "min_tile_col": tile_pyramid.tile_from_xy(
+                    left, top, zoom, on_edge_use="rb"
+                ).col,
+                "max_tile_col": tile_pyramid.tile_from_xy(
+                    right, bottom, zoom, on_edge_use="lt"
+                ).col,
+                "min_tile_row": tile_pyramid.tile_from_xy(
+                    left, top, zoom, on_edge_use="rb"
+                ).row,
+                "max_tile_row": tile_pyramid.tile_from_xy(
+                    right, bottom, zoom, on_edge_use="lt"
+                ).row,
+            }
+            for zoom in range(min_zoom, max_zoom + 1)
+        },
+    }
+
+    stac_extensions = ["tiled-assets"]
+    if "eo:bands" in item_metadata:
+        stac_extensions.append("eo")
+
+    out = {
+        "stac_version": get_stac_version(),
+        "stac_extensions": stac_extensions,
+        "id": item_id,
+        "type": "Feature",
+        "bbox": bounds_4326,
+        "geometry": mapping(geometry_4326),
+        "properties": {
+            **item_metadata.get("properties", {}),
+            "datetime": timestamp,
+            "collection": item_id,
+            "tiles:tile_matrix_links": {tile_matrix_set_identifier: tile_matrix_links},
+            "tiles:tile_matrix_sets": {tile_matrix_set_identifier: tile_matrix_set},
+        },
+        "asset_templates": {"bands": {"href": bands_schema, "type": bands_type}},
+        "assets": {
+            # "thumbnail": {
+            #     "href": thumbnail_href,
+            #     "title": "Overview of the whole tiled dataset.",
+            #     "type": thumbnail_type,
+            # }
+        },
+    }
+    if "eo:bands" in item_metadata:
+        out["asset_templates"]["bands"]["eo:bands"] = item_metadata["eo:bands"]
+        # out["assets"]["thumbnail"]["eo:bands"] = item_metadata["eo:bands"]
+    out["links"] = item_metadata.get("links", [])
+    if item_path:
+        out["links"].extend([{"rel": "self", "href": item_path}])
+
+    return pystac.read_dict(out)
+
+
+def _scale(grid, pixel_x_size, unit_to_meter=1):
+    return UNIT_TO_METER.get(grid, unit_to_meter) * pixel_x_size / OUT_PIXEL_SIZE
+
+
+def _cleanup_datetime(d):
+    """Convert datetime objects in dictionary to strings."""
+    return OrderedDict(
+        (k, _cleanup_datetime(v))
+        if isinstance(v, dict)
+        else (k, str(v))
+        if isinstance(v, datetime.date)
+        else (k, v)
+        for k, v in d.items()
+    )
diff --git a/mapchete/validate.py b/mapchete/validate.py
index ce4289b2..81b9971c 100644
--- a/mapchete/validate.py
+++ b/mapchete/validate.py
@@ -124,7 +124,9 @@ def validate_values(config, values):
         if value not in config:
             raise ValueError("%s not given" % value)
         if not isinstance(config[value], vtype):
-            raise TypeError("%s must be %s" % (value, vtype))
+            if config[value] is None:
+                raise ValueError("%s not given" % value)
+            raise TypeError("%s must be %s, not %s" % (value, vtype, config[value]))
     return True
 
 
diff --git a/requirements.txt b/requirements.txt
index 945401b5..b9802bd6 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -9,7 +9,7 @@ click-plugins>=1.1.1
 click-spinner>=0.1.10
 fiona>=1.8.13.post1
 Flask>=1.1.2
-fsspec[s3]==0.8.4
+fsspec[s3]>=0.8.4
 geobuf>=1.1.1
 importlib-metadata>=1.6.0
 importlib-resources>=1.5.0
@@ -17,9 +17,11 @@ lxml>=4.5.0
 matplotlib>=3.2.1
 numpy>=1.16
 oyaml>=0.9
+pyproj
 retry>=0.9.2
+rasterio>=1.0.28,<1.2.7
 rasterio>=1.0.28
-s3fs==0.5.2
+s3fs>=0.5.2
 Shapely>=1.7.1
 tilematrix>=0.20
 tqdm>=4.46.0
diff --git a/setup.py b/setup.py
index 4d6288bf..fc4d8030 100644
--- a/setup.py
+++ b/setup.py
@@ -32,19 +32,31 @@
     "importlib-resources",
     "numpy>=1.16",
     "oyaml",
+    "pyproj",
     "retry",
-    "rasterio>=1.0.28",
+    "rasterio>=1.0.28,<1.2.7",
     "shapely",
     "tilematrix>=0.20",
     "tqdm",
 ]
 req_contours = ["matplotlib"]
+req_dask = ["dask", "distributed"]
 req_geobuf = ["geobuf"]
 req_http = ["fsspec[http]", "aiohttp", "requests"]
 req_s3 = ["boto3", "fsspec[s3]", "s3fs>=0.5.1"]
 req_serve = ["flask", "werkzeug>=0.15"]
+req_stac = ["pystac"]
 req_vrt = ["lxml"]
-req_complete = req_contours + req_geobuf + req_http + req_s3 + req_serve + req_vrt
+req_complete = (
+    req_contours
+    + req_dask
+    + req_geobuf
+    + req_http
+    + req_s3
+    + req_serve
+    + req_stac
+    + req_vrt
+)
 
 setup(
     name="mapchete",
@@ -66,8 +78,9 @@
             "formats=mapchete.cli.default.formats:formats",
             "index=mapchete.cli.default.index:index",
             "processes=mapchete.cli.default.processes:processes",
-            "rm=mapchete.cli.default.rm:rm_",
+            "rm=mapchete.cli.default.rm:rm",
             "serve=mapchete.cli.default.serve:serve",
+            "stac=mapchete.cli.default.stac:stac",
         ],
         "mapchete.formats.drivers": [
             "flatgeobuf=mapchete.formats.default.flatgeobuf",
@@ -94,9 +107,11 @@
     extras_require={
         "complete": req_complete,
         "contours": req_contours,
+        "dask": req_dask,
         "geobuf": req_geobuf,
         "s3": req_s3,
         "serve": req_serve,
+        "stac": req_stac,
         "vrt": req_vrt,
     },
     tests_require=[


Test patch:
diff --git a/test/conftest.py b/test/conftest.py
index 80860dfd..5f351412 100644
--- a/test/conftest.py
+++ b/test/conftest.py
@@ -73,6 +73,12 @@ def wkt_geom():
     return "Polygon ((2.8125 11.25, 2.8125 14.0625, 0 14.0625, 0 11.25, 2.8125 11.25))"
 
 
+@pytest.fixture
+def wkt_geom_tl():
+    """Example WKT geometry."""
+    return "Polygon ((-176.04949 85.59671, -174.57652 73.86651, -159.98073 74.58961, -161.74829 83.05249, -176.04949 85.59671))"
+
+
 # example files
 @pytest.fixture
 def http_raster():
@@ -355,6 +361,13 @@ def example_mapchete():
     return ExampleConfig(path=path, dict=_dict_from_mapchete(path))
 
 
+@pytest.fixture
+def example_custom_process_mapchete():
+    """Fixture for example.mapchete."""
+    path = os.path.join(TESTDATA_DIR, "example_custom_process.mapchete")
+    return ExampleConfig(path=path, dict=_dict_from_mapchete(path))
+
+
 @pytest.fixture
 def zoom_mapchete():
     """Fixture for zoom.mapchete."""
@@ -383,6 +396,13 @@ def cleantopo_br():
     return ExampleConfig(path=path, dict=_dict_from_mapchete(path))
 
 
+@pytest.fixture
+def cleantopo_br_metatiling_1():
+    """Fixture for cleantopo_br.mapchete."""
+    path = os.path.join(TESTDATA_DIR, "cleantopo_br_metatiling_1.mapchete")
+    return ExampleConfig(path=path, dict=_dict_from_mapchete(path))
+
+
 @pytest.fixture
 def cleantopo_remote():
     """Fixture for cleantopo_remote.mapchete."""
diff --git a/test/requirements.txt b/test/requirements.txt
index be956ca6..6358d18f 100644
--- a/test/requirements.txt
+++ b/test/requirements.txt
@@ -1,6 +1,7 @@
 coverage>=5.1
 coveralls>=2.2.0
 flake8>=3.8.4
+pystac
 pytest>=5.4.1
 pytest-cov>=2.8.1
 pytest-flask>=1.0.0
diff --git a/test/test_cli.py b/test/test_cli.py
index 88697cb7..9c2ee13a 100644
--- a/test/test_cli.py
+++ b/test/test_cli.py
@@ -16,7 +16,7 @@
 
 import mapchete
 from mapchete.cli.main import main as mapchete_cli
-from mapchete.cli import utils
+from mapchete.cli import options
 from mapchete.errors import MapcheteProcessOutputError
 
 
@@ -41,10 +41,15 @@ def version_is_greater_equal(a, b):
 
 
 def run_cli(args, expected_exit_code=0, output_contains=None, raise_exc=True):
-    result = CliRunner(env=dict(MAPCHETE_TEST="TRUE")).invoke(mapchete_cli, args)
+    result = CliRunner(env=dict(MAPCHETE_TEST="TRUE"), mix_stderr=True).invoke(
+        mapchete_cli, args
+    )
     if output_contains:
-        assert output_contains in result.output
+        assert output_contains in result.output or output_contains in str(
+            result.exception
+        )
     if raise_exc and result.exception:
+        print(result.output)
         raise result.exception
     assert result.exit_code == expected_exit_code
 
@@ -69,17 +74,7 @@ def test_main(mp_tmpdir):
     )
 
 
-def test_missing_input_file(mp_tmpdir):
-    """Check if IOError is raised if input-file is invalid."""
-    run_cli(
-        ["execute", "process.mapchete", "--input-file", "invalid.tif"],
-        expected_exit_code=2,
-        output_contains="Path 'process.mapchete' does not exist.",
-        raise_exc=False,
-    )
-
-
-def test_create_and_execute(mp_tmpdir, cleantopo_br_tif):
+def test_create(mp_tmpdir, cleantopo_br_tif):
     """Run mapchete create and execute."""
     temp_mapchete = os.path.join(mp_tmpdir, "temp.mapchete")
     temp_process = os.path.join(mp_tmpdir, "temp.py")
@@ -102,22 +97,6 @@ def test_create_and_execute(mp_tmpdir, cleantopo_br_tif):
         config["output"].update(bands=1, dtype="uint8", path=mp_tmpdir)
     with open(temp_mapchete, "w") as config_file:
         config_file.write(yaml.dump(config, default_flow_style=False))
-    # run process for single tile, this creates an empty output error
-    with pytest.raises(MapcheteProcessOutputError):
-        run_cli(
-            [
-                "execute",
-                temp_mapchete,
-                "--tile",
-                "6",
-                "62",
-                "124",
-                "--input-file",
-                cleantopo_br_tif,
-                "-d",
-            ],
-            expected_exit_code=-1,
-        )
 
 
 def test_create_existing(mp_tmpdir):
@@ -140,50 +119,72 @@ def test_create_existing(mp_tmpdir):
         run_cli(args, expected_exit_code=-1)
 
 
-def test_execute_multiprocessing(mp_tmpdir, cleantopo_br, cleantopo_br_tif):
-    """Run mapchete execute with multiple workers."""
-    temp_mapchete = os.path.join(mp_tmpdir, "temp.mapchete")
-    temp_process = os.path.join(mp_tmpdir, "temp.py")
-    out_format = "GTiff"
-    # create from template
+def test_execute_concurrent_processes(mp_tmpdir, cleantopo_br_metatiling_1):
+    # """Run mapchete execute with multiple workers."""
     run_cli(
         [
-            "create",
-            temp_mapchete,
-            temp_process,
-            out_format,
-            "--pyramid-type",
-            "geodetic",
+            "execute",
+            cleantopo_br_metatiling_1.path,
+            "--zoom",
+            "5",
+            "--workers",
+            "2",
+            "-d",
+            "--concurrency",
+            "processes",
         ]
     )
-    # edit configuration
-    with open(temp_mapchete, "r") as config_file:
-        config = yaml.safe_load(config_file)
-        config["output"].update(bands=1, dtype="uint8", path=mp_tmpdir)
-    with open(temp_mapchete, "w") as config_file:
-        config_file.write(yaml.dump(config, default_flow_style=False))
-    # run process with multiprocessing
-    with pytest.raises(MapcheteProcessOutputError):
-        run_cli(
-            [
-                "execute",
-                temp_mapchete,
-                "--zoom",
-                "5",
-                "--input-file",
-                cleantopo_br_tif,
-                "-m",
-                "2",
-                "-d",
-            ]
-        )
-    # run example process with multiprocessing
-    run_cli(["execute", cleantopo_br.path, "--zoom", "5", "-m", "2", "-d"])
+
+
+def test_execute_concurrent_threads(mp_tmpdir, cleantopo_br_metatiling_1):
+    """Run mapchete execute with multiple workers."""
+    run_cli(
+        [
+            "execute",
+            cleantopo_br_metatiling_1.path,
+            "--zoom",
+            "5",
+            "--workers",
+            "2",
+            "-d",
+            "--concurrency",
+            "threads",
+        ],
+    )
+
+
+def test_execute_concurrent_dask(mp_tmpdir, cleantopo_br_metatiling_1):
+    """Run mapchete execute with multiple workers."""
+    run_cli(
+        [
+            "execute",
+            cleantopo_br_metatiling_1.path,
+            "--zoom",
+            "5",
+            "--workers",
+            "2",
+            "-d",
+            "--concurrency",
+            "dask",
+        ],
+    )
 
 
 def test_execute_debug(mp_tmpdir, example_mapchete):
     """Using debug output."""
-    run_cli(["execute", example_mapchete.path, "-t", "10", "500", "1040", "--debug"])
+    run_cli(
+        [
+            "execute",
+            example_mapchete.path,
+            "-t",
+            "10",
+            "500",
+            "1040",
+            "--debug",
+            "--concurrency",
+            "none",
+        ]
+    )
 
 
 def test_execute_vrt(mp_tmpdir, cleantopo_br):
@@ -196,7 +197,17 @@ def test_execute_vrt(mp_tmpdir, cleantopo_br):
 
     # run again, this time with custom output directory
     run_cli(
-        ["execute", cleantopo_br.path, "-z", "5", "--vrt", "--idx-out-dir", mp_tmpdir]
+        [
+            "execute",
+            cleantopo_br.path,
+            "-z",
+            "5",
+            "--vrt",
+            "--idx-out-dir",
+            mp_tmpdir,
+            "--concurrency",
+            "none",
+        ]
     )
     with mapchete.open(cleantopo_br.dict) as mp:
         vrt_path = os.path.join(mp_tmpdir, "5.vrt")
@@ -204,15 +215,51 @@ def test_execute_vrt(mp_tmpdir, cleantopo_br):
             assert src.read().any()
 
     # run with single tile
-    run_cli(["execute", cleantopo_br.path, "-t", "5", "3", "7", "--vrt"])
+    run_cli(
+        [
+            "execute",
+            cleantopo_br.path,
+            "-t",
+            "5",
+            "3",
+            "7",
+            "--vrt",
+            "--concurrency",
+            "none",
+        ]
+    )
 
     # no new entries
-    run_cli(["execute", cleantopo_br.path, "-t", "5", "0", "0", "--vrt"])
+    run_cli(
+        [
+            "execute",
+            cleantopo_br.path,
+            "-t",
+            "5",
+            "0",
+            "0",
+            "--vrt",
+            "--concurrency",
+            "none",
+        ]
+    )
 
 
 def test_execute_verbose(mp_tmpdir, example_mapchete):
     """Using verbose output."""
-    run_cli(["execute", example_mapchete.path, "-t", "10", "500", "1040", "--verbose"])
+    run_cli(
+        [
+            "execute",
+            example_mapchete.path,
+            "-t",
+            "10",
+            "500",
+            "1040",
+            "--verbose",
+            "--concurrency",
+            "none",
+        ]
+    )
 
 
 def test_execute_logfile(mp_tmpdir, example_mapchete):
@@ -228,6 +275,8 @@ def test_execute_logfile(mp_tmpdir, example_mapchete):
             "1040",
             "--logfile",
             logfile,
+            "--concurrency",
+            "none",
         ]
     )
     assert os.path.isfile(logfile)
@@ -235,14 +284,11 @@ def test_execute_logfile(mp_tmpdir, example_mapchete):
         assert "DEBUG" in log.read()
 
 
-def test_execute_wkt_bounds(mp_tmpdir, example_mapchete, wkt_geom):
-    """Using bounds from WKT."""
-    run_cli(["execute", example_mapchete.path, "--wkt-geometry", wkt_geom])
-
-
 def test_execute_wkt_area(mp_tmpdir, example_mapchete, wkt_geom):
     """Using area from WKT."""
-    run_cli(["execute", example_mapchete.path, "--area", wkt_geom])
+    run_cli(
+        ["execute", example_mapchete.path, "--area", wkt_geom, "--concurrency", "none"]
+    )
 
 
 def test_execute_point(mp_tmpdir, example_mapchete, wkt_geom):
@@ -255,6 +301,8 @@ def test_execute_point(mp_tmpdir, example_mapchete, wkt_geom):
             "--point",
             str(g.centroid.x),
             str(g.centroid.y),
+            "--concurrency",
+            "none",
         ]
     )
 
@@ -274,7 +322,17 @@ def test_formats(capfd):
 
 def test_convert_geodetic(cleantopo_br_tif, mp_tmpdir):
     """Automatic geodetic tile pyramid creation of raster files."""
-    run_cli(["convert", cleantopo_br_tif, mp_tmpdir, "--output-pyramid", "geodetic"])
+    run_cli(
+        [
+            "convert",
+            cleantopo_br_tif,
+            mp_tmpdir,
+            "--output-pyramid",
+            "geodetic",
+            "--concurrency",
+            "none",
+        ]
+    )
     for zoom, row, col in [(4, 15, 31), (3, 7, 15), (2, 3, 7), (1, 1, 3)]:
         out_file = os.path.join(*[mp_tmpdir, str(zoom), str(row), str(col) + ".tif"])
         with rasterio.open(out_file, "r") as src:
@@ -286,7 +344,17 @@ def test_convert_geodetic(cleantopo_br_tif, mp_tmpdir):
 
 def test_convert_mercator(cleantopo_br_tif, mp_tmpdir):
     """Automatic mercator tile pyramid creation of raster files."""
-    run_cli(["convert", cleantopo_br_tif, mp_tmpdir, "--output-pyramid", "mercator"])
+    run_cli(
+        [
+            "convert",
+            cleantopo_br_tif,
+            mp_tmpdir,
+            "--output-pyramid",
+            "mercator",
+            "--concurrency",
+            "none",
+        ]
+    )
     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"])
         with rasterio.open(out_file, "r") as src:
@@ -307,6 +375,8 @@ def test_convert_png(cleantopo_br_tif, mp_tmpdir):
             "mercator",
             "--output-format",
             "PNG",
+            "--concurrency",
+            "none",
         ]
     )
     for zoom, row, col in [(4, 15, 15), (3, 7, 7)]:
@@ -334,6 +404,8 @@ def test_convert_bidx(cleantopo_br_tif, mp_tmpdir):
             "3",
             "--bidx",
             "1",
+            "--concurrency",
+            "none",
         ]
     )
     with rasterio.open(single_gtiff, "r") as src:
@@ -356,6 +428,8 @@ def test_convert_single_gtiff(cleantopo_br_tif, mp_tmpdir):
             "geodetic",
             "-z",
             "3",
+            "--concurrency",
+            "none",
         ]
     )
     with rasterio.open(single_gtiff, "r") as src:
@@ -379,6 +453,8 @@ def test_convert_single_gtiff_cog(cleantopo_br_tif, mp_tmpdir):
             "-z",
             "3",
             "--cog",
+            "--concurrency",
+            "none",
         ]
     )
     with rasterio.open(single_gtiff, "r") as src:
@@ -404,6 +480,10 @@ def test_convert_single_gtiff_overviews(cleantopo_br_tif, mp_tmpdir):
             "--overviews",
             "--overviews-resampling-method",
             "bilinear",
+            "--multi",
+            "1",
+            "--concurrency",
+            "none",
         ]
     )
     with rasterio.open(single_gtiff, "r") as src:
@@ -426,6 +506,8 @@ def test_convert_remote_single_gtiff(http_raster, mp_tmpdir):
             "geodetic",
             "-z",
             "1",
+            "--concurrency",
+            "none",
         ]
     )
     with rasterio.open(single_gtiff, "r") as src:
@@ -446,6 +528,8 @@ def test_convert_dtype(cleantopo_br_tif, mp_tmpdir):
             "mercator",
             "--output-dtype",
             "uint8",
+            "--concurrency",
+            "none",
         ]
     )
     for zoom, row, col in [(4, 15, 15), (3, 7, 7)]:
@@ -470,6 +554,8 @@ def test_convert_scale_ratio(cleantopo_br_tif, mp_tmpdir):
             "uint8",
             "--scale-ratio",
             "0.003",
+            "--concurrency",
+            "none",
         ]
     )
     for zoom, row, col in [(4, 15, 15), (3, 7, 7)]:
@@ -495,6 +581,8 @@ def test_convert_scale_offset(cleantopo_br_tif, mp_tmpdir):
             "uint8",
             "--scale-offset",
             "1",
+            "--concurrency",
+            "none",
         ]
     )
     for zoom, row, col in [(4, 15, 15), (3, 7, 7)]:
@@ -518,7 +606,9 @@ def test_convert_clip(cleantopo_br_tif, mp_tmpdir, landpoly):
             "geodetic",
             "--clip-geometry",
             landpoly,
-            "-d",
+            "-v",
+            "--concurrency",
+            "none",
         ],
         output_contains="Process area is empty",
     )
@@ -535,6 +625,8 @@ def test_convert_zoom(cleantopo_br_tif, mp_tmpdir):
             "mercator",
             "-z",
             "3",
+            "--concurrency",
+            "none",
         ]
     )
     for zoom, row, col in [(4, 15, 15), (2, 3, 0)]:
@@ -553,6 +645,8 @@ def test_convert_zoom_minmax(cleantopo_br_tif, mp_tmpdir):
             "mercator",
             "-z",
             "3,4",
+            "--concurrency",
+            "none",
         ]
     )
     for zoom, row, col in [(2, 3, 0)]:
@@ -571,6 +665,8 @@ def test_convert_zoom_maxmin(cleantopo_br_tif, mp_tmpdir):
             "mercator",
             "-z",
             "4,3",
+            "--concurrency",
+            "none",
         ]
     )
     for zoom, row, col in [(2, 3, 0)]:
@@ -592,6 +688,8 @@ def test_convert_mapchete(cleantopo_br, mp_tmpdir):
             "--output-metatiling",
             "1",
             "-d",
+            "--concurrency",
+            "none",
         ]
     )
     for zoom, row, col in [(4, 15, 31), (3, 7, 15), (2, 3, 7), (1, 1, 3)]:
@@ -621,6 +719,8 @@ def test_convert_tiledir(cleantopo_br, mp_tmpdir):
             "--zoom",
             "1,4",
             "-d",
+            "--concurrency",
+            "none",
         ]
     )
     for zoom, row, col in [(4, 15, 31), (3, 7, 15), (2, 3, 7), (1, 1, 3)]:
@@ -642,6 +742,8 @@ def test_convert_geojson(landpoly, mp_tmpdir):
             "geodetic",
             "--zoom",
             "4",
+            "--concurrency",
+            "none",
         ]
     )
     for (zoom, row, col), control in zip([(4, 0, 7), (4, 1, 7)], [9, 32]):
@@ -668,6 +770,8 @@ def test_convert_geobuf(landpoly, mp_tmpdir):
             "4",
             "--output-format",
             "Geobuf",
+            "--concurrency",
+            "none",
         ]
     )
     for (zoom, row, col), control in zip([(4, 0, 7), (4, 1, 7)], [9, 32]):
@@ -698,6 +802,8 @@ def test_convert_geobuf(landpoly, mp_tmpdir):
             "4",
             "--output-format",
             "GeoJSON",
+            "--concurrency",
+            "none",
         ]
     )
     for (zoom, row, col), control in zip([(4, 0, 7), (4, 1, 7)], [9, zoom9_control]):
@@ -724,6 +830,8 @@ def test_convert_geobuf_multipolygon(landpoly, mp_tmpdir):
             "Geobuf",
             "--output-geometry-type",
             "MultiPolygon",
+            "--concurrency",
+            "none",
         ]
     )
     for (zoom, row, col), control in zip([(4, 0, 7), (4, 1, 7)], [7, 30]):
@@ -740,11 +848,36 @@ def test_convert_geobuf_multipolygon(landpoly, mp_tmpdir):
             assert multipolygons
 
 
+def test_convert_vrt(cleantopo_br_tif, mp_tmpdir):
+    """Automatic geodetic tile pyramid creation of raster files."""
+    run_cli(
+        [
+            "convert",
+            cleantopo_br_tif,
+            mp_tmpdir,
+            "--output-pyramid",
+            "geodetic",
+            "--vrt",
+            "--zoom",
+            "1,4",
+            "--concurrency",
+            "none",
+        ]
+    )
+    for zoom in [4, 3, 2, 1]:
+        out_file = os.path.join(*[mp_tmpdir, str(zoom) + ".vrt"])
+        with rasterio.open(out_file, "r") as src:
+            assert src.meta["driver"] == "VRT"
+            assert src.meta["dtype"] == "uint16"
+            data = src.read(masked=True)
+            assert data.mask.any()
+
+
 def test_convert_errors(s2_band_jp2, mp_tmpdir, s2_band, cleantopo_br, landpoly):
     # output format required
     run_cli(
         ["convert", s2_band_jp2, mp_tmpdir, "--output-pyramid", "geodetic"],
-        expected_exit_code=2,
+        expected_exit_code=1,
         output_contains="Output format required.",
         raise_exc=False,
     )
@@ -752,7 +885,7 @@ def test_convert_errors(s2_band_jp2, mp_tmpdir, s2_band, cleantopo_br, landpoly)
     # output pyramid reqired
     run_cli(
         ["convert", s2_band, mp_tmpdir],
-        expected_exit_code=2,
+        expected_exit_code=1,
         output_contains="Output pyramid required.",
         raise_exc=False,
     )
@@ -772,8 +905,10 @@ def test_convert_errors(s2_band_jp2, mp_tmpdir, s2_band, cleantopo_br, landpoly)
             mp_tmpdir,
             "--output-pyramid",
             "geodetic",
+            "--concurrency",
+            "none",
         ],
-        expected_exit_code=2,
+        expected_exit_code=1,
         output_contains="Zoom levels required.",
         raise_exc=False,
     )
@@ -790,8 +925,10 @@ def test_convert_errors(s2_band_jp2, mp_tmpdir, s2_band, cleantopo_br, landpoly)
             "5",
             "--output-format",
             "GeoJSON",
+            "--concurrency",
+            "none",
         ],
-        expected_exit_code=2,
+        expected_exit_code=1,
         output_contains=(
             "Output format type (vector) is incompatible with input format (raster)."
         ),
@@ -808,8 +945,10 @@ def test_convert_errors(s2_band_jp2, mp_tmpdir, s2_band, cleantopo_br, landpoly)
             "geodetic",
             "--zoom",
             "5",
+            "--concurrency",
+            "none",
         ],
-        expected_exit_code=2,
+        expected_exit_code=1,
         output_contains=("Could not determine output from extension"),
         raise_exc=False,
     )
@@ -838,7 +977,6 @@ def test_serve_cli_params(cleantopo_br, mp_tmpdir):
         ["serve", cleantopo_br.path, "--overwrite"],
         ["serve", cleantopo_br.path, "--readonly"],
         ["serve", cleantopo_br.path, "--memory"],
-        ["serve", cleantopo_br.path, "--input-file", cleantopo_br.path],
     ]:
         run_cli(args)
 
@@ -882,7 +1020,9 @@ def test_serve(client, mp_tmpdir):
 
 def test_index_geojson(mp_tmpdir, cleantopo_br):
     # execute process at zoom 3
-    run_cli(["execute", cleantopo_br.path, "-z", "3", "--debug"])
+    run_cli(
+        ["execute", cleantopo_br.path, "-z", "3", "--debug", "--concurrency", "none"]
+    )
 
     # generate index for zoom 3
     run_cli(["index", cleantopo_br.path, "-z", "3", "--geojson", "--debug"])
@@ -898,7 +1038,9 @@ def test_index_geojson(mp_tmpdir, cleantopo_br):
 
 def test_index_geojson_fieldname(mp_tmpdir, cleantopo_br):
     # execute process at zoom 3
-    run_cli(["execute", cleantopo_br.path, "-z", "3", "--debug"])
+    run_cli(
+        ["execute", cleantopo_br.path, "-z", "3", "--debug", "--concurrency", "none"]
+    )
 
     # index and rename "location" to "new_fieldname"
     run_cli(
@@ -924,7 +1066,9 @@ def test_index_geojson_fieldname(mp_tmpdir, cleantopo_br):
 
 def test_index_geojson_basepath(mp_tmpdir, cleantopo_br):
     # execute process at zoom 3
-    run_cli(["execute", cleantopo_br.path, "-z", "3", "--debug"])
+    run_cli(
+        ["execute", cleantopo_br.path, "-z", "3", "--debug", "--concurrency", "none"]
+    )
 
     basepath = "http://localhost"
     # index and rename "location" to "new_fieldname"
@@ -979,7 +1123,19 @@ def test_index_geojson_for_gdal(mp_tmpdir, cleantopo_br):
 
 def test_index_geojson_tile(mp_tmpdir, cleantopo_tl):
     # execute process for single tile
-    run_cli(["execute", cleantopo_tl.path, "-t", "3", "0", "0", "--debug"])
+    run_cli(
+        [
+            "execute",
+            cleantopo_tl.path,
+            "-t",
+            "3",
+            "0",
+            "0",
+            "--debug",
+            "--concurrency",
+            "none",
+        ]
+    )
     # generate index
     run_cli(["index", cleantopo_tl.path, "-t", "3", "0", "0", "--geojson", "--debug"])
     with mapchete.open(cleantopo_tl.dict) as mp:
@@ -990,25 +1146,20 @@ def test_index_geojson_tile(mp_tmpdir, cleantopo_tl):
         assert len(list(src)) == 1
 
 
-def test_index_geojson_wkt_geom(mp_tmpdir, cleantopo_br, wkt_geom):
+def test_index_geojson_wkt_area(mp_tmpdir, cleantopo_br, wkt_geom):
     # execute process at zoom 3
-    run_cli(["execute", cleantopo_br.path, "--debug", "--wkt-geometry", wkt_geom])
-
-    # generate index for zoom 3
     run_cli(
-        ["index", cleantopo_br.path, "--geojson", "--debug", "--wkt-geometry", wkt_geom]
+        [
+            "execute",
+            cleantopo_br.path,
+            "--debug",
+            "--area",
+            wkt_geom,
+            "--concurrency",
+            "none",
+        ]
     )
 
-    with mapchete.open(cleantopo_br.dict) as mp:
-        files = os.listdir(mp.config.output.path)
-        assert len(files) == 7
-        assert "3.geojson" in files
-
-
-def test_index_geojson_wkt_area(mp_tmpdir, cleantopo_br, wkt_geom):
-    # execute process at zoom 3
-    run_cli(["execute", cleantopo_br.path, "--debug", "--area", wkt_geom])
-
     # generate index for zoom 3
     run_cli(["index", cleantopo_br.path, "--geojson", "--debug", "--area", wkt_geom])
 
@@ -1020,7 +1171,9 @@ def test_index_geojson_wkt_area(mp_tmpdir, cleantopo_br, wkt_geom):
 
 def test_index_gpkg(mp_tmpdir, cleantopo_br):
     # execute process
-    run_cli(["execute", cleantopo_br.path, "-z", "5", "--debug"])
+    run_cli(
+        ["execute", cleantopo_br.path, "-z", "5", "--debug", "--concurrency", "none"]
+    )
 
     # generate index
     run_cli(["index", cleantopo_br.path, "-z", "5", "--gpkg", "--debug"])
@@ -1045,7 +1198,9 @@ def test_index_gpkg(mp_tmpdir, cleantopo_br):
 
 def test_index_shp(mp_tmpdir, cleantopo_br):
     # execute process
-    run_cli(["execute", cleantopo_br.path, "-z", "5", "--debug"])
+    run_cli(
+        ["execute", cleantopo_br.path, "-z", "5", "--debug", "--concurrency", "none"]
+    )
 
     # generate index
     run_cli(["index", cleantopo_br.path, "-z", "5", "--shp", "--debug"])
@@ -1070,7 +1225,9 @@ def test_index_shp(mp_tmpdir, cleantopo_br):
 
 def test_index_text(cleantopo_br):
     # execute process
-    run_cli(["execute", cleantopo_br.path, "-z", "5", "--debug"])
+    run_cli(
+        ["execute", cleantopo_br.path, "-z", "5", "--debug", "--concurrency", "none"]
+    )
 
     # generate index
     run_cli(["index", cleantopo_br.path, "-z", "5", "--txt", "--debug"])
@@ -1096,7 +1253,7 @@ def test_index_text(cleantopo_br):
 
 
 def test_index_errors(mp_tmpdir, cleantopo_br):
-    with pytest.raises(SystemExit):
+    with pytest.raises(ValueError):
         run_cli(["index", cleantopo_br.path, "-z", "5", "--debug"])
 
     with pytest.raises(SystemExit):
@@ -1110,13 +1267,13 @@ def test_processes():
 
 def test_callback_errors(cleantopo_tl):
     run_cli(
-        ["execute", cleantopo_tl.path, "--zoom", "4,5,7"],
+        ["execute", cleantopo_tl.path, "--zoom", "4,5,7", "--concurrency", "none"],
         expected_exit_code=2,
         raise_exc=False,
         output_contains="zooms can be maximum two items",
     )
     run_cli(
-        ["execute", cleantopo_tl.path, "--zoom", "invalid"],
+        ["execute", cleantopo_tl.path, "--zoom", "invalid", "--concurrency", "none"],
         expected_exit_code=2,
         raise_exc=False,
         output_contains="zoom levels must be integer values",
@@ -1137,6 +1294,8 @@ def test_cp(mp_tmpdir, cleantopo_br, wkt_geom):
             "-90",
             "180",
             "-80.18582802550002",
+            "--concurrency",
+            "none",
         ]
     )
     out_path = os.path.join(TESTDATA_DIR, cleantopo_br.dict["output"]["path"])
@@ -1154,6 +1313,8 @@ def test_cp(mp_tmpdir, cleantopo_br, wkt_geom):
             "-90",
             "180",
             "-80.18582802550002",
+            "--concurrency",
+            "none",
         ]
     )
     # copy all tiles
@@ -1164,14 +1325,36 @@ def test_cp(mp_tmpdir, cleantopo_br, wkt_geom):
             os.path.join(mp_tmpdir, "all"),
             "-z",
             "5",
+            "--concurrency",
+            "none",
         ]
     )
     # copy tiles and subset by area
     run_cli(
-        ["cp", out_path, os.path.join(mp_tmpdir, "all"), "-z", "5", "--area", wkt_geom]
+        [
+            "cp",
+            out_path,
+            os.path.join(mp_tmpdir, "all"),
+            "-z",
+            "5",
+            "--area",
+            wkt_geom,
+            "--concurrency",
+            "none",
+        ]
+    )
+    # copy local tiles wit using threads
+    run_cli(
+        [
+            "cp",
+            out_path,
+            os.path.join(mp_tmpdir, "all"),
+            "-z",
+            "5",
+            "--concurrency",
+            "threads",
+        ]
     )
-    # copy local tiles without using threads
-    run_cli(["cp", out_path, os.path.join(mp_tmpdir, "all"), "-z", "5", "--multi", "1"])
 
 
 def test_cp_http(mp_tmpdir, http_tiledir):
@@ -1188,6 +1371,8 @@ def test_cp_http(mp_tmpdir, http_tiledir):
             "1.0",
             "4.0",
             "2.0",
+            "--concurrency",
+            "none",
         ]
     )
 
@@ -1204,6 +1389,8 @@ def test_rm(mp_tmpdir, cleantopo_br):
             "-90",
             "180",
             "-80.18582802550002",
+            "--concurrency",
+            "none",
         ]
     )
     out_path = os.path.join(TESTDATA_DIR, cleantopo_br.dict["output"]["path"])
@@ -1249,7 +1436,7 @@ def test_rm_storage_option_errors(cleantopo_br):
 
 
 def test_fs_opt_extractor():
-    kwargs = utils._cb_key_val(
+    kwargs = options._cb_key_val(
         None,
         None,
         [
@@ -1274,3 +1461,22 @@ def test_fs_opt_extractor():
     assert kwargs["bool4"] is False
     assert kwargs["none"] is None
     assert kwargs["none2"] is None
+
+
+def test_stac_mapchete_file(cleantopo_br):
+    run_cli(["stac", "create-item", cleantopo_br.path, "-z", "5", "--force"])
+
+
+def test_stac_tiledir(http_tiledir, mp_tmpdir):
+    run_cli(
+        [
+            "stac",
+            "create-item",
+            http_tiledir,
+            "-z",
+            "5",
+            "--force",
+            "--item-path",
+            f"{mp_tmpdir}/stac_example.json",
+        ]
+    )
diff --git a/test/test_commands.py b/test/test_commands.py
new file mode 100644
index 00000000..86460357
--- /dev/null
+++ b/test/test_commands.py
@@ -0,0 +1,659 @@
+import fiona
+import geobuf
+import os
+import pytest
+import rasterio
+from rio_cogeo.cogeo import cog_validate
+from shapely.geometry import box, shape
+from tilematrix import TilePyramid
+import warnings
+
+import mapchete
+from mapchete.commands import convert, cp, execute, index, rm
+
+
+SCRIPTDIR = os.path.dirname(os.path.realpath(__file__))
+TESTDATA_DIR = os.path.join(SCRIPTDIR, "testdata")
+
+
+def test_cp(mp_tmpdir, cleantopo_br, wkt_geom):
+    # generate TileDirectory
+    with mapchete.open(
+        cleantopo_br.path, bounds=[169.19251592399996, -90, 180, -80.18582802550002]
+    ) as mp:
+        mp.batch_process(zoom=5)
+    out_path = os.path.join(TESTDATA_DIR, cleantopo_br.dict["output"]["path"])
+
+    # copy tiles and subset by bounds
+    tiles = cp(
+        out_path,
+        os.path.join(mp_tmpdir, "bounds"),
+        zoom=5,
+        bounds=[169.19251592399996, -90, 180, -80.18582802550002],
+    )
+    assert len(tiles)
+
+    # copy all tiles
+    tiles = cp(
+        out_path,
+        os.path.join(mp_tmpdir, "all"),
+        zoom=5,
+    )
+    assert len(tiles)
+
+    # copy tiles and subset by area
+    tiles = cp(out_path, os.path.join(mp_tmpdir, "area"), zoom=5, area=wkt_geom)
+    assert len(tiles)
+
+    # copy local tiles without using threads
+    tiles = cp(out_path, os.path.join(mp_tmpdir, "nothreads"), zoom=5, multi=1)
+    assert len(tiles)
+
+
+def test_cp_http(mp_tmpdir, http_tiledir):
+    # copy tiles and subset by bounds
+    tiles = cp(
+        http_tiledir, os.path.join(mp_tmpdir, "http"), zoom=1, bounds=[3, 1, 4, 2]
+    )
+    assert len(tiles)
+
+
+def test_rm(mp_tmpdir, cleantopo_br):
+    # generate TileDirectory
+    with mapchete.open(
+        cleantopo_br.path, bounds=[169.19251592399996, -90, 180, -80.18582802550002]
+    ) as mp:
+        mp.batch_process(zoom=5)
+    out_path = os.path.join(TESTDATA_DIR, cleantopo_br.dict["output"]["path"])
+
+    # remove tiles
+    tiles = rm(out_path, zoom=5)
+    assert len(tiles) > 0
+
+    # remove tiles but this time they should already have been removed
+    tiles = rm(out_path, zoom=5)
+    assert len(tiles) == 0
+
+
+def test_execute(mp_tmpdir, cleantopo_br, cleantopo_br_tif):
+    zoom = 5
+    config = cleantopo_br.dict
+    config["pyramid"].update(metatiling=1)
+    tp = TilePyramid("geodetic")
+    tiles = list(tp.tiles_from_bounds(rasterio.open(cleantopo_br_tif).bounds, zoom))
+    job = execute(config, zoom=zoom)
+    for t in job:
+        assert t
+    assert len(tiles) == len(job)
+    with mapchete.open(config) as 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, cleantopo_br_tif):
+    zoom = 5
+    config = cleantopo_br.dict
+    config["pyramid"].update(metatiling=1)
+    job = execute(config, zoom=zoom, as_iterator=True)
+    for i, t in enumerate(job):
+        job.cancel()
+        break
+    assert i == 0
+    assert job.status == "cancelled"
+
+
+def test_execute_tile(mp_tmpdir, cleantopo_br):
+    tile = (5, 30, 63)
+
+    config = cleantopo_br.dict
+    config["pyramid"].update(metatiling=1)
+    job = execute(config, tile=tile)
+
+    assert len(job) == 1
+
+    with mapchete.open(config) as mp:
+        with rasterio.open(
+            mp.config.output.get_path(mp.config.output_pyramid.tile(*tile))
+        ) as src:
+            assert not src.read(masked=True).mask.all()
+
+
+def test_execute_point(mp_tmpdir, example_mapchete, dummy2_tif):
+    """Using bounds from WKT."""
+    with rasterio.open(dummy2_tif) as src:
+        g = box(*src.bounds)
+    job = execute(example_mapchete.path, point=[g.centroid.x, g.centroid.y], zoom=10)
+    assert len(job) == 1
+
+
+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")
+    assert len(job)
+    for zoom, row, col in [(4, 15, 31), (3, 7, 15), (2, 3, 7), (1, 1, 3)]:
+        out_file = os.path.join(*[mp_tmpdir, str(zoom), str(row), str(col) + ".tif"])
+        with rasterio.open(out_file, "r") as src:
+            assert src.meta["driver"] == "GTiff"
+            assert src.meta["dtype"] == "uint16"
+            data = src.read(masked=True)
+            assert data.mask.any()
+
+
+def test_convert_mercator(cleantopo_br_tif, mp_tmpdir):
+    """Automatic mercator tile pyramid creation of raster files."""
+    job = convert(cleantopo_br_tif, mp_tmpdir, output_pyramid="mercator")
+    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"])
+        with rasterio.open(out_file, "r") as src:
+            assert src.meta["driver"] == "GTiff"
+            assert src.meta["dtype"] == "uint16"
+            data = src.read(masked=True)
+            assert data.mask.any()
+
+
+def test_convert_png(cleantopo_br_tif, mp_tmpdir):
+    """Automatic PNG tile pyramid creation of raster files."""
+    job = convert(
+        cleantopo_br_tif, mp_tmpdir, output_pyramid="mercator", output_format="PNG"
+    )
+    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) + ".png"])
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            with rasterio.open(out_file, "r") as src:
+                assert src.meta["driver"] == "PNG"
+                assert src.meta["dtype"] == "uint8"
+                data = src.read(masked=True)
+                assert data.mask.any()
+
+
+def test_convert_bidx(cleantopo_br_tif, mp_tmpdir):
+    """Automatic geodetic tile pyramid creation of raster files."""
+    single_gtiff = os.path.join(mp_tmpdir, "single_out_bidx.tif")
+    job = convert(
+        cleantopo_br_tif, single_gtiff, output_pyramid="geodetic", zoom=3, bidx=1
+    )
+    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 src.overviews(1)
+
+
+def test_convert_single_gtiff(cleantopo_br_tif, mp_tmpdir):
+    """Automatic geodetic tile pyramid creation of raster files."""
+    single_gtiff = os.path.join(mp_tmpdir, "single_out.tif")
+    job = convert(cleantopo_br_tif, single_gtiff, output_pyramid="geodetic", zoom=3)
+    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 src.overviews(1)
+
+
+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
+    )
+    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 cog_validate(single_gtiff, strict=True)
+
+
+def test_convert_single_gtiff_overviews(cleantopo_br_tif, mp_tmpdir):
+    """Automatic geodetic tile pyramid creation of raster files."""
+    single_gtiff = os.path.join(mp_tmpdir, "single_out.tif")
+    job = convert(
+        cleantopo_br_tif,
+        single_gtiff,
+        output_pyramid="geodetic",
+        zoom=7,
+        overviews=True,
+        overviews_resampling_method="bilinear",
+        concurrency=None,
+    )
+    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 src.overviews(1)
+
+
+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")
+    job = convert(
+        http_raster, single_gtiff, output_pyramid="geodetic", zoom=1, concurrency=None
+    )
+    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.any()
+
+
+def test_convert_dtype(cleantopo_br_tif, mp_tmpdir):
+    """Automatic tile pyramid creation using dtype scale."""
+    job = convert(
+        cleantopo_br_tif, mp_tmpdir, output_pyramid="mercator", output_dtype="uint8"
+    )
+    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"])
+        with rasterio.open(out_file, "r") as src:
+            assert src.meta["driver"] == "GTiff"
+            assert src.meta["dtype"] == "uint8"
+            data = src.read(masked=True)
+            assert data.mask.any()
+
+
+def test_convert_scale_ratio(cleantopo_br_tif, mp_tmpdir):
+    """Automatic tile pyramid creation cropping data."""
+    job = convert(
+        cleantopo_br_tif,
+        mp_tmpdir,
+        output_pyramid="mercator",
+        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"])
+        with rasterio.open(out_file, "r") as src:
+            assert src.meta["driver"] == "GTiff"
+            assert src.meta["dtype"] == "uint8"
+            data = src.read(masked=True)
+            assert data.mask.any()
+            assert not data.mask.all()
+
+
+def test_convert_scale_offset(cleantopo_br_tif, mp_tmpdir):
+    """Automatic tile pyramid creation cropping data."""
+    job = convert(
+        cleantopo_br_tif,
+        mp_tmpdir,
+        output_pyramid="mercator",
+        output_dtype="uint8",
+        scale_offset=1,
+    )
+    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"])
+        with rasterio.open(out_file, "r") as src:
+            assert src.meta["driver"] == "GTiff"
+            assert src.meta["dtype"] == "uint8"
+            data = src.read(masked=True)
+            assert data.mask.any()
+            assert not data.mask.all()
+
+
+def test_convert_clip(cleantopo_br_tif, mp_tmpdir, landpoly):
+    """Automatic tile pyramid creation cropping data."""
+    job = convert(
+        cleantopo_br_tif, mp_tmpdir, output_pyramid="geodetic", clip_geometry=landpoly
+    )
+    assert len(job) == 0
+
+
+def test_convert_zoom(cleantopo_br_tif, mp_tmpdir):
+    """Automatic tile pyramid creation using a specific zoom."""
+    job = convert(cleantopo_br_tif, mp_tmpdir, output_pyramid="mercator", zoom=3)
+    assert len(job)
+    for zoom, row, col in [(4, 15, 15), (2, 3, 0)]:
+        out_file = os.path.join(*[mp_tmpdir, str(zoom), str(row), str(col) + ".tif"])
+        assert not os.path.isfile(out_file)
+
+
+def test_convert_zoom_minmax(cleantopo_br_tif, mp_tmpdir):
+    """Automatic tile pyramid creation using min max zoom."""
+    job = convert(cleantopo_br_tif, mp_tmpdir, output_pyramid="mercator", zoom=[3, 4])
+    assert len(job)
+    for zoom, row, col in [(2, 3, 0)]:
+        out_file = os.path.join(*[mp_tmpdir, str(zoom), str(row), str(col) + ".tif"])
+        assert not os.path.isfile(out_file)
+
+
+def test_convert_zoom_maxmin(cleantopo_br_tif, mp_tmpdir):
+    """Automatic tile pyramid creation using max min zoom."""
+    job = convert(cleantopo_br_tif, mp_tmpdir, output_pyramid="mercator", zoom=[4, 3])
+    assert len(job)
+    for zoom, row, col in [(2, 3, 0)]:
+        out_file = os.path.join(*[mp_tmpdir, str(zoom), str(row), str(col) + ".tif"])
+        assert not os.path.isfile(out_file)
+
+
+def test_convert_mapchete(cleantopo_br, mp_tmpdir):
+    # prepare data
+    job = execute(cleantopo_br.path, zoom=[1, 4])
+    assert len(job)
+
+    job = convert(
+        cleantopo_br.path, mp_tmpdir, output_pyramid="geodetic", output_metatiling=1
+    )
+    assert len(job)
+    for zoom, row, col in [(4, 15, 31), (3, 7, 15), (2, 3, 7), (1, 1, 3)]:
+        out_file = os.path.join(*[mp_tmpdir, str(zoom), str(row), str(col) + ".tif"])
+        with rasterio.open(out_file, "r") as src:
+            assert src.meta["driver"] == "GTiff"
+            assert src.meta["dtype"] == "uint16"
+            data = src.read(masked=True)
+            assert data.mask.any()
+
+
+def test_convert_tiledir(cleantopo_br, mp_tmpdir):
+    # prepare data
+    with mapchete.open(cleantopo_br.path) as mp:
+        mp.batch_process(zoom=[1, 4])
+    job = convert(
+        os.path.join(
+            cleantopo_br.dict["config_dir"], cleantopo_br.dict["output"]["path"]
+        ),
+        mp_tmpdir,
+        output_pyramid="geodetic",
+        output_metatiling=1,
+        zoom=[1, 4],
+    )
+    assert len(job)
+    for zoom, row, col in [(4, 15, 31), (3, 7, 15), (2, 3, 7), (1, 1, 3)]:
+        out_file = os.path.join(*[mp_tmpdir, str(zoom), str(row), str(col) + ".tif"])
+        with rasterio.open(out_file, "r") as src:
+            assert src.meta["driver"] == "GTiff"
+            assert src.meta["dtype"] == "uint16"
+            data = src.read(masked=True)
+            assert data.mask.any()
+
+
+def test_convert_geojson(landpoly, mp_tmpdir):
+    job = convert(landpoly, mp_tmpdir, output_pyramid="geodetic", zoom=4)
+    assert len(job)
+    for (zoom, row, col), control in zip([(4, 0, 7), (4, 1, 7)], [9, 32]):
+        out_file = os.path.join(
+            *[mp_tmpdir, str(zoom), str(row), str(col) + ".geojson"]
+        )
+        with fiona.open(out_file, "r") as src:
+            assert len(src) == control
+            for f in src:
+                assert shape(f["geometry"]).is_valid
+
+
+def test_convert_geobuf(landpoly, mp_tmpdir):
+    # convert to geobuf
+    geobuf_outdir = os.path.join(mp_tmpdir, "geobuf")
+    job = convert(
+        landpoly,
+        geobuf_outdir,
+        output_pyramid="geodetic",
+        zoom=4,
+        output_format="Geobuf",
+    )
+    assert len(job)
+    for (zoom, row, col), control in zip([(4, 0, 7), (4, 1, 7)], [9, 32]):
+        out_file = os.path.join(
+            *[geobuf_outdir, str(zoom), str(row), str(col) + ".pbf"]
+        )
+        with open(out_file, "rb") as src:
+            features = geobuf.decode(src.read())["features"]
+            assert len(features) == control
+            for f in features:
+                assert f["geometry"]["type"] == "Polygon"
+                assert shape(f["geometry"]).area
+
+
+def test_convert_errors(s2_band_jp2, mp_tmpdir, s2_band, cleantopo_br, landpoly):
+    # output format required
+    with pytest.raises(ValueError):
+        convert(s2_band_jp2, mp_tmpdir, output_pyramid="geodetic")
+
+    # output pyramid reqired
+    with pytest.raises(ValueError):
+        convert(s2_band, mp_tmpdir)
+
+    # prepare data for tiledir input
+    with mapchete.open(cleantopo_br.path) as mp:
+        mp.batch_process(zoom=[1, 4])
+    tiledir_path = os.path.join(
+        cleantopo_br.dict["config_dir"], cleantopo_br.dict["output"]["path"]
+    )
+
+    # zoom level required
+    with pytest.raises(ValueError):
+        convert(tiledir_path, mp_tmpdir, output_pyramid="geodetic")
+
+    # incompatible formats
+    with pytest.raises(ValueError):
+        convert(
+            tiledir_path,
+            mp_tmpdir,
+            output_pyramid="geodetic",
+            zoom=5,
+            output_format="GeoJSON",
+        )
+
+    # unsupported output format extension
+    with pytest.raises(ValueError):
+        convert(s2_band_jp2, "output.jp2", output_pyramid="geodetic", zoom=5)
+
+    # malformed band index
+    with pytest.raises(ValueError):
+        convert(s2_band_jp2, "output.tif", bidx="invalid")
+
+
+def test_index_geojson(mp_tmpdir, cleantopo_br):
+    # execute process at zoom 3
+    job = execute(cleantopo_br.path, zoom=3)
+    assert len(job)
+
+    # generate index for zoom 3
+    job = index(cleantopo_br.path, zoom=3, geojson=True)
+    assert len(job)
+
+    with mapchete.open(cleantopo_br.dict) as mp:
+        files = os.listdir(mp.config.output.path)
+        assert len(files) == 3
+        assert "3.geojson" in files
+    with fiona.open(os.path.join(mp.config.output.path, "3.geojson")) as src:
+        for f in src:
+            assert "location" in f["properties"]
+        assert len(list(src)) == 1
+
+
+def test_index_geojson_fieldname(mp_tmpdir, cleantopo_br):
+    # execute process at zoom 3
+    job = execute(cleantopo_br.path, zoom=3)
+    assert len(job)
+
+    # index and rename "location" to "new_fieldname"
+    job = index(
+        cleantopo_br.path,
+        zoom=3,
+        geojson=True,
+        fieldname="new_fieldname",
+    )
+    assert len(job)
+    with mapchete.open(cleantopo_br.dict) as mp:
+        files = os.listdir(mp.config.output.path)
+        assert "3.geojson" in files
+    with fiona.open(os.path.join(mp.config.output.path, "3.geojson")) as src:
+        for f in src:
+            assert "new_fieldname" in f["properties"]
+        assert len(list(src)) == 1
+
+
+def test_index_geojson_basepath(mp_tmpdir, cleantopo_br):
+    # execute process at zoom 3
+    job = execute(cleantopo_br.path, zoom=3)
+    assert len(job)
+
+    basepath = "http://localhost"
+    # index and rename "location" to "new_fieldname"
+    job = index(cleantopo_br.path, zoom=3, geojson=True, basepath=basepath)
+    assert len(job)
+
+    with mapchete.open(cleantopo_br.dict) as mp:
+        files = os.listdir(mp.config.output.path)
+        assert "3.geojson" in files
+    with fiona.open(os.path.join(mp.config.output.path, "3.geojson")) as src:
+        for f in src:
+            assert f["properties"]["location"].startswith(basepath)
+        assert len(list(src)) == 1
+
+
+def test_index_geojson_for_gdal(mp_tmpdir, cleantopo_br):
+    # execute process at zoom 3
+    job = execute(cleantopo_br.path, zoom=3)
+    assert len(job)
+
+    basepath = "http://localhost"
+    # index and rename "location" to "new_fieldname"
+    job = index(
+        cleantopo_br.path, zoom=3, geojson=True, basepath=basepath, for_gdal=True
+    )
+    assert len(job)
+
+    with mapchete.open(cleantopo_br.dict) as mp:
+        files = os.listdir(mp.config.output.path)
+        assert "3.geojson" in files
+    with fiona.open(os.path.join(mp.config.output.path, "3.geojson")) as src:
+        for f in src:
+            assert f["properties"]["location"].startswith("/vsicurl/" + basepath)
+        assert len(list(src)) == 1
+
+
+def test_index_geojson_tile(mp_tmpdir, cleantopo_tl):
+    # execute process at zoom 3
+    job = execute(cleantopo_tl.path, zoom=3)
+    assert len(job)
+
+    # generate index
+    job = index(cleantopo_tl.path, tile=(3, 0, 0), geojson=True)
+    assert len(job)
+
+    with mapchete.open(cleantopo_tl.dict) as mp:
+        files = os.listdir(mp.config.output.path)
+        assert len(files) == 3
+        assert "3.geojson" in files
+    with fiona.open(os.path.join(mp.config.output.path, "3.geojson")) as src:
+        assert len(list(src)) == 1
+
+
+def test_index_geojson_wkt_area(mp_tmpdir, cleantopo_tl, wkt_geom_tl):
+    # execute process at zoom 3
+    job = execute(cleantopo_tl.path, area=wkt_geom_tl)
+    assert len(job)
+
+    # generate index for zoom 3
+    job = index(cleantopo_tl.path, geojson=True, area=wkt_geom_tl)
+    assert len(job)
+
+    with mapchete.open(cleantopo_tl.dict) as mp:
+        files = os.listdir(mp.config.output.path)
+        assert len(files) == 13
+        assert "3.geojson" in files
+
+
+def test_index_gpkg(mp_tmpdir, cleantopo_br):
+    # execute process
+    job = execute(cleantopo_br.path, zoom=5)
+    assert len(job)
+
+    # generate index
+    job = index(cleantopo_br.path, zoom=5, gpkg=True)
+    assert len(job)
+    with mapchete.open(cleantopo_br.dict) as mp:
+        files = os.listdir(mp.config.output.path)
+        assert "5.gpkg" in files
+    with fiona.open(os.path.join(mp.config.output.path, "5.gpkg")) as src:
+        for f in src:
+            assert "location" in f["properties"]
+        assert len(list(src)) == 1
+
+    # write again and assert there is no new entry because there is already one
+    job = index(cleantopo_br.path, zoom=5, gpkg=True)
+    assert len(job)
+    with mapchete.open(cleantopo_br.dict) as mp:
+        files = os.listdir(mp.config.output.path)
+        assert "5.gpkg" in files
+    with fiona.open(os.path.join(mp.config.output.path, "5.gpkg")) as src:
+        for f in src:
+            assert "location" in f["properties"]
+        assert len(list(src)) == 1
+
+
+def test_index_shp(mp_tmpdir, cleantopo_br):
+    # execute process
+    job = execute(cleantopo_br.path, zoom=5)
+    assert len(job)
+
+    # generate index
+    job = index(cleantopo_br.path, zoom=5, shp=True)
+    assert len(job)
+    with mapchete.open(cleantopo_br.dict) as mp:
+        files = os.listdir(mp.config.output.path)
+        assert "5.shp" in files
+    with fiona.open(os.path.join(mp.config.output.path, "5.shp")) as src:
+        for f in src:
+            assert "location" in f["properties"]
+        assert len(list(src)) == 1
+
+    # write again and assert there is no new entry because there is already one
+    job = index(cleantopo_br.path, zoom=5, shp=True)
+    assert len(job)
+    with mapchete.open(cleantopo_br.dict) as mp:
+        files = os.listdir(mp.config.output.path)
+        assert "5.shp" in files
+    with fiona.open(os.path.join(mp.config.output.path, "5.shp")) as src:
+        for f in src:
+            assert "location" in f["properties"]
+        assert len(list(src)) == 1
+
+
+def test_index_text(cleantopo_br):
+    # execute process
+    job = execute(cleantopo_br.path, zoom=5)
+    assert len(job)
+
+    # generate index
+    job = index(cleantopo_br.path, zoom=5, txt=True)
+    assert len(job)
+    with mapchete.open(cleantopo_br.dict) as mp:
+        files = os.listdir(mp.config.output.path)
+        assert "5.txt" in files
+    with open(os.path.join(mp.config.output.path, "5.txt")) as src:
+        lines = list(src)
+        assert len(lines) == 1
+        for l in lines:
+            assert l.endswith("7.tif\n")
+
+    # write again and assert there is no new entry because there is already one
+    job = index(cleantopo_br.path, zoom=5, txt=True)
+    assert len(job)
+    with mapchete.open(cleantopo_br.dict) as mp:
+        files = os.listdir(mp.config.output.path)
+        assert "5.txt" in files
+    with open(os.path.join(mp.config.output.path, "5.txt")) as src:
+        lines = list(src)
+        assert len(lines) == 1
+        for l in lines:
+            assert l.endswith("7.tif\n")
+
+
+def test_index_errors(mp_tmpdir, cleantopo_br):
+    with pytest.raises(ValueError):
+        index(cleantopo_br.path, zoom=5)
diff --git a/test/test_config.py b/test/test_config.py
index 0e823a5b..35339852 100644
--- a/test/test_config.py
+++ b/test/test_config.py
@@ -442,3 +442,8 @@ def test_init_overrides_config(example_mapchete):
     with mapchete.open(dict(example_mapchete.dict, bounds=process_bounds)) as mp:
         assert mp.config.bounds == process_bounds
         assert mp.config.init_bounds == process_bounds
+
+
+def test_custom_process(example_custom_process_mapchete):
+    with mapchete.open(example_custom_process_mapchete.dict) as mp:
+        assert callable(mp.config.process_func)
diff --git a/test/test_errors.py b/test/test_errors.py
index 52290a9c..51e714c6 100644
--- a/test/test_errors.py
+++ b/test/test_errors.py
@@ -300,7 +300,7 @@ def test_finished_task():
     assert task.exception()
     with pytest.raises(ZeroDivisionError):
         task.result()
-    assert "FinishedTask" in str(task)
+    assert "FakeFuture" in str(task)
 
 
 def test_strip_zoom_error(files_zooms):
diff --git a/test/test_executor.py b/test/test_executor.py
new file mode 100644
index 00000000..a7d749eb
--- /dev/null
+++ b/test/test_executor.py
@@ -0,0 +1,127 @@
+import pytest
+import time
+
+from mapchete import Executor
+from mapchete._executor import FakeFuture
+
+
+def _dummy_process(i, sleep=0):
+    time.sleep(sleep)
+    return i + 1
+
+
+def test_sequential_executor_as_completed():
+    items = 10
+    count = 0
+    with Executor(concurrency=None) as executor:
+        # process all
+        for future in executor.as_completed(_dummy_process, range(items)):
+            count += 1
+            assert future.result()
+        assert items == count
+
+        # abort
+        cancelled = False
+        for future in executor.as_completed(_dummy_process, range(items)):
+            if cancelled:
+                raise RuntimeError()
+            assert future.result()
+            cancelled = True
+            executor.cancel()
+
+
+def test_sequential_executor_map():
+    items = list(range(10))
+    with Executor(concurrency=None) as executor:
+        result = executor.map(_dummy_process, items)
+        assert [i + 1 for i in items] == result
+
+
+def test_concurrent_futures_processes_executor_as_completed():
+    items = 10
+    with Executor(concurrency="processes") as executor:
+        # process all
+        count = 0
+        for future in executor.as_completed(_dummy_process, range(items)):
+            count += 1
+            assert future.result()
+
+
+def test_concurrent_futures_processes_executor_cancel_as_completed():
+    items = 100
+    with Executor(concurrency="processes", max_workers=2) as executor:
+        # abort
+        for future in executor.as_completed(
+            _dummy_process, range(items), fkwargs=dict(sleep=2)
+        ):
+            assert future.result()
+            executor.cancel()
+            break
+
+        assert any([future.cancelled() for future in executor.futures])
+
+
+def test_concurrent_futures_processes_executor_map():
+    items = list(range(10))
+    with Executor(concurrency="processes") as executor:
+        result = executor.map(_dummy_process, items)
+        assert [i + 1 for i in items] == result
+
+
+def test_concurrent_futures_threads_executor_as_completed():
+    items = 100
+    with Executor(concurrency="threads", max_workers=2) as executor:
+        # abort
+        for future in executor.as_completed(
+            _dummy_process, range(items), fkwargs=dict(sleep=2)
+        ):
+            assert future.result()
+            executor.cancel()
+            break
+
+        assert any([future.cancelled() for future in executor.futures])
+
+
+def test_concurrent_futures_threads_executor_map():
+    items = list(range(10))
+    with Executor(concurrency="threads") as executor:
+        result = executor.map(_dummy_process, items)
+        assert [i + 1 for i in items] == result
+
+
+def test_dask_executor_as_completed():
+    items = 100
+    with Executor(concurrency="dask", max_workers=2) as executor:
+        # abort
+        for future in executor.as_completed(
+            _dummy_process, range(items), fkwargs=dict(sleep=2)
+        ):
+            assert future.result()
+            executor.cancel()
+            break
+
+        assert any([future.cancelled() for future in executor.futures])
+
+
+def test_concurrent_futures_dask_executor_map():
+    items = list(range(10))
+    with Executor(concurrency="dask") as executor:
+        result = executor.map(_dummy_process, items)
+        assert [i + 1 for i in items] == result
+
+
+def test_fake_future():
+    def task(*args, **kwargs):
+        return True
+
+    def failing_task(*args, **kwargs):
+        raise RuntimeError()
+
+    future = FakeFuture(task, fargs=[1, True], fkwargs=dict(foo="bar"))
+    assert future.result()
+    assert not future.exception()
+
+    future = FakeFuture(failing_task, fargs=[1, True], fkwargs=dict(foo="bar"))
+    with pytest.raises(RuntimeError):
+        future.result()
+    assert future.exception()
diff --git a/test/test_formats_tiledir_input.py b/test/test_formats_tiledir_input.py
index d4cf8bc5..53815185 100644
--- a/test/test_formats_tiledir_input.py
+++ b/test/test_formats_tiledir_input.py
@@ -157,18 +157,6 @@ def test_no_metadata_json(mp_tmpdir, cleantopo_br_tiledir):
         )
 
 
-@pytest.mark.remote
-def test_read_remote_raster_data(mp_tmpdir, cleantopo_remote):
-    """Read raster data."""
-    with mapchete.open(cleantopo_remote.path) as mp:
-        assert all(
-            [
-                next(iter(mp.config.input.values())).open(tile).read().any()
-                for tile in mp.get_process_tiles(1)
-            ]
-        )
-
-
 def test_parse_errors(geojson_tiledir, cleantopo_br_tiledir):
     """Different configuration exceptions."""
     # without path
diff --git a/test/test_mapchete.py b/test/test_mapchete.py
index 66e3632e..736d00d5 100644
--- a/test/test_mapchete.py
+++ b/test/test_mapchete.py
@@ -11,11 +11,9 @@
 import pkg_resources
 
 try:
-    from cPickle import dumps
+    from cPickle import dumps as pickle_dumps
 except ImportError:
-    from pickle import dumps
-from functools import partial
-from multiprocessing import Pool
+    from pickle import dumps as pickle_dumps
 from shapely.geometry import box, shape
 from shapely.ops import unary_union
 
@@ -387,27 +385,14 @@ def test_processing(mp_tmpdir, cleantopo_br, cleantopo_tl):
                     shutil.rmtree(mp_tmpdir, ignore_errors=True)
 
 
-def test_multiprocessing(mp_tmpdir, cleantopo_tl):
+def test_pickleability(mp_tmpdir, cleantopo_tl):
     """Test parallel tile processing."""
     with mapchete.open(cleantopo_tl.path) as mp:
-        assert dumps(mp)
-        assert dumps(mp.config)
-        assert dumps(mp.config.output)
+        assert pickle_dumps(mp)
+        assert pickle_dumps(mp.config)
+        assert pickle_dumps(mp.config.output)
         for tile in mp.get_process_tiles():
-            assert dumps(tile)
-        f = partial(_worker, mp)
-        try:
-            pool = Pool()
-            for zoom in reversed(mp.config.zoom_levels):
-                for tile, raw_output in pool.imap_unordered(
-                    f, mp.get_process_tiles(zoom), chunksize=8
-                ):
-                    mp.write(tile, raw_output)
-        except KeyboardInterrupt:
-            pool.terminate()
-        finally:
-            pool.close()
-            pool.join()
+            assert pickle_dumps(tile)
 
 
 def _worker(mp, tile):
@@ -500,10 +485,10 @@ def test_skip_tiles(mp_tmpdir, cleantopo_tl):
 def test_custom_grid(mp_tmpdir, custom_grid):
     """Cutom grid processing."""
     # process and save
-    with mapchete.open(custom_grid.dict) as mp:
+    with mapchete.open(custom_grid.path) as mp:
         mp.batch_process()
     # read written output
-    with mapchete.open(custom_grid.dict) as mp:
+    with mapchete.open(custom_grid.path) as mp:
         for tile in mp.get_process_tiles(5):
             data = mp.config.output.read(tile)
             assert data.any()
diff --git a/test/test_preprocessing.py b/test/test_preprocessing.py
new file mode 100644
index 00000000..3f3fd184
--- /dev/null
+++ b/test/test_preprocessing.py
@@ -0,0 +1,89 @@
+import mapchete
+import pytest
+
+
+def _trivial_func(arg, kwarg=None):
+    return arg + kwarg + "bar"
+
+
+def test_preprocessing_empty_tasks_per_input(example_mapchete):
+    with mapchete.open(example_mapchete.dict) as mp:
+        assert len(mp.config.input) == len(mp.config.preprocessing_tasks_per_input())
+
+
+def test_preprocessing_empty_tasks(example_mapchete):
+    with mapchete.open(example_mapchete.dict) as mp:
+        assert len(mp.config.preprocessing_tasks()) == 0
+
+
+def test_preprocessing_empty_tasks_count(example_mapchete):
+    with mapchete.open(example_mapchete.dict) as mp:
+        assert mp.config.preprocessing_tasks_count() == 0
+
+
+def test_preprocessing_set_empty_preprocessing_task_result(example_mapchete):
+    with mapchete.open(example_mapchete.dict) as mp:
+        with pytest.raises(KeyError):
+            mp.config.set_preprocessing_task_result("foo", "bar")
+
+
+def test_add_preprocessing_task(example_mapchete):
+    with mapchete.open(example_mapchete.dict) as mp:
+        # get input object
+        inp = mp.config.input_at_zoom("file1", 10)
+        tasks = mp.count_tasks()
+        # add a preprocessing task and make sure it is registered
+        inp.add_preprocessing_task(_trivial_func, fargs=("foo"))
+        assert len(mp.config.preprocessing_tasks()) == 1
+        assert mp.config.preprocessing_tasks_count() == 1
+        assert mp.count_tasks() == tasks + 1
+
+
+def test_run_preprocessing_task(example_mapchete):
+    with mapchete.open(example_mapchete.dict) as mp:
+        # get input object
+        inp = mp.config.input_at_zoom("file1", 10)
+        with pytest.raises(KeyError):
+            inp.get_preprocessing_task_result("test_task")
+        # add a preprocessing task and make sure it is registered
+        inp.add_preprocessing_task(
+            _trivial_func, key="test_task", fargs="foo", fkwargs={"kwarg": "foo"}
+        )
+        with pytest.raises(ValueError):
+            inp.get_preprocessing_task_result("test_task")
+        mp.batch_preprocess()
+        assert inp.get_preprocessing_task_result("test_task") == "foofoobar"
+
+
+def test_run_preprocessing_tasks(example_mapchete):
+    with mapchete.open(example_mapchete.dict) as mp:
+        inp1 = mp.config.input_at_zoom("file1", 10)
+        inp1.add_preprocessing_task(
+            _trivial_func, key="test_task", fargs="foo", fkwargs={"kwarg": "foo"}
+        )
+        inp1.add_preprocessing_task(
+            _trivial_func, key="test_other_task", fargs="bar", fkwargs={"kwarg": "foo"}
+        )
+        inp2 = mp.config.input_at_zoom("file2", 10)
+        inp2.add_preprocessing_task(
+            _trivial_func, key="test_task", fargs="foo", fkwargs={"kwarg": "foo"}
+        )
+
+        # total number of tasks should be 3
+        assert (
+            sum(
+                [
+                    len(tasks)
+                    for tasks in mp.config.preprocessing_tasks_per_input().values()
+                ]
+            )
+            == 3
+        )
+
+        # number of unique tasks should be 2
+        assert mp.config.preprocessing_tasks_count() == 2
+
+        mp.batch_preprocess()
+        assert inp1.get_preprocessing_task_result("test_task") == "foofoobar"
+        assert inp1.get_preprocessing_task_result("test_other_task") == "barfoobar"
+        assert inp2.get_preprocessing_task_result("test_task") == "foofoobar"
diff --git a/test/test_stac.py b/test/test_stac.py
new file mode 100644
index 00000000..1b8bdba6
--- /dev/null
+++ b/test/test_stac.py
@@ -0,0 +1,160 @@
+import pystac
+import pytest
+from shapely.geometry import box, shape
+
+from mapchete.io.vector import reproject_geometry
+from mapchete.stac import tile_directory_stac_item
+from mapchete.tile import BufferedTilePyramid
+
+
+def test_wkss_geodetic():
+    tp = BufferedTilePyramid("geodetic")
+    item = tile_directory_stac_item(
+        item_id="foo", item_path="foo/bar.json", tile_pyramid=tp, max_zoom=5
+    )
+    assert item.id == "foo"
+    assert shape(item.geometry).difference(box(*tp.bounds)).is_empty
+    assert item.bbox == list(tp.bounds)
+    assert item.datetime
+    assert "tiled-assets" in item.stac_extensions
+    assert "bands" in item.extra_fields["asset_templates"]
+    assert "tiles:tile_matrix_links" in item.properties
+    assert "tiles:tile_matrix_sets" in item.properties
+    assert "WorldCRS84Quad" in item.properties["tiles:tile_matrix_sets"]
+
+
+def test_wkss_mercator():
+    tp = BufferedTilePyramid("mercator")
+    item = tile_directory_stac_item(
+        item_id="foo", item_path="foo/bar.json", tile_pyramid=tp, max_zoom=5
+    )
+    assert item.id == "foo"
+    item_geometry = reproject_geometry(
+        shape(item.geometry), src_crs="EPSG:4326", dst_crs=tp.crs
+    )
+    assert item_geometry.difference(box(*tp.bounds)).is_empty
+    assert item.datetime
+    assert "tiled-assets" in item.stac_extensions
+    assert "bands" in item.extra_fields["asset_templates"]
+    assert "tiles:tile_matrix_links" in item.properties
+    assert "tiles:tile_matrix_sets" in item.properties
+    assert "WebMercatorQuad" in item.properties["tiles:tile_matrix_sets"]
+
+
+def test_custom_datetime():
+    item = tile_directory_stac_item(
+        item_id="foo",
+        item_path="foo/bar.json",
+        tile_pyramid=BufferedTilePyramid("geodetic"),
+        max_zoom=5,
+        item_metadata=dict(properties=dict(start_datetime="2021-01-01 00:00:00")),
+    )
+    assert str(item.datetime) == "2021-01-01 00:00:00"
+
+
+def test_custom_tilematrix():
+    tp = BufferedTilePyramid(
+        grid=dict(
+            shape=[117, 9],
+            bounds=[145980, 0, 883260, 9584640.0],
+            is_global=False,
+            epsg=32630,
+        ),
+        metatiling=4,
+    )
+    item = tile_directory_stac_item(
+        item_id="foo",
+        item_path="foo/bar.json",
+        tile_pyramid=tp,
+        max_zoom=5,
+        item_metadata=dict(properties=dict(start_datetime="2021-01-01 00:00:00")),
+    )
+    assert str(item.datetime) == "2021-01-01 00:00:00"
+    assert "custom" in item.properties["tiles:tile_matrix_sets"]
+
+
+def test_tiled_asset_path():
+    # default: create absolute path from item basepath
+    item = tile_directory_stac_item(
+        item_id="foo",
+        item_path="foo/bar.json",
+        tile_pyramid=BufferedTilePyramid("geodetic"),
+        max_zoom=0,
+    )
+    basepath = item.to_dict()["asset_templates"]["bands"]["href"]
+    assert basepath.startswith("foo/")
+
+    # use alternative asset basepath
+    item = tile_directory_stac_item(
+        item_id="foo",
+        item_path="foo/bar.json",
+        asset_basepath="s3://bar/",
+        tile_pyramid=BufferedTilePyramid("geodetic"),
+        max_zoom=0,
+    )
+    basepath = item.to_dict()["asset_templates"]["bands"]["href"]
+    assert basepath.startswith("s3://bar/")
+
+    # create relative path
+    item = tile_directory_stac_item(
+        item_id="foo",
+        relative_paths=True,
+        tile_pyramid=BufferedTilePyramid("geodetic"),
+        max_zoom=0,
+    )
+    basepath = item.to_dict()["asset_templates"]["bands"]["href"]
+    assert basepath.startswith("{TileMatrix}/{TileRow}")
+
+
+def test_tiled_asset_eo_bands_metadata():
+    item = tile_directory_stac_item(
+        item_id="foo",
+        item_path="foo/bar.json",
+        tile_pyramid=BufferedTilePyramid("geodetic"),
+        max_zoom=0,
+        item_metadata={"eo:bands": {"foo": "bar"}},
+    )
+    assert "eo" in item.to_dict()["stac_extensions"]
+    assert "eo:bands" in item.to_dict()["asset_templates"]["bands"]
+
+
+def test_create_stac_item_errors():
+    tp = BufferedTilePyramid("geodetic")
+    # no item_id
+    with pytest.raises(ValueError):
+        tile_directory_stac_item(
+            item_path="foo/bar.json",
+            tile_pyramid=tp,
+            max_zoom=5,
+        )
+
+    # no max_zoom
+    with pytest.raises(ValueError):
+        tile_directory_stac_item(
+            item_id="foo",
+            item_path="foo/bar.json",
+            tile_pyramid=tp,
+        )
+
+    # no tile_pyramid
+    with pytest.raises(ValueError):
+        tile_directory_stac_item(
+            item_id="foo",
+            item_path="foo/bar.json",
+            max_zoom=5,
+        )
+
+    # no item_path or asset_basepath
+    with pytest.raises(ValueError):
+        tile_directory_stac_item(
+            item_id="foo",
+            tile_pyramid=tp,
+            max_zoom=5,
+        )
+
+
+# def test_single_file():
+#     tile_directory_stac_item(
+#         item_id="foo",
+
+#     )
diff --git a/test/testdata/cleantopo_br_metatiling_1.mapchete b/test/testdata/cleantopo_br_metatiling_1.mapchete
new file mode 100644
index 00000000..43e651d5
--- /dev/null
+++ b/test/testdata/cleantopo_br_metatiling_1.mapchete
@@ -0,0 +1,15 @@
+process: ../example_process.py
+zoom_levels:
+    min: 0
+    max: 5
+pyramid:
+    grid: geodetic
+    pixelbuffer: 20
+input:
+    file1: cleantopo_br.tif
+output:
+    dtype: uint16
+    bands: 1
+    format: GTiff
+    path: tmp/cleantopo_br_metatiling_1
+    pixelbuffer: 20
diff --git a/test/testdata/example_custom_process.mapchete b/test/testdata/example_custom_process.mapchete
new file mode 100644
index 00000000..782c1181
--- /dev/null
+++ b/test/testdata/example_custom_process.mapchete
@@ -0,0 +1,54 @@
+# mandatory parameters
+######################
+# this is the location of user python code:
+process:
+    - '"""Example process file."""'
+    - ''
+    - ''
+    - 'def execute(mp):'
+    - '    """User defined process."""'
+    - '    # Reading and writing data works like this:'
+    - '    with mp.open("file1") as raster_file:'
+    - '        if raster_file.is_empty():'
+    - '            return "empty"'
+    - '            # This assures a transparent tile instead of a pink error tile'
+    - '            # is returned when using mapchete serve.'
+    - '        dem = raster_file.read(resampling="bilinear")'
+    - '    return dem'
+
+
+# zoom level range:
+zoom_levels:
+    min: 7
+    max: 11
+# or define single zoom level
+# zoom_levels: 5
+
+# geographical subset:
+# bounds: [1.0, 2.0, 3.0, 4.0]
+
+# output pyramid definition
+
+pyramid:
+    grid: geodetic
+    metatiling: 1 # can be 1, 2, 4, 8, 16 (default 1)
+
+
+input:
+    file1:
+        zoom>=10: dummy1.tif
+    file2: dummy2.tif
+output:
+    path: tmp/example_custom_process
+    format: GTiff
+    dtype: float32
+    bands: 3
+
+# free parameters
+#################
+some_integer_parameter: 12
+some_float_parameter: 5.3
+some_string_parameter:
+    zoom<=7: string1
+    zoom>7: string2
+some_bool_parameter: true
diff --git a/test/testdata/metadata.json b/test/testdata/metadata.json
new file mode 100644
index 00000000..9975596a
--- /dev/null
+++ b/test/testdata/metadata.json
@@ -0,0 +1,64 @@
+{
+    "driver": {
+        "bands": 1,
+        "compress": "deflate",
+        "delimiters": {
+            "bounds": [
+                -180.0,
+                -90.0,
+                180.0,
+                90.0
+            ],
+            "effective_bounds": [
+                -180.17578125,
+                -90.0,
+                180.17578125,
+                90.0
+            ],
+            "process_bounds": [
+                -180.0,
+                -90.0,
+                180.0,
+                90.0
+            ],
+            "zoom": [
+                0,
+                1,
+                2,
+                3,
+                4,
+                5,
+                6,
+                7,
+                8
+            ]
+        },
+        "dtype": "float32",
+        "format": "GTiff",
+        "mode": "continue",
+        "nodata": -32768.0,
+        "predictor": 3
+    },
+    "pyramid": {
+        "grid": {
+            "bounds": [
+                -180.0,
+                -90.0,
+                180.0,
+                90.0
+            ],
+            "is_global": true,
+            "shape": [
+                1,
+                2
+            ],
+            "srs": {
+                "wkt": "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]"
+            },
+            "type": "geodetic"
+        },
+        "metatiling": 8,
+        "pixelbuffer": 64,
+        "tile_size": 256
+    }
+}
\ No newline at end of file