Repository: ungarj/mapchete (language: python)

Issue:
commands package

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:

Fix patch:
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/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..8c50305f 100644
--- a/mapchete/__init__.py
+++ b/mapchete/__init__.py
@@ -1,12 +1,13 @@
 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"]
+__all__ = ["open", "count_tiles", "Mapchete", "MapcheteProcess", "ProcessInfo", "Timer"]
 __version__ = "0.40"
 
 
diff --git a/mapchete/_core.py b/mapchete/_core.py
index 7d4719fc..13c15e91 100644
--- a/mapchete/_core.py
+++ b/mapchete/_core.py
@@ -6,7 +6,7 @@
 import os
 import threading
 
-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
@@ -76,15 +76,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 +155,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
 
@@ -188,11 +188,14 @@ def batch_process(
         self,
         zoom=None,
         tile=None,
+        distributed=False,
+        dask_scheduler=None,
         multi=None,
         max_chunksize=1,
         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.
@@ -217,7 +220,7 @@ def batch_process(
             (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;
@@ -226,11 +229,14 @@ def batch_process(
             self.batch_processor(
                 zoom=zoom,
                 tile=tile,
+                distributed=distributed,
+                dask_scheduler=dask_scheduler,
                 multi=multi or multiprocessing.cpu_count(),
                 max_chunksize=max_chunksize,
                 multiprocessing_module=multiprocessing_module or multiprocessing,
                 multiprocessing_start_method=multiprocessing_start_method,
                 skip_output_check=skip_output_check,
+                executor=executor,
             )
         )
 
@@ -238,11 +244,14 @@ def batch_processor(
         self,
         zoom=None,
         tile=None,
+        distributed=False,
+        dask_scheduler=None,
         multi=None,
         max_chunksize=1,
         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.
@@ -265,7 +274,7 @@ def batch_processor(
             (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;
@@ -276,22 +285,27 @@ def batch_processor(
         # 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)),
+                distributed=distributed,
+                dask_scheduler=dask_scheduler,
                 multi=multi or multiprocessing.cpu_count(),
                 max_chunksize=max_chunksize,
                 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_tiles(self, minzoom=None, maxzoom=None, init_zoom=0):
         """
         Count number of tiles intersecting with process area at given zoom levels.
 
@@ -307,13 +321,15 @@ def count_tiles(self, minzoom, maxzoom, init_zoom=0):
         -------
         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)]
 
@@ -342,8 +358,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 +558,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..c6817274
--- /dev/null
+++ b/mapchete/_executor.py
@@ -0,0 +1,262 @@
+from cached_property import cached_property
+import concurrent.futures
+from functools import partial
+from itertools import chain
+import logging
+import multiprocessing
+import os
+
+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):
+            yield future
+
+    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 _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,
+            initializer=set_log_level,
+            initargs=(logger.getEffectiveLevel(),),
+        )
+        if concurrency == "processes":
+            self._executor_cls = concurrent.futures.ProcessPoolExecutor
+            start_method = (
+                kwargs.get("multiprocessing_start_method")
+                or MULTIPROCESSING_DEFAULT_START_METHOD
+            )
+            self._executor_kwargs.update(
+                mp_context=multiprocessing.get_context(method=start_method)
+            )
+        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
+
+
+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 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..8f2e1ff9 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,109 +408,12 @@ 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."""
-
-    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 __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)
-
-
-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)
-
-
-def _exception_wrapper(func, fargs, fkwargs, i):
-    """Wrap function around FinishedTask object."""
-    return FinishedTask(func, list(chain([i], fargs)), fkwargs)
-
-
 ###########################
 # batch execution options #
 ###########################
 
 
-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,8 +429,11 @@ def _run_on_single_tile(process=None, tile=None):
 
 
 def _run_area(
+    executor=None,
     process=None,
     zoom_levels=None,
+    distributed=False,
+    dask_scheduler=None,
     multi=None,
     max_chunksize=None,
     multiprocessing_module=None,
@@ -475,9 +446,12 @@ 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,
+            distributed=distributed,
+            dask_scheduler=dask_scheduler,
             multi=multi,
             multiprocessing_start_method=multiprocessing_start_method,
             multiprocessing_module=multiprocessing_module,
@@ -490,10 +464,13 @@ 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,
+            dask_scheduler=dask_scheduler,
+            distributed=distributed,
             multi=multi,
             multiprocessing_start_method=multiprocessing_start_method,
             multiprocessing_module=multiprocessing_module,
@@ -520,9 +497,12 @@ def _filter_skipable(process=None, tiles=None, todo=None, target_set=None):
 
 
 def _run_multi(
+    executor=None,
     func=None,
     zoom_levels=None,
     process=None,
+    distributed=False,
+    dask_scheduler=None,
     multi=None,
     multiprocessing_start_method=None,
     multiprocessing_module=None,
@@ -534,38 +514,102 @@ def _run_multi(
     total_tiles = process.count_tiles(min(zoom_levels), max(zoom_levels))
     workers = min([multi, 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 +619,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 +631,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 +665,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 +697,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..2cb32d21 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="tile",
+                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..b85996a4
--- a/mapchete/cli/default/execute.py
+++ b/mapchete/cli/default/execute.py
@@ -1,108 +1,84 @@
-"""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_max_chunksize
+@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="tile",
+                    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/utils.py b/mapchete/cli/options.py
similarity index 58%
rename from mapchete/cli/utils.py
rename to mapchete/cli/options.py
index 7e5ebd14..6422c89d 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.",
 )
 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."
@@ -257,18 +245,17 @@ def _cb_key_val(ctx, param, value):
 )
 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)",
+    help="Maximum number of tasks 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 +343,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..e6c411cd
--- /dev/null
+++ b/mapchete/commands/_execute.py
@@ -0,0 +1,184 @@
+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,
+    max_chunksize: 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.
+    max_chunksize : int
+        Maximum number of process tiles to be queued for each  worker. (default: 1)
+    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:
+        tiles_count = mp.count_tiles()
+        if tile:
+            msg_callback("processing 1 tile")
+        else:
+            msg_callback(f"processing {tiles_count} tile(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 tiles_count == 1 or workers == 1:
+            logger.debug(
+                f"using sequential Executor because there is only one {'tile' if tiles_count == 1 else 'worker'}"
+            )
+            concurrency = None
+        return mapchete.Job(
+            _msg_wrapper,
+            fargs=(
+                msg_callback,
+                mp,
+            ),
+            fkwargs=dict(
+                tile=tile,
+                multi=workers,
+                zoom=None if tile else zoom,
+            ),
+            executor_concurrency=concurrency,
+            executor_kwargs=dict(
+                dask_scheduler=dask_scheduler,
+                dask_client=dask_client,
+                max_chunksize=max_chunksize,
+                multiprocessing_start_method=multiprocessing_start_method,
+            ),
+            as_iterator=as_iterator,
+            total=1 if tile else tiles_count,
+        )
+    # explicitly exit the mp object on failure
+    except Exception:  # pragma: no cover
+        mp.__exit__(None, None, None)
+        raise
+
+
+def _msg_wrapper(msg_callback, mp, executor=None, **kwargs):
+    try:
+        for process_info in mp.batch_processor(executor=executor, **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..4525940a 100644
--- a/mapchete/config.py
+++ b/mapchete/config.py
@@ -26,6 +26,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 +59,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
@@ -93,6 +94,8 @@
     "mode",
 ]
 
+MULTIPROCESSING_DEFAULT_START_METHOD = "spawn"
+
 
 class MapcheteConfig(object):
     """
@@ -206,7 +209,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
@@ -548,7 +551,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 +593,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 +635,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):
@@ -991,11 +998,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 +1013,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/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/_path.py b/mapchete/io/_path.py
index 1730e803..96053dea 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.
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/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/setup.py b/setup.py
index 4d6288bf..22d8b59a 100644
--- a/setup.py
+++ b/setup.py
@@ -39,12 +39,15 @@
     "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_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_vrt
+)
 
 setup(
     name="mapchete",
@@ -66,7 +69,7 @@
             "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",
         ],
         "mapchete.formats.drivers": [
@@ -94,6 +97,7 @@
     extras_require={
         "complete": req_complete,
         "contours": req_contours,
+        "dask": req_dask,
         "geobuf": req_geobuf,
         "s3": req_s3,
         "serve": req_serve,


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/test_cli.py b/test/test_cli.py
index 88697cb7..cf28e392 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,9 +41,13 @@ 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:
         raise result.exception
     assert result.exit_code == expected_exit_code
@@ -69,17 +73,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 +96,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 +118,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 +196,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 +214,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 +274,8 @@ def test_execute_logfile(mp_tmpdir, example_mapchete):
             "1040",
             "--logfile",
             logfile,
+            "--concurrency",
+            "none",
         ]
     )
     assert os.path.isfile(logfile)
@@ -235,14 +283,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 +300,8 @@ def test_execute_point(mp_tmpdir, example_mapchete, wkt_geom):
             "--point",
             str(g.centroid.x),
             str(g.centroid.y),
+            "--concurrency",
+            "none",
         ]
     )
 
@@ -274,7 +321,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 +343,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 +374,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 +403,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 +427,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 +452,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 +479,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 +505,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 +527,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 +553,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 +580,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 +605,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 +624,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 +644,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 +664,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 +687,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 +718,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 +741,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 +769,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 +801,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 +829,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 +847,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 +884,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 +904,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 +924,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 +944,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 +976,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 +1019,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 +1037,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 +1065,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 +1122,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 +1145,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 +1170,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 +1197,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 +1224,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 +1252,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 +1266,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 +1293,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 +1312,8 @@ def test_cp(mp_tmpdir, cleantopo_br, wkt_geom):
             "-90",
             "180",
             "-80.18582802550002",
+            "--concurrency",
+            "none",
         ]
     )
     # copy all tiles
@@ -1164,14 +1324,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 +1370,8 @@ def test_cp_http(mp_tmpdir, http_tiledir):
             "1.0",
             "4.0",
             "2.0",
+            "--concurrency",
+            "none",
         ]
     )
 
@@ -1204,6 +1388,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 +1435,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,
         [
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..9c33b4d2
--- /dev/null
+++ b/test/test_executor.py
@@ -0,0 +1,99 @@
+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():
+    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_concurrent_futures_processes_executor():
+    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():
+    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_threads_executor():
+    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_dask_executor():
+    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_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..56a31072 100644
--- a/test/test_mapchete.py
+++ b/test/test_mapchete.py
@@ -500,10 +500,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/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