Repository: theroggy/cropclassification (language: python)

Issue:
Add support to use openeo for image retrieval/calculation
Images can be created as periodic mosaics using openeo like it worked using gee in the (distant) past.

Todo-list:
- [x] generate openeo "s2-landcover" and "s2-ndvi" weekly images for the period defined during marker calculation
- [x] generate openeo "s2-agri" weekly images for the period defined in the marker during marker calculation
- [x] generate openeo "s1-grd" weekly images for the period defined in the marker during marker calculation
- [x] generate openeo "s1-coh" weekly images for the period defined in the marker during marker calculation
- [x] possibility to only generate openeo images in a seperate task
- [x] optimize cropgroup/landuse markers accuracy for these images



Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index 040a9b3c..2c2f7cd6 100644
--- a/README.md
+++ b/README.md
@@ -23,7 +23,7 @@ Once you have anaconda installed, you can open an anaconda terminal window and f
       ```
       2. Install the dependencies for the crop classification scripts:
       ```
-      conda install python=3.9 geopandas geofileops "h5py<3" psutil rasterio "rasterstats<0.16.0" scikit-learn
+      conda install python=3.9 geopandas geofileops "h5py<3" openeo psutil rasterio rasterstats scikit-learn
       ```
       3. If it was the first time you installed anaconda/geopandas, you might have to restart your computer to proceed
       4. Start the anaconda terminal window again and activate the environment
diff --git a/bin_util/calc_dias_weekly.py b/bin_util/calc_dias_weekly.py
index 88f79a80..965c1a0c 100644
--- a/bin_util/calc_dias_weekly.py
+++ b/bin_util/calc_dias_weekly.py
@@ -46,12 +46,11 @@ def main():
 
     # Calculate!
     input_parcel_path = conf.dirs.getpath("input_dir") / input_features_filename
-    ts_helper.calculate_periodic_data(
-        input_parcel_path=input_parcel_path,
-        input_base_dir=timeseries_per_image_dir,
-        start_date_str=f"{year}-03-15",
-        end_date_str=f"{year}-08-15",
-        # sensordata_to_get=conf.marker.getlist('sensordata_to_use'),
+    ts_helper.calculate_periodic_timeseries(
+        parcel_path=input_parcel_path,
+        timeseries_per_image_dir=timeseries_per_image_dir,
+        start_date=datetime.fromisoformat(f"{year}-03-15"),
+        end_date=datetime.fromisoformat(f"{year}-08-15"),
         sensordata_to_get=["SENSORDATA_S1_COHERENCE"],
         dest_data_dir=timeseries_periodic_dir,
         force=False,
diff --git a/cropclassification/calc_marker.py b/cropclassification/calc_marker.py
index 5be0301b..35cb7765 100644
--- a/cropclassification/calc_marker.py
+++ b/cropclassification/calc_marker.py
@@ -102,7 +102,6 @@ def calc_marker_task(config_paths: List[Path], default_basedir: Path):
     # Read the info about the run
     input_parcel_filename = conf.calc_marker_params.getpath("input_parcel_filename")
     input_parcel_filetype = conf.calc_marker_params["input_parcel_filetype"]
-    country_code = conf.calc_marker_params["country_code"]
     classes_refe_filename = conf.calc_marker_params.getpath("classes_refe_filename")
     input_groundtruth_filename = conf.calc_marker_params.getpath(
         "input_groundtruth_filename"
@@ -174,24 +173,22 @@ def calc_marker_task(config_paths: List[Path], default_basedir: Path):
 
     # STEP 2: Get the timeseries data needed for the classification
     # -------------------------------------------------------------
-    # Get the time series data (S1 and S2) to be used for the classification
+    # Get the time series data (eg. S1, S2,...) to be used for the classification
     # Result: data is put in files in timeseries_periodic_dir, in one file per
     #         date/period
     timeseries_periodic_dir = conf.dirs.getpath("timeseries_periodic_dir")
+    timeseries_periodic_dir = timeseries_periodic_dir / imagedata_input_parcel_path.stem
     start_date_str = conf.marker["start_date_str"]
     end_date_str = conf.marker["end_date_str"]
-    sensordata_to_use = conf.marker.getlist("sensordata_to_use")
+    sensordata_to_use = conf.parse_sensordata_to_use(conf.marker["sensordata_to_use"])
     parceldata_aggregations_to_use = conf.marker.getlist(
         "parceldata_aggregations_to_use"
     )
-    base_filename = f"{input_parcel_filename.stem}_bufm{buffer}_weekly"
     ts.calc_timeseries_data(
         input_parcel_path=imagedata_input_parcel_path,
-        input_country_code=country_code,
         start_date_str=start_date_str,
         end_date_str=end_date_str,
         sensordata_to_get=sensordata_to_use,
-        base_filename=base_filename,
         dest_data_dir=timeseries_periodic_dir,
     )
 
@@ -211,13 +208,13 @@ def calc_marker_task(config_paths: List[Path], default_basedir: Path):
     #             Is -1 if the parcel doesn't have any S1/S2 data.
     classtype_to_prepare = conf.preprocess["classtype_to_prepare"]
     parcel_path = run_dir / f"{input_parcel_filename.stem}_parcel{data_ext}"
-    parcel_pixcount_path = (
-        timeseries_periodic_dir / f"{base_filename}_pixcount{data_ext}"
-    )
+    base_filename = f"{input_parcel_filename.stem}_bufm{buffer}_weekly"
     class_pre.prepare_input(
         input_parcel_path=input_parcel_nogeo_path,
         input_parcel_filetype=input_parcel_filetype,
-        input_parcel_pixcount_path=parcel_pixcount_path,
+        timeseries_periodic_dir=timeseries_periodic_dir,
+        base_filename=base_filename,
+        data_ext=data_ext,
         classtype_to_prepare=classtype_to_prepare,
         classes_refe_path=classes_refe_path,
         output_parcel_path=parcel_path,
@@ -352,7 +349,9 @@ def calc_marker_task(config_paths: List[Path], default_basedir: Path):
         class_pre.prepare_input(
             input_parcel_path=input_groundtruth_path,
             input_parcel_filetype=input_parcel_filetype,
-            input_parcel_pixcount_path=parcel_pixcount_path,
+            timeseries_periodic_dir=timeseries_periodic_dir,
+            base_filename=base_filename,
+            data_ext=data_ext,
             classtype_to_prepare=conf.preprocess["classtype_to_prepare_groundtruth"],
             classes_refe_path=classes_refe_path,
             output_parcel_path=groundtruth_path,
diff --git a/cropclassification/calc_timeseries.py b/cropclassification/calc_timeseries.py
index e1792752..40376ebe 100644
--- a/cropclassification/calc_timeseries.py
+++ b/cropclassification/calc_timeseries.py
@@ -16,8 +16,7 @@
 
 from cropclassification.helpers import config_helper as conf
 from cropclassification.helpers import log_helper
-from cropclassification.helpers import raster_helper
-from cropclassification.preprocess import timeseries_calc_dias_onda_per_image as calc_ts
+from cropclassification.util import zonal_stats_bulk
 from cropclassification.preprocess import _timeseries_helper as ts_helper
 
 
@@ -130,7 +129,7 @@ def calc_timeseries_task(config_paths: List[Path], default_basedir: Path):
         input_image_path = Path(input_image_path)
         # Get more detailed info about the image
         try:
-            image_info = raster_helper.get_image_info(input_image_path)
+            image_info = zonal_stats_bulk.get_image_info(input_image_path)
         except Exception:
             # If not possible to get info for image, log and skip it
             logger.exception(f"SKIP image: error getting info for {input_image_path}")
@@ -139,7 +138,9 @@ def calc_timeseries_task(config_paths: List[Path], default_basedir: Path):
         # Fast-24h <- 2021-02-23 -> NRT-3H
         productTimelinessCategory = (
             "Fast-24h"
-            if dateutil.parser.isoparse(image_info["acquisition_date"]).timestamp()
+            if dateutil.parser.isoparse(
+                image_info.extra["acquisition_date"]
+            ).timestamp()
             < esaSwitchedProcessingMethod
             else "NRT-3h"
         )
@@ -147,13 +148,14 @@ def calc_timeseries_task(config_paths: List[Path], default_basedir: Path):
         # If sentinel1 and wrong productTimelinessCategory, skip: we only
         # want 1 type to evade images used twice
         if (
-            image_info["satellite"].startswith("S1")
-            and image_info["productTimelinessCategory"] != productTimelinessCategory
+            image_info.imagetype.startswith("S1")
+            and image_info.extra["productTimelinessCategory"]
+            != productTimelinessCategory
         ):
             logger.info(
                 f"SKIP image, productTimelinessCategory should be "
                 f"'{productTimelinessCategory}', but is: "
-                f"{image_info['productTimelinessCategory']} for "
+                f"{image_info.extra['productTimelinessCategory']} for "
                 f"{input_image_path}"
             )
             continue
@@ -171,11 +173,11 @@ def calc_timeseries_task(config_paths: List[Path], default_basedir: Path):
         )
 
     try:
-        calc_ts.calc_stats_per_image(
+        images_bands = [(path, ["VV", "VH"]) for path in input_image_paths]
+        zonal_stats_bulk.zonal_stats(
             features_path=input_features_path,
             id_column=conf.columns["id"],
-            image_paths=input_image_paths,
-            bands=["VV", "VH"],
+            images_bands=images_bands,
             output_dir=output_dir,
             temp_dir=temp_dir,
             log_dir=log_dir,
@@ -221,7 +223,7 @@ def calc_timeseries_task(config_paths: List[Path], default_basedir: Path):
 
         # Get more detailed info about the image
         try:
-            image_info = raster_helper.get_image_info(input_image_path)
+            image_info = zonal_stats_bulk.get_image_info(input_image_path)
         except Exception:
             # If not possible to get info for image, log and skip it
             logger.exception(f"SKIP image: error getting info for {input_image_path}")
@@ -230,12 +232,12 @@ def calc_timeseries_task(config_paths: List[Path], default_basedir: Path):
         # If sentinel2 and cloud coverage too high... skip
         if (
             max_cloudcover_pct >= 0
-            and image_info["satellite"].startswith("S2")
-            and image_info["Cloud_Coverage_Assessment"] > max_cloudcover_pct
+            and image_info.imagetype.startswith("S2")
+            and image_info.extra["Cloud_Coverage_Assessment"] > max_cloudcover_pct
         ):
             logger.info(
                 "SKIP image, Cloud_Coverage_Assessment: "
-                f"{image_info['Cloud_Coverage_Assessment']:0.2f} > "
+                f"{image_info.extra['Cloud_Coverage_Assessment']:0.2f} > "
                 f"{max_cloudcover_pct} for {input_image_path}"
             )
             continue
@@ -243,11 +245,12 @@ def calc_timeseries_task(config_paths: List[Path], default_basedir: Path):
         tmp_input_image_paths.append(input_image_path)
 
     try:
-        calc_ts.calc_stats_per_image(
+        bands = conf.timeseries.getlist("s2bands")
+        images_bands = [(path, bands) for path in input_image_paths]
+        zonal_stats_bulk.zonal_stats(
             features_path=input_features_path,
             id_column=conf.columns["id"],
-            image_paths=input_image_paths,
-            bands=conf.timeseries.getlist("s2bands"),
+            images_bands=images_bands,
             output_dir=output_dir,
             temp_dir=temp_dir,
             log_dir=log_dir,
@@ -287,11 +290,11 @@ def calc_timeseries_task(config_paths: List[Path], default_basedir: Path):
     ]
 
     try:
-        calc_ts.calc_stats_per_image(
+        images_bands = [(path, ["VV", "VH"]) for path in input_image_paths]
+        zonal_stats_bulk.zonal_stats(
             features_path=input_features_path,
             id_column=conf.columns["id"],
-            image_paths=input_image_paths,
-            bands=["VV", "VH"],
+            images_bands=images_bands,
             output_dir=output_dir,
             temp_dir=temp_dir,
             log_dir=log_dir,
diff --git a/cropclassification/general.ini b/cropclassification/general.ini
index 6ac3c268..6eba7c05 100644
--- a/cropclassification/general.ini
+++ b/cropclassification/general.ini
@@ -52,31 +52,6 @@ geofile_ext = .gpkg
 # The log level to use
 log_level = INFO
 
-# Constants for types of aggregation to use
-# Mean value of the pixels values in a parcel.
-PARCELDATA_AGGRAGATION_MEAN = mean
-# Median value of the pixels values in a parcel.
-PARCELDATA_AGGRAGATION_MEDIAN = median
-# std dev of the values of the pixels in a parcel
-PARCELDATA_AGGRAGATION_STDDEV = std
-
-# Constants for types of sensor data
-SENSORDATA_S1 = S1
-# Sentinel 1 data
-# Sentinel 1 data, in dB
-SENSORDATA_S1DB = S1dB
-# Sentinel 1 data, divided in Ascending and Descending passes
-SENSORDATA_S1_ASCDESC = S1AscDesc
-# Sentinel 1 data, in dB, divided in Ascending and Descending passes
-SENSORDATA_S1DB_ASCDESC = S1dBAscDesc
-# Sentinel 2 data
-SENSORDATA_S2 = S2
-# Sentinel 2 data (B2,B3,B4,B8) IF available for 95% or area
-SENSORDATA_S2gt95 = S2gt95
-
-#SENSORDATA_S1_COHERENCE
-SENSORDATA_S1_COHERENCE = S1Coh
-
 [marker]
 # markertype, must be overriden in child ini files
 markertype = MUST_OVERRIDE
@@ -123,29 +98,40 @@ classes_to_ignore_unimportant_default = IGNORE_UNIMPORTANT
 #   * BALANCING_STRATEGY_EQUAL: for each input class, the same amount of samples is used as training. For classes with few samples, (samples will be duplicated if needed)
 balancing_strategy = BALANCING_STRATEGY_MEDIUM2
 
-# The sensor data to be used for this marker
-sensordata_to_use = ${general:SENSORDATA_S1_ASCDESC}, ${general:SENSORDATA_S2gt95}
-# The aggregation type to use on parcel level
-parceldata_aggregations_to_use = ${general:PARCELDATA_AGGRAGATION_MEAN}
+# Sensordata to use for the markers. This is a list of ImageProfile names:
+#   - S1dB: Sentinel 1 data, in dB
+#   - S1ascDesc: Sentinel 1 data, divided in Ascending and Descending passes
+#   - S1dBAscDesc: Sentinel 1 data, in dB, divided in Ascending and Descending passes
+#   - S1Coh: Sentinel 1 coherence
+#   - S2gt95: Sentinel 2 data (B2,B3,B4,B8) IF available for 95% or area
+#
+#   - s2-landcover: Sentinel 2 landcover, based on SCL mask
+#   - s2-ndvi: Sentinel 2 ndvi
+#   - s2-agri: Sentinel 2 bands relevant for agriculture (B02, B03, B04, B08, B11, B12)
+#   - s1-grd-sigma0-asc: Sentinel 1 GRD, ascending orbit (VV, VH, angle)
+#   - s1-grd-sigma0-desc: Sentinel 1 GRD, descending orbit (VV, VH, angle)
+#   - s1-coh: Sentinel 1, 6 day or 12 day coherence, depending on availability (VV, VH)
+#
+# Optionally, it is possible to specify the bands of the ImageProfile to be used like this:
+#   eg. [..., {"s2-grd-sigma0-asc": ["VV", "VH"]}]
+sensordata_to_use = ["s2-agri", {"s1-grd-sigma0-asc": ["VV", "VH"]}, {"s1-grd-sigma0-desc": ["VV", "VH"]}]
+
+# The aggregation type to use on parcel level. Following options are avalable:
+#    - mean: the mean of the pixel values in a parcel
+#    - median: the median of the pixel values in a parcel
+#    - std: the standard deviation of the pixel values in a parcel
+parceldata_aggregations_to_use = median
 
 # Postprocess...
 postprocess_to_groups = 
 
 [timeseries]
-# The way timeseries are calculated
-#     * gee: google earth engine
-#     * onda: ONDA DIAS
-timeseries_calc_type = onda
-
 # The maximum percentage cloudcover an (S2) image can have to be used
 max_cloudcover_pct = 15
 
 # The min percentage of parcels that need to have valid data for a time+sensor to use it 
 min_parcels_with_data_pct = 90
 
-# Bands to extract for S2 images 
-s2bands = B02-10m, B03-10m, B04-10m, B08-10m, B11-20m, B12-20m, SCL-20m
-
 [preprocess]
 # 
 dedicated_data_columns = ${columns:id}, ${columns:class}, ${columns:class}_orig, ${columns:class_declared}, ${columns:class_declared}_orig, ${columns:pixcount_s1s2}
@@ -242,11 +228,13 @@ prediction_conclusion_cons = pred_conclusion_cons
 # data_dir can be specified as a relative path which will be resolved to
 # the parent dir of the tasks dir 
 data_dir = .
-temp_dir = ${dirs:data_dir}/tmp
+# temp_dir can contain placeholder {tmp_dir}: resolves to the default system tempdir
+temp_dir = {tmp_dir}/cropclassification
 # marker_basedir can be specified as a relative path which will be resolved to
 # the parent dir of the tasks dir
 marker_basedir = .
 marker_dir = ${dirs:marker_basedir}/${marker:year}_${marker:markertype}
+images_periodic_dir = ${dirs:data_dir}/_images_periodic/${calc_marker_params:country_code}
 timeseries_periodic_dir = ${dirs:data_dir}/_ts_periodic
 timeseries_per_image_dir = ${dirs:data_dir}/_ts_per_image
 input_dir = ${dirs:data_dir}/_inputdata
diff --git a/cropclassification/helpers/config_helper.py b/cropclassification/helpers/config_helper.py
index 415f8c5b..4090cd42 100644
--- a/cropclassification/helpers/config_helper.py
+++ b/cropclassification/helpers/config_helper.py
@@ -7,13 +7,34 @@
 import json
 from pathlib import Path
 import pprint
-from typing import List, Optional
+import tempfile
+from typing import Dict, List, Optional
+
+from cropclassification.util.openeo_util import ImageProfile
 
 # -------------------------------------------------------------
 # The real work
 # -------------------------------------------------------------
 
 
+class SensorData:
+    def __init__(
+        self,
+        imageprofile_name: str,
+        imageprofile: Optional[ImageProfile] = None,
+        bands: Optional[List[str]] = None,
+    ):
+        self.imageprofile_name = imageprofile_name
+        if imageprofile is not None:
+            self.imageprofile = imageprofile
+        else:
+            self.imageprofile = image_profiles[imageprofile_name]
+        if bands is not None:
+            self.bands = bands
+        else:
+            self.bands = self.imageprofile.bands
+
+
 def read_config(config_paths: List[Path], default_basedir: Optional[Path] = None):
     # Read the configuration
     global config
@@ -68,6 +89,10 @@ def read_config(config_paths: List[Path], default_basedir: Optional[Path] = None
         )
         config["dirs"]["marker_basedir"] = marker_basedir_absolute.as_posix()
 
+    # Fill out placeholder in the temp_dir (if it is there)
+    tmp_dir_str = tempfile.gettempdir()
+    config["dirs"]["temp_dir"] = config["dirs"]["temp_dir"].format(tmp_dir=tmp_dir_str)
+
     global config_paths_used
     config_paths_used = config_paths
 
@@ -93,6 +118,96 @@ def read_config(config_paths: List[Path], default_basedir: Optional[Path] = None
     global dirs
     dirs = config["dirs"]
 
+    global image_profiles
+    image_profiles = _get_raster_profiles()
+
+
+def parse_sensordata_to_use(input) -> Dict[str, SensorData]:
+    result = None
+    sensordata_parsed = None
+    try:
+        sensordata_parsed = json.loads(input)
+    except Exception:
+        pass
+
+    if sensordata_parsed is not None:
+        # It was a json object, so parse as such
+        result = {}
+        for imageprofile in sensordata_parsed:
+            if isinstance(imageprofile, str):
+                result[imageprofile] = SensorData(imageprofile)
+            elif isinstance(imageprofile, dict):
+                if len(imageprofile) != 1:
+                    raise ValueError(
+                        "invalid sensordata_to_use: this should be a single key dict: "
+                        f"{imageprofile}"
+                    )
+                imageprofile_name = list(imageprofile.keys())[0]
+                bands = list(imageprofile.values())[0]
+                result[imageprofile_name] = SensorData(imageprofile_name, bands=bands)
+            else:
+                raise ValueError(
+                    "invalid sensordata_to_use: only str or dict elements allowed, "
+                    f"not: {imageprofile}"
+                )
+    else:
+        # It was no json object, so it must be a list
+        result = {i.strip(): SensorData(i.strip()) for i in input.split(",")}
+
+    return result
+
+
+def _get_raster_profiles() -> Dict[str, ImageProfile]:
+    # TODO: this should move to a config file
+    profiles = {}
+    profiles["s2-agri"] = ImageProfile(
+        name="s2-agri",
+        satellite="s2",
+        collection="TERRASCOPE_S2_TOC_V2",
+        bands=["B02", "B03", "B04", "B08", "B11", "B12"],
+        # Use the "min" reducer filters out "lightly clouded areas"
+        process_options={"time_dimension_reducer": "min", "cloud_filter_band": "SCL"},
+    )
+    profiles["s2-ndvi"] = ImageProfile(
+        name="s2-ndvi",
+        satellite="s2",
+        collection="TERRASCOPE_S2_NDVI_V2",
+        bands=["NDVI"],
+        process_options={
+            "time_dimension_reducer": "max",
+            "cloud_filter_band": "SCENECLASSIFICATION_20M",
+        },
+    )
+    profiles["s1-grd-sigma0-asc"] = ImageProfile(
+        name="s1-grd-sigma0-asc",
+        satellite="s1",
+        collection="S1_GRD_SIGMA0_ASCENDING",
+        bands=["VV", "VH", "angle"],
+        process_options={
+            "time_dimension_reducer": "min",
+        },
+    )
+    profiles["s1-grd-sigma0-desc"] = ImageProfile(
+        name="s1-grd-sigma0-desc",
+        satellite="s1",
+        collection="S1_GRD_SIGMA0_DESCENDING",
+        bands=["VV", "VH", "angle"],
+        process_options={
+            "time_dimension_reducer": "min",
+        },
+    )
+    profiles["s1-coh"] = ImageProfile(
+        name="s1-coh",
+        satellite="s1",
+        collection="TERRASCOPE_S1_SLC_COHERENCE_V1",
+        bands=["VV", "VH"],
+        process_options={
+            "time_dimension_reducer": "min",
+        },
+    )
+
+    return profiles
+
 
 def pformat_config():
     message = f"Config files used: {pprint.pformat(config_paths_used)} \n"
diff --git a/cropclassification/helpers/pandas_helper.py b/cropclassification/helpers/pandas_helper.py
index 40631b5a..a89b7613 100644
--- a/cropclassification/helpers/pandas_helper.py
+++ b/cropclassification/helpers/pandas_helper.py
@@ -5,7 +5,7 @@
 
 from pathlib import Path
 import os
-from typing import List, Optional
+from typing import List, Optional, Union
 
 import pandas as pd
 import sqlite3
@@ -78,7 +78,7 @@ def read_file(
 
 
 def to_file(
-    df: pd.DataFrame,
+    df: Union[pd.DataFrame, pd.Series],
     path: Path,
     table_name: str = "info",
     index: bool = True,
diff --git a/cropclassification/helpers/raster_helper.py b/cropclassification/helpers/raster_helper.py
deleted file mode 100644
index 4ac9448f..00000000
--- a/cropclassification/helpers/raster_helper.py
+++ /dev/null
@@ -1,656 +0,0 @@
-import logging
-import os
-from pathlib import Path
-import time
-from typing import List, Tuple
-import xml.etree.ElementTree as ET
-
-import rasterio
-from shapely.geometry import polygon as sh_polygon
-
-from cropclassification.util import io_util
-
-# General init
-logger = logging.getLogger(__name__)
-
-
-def get_image_data(
-    image_path: Path,
-    bounds: Tuple[float, float, float, float],
-    bands: List[str],
-    pixel_buffer: int = 0,
-) -> dict:
-    """
-    Reads the data from the image.
-
-    Adds a small buffer around the bounds asked to evade possible rounding issues.
-
-    Returns a dict of the following structure:
-    imagedata[band]['data']: the data read as numpy array
-                   ['transform']: the Affine transform of the band read
-
-    Args
-        image_path: the path the image
-        bounds: the bounds to be read, in coordinates in the projection of the image
-        bands: list of bands to be read, eg. "VV", "VH",...
-        pixel_buffer: number to pixels to take as buffer around the bounds provided in
-            pixels
-    """
-    logger.info("get_image_data start")
-    # Get info about the image
-    image_info = get_image_info(image_path)
-
-    # Now get the data
-    image_data = {}  # Dict for the transforms and the data per band
-    if image_info["image_type"] in ("CARD", "SAFE"):
-        # Loop over bands and get data
-        for band in bands:
-            band_relative_path = image_info["bands"][band]["relative_path"]
-            image_band_path = image_path / band_relative_path
-            band_index = image_info["bands"][band]["bandindex"]
-            logger.info(
-                f"Read image data from {image_band_path}, with band_index: {band_index}"
-            )
-            image_data[band] = {}
-            with rasterio.open(str(image_band_path)) as src:
-                # Determine the window we need to read from the image:
-                window_to_read = projected_bounds_to_window(
-                    bounds, src.transform, src.width, src.height, pixel_buffer
-                )
-                image_data[band][
-                    "transform"
-                ] = rasterio.windows.transform(  # type: ignore
-                    window_to_read, src.transform
-                )
-                # Read!
-                # Remark: bandindex in rasterio is 1-based instead of 0-based -> +1
-                logger.debug(f"Read image data from {image_band_path}")
-                image_data[band]["data"] = src.read(
-                    band_index + 1, window=window_to_read
-                )
-                logger.info("Image data read")
-    else:
-        message = f"Format currently not supported: {image_path}"
-        logger.error(message)
-        raise NotImplementedError(message)
-
-    logger.info("get_image_data ready")
-
-    return image_data
-
-
-def get_image_info(image_path: Path) -> dict:
-    image_info = {}
-
-    # First determine general info of image
-    image_path_noext, image_ext = os.path.splitext(image_path)
-    _, image_basename_noext = os.path.split(image_path_noext)
-    image_basename_noext_split = image_basename_noext.split("_")
-    image_info["satellite"] = image_basename_noext_split[0]
-
-    # Specific code per image type
-    if image_ext.upper() == ".CARD":
-        # This is a sentinel 1 image (GRD or coherence)
-        # First extract and fill out some basic info
-        image_info["image_type"] = "CARD"
-        image_info["image_id"] = image_basename_noext
-
-        # Read info from the metadata file
-        metadata_xml_path = image_path / "metadata.xml"
-        metadata = ET.parse(str(metadata_xml_path))
-        metadata_root = metadata.getroot()
-
-        logger.debug(f"Parse metadata info from {metadata_xml_path}")
-
-        try:
-            # Get the filename
-            image_info["filename"] = metadata_root.find(
-                "productFilename"
-            ).text  # type: ignore
-
-            # Get the footprint
-            image_info["footprint"] = {}
-            footprint = metadata_root.find("footprint")
-            assert footprint is not None
-            poly = footprint.find("{http://www.opengis.net/gml}Polygon")
-            assert poly is not None
-            image_info["footprint"]["srsname"] = poly.attrib.get("srsName")
-            linear_ring = []
-            # for coord in poly.findall(
-            #   "{http://www.opengis.net/gml}outerBoundaryIs/"
-            #   "{http://www.opengis.net/gml}LinearRing/"
-            #   "{http://www.opengis.net/gml}coord"
-            # ):
-            #    linear_ring.append(
-            #        (float(coord.findtext(
-            #            "{http://www.opengis.net/gml}X"
-            #        )),
-            #        float(coord.findtext("{http://www.opengis.net/gml}Y")))
-            #    )
-            coord_str = poly.find(
-                "{http://www.opengis.net/gml}outerBoundaryIs/"
-                "{http://www.opengis.net/gml}LinearRing/"
-                "{http://www.opengis.net/gml}coordinates"
-            )
-            assert coord_str is not None and coord_str.text is not None
-            logger.debug(f"coord_str: {coord_str}, coord_str.text: {coord_str.text}")
-            coord_list = coord_str.text.split(" ")
-            for coord in coord_list:
-                # Watch out, latitude (~y) is first, than longitude (~x)
-                # TODO: add check if the projection is in degrees (latlon) or
-                # coordinates (xy) instead of hardcoded latlon
-                y, x = coord.split(",")
-                linear_ring.append((float(x), float(y)))
-            image_info["footprint"]["shape"] = sh_polygon.Polygon(linear_ring)
-
-            # get epsg
-            epsg = metadata_root.find("imageProjection/EPSG")
-            assert epsg is not None and epsg.text is not None
-            image_info["image_epsg"] = float(epsg.text)
-            image_info["image_crs"] = f"EPSG:{epsg.text}"
-        except Exception as ex:
-            raise Exception(
-                f"Exception extracting info from {metadata_xml_path}"
-            ) from ex
-
-        # Read info from the manifest.safe file
-        manifest_xml_searchstring = "*_manifest.safe"
-        manifest_xml_paths = list(image_path.glob(manifest_xml_searchstring))
-
-        # The number of .safe indicates whether it is a GRD or a Coherence image
-        nb_safefiles = len(manifest_xml_paths)
-        if nb_safefiles == 1:
-            # Now parse the .safe file
-            manifest_xml_path = manifest_xml_paths[0]
-
-            try:
-                manifest = ET.parse(str(manifest_xml_path))
-                manifest_root = manifest.getroot()
-
-                # Define namespaces...
-                ns = {
-                    "safe": "http://www.esa.int/safe/sentinel-1.0",
-                    "s1": "http://www.esa.int/safe/sentinel-1.0/sentinel-1",
-                    "s1sarl1": (
-                        "http://www.esa.int/safe/sentinel-1.0/sentinel-1/sar/level-1"
-                    ),
-                }
-
-                logger.debug(f"Parse manifest info from {metadata_xml_path}")
-                image_info["transmitter_receiver_polarisation"] = []
-                for polarisation in manifest_root.findall(
-                    "metadataSection/metadataObject/metadataWrap/xmlData/"
-                    "s1sarl1:standAloneProductInformation/"
-                    "s1sarl1:transmitterReceiverPolarisation",
-                    ns,
-                ):
-                    image_info["transmitter_receiver_polarisation"].append(
-                        polarisation.text
-                    )
-                image_info["productTimelinessCategory"] = manifest_root.find(
-                    "metadataSection/metadataObject/metadataWrap/xmlData/"
-                    "s1sarl1:standAloneProductInformation/"
-                    "s1sarl1:productTimelinessCategory",
-                    ns,
-                ).text  # type: ignore
-
-                image_info["instrument_mode"] = manifest_root.find(
-                    "metadataSection/metadataObject/metadataWrap/xmlData/"
-                    "safe:platform/safe:instrument/safe:extension/"
-                    "s1sarl1:instrumentMode/s1sarl1:mode",
-                    ns,
-                ).text  # type: ignore
-                image_info["orbit_properties_pass"] = manifest_root.find(
-                    "metadataSection/metadataObject/metadataWrap/xmlData/"
-                    "safe:orbitReference/safe:extension/s1:orbitProperties/s1:pass",
-                    ns,
-                ).text  # type: ignore
-                image_info["acquisition_date"] = manifest_root.find(
-                    "metadataSection/metadataObject/metadataWrap/xmlData/"
-                    "safe:acquisitionPeriod/safe:startTime",
-                    ns,
-                ).text  # type: ignore
-
-                # Now have a look in the files themselves to get band info,...
-                # TODO: probably cleaner/easier to read from metadata files...
-                image_basename_noext_nodate = image_basename_noext[
-                    : image_basename_noext.rfind("_", 0)
-                ]
-                image_datadirname = f"{image_basename_noext_nodate}.data"
-                image_datadir = image_path / image_datadirname
-                band_paths = list(image_datadir.glob("*.img"))
-
-                # If no files were found, check if the images are in .tif format!
-                if len(band_paths) == 0:
-                    tif_path = image_path / "Gamma0_VH.tif"
-                    if tif_path.exists():
-                        band_paths.append(tif_path)
-                    tif_path = image_path / "Gamma0_VV.tif"
-                    if tif_path.exists():
-                        band_paths.append(tif_path)
-
-                # If no files were found, error
-                if len(band_paths) == 0:
-                    message = (
-                        f"No image files found in {image_datadir} or {image_path}!"
-                    )
-                    logger.error(message)
-                    raise Exception(message)
-
-                image_info["bands"] = {}
-                for i, band_path in enumerate(band_paths):
-                    # Extract bound,... info from the first file only
-                    # (they are all the same)
-                    if i == 0:
-                        # logger.debug(f"Read image metadata from {band_path}")
-                        with rasterio.open(str(band_path)) as src:
-                            image_info["image_bounds"] = src.bounds
-                            image_info["image_affine"] = src.transform
-                            if src.crs is not None:
-                                image_info["image_crs"] = src.crs.to_string()
-                                image_info["image_epsg"] = src.crs.to_epsg()
-                        # logger.debug(f"Image metadata read: {image_info}")
-                    if band_path.stem == "Gamma0_VH":
-                        band = "VH"
-                    elif band_path.stem == "Gamma0_VV":
-                        band = "VV"
-                    else:
-                        raise NotImplementedError(
-                            f"Filename not supported: {band_path}"
-                        )
-
-                    # Add specific info about the band
-                    image_info["bands"][band] = {}
-                    image_info["bands"][band]["path"] = band_path
-                    image_info["bands"][band]["relative_path"] = str(band_path).replace(
-                        str(image_path) + os.pathsep, ""
-                    )
-                    image_info["bands"][band]["filename"] = band_path.name
-                    image_info["bands"][band]["bandindex"] = 0
-
-            except Exception as ex:
-                raise Exception(
-                    f"Exception extracting info from {manifest_xml_path}"
-                ) from ex
-
-        elif nb_safefiles == 2 or nb_safefiles == 3:
-            # 2 safe files -> coherence
-            # Now parse the first .safe file
-            # TODO: maybe check if the info in all safe files  are the same or?.?
-            manifest_xml_path = manifest_xml_paths[0]
-
-            try:
-                manifest = ET.parse(str(manifest_xml_path))
-                manifest_root = manifest.getroot()
-
-                # Define namespaces...
-                ns = {
-                    "safe": "http://www.esa.int/safe/sentinel-1.0",
-                    "s1": "http://www.esa.int/safe/sentinel-1.0/sentinel-1",
-                    "s1sarl1": (
-                        "http://www.esa.int/safe/sentinel-1.0/sentinel-1/sar/level-1"
-                    ),
-                }
-
-                # logger.debug(f"Parse manifest info from {metadata_xml_path}")
-                image_info["transmitter_receiver_polarisation"] = []
-                for polarisation in manifest_root.findall(
-                    "metadataSection/metadataObject/metadataWrap/xmlData/"
-                    "s1sarl1:standAloneProductInformation/"
-                    "s1sarl1:transmitterReceiverPolarisation",
-                    ns,
-                ):
-                    image_info["transmitter_receiver_polarisation"].append(
-                        polarisation.text
-                    )
-                image_info["productTimelinessCategory"] = manifest_root.find(
-                    "metadataSection/metadataObject/metadataWrap/xmlData/"
-                    "s1sarl1:standAloneProductInformation/"
-                    "s1sarl1:productTimelinessCategory",
-                    ns,
-                ).text  # type: ignore
-                image_info["instrument_mode"] = manifest_root.find(
-                    "metadataSection/metadataObject/metadataWrap/xmlData/"
-                    "safe:platform/safe:instrument/safe:extension/"
-                    "s1sarl1:instrumentMode/s1sarl1:mode",
-                    ns,
-                ).text  # type: ignore
-                image_info["orbit_properties_pass"] = manifest_root.find(
-                    "metadataSection/metadataObject/metadataWrap/xmlData/"
-                    "safe:orbitReference/safe:extension/s1:orbitProperties/s1:pass",
-                    ns,
-                ).text  # type: ignore
-                image_info["acquisition_date"] = manifest_root.find(
-                    "metadataSection/metadataObject/metadataWrap/xmlData/"
-                    "safe:acquisitionPeriod/safe:startTime",
-                    ns,
-                ).text  # type: ignore
-
-                # Now have a look in the files themselves to get band info,...
-                # TODO: probably cleaner/easier to read from metadata files...
-                # For coherence filename, remove extra .LCO1 extension + date
-                # image_datafilebasename, _ = os.path.splitext(image_basename_noext)
-                # image_datafilebasename = image_datafilebasename[
-                #     :image_datafilebasename.rindex('_')
-                # ]
-
-                metadata = ET.parse(str(image_path / "metadata.xml"))
-                metadata_root = metadata.getroot()
-
-                image_basename_noext_nodate = image_basename_noext[
-                    : image_basename_noext.rfind("_", 0)
-                ]
-                image_datafilename = f"{image_basename_noext_nodate}_byte.tif"
-                image_datapath = image_path / image_datafilename
-
-                # logger.debug(f"Read image metadata from {image_datapath}")
-                with rasterio.open(str(image_datapath)) as src:
-                    image_info["image_bounds"] = src.bounds
-                    image_info["image_affine"] = src.transform
-                    if src.crs is not None:
-                        image_info["image_crs"] = src.crs.to_string()
-                        image_info["image_epsg"] = src.crs.to_epsg()
-
-                # Add specific info about the bands
-                image_info["bands"] = {}
-                band = "VH"
-                image_info["bands"][band] = {}
-                image_info["bands"][band]["path"] = image_datapath
-                image_info["bands"][band]["relative_path"] = image_datafilename
-                image_info["bands"][band]["filename"] = image_datafilename
-                image_info["bands"][band]["bandindex"] = 0
-
-                band = "VV"
-                image_info["bands"][band] = {}
-                image_info["bands"][band]["path"] = image_datapath
-                image_info["bands"][band]["relative_path"] = image_datafilename
-                image_info["bands"][band]["filename"] = image_datafilename
-                image_info["bands"][band]["bandindex"] = 1
-
-            except Exception as ex:
-                raise Exception(
-                    f"Exception extracting info from {manifest_xml_path}"
-                ) from ex
-        else:
-            message = (
-                f"Error: found {nb_safefiles} .safe files doing glob in {image_path} "
-                f"with {manifest_xml_searchstring}"
-            )
-            logger.error(message)
-            raise Exception(message)
-
-    elif image_ext.upper() == ".SAFE":
-        # This is a level 2 sentinel 2 file
-
-        # First extract and fill out some basic info
-        image_info["image_type"] = "SAFE"
-        _, image_basename_noext = os.path.split(image_path_noext)
-        image_info["image_id"] = image_basename_noext
-
-        # Read info from the manifest.safe file
-        """
-        manifest_xml_searchstring = os.path.join(image_path, f"*_manifest.safe")
-        manifest_xml_paths = glob.glob(manifest_xml_searchstring)
-
-        # We should only find one .safe file (at the moment)
-        nb_safefiles = len(manifest_xml_paths)
-        if nb_safefiles == 1:
-
-        # Now parse the metadata xml file
-        manifest_xml_path = manifest_xml_paths[0]
-        """
-        metadata_xml_path = image_path / "MTD_MSIL2A.xml"
-        try:
-            metadata = ET.parse(str(metadata_xml_path))
-            metadata_root = metadata.getroot()
-
-            # Define namespaces...
-            # xsi:schemaLocation=
-            #     "https://psd-14.sentinel2.eo.esa.int/PSD/User_Product_Level-2A.xsd"
-            ns = {
-                "n1": (
-                    "https://psd-14.sentinel2.eo.esa.int/PSD/User_Product_Level-2A.xsd"
-                ),
-                "xsi": "http://www.w3.org/2001/XMLSchema-instance",
-            }
-
-            # logger.debug(f"Parse metadata info from {metadata_xml_path}")
-            image_info["Cloud_Coverage_Assessment"] = float(
-                metadata_root.find(
-                    "n1:Quality_Indicators_Info/Cloud_Coverage_Assessment", ns
-                ).text  # type: ignore
-            )
-            image_info["acquisition_date"] = metadata_root.find(
-                "n1:General_Info/Product_Info/PRODUCT_START_TIME", ns
-            ).text  # type: ignore
-
-        except Exception as ex:
-            raise Exception(
-                f"Exception extracting info from {metadata_xml_path}"
-            ) from ex
-
-        # Now have a look in the files themselves to get band info,...
-        image_datadir = image_path / "GRANULE"
-        band_paths = list(image_datadir.rglob("*.jp2"))
-
-        metadata_path = list(image_datadir.rglob("MTD_TL.xml"))[0]
-        try:
-            # read epsg
-            metadata = ET.parse(str(metadata_path))
-            metadata_root = metadata.getroot()
-
-            ns = {
-                "n1": (
-                    "https://psd-14.sentinel2.eo.esa.int/PSD/"
-                    "S2_PDI_Level-2A_Tile_Metadata.xsd"
-                ),
-                "xsi": "http://www.w3.org/2001/XMLSchema-instance",
-            }
-
-            image_info["image_crs"] = metadata_root.find(
-                "n1:Geometric_Info/Tile_Geocoding/HORIZONTAL_CS_CODE", ns
-            ).text  # type: ignore
-            image_info["image_epsg"] = float(
-                image_info["image_crs"].replace("EPSG:", "")  # type: ignore
-            )
-        except Exception as ex:
-            raise Exception(f"Exception extracting info from {metadata_path}") from ex
-
-        # If no files were found, error!
-        if len(band_paths) == 0:
-            message = f"No image files found in {image_datadir}!"
-            logger.error(message)
-            raise Exception(message)
-
-        image_info["bands"] = {}
-        for i, band_path in enumerate(band_paths):
-            band_filename = os.path.basename(band_path)
-            band_filename_noext, _ = os.path.splitext(band_filename)
-            band_filename_noext_split = band_filename_noext.split("_")
-
-            if len(band_filename_noext_split) == 4:
-                # IMG_DATA files
-                band = f"{band_filename_noext_split[2]}-{band_filename_noext_split[3]}"
-            elif len(band_filename_noext_split) == 5:
-                # IMG_DATA files
-                band = f"{band_filename_noext_split[3]}-{band_filename_noext_split[4]}"
-            elif len(band_filename_noext_split) == 3:
-                # QI_DATA files
-                band = band_filename_noext
-            else:
-                message = f"Filename of band doesn't have supported format: {band_path}"
-                logger.error(message)
-                raise Exception(message)
-
-            # Add specific info about the band
-            image_info["bands"][band] = {}
-            image_info["bands"][band]["path"] = band_path
-            image_info["bands"][band]["relative_path"] = str(band_path).replace(
-                str(image_path), ""
-            )[1:]
-            image_info["bands"][band]["filename"] = band_filename
-            image_info["bands"][band]["bandindex"] = 0
-
-            # Extract bound,... info
-            logger.debug(f"Read image metadata from {band_path}")
-            with rasterio.open(str(band_path)) as src:
-                image_info["bands"][band]["bounds"] = src.bounds
-                image_info["bands"][band]["affine"] = src.transform
-                if src.crs is None:
-                    image_info["bands"][band]["crs"] = image_info["image_crs"]
-                    image_info["bands"][band]["epsg"] = image_info["image_epsg"]
-                else:
-                    image_info["bands"][band]["crs"] = src.crs.to_string()
-                    image_info["bands"][band]["epsg"] = src.crs.to_epsg()
-
-            # Store the crs also on image level + check if all bands have the same crs
-            if i == 0:
-                image_info["image_bounds"] = image_info["bands"][band]["bounds"]
-                # image_info['image_crs'] = image_info['bands'][band]['crs']
-                # image_info['image_epsg'] = image_info['bands'][band]['epsg']
-            else:
-                if (
-                    image_info["image_crs"] != image_info["bands"][band]["crs"]
-                    or image_info["image_epsg"] != image_info["bands"][band]["epsg"]
-                ):
-                    message = f"Not all bands have the same crs for {image_info}"
-                    logger.error(message)
-                    raise Exception(message)
-
-        logger.debug(f"Image metadata read: {image_info}")
-
-    else:
-        message = f"Not a supported image format: {image_path}"
-        logger.error(message)
-        raise NotImplementedError(message)
-
-    return image_info
-
-
-def projected_bounds_to_window(
-    projected_bounds,
-    image_transform,
-    image_pixel_width: int,
-    image_pixel_height: int,
-    pixel_buffer: int = 0,
-):
-    """
-    Returns a rasterio.windows.Window to be used in rasterio to read the part of the
-    image specified.
-
-    Args
-        projected_bounds: bounds to created the window from, in projected coordinates
-        image_transform: Affine transform of the image you want to create the pixel
-            window for
-        image_pixel_width: total width of the image you want to create the pixel window
-            for, in pixels
-        image_pixel_height: total height of the image you want to create the pixel
-            window for, in pixels
-        pixel_buffer: number to pixels to take as buffer around the bounds provided in
-            pixels
-    """
-    # Take bounds of the features + convert to image pixels
-    xmin, ymin, xmax, ymax = projected_bounds
-    window_to_read_raw = rasterio.windows.from_bounds(  # type: ignore
-        xmin, ymin, xmax, ymax, image_transform
-    )
-
-    # Round so it only increases window size
-    window_to_read = window_to_read_raw.round_offsets("floor").round_lengths("ceil")
-
-    # Now some general math on window properties, but as they are readonly properties,
-    # work on copy
-    col_off, row_off, width, height = window_to_read.flatten()
-    # Add buffer of 1 pixel extra around
-    col_off -= pixel_buffer
-    row_off -= pixel_buffer
-    width += 2 * pixel_buffer
-    height += 2 * pixel_buffer
-
-    # Make sure the offsets aren't negative, as the pixels that are 'read' there acually
-    # get some value instead of eg. nan...!
-    if col_off < 0:
-        width -= abs(col_off)
-        col_off = 0
-    if row_off < 0:
-        height -= abs(row_off)
-        row_off = 0
-
-    # Make sure there won't be extra pixels to the top and right that will be read
-    if (col_off + width) > image_pixel_width:
-        width = image_pixel_width - col_off
-    if (row_off + height) > image_pixel_height:
-        height = image_pixel_height - row_off
-
-    # Ready... prepare to return...
-    window_to_read = rasterio.windows.Window(  # type: ignore
-        col_off, row_off, width, height
-    )
-
-    # Debug info
-    """
-    bounds_to_read = rasterio.windows.bounds(window_to_read, image_transform)
-    logger.debug(f"projected_bounds: {projected_bounds}, "
-        f"window_to_read_raw: {window_to_read_raw}, window_to_read: {window_to_read}, "
-        f"image_pixel_width: {image_pixel_width}, "
-        f"image_pixel_height: {image_pixel_height}, "
-        f"file transform: {image_transform}, bounds_to_read: {bounds_to_read}"
-    )
-    """
-
-    return window_to_read
-
-
-def prepare_image(image_path: Path, temp_dir: Path) -> Path:
-    """
-    Prepares the input image for usages.
-
-    In case of a zip file, the file is unzipped to the temp dir specified.
-
-    Returns the path to the prepared file/directory.
-    """
-
-    # If the input path is not a zip file, don't make local copy + just return path
-    if image_path.suffix.lower() != ".zip":
-        return image_path
-    else:
-        # It is a zip file, so it needs to be unzipped first...
-        # Create destination file path
-        image_basename_withzipext = os.path.basename(image_path)
-        image_basename = os.path.splitext(image_basename_withzipext)[0]
-        image_unzipped_path = temp_dir / image_basename
-
-        image_unzipped_path_busy = Path(f"{image_unzipped_path}_busy")
-        # If the input is a zip file, unzip file to temp local location if it doesn't
-        # exist yet. If the file doesn't exist yet in right projection, read original
-        # input file to reproject/write to new file with correct epsg
-        if not (image_unzipped_path_busy.exists() or image_unzipped_path.exists()):
-            # Create temp dir if it doesn't exist yet
-            os.makedirs(temp_dir, exist_ok=True)
-
-            # Create lock file in an atomic way, so we are sure we are the only process
-            # working on it. If function returns true, there isn't any other
-            # thread/process already working on it
-            if io_util.create_file_atomic(image_unzipped_path_busy):
-                try:
-                    logger.info(
-                        f"Unzip image {image_path} to local location {temp_dir}"
-                    )
-
-                    # Create the dest dir where the file will be unzipped to + unzip!
-                    if not image_unzipped_path.exists():
-                        import zipfile
-
-                        with zipfile.ZipFile(image_path, "r") as zippedfile:
-                            zippedfile.extractall(temp_dir)
-                finally:
-                    # Remove lock file when we are ready
-                    os.remove(image_unzipped_path_busy)
-
-        # If a "busy file" still exists, the file isn't ready yet, but another process
-        # is working on it, so wait till it disappears
-        while image_unzipped_path_busy.exists():
-            time.sleep(1)
-
-        # Now we are ready to return the path...
-        return image_unzipped_path
diff --git a/cropclassification/predict/classification.py b/cropclassification/predict/classification.py
index 16c0d32d..4b801b37 100644
--- a/cropclassification/predict/classification.py
+++ b/cropclassification/predict/classification.py
@@ -82,6 +82,10 @@ def train_test_predict(
     if input_parcel_classification_data_df.index.name != conf.columns["id"]:
         input_parcel_classification_data_df.set_index(conf.columns["id"], inplace=True)
     logger.debug("Read classification data file ready")
+    if len(input_parcel_classification_data_df.columns) == 0:
+        raise ValueError(
+            f"no classification data found in {input_parcel_classification_data_path}"
+        )
 
     # Train the classification
     output_classifier_path = train(
diff --git a/cropclassification/preprocess/_timeseries_calc_openeo.py b/cropclassification/preprocess/_timeseries_calc_openeo.py
new file mode 100644
index 00000000..b54d97c0
--- /dev/null
+++ b/cropclassification/preprocess/_timeseries_calc_openeo.py
@@ -0,0 +1,77 @@
+from datetime import datetime
+import logging
+from pathlib import Path
+import tempfile
+from typing import List
+
+from cropclassification.helpers import config_helper as conf
+from cropclassification.util import openeo_util
+from cropclassification.util.openeo_util import ImageProfile
+from cropclassification.util import zonal_stats_bulk
+
+# First define/init some general variables/constants
+# -------------------------------------------------------------
+
+# Get a logger...
+logger = logging.getLogger(__name__)
+
+# The real work
+# -------------------------------------------------------------
+
+
+def calculate_periodic_timeseries(
+    input_parcel_path: Path,
+    start_date: datetime,
+    end_date: datetime,
+    sensordata_to_get: List[ImageProfile],
+    dest_image_data_dir: Path,
+    dest_data_dir: Path,
+):
+    """
+    Calculate timeseries data for the input parcels.
+
+    args
+        data_to_get: an array with data you want to be calculated: check out the
+            constants starting with DATA_TO_GET... for the options.
+    """
+    # As we want a weekly calculation, get nearest monday for start and stop day
+    days_per_period = 7
+    periodic_images_result = openeo_util.calc_periodic_mosaic(
+        roi_path=input_parcel_path,
+        start_date=start_date,
+        end_date=end_date,
+        days_per_period=days_per_period,
+        output_dir=dest_image_data_dir,
+        images_to_get=sensordata_to_get,
+        force=False,
+    )
+
+    # Now calculate the timeseries
+    images_bands = [
+        (image_path, imageprofile.bands)
+        for image_path, imageprofile in periodic_images_result
+    ]
+    temp_dir = conf.dirs.getpath("temp_dir")
+    if temp_dir == "None":
+        temp_dir = Path(tempfile.gettempdir())
+    zonal_stats_bulk.zonal_stats(
+        features_path=input_parcel_path,
+        id_column=conf.columns["id"],
+        images_bands=images_bands,
+        output_dir=dest_data_dir,
+        temp_dir=temp_dir,
+        log_dir=conf.dirs.getpath("log_dir"),
+        log_level=conf.general.get("log_level"),
+    )
+
+    """
+    for image_path in periodic_images:
+        output_path = (
+            dest_data_dir / f"{input_parcel_path.stem}__{image_path.stem}.gpkg"
+        )
+        geoops_util.zonal_stats(
+            input_vector_path=input_parcel_path,
+            input_raster_path=image_path,
+            output_path=output_path,
+        )
+    """
diff --git a/cropclassification/preprocess/_timeseries_helper.py b/cropclassification/preprocess/_timeseries_helper.py
index 46342514..f518fb68 100644
--- a/cropclassification/preprocess/_timeseries_helper.py
+++ b/cropclassification/preprocess/_timeseries_helper.py
@@ -3,12 +3,12 @@
 Calculates periodic timeseries for input parcels.
 """
 
-from datetime import datetime
+from datetime import datetime, timedelta
 import logging
 import gc
 import os
 from pathlib import Path
-from typing import List, Optional
+from typing import List, Optional, Union
 
 import geofileops as gfo
 import numpy as np
@@ -34,7 +34,7 @@ def prepare_input(
     output_imagedata_parcel_input_path: Path,
     output_parcel_nogeo_path: Optional[Path] = None,
     force: bool = False,
-):
+) -> bool:
     """
     This function creates a file that is preprocessed to be a good input file for
     timeseries extraction of sentinel images.
@@ -44,7 +44,6 @@ def prepare_input(
         output_imagedata_parcel_input_path (Path): prepared output file
         output_parcel_nogeo_path (Path): output file with a copy of the non-geo data
         force: force creation, even if output file(s) exist already
-
     """
     # Check if parameters are OK and init some extra params
     if not input_parcel_path.exists():
@@ -66,11 +65,11 @@ def prepare_input(
         and (output_parcel_nogeo_path is None or output_parcel_nogeo_path.exists())
     ):
         logger.warning(
-            "prepare_input: force == False and output files exist, so stop: "
+            "prepare_input: force is False and output files exist, so stop: "
             + f"{output_imagedata_parcel_input_path}, "
             + f"{output_parcel_nogeo_path}"
         )
-        return
+        return False
 
     logger.info(f"Process input file {input_parcel_path}")
 
@@ -107,10 +106,10 @@ def prepare_input(
     # If force == False Check and the output file exists already, stop.
     if force is False and output_imagedata_parcel_input_path.exists():
         logger.warning(
-            "prepare_input: force == False and output files exist, so stop: "
+            "prepare_input: force is False and output files exist, so stop: "
             f"{output_imagedata_parcel_input_path}"
         )
-        return
+        return False
 
     logger.info("Apply buffer on parcel")
     parceldata_buf_gdf = parceldata_gdf.copy()
@@ -166,36 +165,29 @@ def prepare_input(
     gfo.to_file(parceldata_buf_poly_gdf, output_imagedata_parcel_input_path)
     logger.info(parceldata_buf_poly_gdf)
 
-    message = (
-        "The buffered file just has been prepared, so probably you now you probably "
-        "need to sync it to the DIAS and start the timeseries data extraction before "
-        "proceding!"
-    )
-    logger.warning(message)
-    raise Exception(message)
+    return True
 
 
-def calculate_periodic_data(
-    input_parcel_path: Path,
-    input_base_dir: Path,
-    start_date_str: str,
-    end_date_str: str,
+def calculate_periodic_timeseries(
+    parcel_path: Path,
+    timeseries_per_image_dir: Path,
+    start_date: datetime,
+    end_date: datetime,
     sensordata_to_get: List[str],
     dest_data_dir: Path,
     force: bool = False,
 ):
     """
-    This function creates a file that is a weekly summarize of timeseries images from
-    DIAS.
+    This function creates a file that is a weekly aggregation of timeseries files.
 
     TODO: add possibility to choose which values to extract (mean, min, max,...)?
 
     Args:
-        input_parcel_path (Path): [description]
-        input_base_dir (Path): [description]
-        start_date_str (str): Start date in format %Y-%m-%d. Needs to be aligned
+        parcel_path (Path): [description]
+        timeseries_per_image_dir (Path): [description]
+        start_date (datetime): Start date. Needs to be aligned
             already on the periods wanted + data on this date is included.
-        end_date_str (str): End date in format %Y-%m-%d. Needs to be aligned already on
+        end_date (datetime): End date. Needs to be aligned already on
             the periods wanted + data on this date is excluded.
         sensordata_to_get ([]):
         dest_data_dir (Path): [description]
@@ -204,30 +196,22 @@ def calculate_periodic_data(
     logger.info("calculate_periodic_data")
 
     # Init
-    input_dir = input_base_dir / input_parcel_path.stem
-
     # TODO: in config?
     input_ext = ".sqlite"
     output_ext = ".sqlite"
 
-    start_date = datetime.strptime(start_date_str, "%Y-%m-%d")
-    end_date = datetime.strptime(end_date_str, "%Y-%m-%d")
-    year = start_date_str.split("-")[0]
+    year = start_date.year
 
     # Prepare output dir
-    test = False
-    if test:
-        dest_data_dir = Path(f"{str(dest_data_dir)}_test")
-    if not dest_data_dir.exists():
-        os.mkdir(dest_data_dir)
+    dest_data_dir.mkdir(parents=True, exist_ok=True)
 
     # Create Dataframe with all files with their info
     logger.debug("Create Dataframe with all files and their properties")
     file_info_list = []
-    for filename in os.listdir(input_dir):
+    for filename in os.listdir(timeseries_per_image_dir):
         if filename.endswith(input_ext):
             # Get seperate filename parts
-            file_info = get_file_info(input_dir / filename)
+            file_info = get_fileinfo_timeseries(input_dir / filename)
             file_info_list.append(file_info)
 
     all_inputfiles_df = pd.DataFrame(file_info_list)
@@ -239,7 +223,7 @@ def calculate_periodic_data(
             "Get files we need based on start- & stopdates, sensordata_to_get,..."
         )
         orbits = [None]
-        if sensordata_type == conf.general["SENSORDATA_S1_ASCDESC"]:
+        if sensordata_type == "S1AscDesc":
             # Filter files to the ones we need
             # satellitetype = "S1"
             imagetype = IMAGETYPE_S1_GRD
@@ -252,8 +236,19 @@ def calculate_periodic_data(
                 & (all_inputfiles_df.band.isin(bands))
                 & (all_inputfiles_df.orbit.isin(orbits))
             ]
-        elif sensordata_type == conf.general["SENSORDATA_S2gt95"]:
-            # satellitetype = "S2"
+        elif sensordata_type == "S1Coh":
+            satellitetype = "S1"
+            imagetype = IMAGETYPE_S1_COHERENCE
+            bands = ["VV", "VH"]
+            orbits = ["ASC", "DESC"]
+            needed_inputfiles_df = all_inputfiles_df.loc[
+                (all_inputfiles_df.date >= start_date)
+                & (all_inputfiles_df.date < end_date)
+                & (all_inputfiles_df.imagetype == imagetype)
+                & (all_inputfiles_df.band.isin(bands))
+            ]
+        elif sensordata_type == "S2gt95":
+            satellitetype = "S2"
             imagetype = IMAGETYPE_S2_L2A
             bands = ["B02-10m", "B03-10m", "B04-10m", "B08-10m", "B11-20m", "B12-20m"]
             needed_inputfiles_df = all_inputfiles_df.loc[
@@ -262,11 +257,20 @@ def calculate_periodic_data(
                 & (all_inputfiles_df.imagetype == imagetype)
                 & (all_inputfiles_df.band.isin(bands))
             ]
-        elif sensordata_type == conf.general["SENSORDATA_S1_COHERENCE"]:
-            # satellitetype = "S1"
-            imagetype = IMAGETYPE_S1_COHERENCE
-            bands = ["VV", "VH"]
-            orbits = ["ASC", "DESC"]
+        elif sensordata_type == "S2-landcover":
+            satellitetype = "S2"
+            imagetype = IMAGETYPE_S2_L2A
+            bands = ["landcover"]
+            needed_inputfiles_df = all_inputfiles_df.loc[
+                (all_inputfiles_df.date >= start_date)
+                & (all_inputfiles_df.date < end_date)
+                & (all_inputfiles_df.imagetype == imagetype)
+                & (all_inputfiles_df.band.isin(bands))
+            ]
+        elif sensordata_type == "S2-ndvi":
+            satellitetype = "S2"
+            imagetype = IMAGETYPE_S2_L2A
+            bands = ["ndvi"]
             needed_inputfiles_df = all_inputfiles_df.loc[
                 (all_inputfiles_df.date >= start_date)
                 & (all_inputfiles_df.date < end_date)
@@ -274,10 +278,10 @@ def calculate_periodic_data(
                 & (all_inputfiles_df.band.isin(bands))
             ]
         else:
-            raise Exception(f"Unsupported sensordata_type: {sensordata_type}")
+            raise ValueError(f"Unsupported sensordata_type: {sensordata_type}")
 
         # There should also be one pixcount file
-        pixcount_filename = f"{input_parcel_path.stem}_weekly_pixcount{output_ext}"
+        pixcount_filename = f"{parcel_path.stem}_weekly_pixcount{output_ext}"
         pixcount_path = dest_data_dir / pixcount_filename
 
         # For each week
@@ -293,8 +297,8 @@ def calculate_periodic_data(
             # New file name
             period_date_str_long = period_date.strftime("%Y-%m-%d")
             period_data_filename = (
-                f"{input_parcel_path.stem}_weekly_{period_date_str_long}_"
-                f"{sensordata_type}{output_ext}"
+                f"{parcel_path.stem}_weekly_{period_date_str_long}_{sensordata_type}"
+                f"{output_ext}"
             )
             period_data_path = dest_data_dir / period_data_filename
 
@@ -486,7 +490,7 @@ def calculate_periodic_data(
                     pdh.to_file(pixcount_df, pixcount_path)
 
 
-def get_file_info(path: Path) -> dict:
+def get_fileinfo_timeseries(path: Path) -> dict:
     """
     This function gets info of a timeseries data file.
 
@@ -500,45 +504,70 @@ def get_file_info(path: Path) -> dict:
     try:
         # Split name on parcelinfo versus imageinfo
         filename_splitted = path.stem.split("__")
-        filename_imageinfo = filename_splitted[1]
+        parcel_part = filename_splitted[0]
+        imageinfo_part = filename_splitted[1]
 
         # Extract imageinfo
-        imageinfo_values = filename_imageinfo.split("_")
-
-        # Satellite
-        satellite = imageinfo_values[0]
-
-        # Get the date taken from the filename, depending on the satellite type
-        # Remark: the datetime is in this format: '20180101T055812'
-        imagetype = None
-        filedatetime = None
-        if satellite.startswith("S1"):
-            # Check if it is a GRDH image
-            if imageinfo_values[2] == "GRDH":
-                imagetype = IMAGETYPE_S1_GRD
-                filedatetime = imageinfo_values[4]
-            elif imageinfo_values[1].startswith("S1"):
-                imagetype = IMAGETYPE_S1_COHERENCE
-                filedatetime = imageinfo_values[2]
-        elif satellite.startswith("S2"):
-            imagetype = IMAGETYPE_S2_L2A
-            filedatetime = imageinfo_values[2]
+        imageinfo_values = imageinfo_part.split("_")
+        orbit = None
+        image_profile = None
+        if "-" in imageinfo_values[0]:
+            # OpenEO mosaic filename format
+            image_profile = imageinfo_values[0]
+            imageprofile_parts = image_profile.split("-")
+            satellite = imageprofile_parts[0]
+            product = imageprofile_parts[1]
+            if satellite == "s1":
+                if product in ["asc", "desc"]:
+                    imagetype = IMAGETYPE_S1_GRD
+                    orbit = product
+                else:
+                    imagetype = IMAGETYPE_S1_COHERENCE
+            elif satellite == "s2":
+                imagetype = IMAGETYPE_S2_L2A
+            else:
+                raise ValueError(f"invalid imageprofile in {path}")
 
-        if imagetype is None or filedatetime is None:
-            raise Exception(f"Unsupported file: {path}")
+            start_date = datetime.fromisoformat(imageinfo_values[1])
+            end_date = datetime.fromisoformat(imageinfo_values[1])
+            band = imageinfo_values[-1]  # =last value
 
-        filedate = filedatetime.split("T")[0]
-        parseddate = datetime.strptime(filedate, "%Y%m%d")
-        fileweek = int(parseddate.strftime("%W"))
+        else:
+            # ONDA/ESA filename format
+            # Satellite
+            satellite = imageinfo_values[0].lower()
+            # Get the date taken from the filename, depending on the satellite type
+            # Remark: the datetime is in this format: '20180101T055812'
+            if satellite.startswith("s1"):
+                # Check if it is a GRDH image
+                if imageinfo_values[2] == "GRDH":
+                    imagetype = IMAGETYPE_S1_GRD
+                    filedatetime = imageinfo_values[4]
+                elif imageinfo_values[1].startswith("S1"):
+                    imagetype = IMAGETYPE_S1_COHERENCE
+                    filedatetime = imageinfo_values[2]
+                else:
+                    raise ValueError(f"Unsupported file: {path}")
+                # Also get the orbit
+                orbit = imageinfo_values[-2].lower()  # =2nd last value
 
-        # Get the band
-        fileband = imageinfo_values[-1]  # =last value
+            elif satellite.startswith("s2"):
+                imagetype = IMAGETYPE_S2_L2A
+                filedatetime = imageinfo_values[2]
+            else:
+                raise ValueError(f"Unsupported file: {path}")
 
-        # For S1 images, get the orbit
-        if satellite.startswith("S1"):
-            fileorbit = imageinfo_values[-2]  # =2nd last value
-        else:
-            fileorbit = None
+            # Parse the data found
+            filedate = filedatetime.split("T")[0]
+            parseddate = datetime.strptime(filedate, "%Y%m%d")
+            start_date = parseddate
+            end_date = start_date
+
+            # Get the band
+            band = imageinfo_values[-1]  # =last value
+
+        # Week
+        fileweek = int(start_date.strftime("%W"))
 
         # The file paths of these files sometimes are longer than 256
         # characters, so use trick on windows to support this anyway
@@ -549,38 +578,72 @@ def get_file_info(path: Path) -> dict:
             else:
                 path_safe = f"//?/{path_safe}"
 
-        filenameparts = {
+        image_metadata = {
             "path": path_safe,
+            "parcel_stem": parcel_part,
             "imagetype": imagetype,
             "filestem": path.stem,
-            "date": parseddate,
+            "start_date": start_date,
+            "end_date": end_date,
             "week": fileweek,
-            "band": fileband,
-            "orbit": fileorbit,
-        }  # ASC/DESC
+            "band": band,
+            "orbit": orbit,  # ASC/DESC
+        }
+        if image_profile is not None:
+            image_metadata["image_profile"] = image_profile
 
     except Exception as ex:
         message = f"Error extracting info from filename {path}"
         logger.exception(message)
         raise Exception(message) from ex
 
-    return filenameparts
+    return image_metadata
+
+
+def get_fileinfo_timeseries_periods(path: Path) -> dict:
+    """
+    This function gets info of a period timeseries data file.
 
+    Args:
+        path (Path): The path to the file to get info about.
 
-def get_monday(input_date: str) -> datetime:
+    Returns:
+        dict: a dict containing info about the file
+    """
+    # If there is no double underscore in name: "old file type"
+    if "__" not in path.stem:
+        stem_parts = path.stem.split("_")
+        if len(stem_parts) < 4:
+            raise ValueError(f"stem doesn't contain enough parts: {path.name}")
+        period_name = stem_parts[-3]
+        start_date = datetime.fromisoformat(stem_parts[-2])
+        if period_name.lower() == "weekly":
+            end_date = start_date + timedelta(days=6)
+        else:
+            raise ValueError(f"unsupported period_name in {path}")
+        return {
+            "path": path,
+            "parcel_stem": "_".join(stem_parts[0:-3]),
+            "period_name": period_name,
+            "start_date": start_date,
+            "end_date": end_date,
+            "image_profile": stem_parts[-1],
+        }
+
+    return get_fileinfo_timeseries(path)
+
+
+def get_monday(date: Union[str, datetime]) -> datetime:
     """
     This function gets the first monday before the date provided.
-    She is being used to adapt start_date and end_date so they are mondays, so it
+    It is being used to adapt start_date and end_date so they are mondays, so it
     becomes easier to reuse timeseries data
        - inputformat:  %Y-%m-%d
        - outputformat: datetime
     """
-    parseddate = datetime.strptime(input_date, "%Y-%m-%d")
-    year_week = parseddate.strftime("%Y_%W")
+    if isinstance(date, str):
+        date = datetime.strptime(date, "%Y-%m-%d")
+
+    year_week = date.strftime("%Y_%W")
     year_week_monday = datetime.strptime(year_week + "_1", "%Y_%W_%w")
     return year_week_monday
-
-
-# If the script is run directly...
-if __name__ == "__main__":
-    raise Exception("Not implemented")
diff --git a/cropclassification/preprocess/classification_preprocess.py b/cropclassification/preprocess/classification_preprocess.py
index 6ae208aa..084a83be 100644
--- a/cropclassification/preprocess/classification_preprocess.py
+++ b/cropclassification/preprocess/classification_preprocess.py
@@ -29,7 +29,9 @@
 def prepare_input(
     input_parcel_path: Path,
     input_parcel_filetype: str,
-    input_parcel_pixcount_path: Path,
+    timeseries_periodic_dir: Path,
+    base_filename: str,
+    data_ext: str,
     classtype_to_prepare: str,
     classes_refe_path: Path,
     output_parcel_path: Path,
@@ -41,7 +43,7 @@ def prepare_input(
     """
 
     # If force == False Check and the output file exists already, stop.
-    if force is False and os.path.exists(output_parcel_path) is True:
+    if not force and output_parcel_path.exists():
         logger.warning(
             f"prepare_input: output file exists + force is False: {output_parcel_path}"
         )
@@ -69,9 +71,20 @@ def prepare_input(
         raise Exception(message)
 
     # Load pixcount data and join it
-    logger.info(f"Read pixcount file {input_parcel_pixcount_path}")
-    df_pixcount = pdh.read_file(input_parcel_pixcount_path)
-    logger.debug(f"Read pixcount file ready, shape: {df_pixcount.shape}")
+    parcel_pixcount_path = (
+        timeseries_periodic_dir / f"{base_filename}_pixcount{data_ext}"
+    )
+    if parcel_pixcount_path.exists():
+        logger.info(f"Read pixcount file {parcel_pixcount_path}")
+        df_pixcount = pdh.read_file(parcel_pixcount_path)
+        logger.debug(f"Read pixcount file ready, shape: {df_pixcount.shape}")
+    else:
+        parcel_pixcount_path = next(timeseries_periodic_dir.glob("*s2*.sqlite"))
+        df_pixcount = pdh.read_file(parcel_pixcount_path)
+        df_pixcount = df_pixcount[[conf.columns["id"], "count"]].rename(
+            columns={"count": conf.columns["pixcount_s1s2"]}
+        )
+
     if df_pixcount.index.name != conf.columns["id"]:
         df_pixcount.set_index(conf.columns["id"], inplace=True)
 
diff --git a/cropclassification/preprocess/timeseries.py b/cropclassification/preprocess/timeseries.py
index cd5e9597..9fff9421 100644
--- a/cropclassification/preprocess/timeseries.py
+++ b/cropclassification/preprocess/timeseries.py
@@ -3,33 +3,31 @@
 This module contains general functions that apply to timeseries data...
 """
 
+from datetime import datetime
 import logging
 import os
 from pathlib import Path
-from typing import List
+from typing import Dict, List
 
 import cropclassification.helpers.config_helper as conf
 import cropclassification.helpers.pandas_helper as pdh
 import cropclassification.preprocess._timeseries_helper as ts_helper
 
-# -------------------------------------------------------------
 # First define/init some general variables/constants
 # -------------------------------------------------------------
+
 # Get a logger...
 logger = logging.getLogger(__name__)
 
-# -------------------------------------------------------------
 # The real work
 # -------------------------------------------------------------
 
 
 def calc_timeseries_data(
     input_parcel_path: Path,
-    input_country_code: str,
     start_date_str: str,
     end_date_str: str,
-    sensordata_to_get: List[str],
-    base_filename: str,
+    sensordata_to_get: Dict[str, conf.SensorData],
     dest_data_dir: Path,
 ):
     """
@@ -37,19 +35,17 @@ def calc_timeseries_data(
 
     Args:
         input_parcel_path (str): [description]
-        input_country_code (str): [description]
         start_date_str (str): [description]
         end_date_str (str): [description]
         sensordata_to_get (List[str]): an array with data you want to be calculated:
-                check out the constants starting with DATA_TO_GET... for the options.
-        base_filename (str): [description]
+            check out the constants starting with DATA_TO_GET... for the options.
         dest_data_dir (str): [description]
     """
     # Check some variables...
     if sensordata_to_get is None:
         raise Exception("sensordata_to_get cannot be None")
     if not dest_data_dir.exists():
-        os.makedirs(dest_data_dir)
+        dest_data_dir.mkdir(parents=True, exist_ok=True)
 
     # As we want a weekly calculation, get nearest monday for start and stop day
     start_date = ts_helper.get_monday(
@@ -61,27 +57,44 @@ def calc_timeseries_data(
         f"Start date {start_date_str} converted to monday before: {start_date}, end "
         f"date {end_date_str} as well: {end_date}"
     )
-    timeseries_calc_type = conf.timeseries["timeseries_calc_type"]
-    if timeseries_calc_type == "onda":
+    sensordata_to_get_onda = [
+        sensor for sensor in sensordata_to_get if sensor not in conf.image_profiles
+    ]
+    sensordata_to_get_openeo = [
+        sensor for sensor in sensordata_to_get if sensor in conf.image_profiles
+    ]
+
+    if len(sensordata_to_get_onda) > 0:
         # Start!
         # TODO: start calculation of per image data on DIAS
         # import cropclassification.preprocess.timeseries_calc_dias_onda_per_image as
         # ts_calc
-        timeseries_per_image_dir = conf.dirs.getpath("timeseries_per_image_dir")
 
         # Now all image data is available per image, calculate periodic data
-        return ts_helper.calculate_periodic_data(
+        ts_helper.calculate_periodic_timeseries(
+            parcel_path=input_parcel_path,
+            timeseries_per_image_dir=conf.dirs.getpath("timeseries_per_image_dir"),
+            start_date=start_date,
+            end_date=end_date,
+            sensordata_to_get=sensordata_to_get_onda,
+            dest_data_dir=dest_data_dir,
+        )
+    if len(sensordata_to_get_openeo) > 0:
+        # Pepare periodic images + calculate base timeseries on them
+        import cropclassification.preprocess._timeseries_calc_openeo as ts_calc_openeo
+
+        sensordata_to_get_info_openeo = [
+            conf.image_profiles[sensordatatype]
+            for sensordatatype in sensordata_to_get_openeo
+        ]
+        ts_calc_openeo.calculate_periodic_timeseries(
             input_parcel_path=input_parcel_path,
-            input_base_dir=timeseries_per_image_dir,
-            start_date_str=start_date_str,
-            end_date_str=end_date_str,
-            sensordata_to_get=sensordata_to_get,
+            start_date=start_date,
+            end_date=end_date,
+            sensordata_to_get=sensordata_to_get_info_openeo,
+            dest_image_data_dir=conf.dirs.getpath("images_periodic_dir"),
             dest_data_dir=dest_data_dir,
         )
-    else:
-        message = f"Unsupported timeseries calculation type: {timeseries_calc_type}"
-        logger.error(message)
-        raise Exception(message)
 
 
 def collect_and_prepare_timeseries_data(
@@ -91,7 +104,7 @@ def collect_and_prepare_timeseries_data(
     output_path: Path,
     start_date_str: str,
     end_date_str: str,
-    sensordata_to_use: List[str],
+    sensordata_to_use: Dict[str, conf.SensorData],
     parceldata_aggregations_to_use: List[str],
     force: bool = False,
 ):
@@ -99,8 +112,6 @@ def collect_and_prepare_timeseries_data(
     Collect all timeseries data to use for the classification and prepare it by applying
     scaling,... as needed.
     """
-    SENSORDATA_S2 = conf.general["SENSORDATA_S2"]  # Sentinel 2 data
-    SENSORDATA_S1_COHERENCE = conf.general["SENSORDATA_S1_COHERENCE"]
 
     # If force == False Check and the output file exists already, stop.
     if force is False and output_path.exists() is True:
@@ -108,38 +119,39 @@ def collect_and_prepare_timeseries_data(
             f"Output file already exists and force == False, so stop: {output_path}"
         )
         return
+    start_date = datetime.fromisoformat(start_date_str)
+    end_date = datetime.fromisoformat(end_date_str)
 
     # Init the result with the id's of the parcels we want to treat
     result_df = pdh.read_file(input_parcel_path, columns=[conf.columns["id"]])
     if result_df.index.name != conf.columns["id"]:
         result_df.set_index(conf.columns["id"], inplace=True)
     nb_input_parcels = len(result_df.index)
-    logger.info(
-        "Parceldata aggregations that need to be used: "
-        f"{parceldata_aggregations_to_use}"
-    )
+    logger.info(f"Parceldata aggregations to use: {parceldata_aggregations_to_use}")
     logger.setLevel(logging.DEBUG)
 
     # Loop over all input timeseries data to find the data we really need
-    data_ext = conf.general["data_ext"]
-    path_start = timeseries_dir / f"{base_filename}_{start_date_str}{data_ext}"
-    path_end = timeseries_dir / f"{base_filename}_{end_date_str}{data_ext}"
-    logger.debug(f"path_start_date: {path_start}")
-    logger.debug(f"path_end_date: {path_end}")
-
-    ts_data_files = timeseries_dir.glob(f"{base_filename}_*{data_ext}")
-    for curr_path in sorted(ts_data_files):
+    glob_pattern = f"*{conf.general['data_ext']}"
+    ts_data_paths = list(timeseries_dir.glob(glob_pattern))
+    if len(ts_data_paths) == 0:
+        raise ValueError(f"No timeseries data found for pattern {glob_pattern}")
+
+    for curr_path in sorted(ts_data_paths):
         # Only process data that is of the right sensor types
-        sensor_type = curr_path.stem.split("_")[-1]
-        if sensor_type not in sensordata_to_use:
+        fileinfo = ts_helper.get_fileinfo_timeseries_periods(curr_path)
+        image_profile = fileinfo["image_profile"].lower()
+        if image_profile not in sensordata_to_use:
             logger.debug(
-                f"SKIP: file is not in sensor types asked ({sensordata_to_use}): "
-                f"{curr_path}"
+                f"SKIP: file not needed (only {sensordata_to_use}): {curr_path}"
             )
             continue
         # The only data we want to process is the data in the range of dates
-        if (str(curr_path) < str(path_start)) or (str(curr_path) >= str(path_end)):
-            logger.debug(f"SKIP: File is not in date range asked: {curr_path}")
+        if fileinfo["start_date"] < start_date or fileinfo["end_date"] >= end_date:
+            logger.debug(f"SKIP: file doesn't match the period asked: {curr_path}")
+            continue
+        band = fileinfo["band"]
+        if band not in sensordata_to_use[image_profile].bands:
+            logger.debug(f"SKIP: file doesn't match the bands asked: {curr_path}")
             continue
         # An empty file signifies that there wasn't any valable data for that
         # period/sensor/...
@@ -167,6 +179,7 @@ def collect_and_prepare_timeseries_data(
             data_read_df.set_index(conf.columns["id"], inplace=True)
 
         # Loop over columns to check if there are columns that need to be dropped.
+        columns_to_rename = {}
         for column in data_read_df.columns:
             # If it is the id column, continue
             if column == conf.columns["id"]:
@@ -177,7 +190,13 @@ def collect_and_prepare_timeseries_data(
             for parceldata_aggregation in parceldata_aggregations_to_use:
                 if column.endswith("_" + parceldata_aggregation):
                     column_ok = True
-            if column_ok is False:
+                elif column == parceldata_aggregation:
+                    curr_start_date_str = fileinfo["start_date"].strftime("%Y%m%d")
+                    columns_to_rename[
+                        column
+                    ] = f"{image_profile}_{curr_start_date_str}_{band}_{column}"
+                    column_ok = True
+            if not column_ok:
                 # Drop column if it doesn't end with something in
                 # parcel_data_aggregations_to_use
                 logger.debug(
@@ -199,27 +218,33 @@ def collect_and_prepare_timeseries_data(
                 )
                 data_read_df.drop(column, axis=1, inplace=True)
 
+        # If there are columns that need renaming, do so
+        if len(columns_to_rename) > 0:
+            data_read_df = data_read_df.rename(columns=columns_to_rename)
+
         # If S2, rescale data
-        if sensor_type.startswith(SENSORDATA_S2):
+        if image_profile.startswith("s2"):
             for column in data_read_df.columns:
                 logger.info(
-                    "Column contains S2 data, so scale it by dividing by 10.000: "
-                    f"{column}"
+                    f"Column with s2 data: scale it by dividing by 10.000: {column}"
                 )
                 data_read_df[column] = data_read_df[column] / 10000
 
-        # If S1 coherence, rescale data
-        if sensor_type == SENSORDATA_S1_COHERENCE:
+        # If s1 coherence, rescale data
+        if image_profile in ["s1coh", "s1-coh"]:
             for column in data_read_df.columns:
                 logger.info(
-                    "Column contains S1 Coherence data, so scale it by dividing by "
-                    f"300: {column}"
+                    f"Column with s1 coherence: scale it by dividing by 300: {column}"
                 )
                 data_read_df[column] = data_read_df[column] / 300
 
         # Join the data to the result...
         result_df = result_df.join(data_read_df, how="left")
 
+    # No timeseries data was found, so stop
+    if len(result_df.columns) == 0:
+        raise ValueError("data collection resulted in 0 columns")
+
     # Remove rows with many null values from result
     max_number_null = int(0.6 * len(result_df.columns))
     parcel_many_null_df = result_df[result_df.isnull().sum(axis=1) > max_number_null]
diff --git a/cropclassification/util/geoops_util.py b/cropclassification/util/geoops_util.py
index 7b850ba0..c7cef2a5 100644
--- a/cropclassification/util/geoops_util.py
+++ b/cropclassification/util/geoops_util.py
@@ -1,4 +1,5 @@
 from concurrent import futures
+import datetime
 import logging
 import math
 import multiprocessing
@@ -13,6 +14,11 @@
 from geofileops.util import _io_util
 import psutil
 from psutil._common import bytes2human
+from osgeo import gdal
+
+# ... and suppress errors
+gdal.PushErrorHandler("CPLQuietErrorHandler")
+import rasterio
 import rasterstats
 
 logger = logging.getLogger()
@@ -37,6 +43,11 @@ def zonal_stats(
     boundless: bool = True,
     force: bool = False,
 ):
+    # TODO: GDAL logging was disabled to get rid of this warning, fix it.
+    # Warning 1: TIFFReadDirectory:Sum of Photometric type-related color channels and
+    # ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as
+    # ExtraSamples.
+
     # Init
     if output_path.exists():
         if force is True:
@@ -45,26 +56,29 @@ def zonal_stats(
             logger.info(f"Output file exists already and force is False: {output_path}")
             return
 
+    start_time = datetime.datetime.now()
     tmp_dir = _io_util.create_tempdir("zonal_stats")
     output_temp_path = tmp_dir / f"{output_path.name}"
 
     info = gfo.get_layerinfo(input_vector_path)
+
     max_nb_per_batch = 10000
     nb_batches = math.ceil(info.featurecount / max_nb_per_batch)
     nb_per_batch = math.floor(info.featurecount / nb_batches)
+    nb_batches_done = 0
 
     nb_parallel = multiprocessing.cpu_count()
     future_infos = {}
 
     logger.info(
-        f"calculate zonal_stats for {input_vector_path.name}, in {nb_batches} batches, "
-        f"with {nb_parallel} workers"
+        f"calculate zonal_stats in {nb_batches} batches, with {nb_parallel} workers, "
+        f"for {input_vector_path.name} on {input_raster_path.name}"
     )
     try:
         with futures.ProcessPoolExecutor(
             max_workers=nb_parallel, initializer=setprocessnice(15)
         ) as calculate_pool:
-            # calculate per batch  to manage memory usage
+            # calculate per batch to manage memory usage
             batch_start = 0
             for batch_id in range(nb_batches):
                 # Prepare the slice for this batch
@@ -118,7 +132,6 @@ def zonal_stats(
                         output_temp_partial_path.exists()
                         and output_temp_partial_path.stat().st_size > 0
                     ):
-                        logger.info(f"Write result for batch {future_info['batch_id']}")
                         gfo.append_to(
                             src=output_temp_partial_path,
                             dst=output_temp_path,
@@ -130,6 +143,15 @@ def zonal_stats(
                     logger.exception(f"Error executing {batch_id}")
                     raise
 
+                nb_batches_done += 1
+                progress = format_progress(
+                    start_time=start_time,
+                    nb_done=nb_batches_done,
+                    nb_todo=nb_batches,
+                    operation="zonal_stats",
+                    nb_parallel=nb_parallel,
+                )
+                print(progress)
                 logger.info(
                     "Available virtual memory: "
                     f"{bytes2human(psutil.virtual_memory().available)}, "
@@ -213,7 +235,52 @@ def zonal_stats_ext(
         geojson_out=geojson_out,
         boundless=boundless,
     )
-    result_gdf = gpd.GeoDataFrame.from_features(result)
+    result_gdf = gpd.GeoDataFrame.from_features(result, crs=vector_gdf.crs)
     gfo.to_file(result_gdf, output_path)
 
     return True
+
+
+def format_progress(
+    start_time: datetime.datetime,
+    nb_done: int,
+    nb_todo: int,
+    operation: Optional[str] = None,
+    nb_parallel: int = 1,
+) -> Optional[str]:
+    # Init
+    time_passed = (datetime.datetime.now() - start_time).total_seconds()
+    pct_progress = 100.0 - (nb_todo - nb_done) * 100 / nb_todo
+    nb_todo_str = f"{nb_todo:n}"
+    nb_decimal = len(nb_todo_str)
+
+    # If we haven't really started yet, don't report time estimate yet
+    if nb_done == 0:
+        return (
+            f" ?: ?: ? left, {operation} done on {nb_done:{nb_decimal}n} of "
+            f"{nb_todo:{nb_decimal}n} ({pct_progress:3.2f}%)    "
+        )
+    else:
+        pct_progress = 100.0 - (nb_todo - nb_done) * 100 / nb_todo
+        if time_passed > 0:
+            # Else, report progress properly...
+            processed_per_hour = (nb_done / time_passed) * 3600
+            # Correct the nb processed per hour if running parallel
+            if nb_done < nb_parallel:
+                processed_per_hour = round(processed_per_hour * nb_parallel / nb_done)
+            hours_to_go = (int)((nb_todo - nb_done) / processed_per_hour)
+            min_to_go = (int)((((nb_todo - nb_done) / processed_per_hour) % 1) * 60)
+            secs_to_go = (int)(
+                ((((nb_todo - nb_done) / processed_per_hour) % 1) * 3600) % 60
+            )
+            time_left_str = f"{hours_to_go:02d}:{min_to_go:02d}:{secs_to_go:02d}"
+            nb_left_str = f"{nb_done:{nb_decimal}n} of {nb_todo:{nb_decimal}n}"
+            pct_str = f"({pct_progress:3.2f}%)    "
+        elif pct_progress >= 100:
+            time_left_str = "00:00:00"
+            nb_left_str = f"{nb_done:{nb_decimal}n} of {nb_todo:{nb_decimal}n}"
+            pct_str = f"({pct_progress:3.2f}%)    "
+        else:
+            return None
+        message = f"{time_left_str} left, {operation} done on {nb_left_str} {pct_str}"
+        return message
diff --git a/cropclassification/util/io_util.py b/cropclassification/util/io_util.py
index 7eba1a3d..4feef6b6 100644
--- a/cropclassification/util/io_util.py
+++ b/cropclassification/util/io_util.py
@@ -1,5 +1,7 @@
 import os
 from pathlib import Path
+import tempfile
+from typing import Optional
 
 
 def create_file_atomic(path: Path):
@@ -18,3 +20,46 @@ def create_file_atomic(path: Path):
     finally:
         if fd is not None:
             os.close(fd)
+
+
+def create_tempdir(base_dirname: str, parent_dir: Optional[Path] = None) -> Path:
+    """
+    Creates a new tempdir in the default temp location.
+
+    Remark: the temp dir won't be cleaned up automatically!
+
+    Examples:
+
+        - base_dirname="foo" -> /tmp/foo_000001
+        - base_dirname="foo/bar" -> /tmp/foo/bar_000001
+
+    Args:
+        base_dirname (str): The name the tempdir will start with. The name will be
+            suffixed with a number to make the directory name unique. If a "/" is part
+            of the base_dirname a subdirectory will be created: e.g. "foo/bar".
+        parent_dir (Path, optional): The dir to create the tempdir in. If None, the
+            system temp dir is used. Defaults to None.
+
+    Raises:
+        Exception: if it wasn't possible to create the temp dir because there
+            wasn't found a unique directory name.
+
+    Returns:
+        Path: the path to the temp dir created.
+    """
+
+    if parent_dir is None:
+        parent_dir = Path(tempfile.gettempdir())
+
+    for i in range(1, 999999):
+        try:
+            tempdir = parent_dir / f"{base_dirname}_{i:06d}"
+            tempdir.mkdir(parents=True)
+            return tempdir
+        except FileExistsError:
+            continue
+
+    raise Exception(
+        f"Wasn't able to create a temporary dir with basedir: "
+        f"{parent_dir / base_dirname}"
+    )
diff --git a/cropclassification/util/openeo_util.py b/cropclassification/util/openeo_util.py
index 5bd453fb..5daff50d 100644
--- a/cropclassification/util/openeo_util.py
+++ b/cropclassification/util/openeo_util.py
@@ -1,24 +1,44 @@
 from collections import defaultdict
+from copy import deepcopy
 from datetime import datetime, timedelta
+import json
 import logging
+
 from pathlib import Path
 import pprint
 import time
-from typing import List, Optional
+from typing import List, Optional, Tuple
 
 import geofileops as gfo
 import openeo
 import openeo.rest.job
 import pyproj
+
+from osgeo import gdal
+
+# ... and suppress errors
+gdal.PushErrorHandler("CPLQuietErrorHandler")
 import rasterio
 
-# First define/init some general variables/constants
-# -------------------------------------------------------------
+
 # Get a logger...
 logger = logging.getLogger(__name__)
 
-# The real work
-# -------------------------------------------------------------
+
+class ImageProfile:
+    def __init__(
+        self,
+        name: str,
+        satellite: str,
+        collection: str,
+        bands: List[str],
+        process_options: dict,
+    ):
+        self.name = name
+        self.satellite = satellite
+        self.collection = collection
+        self.bands = bands
+        self.process_options = process_options
 
 
 def calc_periodic_mosaic(
@@ -26,13 +46,21 @@ def calc_periodic_mosaic(
     start_date: datetime,
     end_date: datetime,
     days_per_period: int,
-    sensordata_to_get: List[str],
+    images_to_get: List[ImageProfile],
     output_dir: Path,
     period_name: Optional[str] = None,
+    delete_existing_openeo_jobs: bool = True,
     raise_errors: bool = True,
     force: bool = False,
-) -> List[Path]:
-    result_paths = []
+) -> List[Tuple[Path, ImageProfile]]:
+    # Validate time_dimension_reducer
+    for imageprofile in images_to_get:
+        reducer = imageprofile.process_options.get("time_dimension_reducer", None)
+        if reducer is None:
+            raise ValueError(
+                "process_options must contain a time_dimension_reducer for "
+                "calc_periodic_mosaic"
+            )
 
     # Prepare period_name:
     if period_name is None:
@@ -65,8 +93,15 @@ def calc_periodic_mosaic(
         "north": roi_bounds[3],
     }
 
-    # Connect with VITO openeo backend
-    conn = openeo.connect("https://openeo.vito.be").authenticate_oidc(provider_id="egi")
+    # Connect with openeo backend as configured in file specified in the
+    # OPENEO_CLIENT_CONFIG environment variable.
+    #
+    # More info about the openeo authentication options can be found here:
+    # https://open-eo.github.io/openeo-python-client/auth.html#default-url-and-auto-auth
+    #
+    # More info on the configuration file format can be found here:
+    # https://open-eo.github.io/openeo-python-client/configuration.html#configuration-files  # noqa: E501
+    conn = openeo.connect()
 
     # Save periodic aggregate of asked images
     logger.info(
@@ -74,168 +109,202 @@ def calc_periodic_mosaic(
         f"to {output_dir}"
     )
 
-    # Delete all jobs that are running already
-    jobs = conn.list_jobs()
-    if len(jobs) > 0:
-        for job in jobs:
-            batch_job = openeo.rest.job.BatchJob(job["id"], conn)
-            logger.info(f"delete job '{job['title']}', id: {job['id']}")
-            batch_job.delete_job()
+    if delete_existing_openeo_jobs:
+        # Delete all jobs that are running already
+        jobs = conn.list_jobs()
+        if len(jobs) > 0:
+            for job in jobs:
+                batch_job = openeo.rest.job.BatchJob(job["id"], conn)
+                logger.info(f"delete job '{job['title']}', id: {job['id']}")
+                batch_job.delete_job()
+    else:
+        # Process jobs that are still on the server
+        image_paths, job_errors = get_job_results(conn, output_dir)
+        for image_path in image_paths:
+            add_overviews(image_path)
+        if raise_errors and len(job_errors) > 0:
+            raise RuntimeError(f"Errors occured: {pprint.pformat(job_errors)}")
 
     period_start_date = start_date
+    job_options = {
+        "executor-memory": "4G",
+        "executor-memoryOverhead": "2G",
+        "executor-cores": "2",
+    }
+    result = []
     while period_start_date <= (end_date - timedelta(days=days_per_period)):
-        for sensordata_type in sensordata_to_get:
-            # Period in openeo is inclusive for startdate and excludes enddate
-            period_end_date = period_start_date + timedelta(days=days_per_period)
-            period = [
-                period_start_date.strftime("%Y-%m-%d"),
-                period_end_date.strftime("%Y-%m-%d"),
-            ]
-            satellite, band = sensordata_type.split("-")
-
-            output_basestem = (
-                f"{satellite}_mosaic_{period_name}_{period[0]}_{period[1]}"
+        # Period in openeo is inclusive for startdate and excludes enddate
+        period_end_date = period_start_date + timedelta(days=days_per_period)
+        for imageprofile in images_to_get:
+            # Prepare path. Remark: in the image path the end_date should be inclusive.
+            end_date_incl = period_end_date
+            if days_per_period > 1:
+                end_date_incl = period_end_date - timedelta(days=1)
+            image_path = prepare_image_path(
+                imageprofile.name,
+                start_date=period_start_date,
+                end_date=end_date_incl,
+                bands=imageprofile.bands,
+                dir=output_dir,
             )
+            if not file_exists(image_path, force):
+                process_options = deepcopy(imageprofile.process_options)
+                reducer = process_options["time_dimension_reducer"]
+                del process_options["time_dimension_reducer"]
+
+                create_mosaic_job(
+                    conn=conn,
+                    collection=imageprofile.collection,
+                    spatial_extent=roi_extent,
+                    start_date=period_start_date,
+                    end_date=period_end_date,
+                    bands=imageprofile.bands,
+                    time_dimension_reducer=reducer,
+                    output_name=image_path.name,
+                    job_options=job_options,
+                    process_options=process_options,
+                ).start_job()
+            result.append((image_path, imageprofile))
+
+        period_start_date = period_end_date
 
-            logger.info(f"Get images for period {period}")
+    # Get the results
+    image_paths, job_errors = get_job_results(conn, output_dir)
 
-            if sensordata_type == "S2-landcover":
-                save_rgb = False
-                save_scl = False
-                save_landcover = True
+    # TODO: ideally, all mosaic_paths would be checked if they have overviews.
+    for image_path in image_paths:
+        add_overviews(image_path)
+    if raise_errors and len(job_errors) > 0:
+        raise RuntimeError(f"Errors occured: {pprint.pformat(job_errors)}")
 
-                # Load cube of relevant S2 images
-                # You can look which layers are available here:
-                # https://openeo.cloud/data-collections/
-                s2_cube = conn.load_collection(
-                    "TERRASCOPE_S2_TOC_V2",
-                    spatial_extent=roi_extent,
-                    temporal_extent=period,  # type: ignore
-                    bands=["B04", "B03", "B02", "SCL"],
-                )
-
-                # Save RGB image
-                if save_rgb:
-                    job_title = f"{output_basestem}_rgb.tif"
-                    output_path = output_dir / job_title
-                    if check_output_file(output_path, force):
-                        # Use mask_scl_dilation for "aggressive" cloud mask based on SCL
-                        # Use the "min" reducer filters out "lightly clouded areas"
-                        _ = (
-                            s2_cube.process(
-                                "mask_scl_dilation",
-                                data=s2_cube,
-                                scl_band_name="SCL",
-                            )
-                            .filter_bands(bands=["B04", "B03", "B02"])
-                            .reduce_dimension(dimension="t", reducer="min")
-                            .create_job(out_format="GTiff", title=job_title)
-                            .start_job()
-                        )
-
-                # Save scl layer if available
-                # SCL mask values:
-                #   * 0: NO_DATA
-                #   * 1: SATURATED_OR_DEFECTIVE
-                #   * 2: DARK_AREA_PIXELS
-                #   * 3: CLOUD_SHADOWS
-                #   * 4: VEGETATION
-                #   * 5: NOT_VEGETATED
-                #   * 6: WATER
-                #   * 7: UNCLASSIFIED
-                #   * 8: CLOUD_MEDIUM_PROBABILITY
-                #   * 9: CLOUD_HIGH_PROBABILITY
-                #   * 10: THIN_CIRRUS
-                #   * 11: SNOW
-                if save_scl:
-                    job_title = f"{output_basestem}_scl.tif"
-                    output_path = output_dir / job_title
-                    if check_output_file(output_path, force):
-                        _ = (
-                            s2_cube.band("SCL")
-                            .reduce_dimension(dimension="t", reducer="max")
-                            .create_job(out_format="GTiff", title=job_title)
-                            .start_job()
-                        )
-
-                if save_landcover:
-                    job_title = f"{output_basestem}_{band}.tif"
-                    if (
-                        job_title
-                        == "S2_mosaic_weekly_2021-11-29_2021-12-06_landcover.tif"
-                    ):
-                        period_start_date = period_end_date
-                        continue
-
-                    output_path = output_dir / job_title
-                    result_paths.append(output_path)
-                    if check_output_file(output_path, force):
-                        # Use mask_scl_dilation for "aggressive" cloud mask based on SCL
-                        s2_cube_scl = s2_cube.filter_bands(bands=["SCL"])
-                        s2_cube_scl_data = s2_cube_scl.rename_labels(
-                            dimension="bands", target=["SCL_DATA"], source=["SCL"]
-                        )
-                        s2_cube_scl_merged = s2_cube_scl_data.merge_cubes(s2_cube_scl)
-                        _ = (
-                            s2_cube_scl_merged.process(
-                                "mask_scl_dilation",
-                                data=s2_cube_scl_merged,
-                                scl_band_name="SCL",
-                            )
-                            .filter_bands(bands=["SCL_DATA"])
-                            .reduce_dimension(dimension="t", reducer="max")
-                            .create_job(out_format="GTiff", title=job_title)
-                            .start_job()
-                        )
-
-            elif sensordata_type == "S2-ndvi":
-                job_title = f"{output_basestem}_{band}.tif"
-                output_path = output_dir / job_title
-                result_paths.append(output_path)
-                if check_output_file(output_path, force):
-                    # Load cube of relevant S2 images
-                    # You can look which layers are available here:
-                    #    https://openeo.cloud/data-collections/
-                    s2_ndvi_cube = conn.load_collection(
-                        "TERRASCOPE_S2_NDVI_V2",
-                        spatial_extent=roi_extent,
-                        temporal_extent=period,  # type: ignore
-                        bands=["NDVI_10M", "SCENECLASSIFICATION_20M"],
-                    )
-
-                    # Create + start job
-                    # Use "aggressive" algorythm to mask more clouds based on SCL
-                    _ = (
-                        s2_ndvi_cube.process(
-                            "mask_scl_dilation",
-                            data=s2_ndvi_cube,
-                            scl_band_name="SCENECLASSIFICATION_20M",
-                        )
-                        .reduce_dimension(dimension="t", reducer="max")
-                        .create_job(out_format="GTiff", title=job_title)
-                        .start_job()
-                    )
-            else:
-                raise ValueError(
-                    f"invalid value for sensordata_to_get: {sensordata_type}"
-                )
+    return result
 
-            period_start_date = period_end_date
 
-    # Get the results
-    output_paths = get_job_results(conn, output_dir, raise_errors=raise_errors)
-    for output_path in output_paths:
-        add_overviews(output_path)
+def prepare_image_path(
+    imageprofile: str,
+    start_date: datetime,
+    end_date: datetime,
+    bands: List[str],
+    dir: Path,
+) -> Path:
+    """
+    Returns an image_path + saves a metadata file for the image as f"{image_path}.json".
 
-    return result_paths
+    Args:
+        imageprofile (str): _description_
+        start_date (datetime): _description_
+        end_date (datetime):
+        bands (List[str]): _description_
+        dir (Path): _description_
 
+    Raises:
+        ValueError: _description_
 
-def check_output_file(path: Path, force: bool) -> bool:
+    Returns:
+        Path: _description_
+    """
+    start_date_str = start_date.strftime("%Y-%m-%d")
+    end_date_str = end_date.strftime("%Y-%m-%d")
+    weeks = list(
+        range(int(start_date.strftime("%W")), int(start_date.strftime("%W")) + 1)
+    )
+    bands_str = "-".join(bands)
+    name = f"{imageprofile}_{start_date_str}_{end_date_str}_{bands_str}.tif"
+
+    # If the image metadata file doesn't exist, create it
+    image_path = dir / name
+    imagemeta_path = dir / f"{name}.json"
+    if not imagemeta_path.exists():
+        imageprofile_parts = imageprofile.split("-")
+        satellite = imageprofile_parts[0].lower()
+        metadata = {
+            "imageprofile": imageprofile,
+            "satellite": satellite,
+            "start_date": start_date_str,
+            "end_date": end_date_str,
+            "weeks": weeks,
+            "bands": bands,
+            "path": image_path.as_posix(),
+        }
+        if imageprofile.lower() == "s1-grd-sigma0-asc":
+            metadata["orbit"] = "asc"
+        elif imageprofile.lower() == "s1-grd-sigma0-desc":
+            metadata["orbit"] = "desc"
+
+        # Write to file
+        with open(imagemeta_path, "w") as outfile:
+            outfile.write(json.dumps(metadata, indent=4))
+
+    return image_path
+
+
+def create_mosaic_job(
+    conn: openeo.Connection,
+    collection: str,
+    spatial_extent,
+    start_date: datetime,
+    end_date: datetime,
+    bands: List[str],
+    output_name: str,
+    time_dimension_reducer: str,
+    job_options: dict,
+    process_options: dict,
+) -> openeo.BatchJob:
+    logger.info(f"schedule {output_name}")
+
+    # Check and process input params
+    period = [
+        start_date.strftime("%Y-%m-%d"),
+        end_date.strftime("%Y-%m-%d"),
+    ]
+
+    bands_to_load = list(bands)
+    cloud_filter_band = None
+    for option in process_options:
+        if option == "cloud_filter_band":
+            cloud_filter_band = process_options[option]
+            bands_to_load.append(cloud_filter_band)
+        else:
+            raise ValueError(f"unknown processing option: {option}")
+
+    # Load cube of relevant images
+    # You can look which layers are available here:
+    # https://openeo.cloud/data-collections/
+    cube = conn.load_collection(
+        collection,
+        spatial_extent=spatial_extent,
+        temporal_extent=period,  # type: ignore
+        bands=bands_to_load,
+    )
+
+    # Use mask_scl_dilation for "aggressive" cloud mask based on SCL
+    if cloud_filter_band is not None:
+        cube = cube.process(
+            "mask_scl_dilation",
+            data=cube,
+            scl_band_name=cloud_filter_band,
+        )
+
+    cube = cube.filter_bands(bands=bands)
+    cube = cube.reduce_dimension(dimension="t", reducer=time_dimension_reducer)
+
+    # The NDVI collection needs to be recalculated
+    if collection == "TERRASCOPE_S2_NDVI_V2":
+        cube = cube.apply(lambda x: 0.004 * x - 0.08)
+
+    return cube.create_job(
+        out_format="GTiff",
+        title=output_name,
+        job_options=job_options,
+    )
+
+
+def file_exists(path: Path, force: bool) -> bool:
     """
     Check if the output file exists already. If force is True, the file is removed.
 
     Args:
-        path (Path): the file to test.
+        path (Path): the file to check for.
         force (bool): If True, remove the file if it exists.
 
     Returns:
@@ -246,20 +315,24 @@ def check_output_file(path: Path, force: bool) -> bool:
             path.unlink()
         else:
             logger.info(f"{path.name} exists, so skip")
-            return False
+            return True
 
-    logger.info(f"{path.name} needs to be processed")
-    return True
+    logger.info(f"{path.name} doesn't exist")
+    return False
 
 
 def get_job_results(
-    conn: openeo.Connection, output_dir: Path, raise_errors: bool
-) -> List[Path]:
+    conn: openeo.Connection, output_dir: Path, ignore_errors: bool = False
+) -> Tuple[List[Path], List[str]]:
     """Get results of the completed jobs."""
 
     output_paths = []
+    errors = []
     while True:
         jobs = conn.list_jobs()
+        if len(jobs) == 0:
+            break
+
         jobs_per_status = defaultdict(list)
         for job in jobs:
             jobs_per_status[job["status"]].append(job)
@@ -273,39 +346,39 @@ def get_job_results(
 
             # If there is no title specified, delete job
             if "title" not in job:
-                logger.info(f"job {job} didn't have a title, so just delete it")
+                logger.info(f"job {job} doesn't have a title, so just delete it")
                 batch_job.delete_job()
             else:
                 # Download results + delete job
                 output_path = output_dir / job["title"]
+                logger.info(f"job {job} finished, so download results")
                 batch_job.get_results().download_file(target=output_path)
                 batch_job.delete_job()
-                output_paths += output_paths
+                output_paths.append(output_path)
 
         # As long as processing is needed, keep polling
         if "queued" in jobs_per_status or "running" in jobs_per_status:
             time.sleep(30)
         else:
             # We are ready, so deal with error jobs and break
-            errors_found = []
             for job in jobs_per_status["error"]:
                 batch_job = openeo.rest.job.BatchJob(job["id"], conn)
                 logger.error(f"Error processing job '{job['title']}', id: {job['id']}")
                 errorlog = pprint.pformat(batch_job.logs())
                 logger.error(errorlog)
                 batch_job.delete_job()
-                errors_found.append(f"Error for job '{job['title']}', id: {job['id']}")
+                errors.append(f"Error for job '{job['title']}', id: {job['id']}")
             for job in jobs_per_status["created"]:
                 batch_job = openeo.rest.job.BatchJob(job["id"], conn)
                 logger.info(f"delete created job '{job['title']}', id: {job['id']}")
                 batch_job.delete_job()
 
-            if errors_found and raise_errors:
-                raise RuntimeError(f"Errors occured: {pprint.pformat(errors_found)}")
-
             break
 
-    return output_paths
+    if not ignore_errors and len(errors) > 0:
+        raise RuntimeError(f"openeo processing errors: {pprint.pformat(errors)}")
+
+    return output_paths, errors
 
 
 def add_overviews(path: Path):
@@ -317,7 +390,6 @@ def add_overviews(path: Path):
             if dst.width / factor < 256 or dst.height / factor < 256:
                 break
             factors.append(factor)
-
         if len(factors) > 0:
             dst.build_overviews(factors, rasterio.enums.Resampling.average)
             dst.update_tags(ns="rio_overview", resampling="average")
diff --git a/cropclassification/util/zonal_stats_bulk/__init__.py b/cropclassification/util/zonal_stats_bulk/__init__.py
new file mode 100644
index 00000000..ae07cd86
--- /dev/null
+++ b/cropclassification/util/zonal_stats_bulk/__init__.py
@@ -0,0 +1,2 @@
+from ._calc_zonal_stats_bulk import *
+from ._image_helper import *
diff --git a/cropclassification/preprocess/timeseries_calc_dias_onda_per_image.py b/cropclassification/util/zonal_stats_bulk/_calc_zonal_stats_bulk.py
similarity index 53%
rename from cropclassification/preprocess/timeseries_calc_dias_onda_per_image.py
rename to cropclassification/util/zonal_stats_bulk/_calc_zonal_stats_bulk.py
index c79cfb8a..a6c58fcb 100644
--- a/cropclassification/preprocess/timeseries_calc_dias_onda_per_image.py
+++ b/cropclassification/util/zonal_stats_bulk/_calc_zonal_stats_bulk.py
@@ -14,37 +14,68 @@
 import signal  # To catch CTRL-C explicitly and kill children
 import sys
 import time
-from typing import List, Optional, Union
+from typing import List, Optional, Tuple, Union
+
+from osgeo import gdal
+
+# Suppress gdal warnings/errors
+gdal.PushErrorHandler("CPLQuietErrorHandler")
 
 from affine import Affine
 import geopandas as gpd
 import geofileops as gfo
 import pandas as pd
 import psutil  # To catch CTRL-C explicitly and kill children
-
-from rasterstats import zonal_stats
+import rasterstats
 
 from cropclassification.helpers import pandas_helper as pdh
-from cropclassification.helpers import raster_helper
+from . import _image_helper
 from cropclassification.util import io_util
 
 # General init
 logger = logging.getLogger(__name__)
 
 
-def calc_stats_per_image(
+# TODO: performantie van rasterstats vergelijken met wat andere bibs:
+#   - rastrio manual: https://pysal.org/scipy2019-intermediate-gds/deterministic/gds2-rasters.html  # noqa: E501
+#   - pygeoprocessing (open?): https://pygeoprocessing.readthedocs.io/en/latest/api/pygeoprocessing.html#pygeoprocessing.zonal_statistics  # noqa: E501
+#   - pyjeo (GPL): https://pyjeo.readthedocs.io/en/latest/2_tutorial.html#tutorial-on-extract-calculating-regional-statistics  # noqa: E501
+#   - pyQGIS (GPL): https://gis.stackexchange.com/questions/421556/performance-problem-with-getting-average-pixel-values-within-buffered-circles  # noqa: E501
+#       -> is snel, maar lastig om te installeren op een server.
+#       -> src: https://github.com/qgis/QGIS/blob/d5626d92360efffb4b8085389c8d64072ef65833/src/analysis/vector/qgszonalstatistics.cpp#L266  # noqa: E501
+#   - https://pygis.io/docs/f_rs_extraction.html
+#
+# Dit zijn bibs die ook gedeeltelijke pixels meenemen op basis van opp overlapping,
+# wat niet ideaal is hiervoor...
+#   - xagg: https://github.com/ks905383/xagg
+#   - exactextract (apache2): https://github.com/isciences/exactextract
+def zonal_stats(
     features_path: Path,
     id_column: str,
-    image_paths: List[Path],
-    bands: List[str],
+    images_bands: List[Tuple[Path, List[str]]],
     output_dir: Path,
     temp_dir: Path,
     log_dir: Path,
     log_level: Union[str, int],
+    cloud_filter_band: Optional[str] = None,
+    calc_bands_parallel: bool = True,
     force: bool = False,
 ):
     """
     Calculate the statistics.
+
+    Args:
+        features_path (Path): _description_
+        id_column (str): _description_
+        images_bands (List[Tuple[Path, List[str]]]): _description_
+        output_dir (Path): _description_
+        temp_dir (Path): _description_
+        log_dir (Path): _description_
+        log_level (Union[str, int]): _description_
+        force (bool, optional): _description_. Defaults to False.
+
+    Raises:
+        Exception: _description_
     """
 
     # TODO: probably need to apply some object oriented approach here for "image",
@@ -54,393 +85,381 @@ def calc_stats_per_image(
     # after calculation together with the processing logic
 
     # Some checks on the input parameters
-    if len(image_paths) == 0:
+    nb_todo = len(images_bands)
+    if nb_todo == 0:
         logger.info("No image paths... so nothing to do, so return")
         return
 
     # General init
     output_dir.mkdir(parents=True, exist_ok=True)
-    temp_dir.mkdir(parents=True, exist_ok=True)
     log_dir.mkdir(parents=True, exist_ok=True)
+    temp_dir = io_util.create_tempdir("calc_stats_per_image", parent_dir=temp_dir)
 
     # Create process pool for parallelisation...
     nb_parallel_max = multiprocessing.cpu_count()
     nb_parallel = nb_parallel_max
-    with futures.ProcessPoolExecutor(nb_parallel) as pool:
-        # Loop over all images to start the data preparation for each of them in
-        # parallel...
-        image_paths.sort()
-        start_time = datetime.now()
-        nb_todo = len(image_paths)
-        nb_errors_max = 10
-        nb_errors = 0
-
-        image_dict = {}
-        calc_stats_batch_dict = {}
-        nb_done_total = 0
-        i = 0
 
-        try:
-            # Keep looping. At the end of the loop there are checks when to
-            # break out of the loop...
-            while True:
-                # Start preparation for calculation on next image + features combo
-
-                # If not all images aren't prepared for processing yet, and there
-                # aren't too many busy yet, prepare the next one
-                if (
-                    i < len(image_paths)
-                    and len(filter_on_status(image_dict, "IMAGE_PREPARE_CALC_BUSY"))
-                    < nb_parallel_max
-                ):
-                    # Not too many busy preparing, so get next image_path to start
-                    # prepare on
-                    image_path_str = image_paths[i]
-                    i += 1
-
-                    image_info = raster_helper.get_image_info(Path(image_path_str))
-                    # Create base output filename
-                    orbit = None
-                    if image_info["satellite"].startswith("S1"):
-                        orbit = image_info["orbit_properties_pass"]
-                    output_base_path = format_output_path(
+    # Loop over all images to start the data preparation for each of them in
+    # parallel...
+    start_time = datetime.now()
+    nb_errors_max = 10
+    nb_errors = 0
+
+    image_dict = {}
+    calc_stats_batch_dict = {}
+    nb_done_total = 0
+    image_idx = 0
+
+    pool = futures.ProcessPoolExecutor(nb_parallel)
+    try:
+        # Keep looping. At the end of the loop there are checks when to
+        # break out of the loop...
+        while True:
+            # Start preparation for calculation on next image + features combo
+            # If not all images aren't prepared for processing yet, and there
+            # aren't too many busy yet, prepare the next one
+            if (
+                image_idx < len(images_bands)
+                and len(_filter_on_status(image_dict, "IMAGE_PREPARE_CALC_BUSY"))
+                < nb_parallel_max
+            ):
+                # Not too many busy preparing, so get next image_path to start
+                # prepare on
+                image_path, bands = images_bands[image_idx]
+                image_path_str = str(image_path)
+                image_idx += 1
+
+                image_info = _image_helper.get_image_info(image_path)
+                # Create base output filename
+                # TODO: hoort hier niet echt thuis
+                orbit = None
+                if image_info.imagetype.lower() == "s1-grd-sigma0-asc":
+                    orbit = "ASCENDING"
+                elif image_info.imagetype.lower() == "s1-grd-sigma0-desc":
+                    orbit = "DESCENDING"
+                output_base_path = _format_output_path(
+                    features_path,
+                    image_path,
+                    output_dir,
+                    orbit,
+                    band=None,
+                )
+                output_base_busy_path = temp_dir / f"BUSY_{output_base_path.name}"
+
+                # Check for which bands there is a valid output file already
+                if bands is None:
+                    bands = image_info.bands
+                    assert bands is not None
+                bands_todo = {}
+                for band in bands:
+                    # Prepare the output paths...
+                    output_band_path = _format_output_path(
                         features_path,
-                        Path(image_path_str),
+                        image_path,
                         output_dir,
                         orbit,
-                        band=None,
+                        band,
                     )
-                    output_base_dir, output_base_filename = os.path.split(
-                        output_base_path
-                    )
-                    output_base_busy_path = (
-                        f"{output_base_dir}{os.sep}BUSY_{output_base_filename}"
+                    output_band_busy_path = temp_dir / f"BUSY_{output_band_path.name}"
+
+                    # If a busy output file exists, remove it, otherwise we can get
+                    # double data in it...
+                    if output_band_busy_path.exists():
+                        output_band_busy_path.unlink()
+
+                    # Check if the output file exists already
+                    if output_band_path.exists():
+                        if not force:
+                            logger.debug(
+                                f"Output file for band exists {output_band_path}"
+                            )
+                            continue
+                        else:
+                            output_band_path.unlink()
+                    bands_todo[band] = output_band_path
+
+                # If all bands already processed, skip image...
+                if len(bands_todo) == 0:
+                    logger.info(
+                        "SKIP image: output files for all bands exist already for "
+                        f"{output_base_path}"
                     )
+                    nb_todo -= 1
+                    continue
 
-                    # Check for which bands there is a valid output file already
-                    bands_done = 0
-                    for band in bands:
-                        # Prepare the output paths...
-                        output_band_path = format_output_path(
-                            features_path,
-                            Path(image_path_str),
-                            output_dir,
-                            orbit,
-                            band,
-                        )
-                        output_band_dir, output_band_filename = os.path.split(
-                            output_band_path
-                        )
-                        output_band_busy_path = (
-                            f"{output_band_dir}{os.sep}BUSY_{output_band_filename}"
-                        )
-
-                        # If a busy output file exists, remove it, otherwise we can get
-                        # double data in it...
-                        if os.path.exists(output_band_busy_path):
-                            os.remove(output_band_busy_path)
-
-                        # Check if the output file exists already
-                        if os.path.exists(output_band_path):
-                            if not force:
-                                logger.debug(
-                                    f"Output file for band exists {output_band_path}"
-                                )
-                                bands_done += 1
-                            else:
-                                os.remove(output_band_path)
-
-                    # If all bands already processed, skip image...
-                    if len(bands) == bands_done:
-                        logger.info(
-                            "SKIP image: output files for all bands exist already for "
-                            f"{output_base_path}"
-                        )
-                        nb_todo -= 1
-                        continue
-
-                    # Start the prepare processing assync
-                    # TODO: possibly it is cleaner to do this per band...
-                    future = pool.submit(
-                        prepare_calc,
-                        features_path,
-                        id_column,
-                        Path(image_path_str),
-                        output_base_path,
-                        temp_dir,
-                        log_dir,
-                        log_level,
-                        nb_parallel_max,
-                    )
-                    image_dict[image_path_str] = {
-                        "features_path": features_path,
-                        "prepare_calc_future": future,
-                        "image_info": image_info,
-                        "prepare_calc_starttime": datetime.now(),
-                        "output_base_path": output_base_path,
-                        "output_base_busy_path": output_base_busy_path,
-                        "status": "IMAGE_PREPARE_CALC_BUSY",
-                    }
-
-                    # Jump to next image to start the prepare_calc for it...
+                # Start the prepare processing assync
+                # TODO: possibly it is cleaner to do this per band...
+                logger.info(f"Start calculation for image {output_base_path}")
+                future = pool.submit(
+                    _prepare_calc,
+                    features_path,
+                    id_column,
+                    image_path,
+                    temp_dir,
+                    log_dir,
+                    log_level,
+                    nb_parallel_max,
+                )
+                image_dict[image_path_str] = {
+                    "features_path": features_path,
+                    "image_path": image_path,
+                    "bands": bands,
+                    "bands_todo": bands_todo,
+                    "prepare_calc_future": future,
+                    "image_info": image_info,
+                    "prepare_calc_starttime": datetime.now(),
+                    "output_base_path": output_base_path,
+                    "output_base_busy_path": output_base_busy_path,
+                    "status": "IMAGE_PREPARE_CALC_BUSY",
+                }
+
+                # Jump to next image to start the prepare_calc for it...
+                continue
+
+            # For images busy preparing: if ready: start real calculation in batches
+
+            # Loop through the busy preparations
+            for image_path_str in _filter_on_status(
+                image_dict, "IMAGE_PREPARE_CALC_BUSY"
+            ):
+                # If still running, go to next
+                image = image_dict[image_path_str]
+                if image["prepare_calc_future"].running():
                     continue
 
-                # For images busy preparing: if ready: start real calculation in batches
+                # Extract the result from the preparation
+                try:
+                    # Get the result from the completed prepare_inputs
+                    prepare_calc_result = image["prepare_calc_future"].result()
+
+                    # If nb_features to be treated is 0... create (empty) output
+                    # files and continue with next...
+                    if prepare_calc_result["nb_features_to_calc_total"] == 0:
+                        for band in image["bands_todo"]:
+                            # Create output file
+                            output_band_path = image["bands_todo"][band]
+                            logger.info(
+                                "No features found in image: create done file: "
+                                f"{output_band_path}"
+                            )
+                            io_util.create_file_atomic(output_band_path)
+                            image["status"] = "IMAGE_CALC_DONE"
 
-                # Loop through the busy preparations
-                for image_path_str in filter_on_status(
-                    image_dict, "IMAGE_PREPARE_CALC_BUSY"
-                ):
-                    # If still running, go to next
-                    image = image_dict[image_path_str]
-                    if image["prepare_calc_future"].running():
+                        # Jump to next image
                         continue
 
-                    # Extract the result from the preparation
-                    try:
-                        # Get the result from the completed  prepare_inputs
-                        prepare_calc_result = image["prepare_calc_future"].result()
-
-                        # If nb_features to be treated is 0... create (empty) output
-                        # files and continue with next...
-                        if prepare_calc_result["nb_features_to_calc_total"] == 0:
-                            for band in bands:
-                                # Prepare the output path...
-                                orbit = None
-                                if image["image_info"]["satellite"].startswith("S1"):
-                                    orbit = image["image_info"]["orbit_properties_pass"]
-                                output_band_path = format_output_path(
-                                    features_path,
-                                    Path(image_path_str),
-                                    output_dir,
-                                    orbit,
-                                    band,
-                                )
-
-                                # Create output file
-                                logger.info(
-                                    "No features found in image: create done file: "
-                                    f"{output_band_path}"
-                                )
-                                io_util.create_file_atomic(output_band_path)
-                                image["status"] = "IMAGE_CALC_DONE"
-
-                            # Jump to next image
-                            continue
+                    # Add info about the result of the prepare_calc to image info
+                    image["image_prepared_path"] = prepare_calc_result[
+                        "image_prepared_path"
+                    ]
+                    image["feature_batches"] = prepare_calc_result["feature_batches"]
+                    image["nb_features_to_calc_total"] = prepare_calc_result[
+                        "nb_features_to_calc_total"
+                    ]
+                    image["temp_features_dir"] = prepare_calc_result[
+                        "temp_features_dir"
+                    ]
+
+                    # Set status to calc_is_busy so we know calculation is busy...
+                    image["status"] = "IMAGE_CALC_BUSY"
+                    image["calc_starttime"] = datetime.now()
+
+                    # Now loop through all prepared feature batches to start the
+                    # statistics calculation for each
+                    logger.info(f"Start statistics calculation for {image_path_str}")
+                    bands = image["bands_todo"]
+                    # If bands shouldn't be calculated in parallel...
+                    if not calc_bands_parallel:
+                        bands = [bands]
 
-                        # Add info about the result of the prepare_calc to image info
-                        image["image_prepared_path"] = prepare_calc_result[
-                            "image_prepared_path"
-                        ]
-                        image["feature_batches"] = prepare_calc_result[
-                            "feature_batches"
-                        ]
-                        image["nb_features_to_calc_total"] = prepare_calc_result[
-                            "nb_features_to_calc_total"
-                        ]
-                        image["temp_features_dir"] = prepare_calc_result[
-                            "temp_features_dir"
-                        ]
-
-                        # Set status to calc_is_busy so we know calculation is busy...
-                        image["status"] = "IMAGE_CALC_BUSY"
-                        image["calc_starttime"] = datetime.now()
-
-                        # Now loop through all prepared feature batches to start the
-                        # statistics calculation for each
-                        logger.info(
-                            f"Start statistics calculation for {image_path_str}"
-                        )
+                    for band in bands:
                         for features_batch in image["feature_batches"]:
                             start_time_batch = datetime.now()
                             future = pool.submit(
-                                calc_stats_image_gdf,
+                                _zonal_stats_image_gdf,
                                 features_batch["path"],
                                 id_column,
                                 image["image_prepared_path"],
-                                bands,
+                                band,
                                 image["output_base_busy_path"],
                                 log_dir,
                                 log_level,
                                 start_time_batch,
+                                cloud_filter_band,
                             )
-                            calc_stats_batch_dict[features_batch["path"]] = {
+                            band_id = band if isinstance(band, str) else "-".join(band)
+                            batch_id = f"{features_batch['id']}_{band_id}"
+                            calc_stats_batch_dict[batch_id] = {
                                 "calc_stats_future": future,
-                                "image_path": image_path_str,
+                                "image_path": image["image_path"],
                                 "image_prepared_path": image["image_prepared_path"],
                                 "start_time_batch": start_time_batch,
                                 "nb_items_batch": features_batch["nb_items"],
                                 "status": "BATCH_CALC_BUSY",
                             }
 
-                    except Exception as ex:
-                        message = (
-                            f"Exception getting result of prepare_calc for {image}"
-                        )
-                        logger.exception(message)
-                        nb_errors += 1
-                        if nb_errors > nb_errors_max:
-                            raise Exception(message) from ex
+                except Exception as ex:
+                    message = f"Exception getting result of prepare_calc for {image}"
+                    logger.exception(message)
+                    nb_errors += 1
+                    if nb_errors > nb_errors_max:
+                        raise Exception(message) from ex
+
+            # For batches with calc busy, check if there are ready ones
+            for calc_stats_batch_id in _filter_on_status(
+                calc_stats_batch_dict, "BATCH_CALC_BUSY"
+            ):
+                # If it is not done yet, continue
+                calc_stats_batch_info = calc_stats_batch_dict[calc_stats_batch_id]
+                if calc_stats_batch_info["calc_stats_future"].done() is False:
+                    continue
 
-                # For batches with calc busy, check if ready
+                try:
+                    # Get the result
+                    result = calc_stats_batch_info["calc_stats_future"].result()
+                    if not result:
+                        raise Exception("Returned False?")
+
+                    # Set the processed flag to True
+                    calc_stats_batch_info["status"] = "BATCH_CALC_DONE"
+                    logger.debug(
+                        "Ready processing batch of "
+                        f"{calc_stats_batch_info['nb_items_batch']} for features "
+                        f"image: {calc_stats_batch_info['image_path']}"
+                    )
 
-                # Loop through the busy batche calculations
-                for calc_stats_batch_id in filter_on_status(
+                except Exception as ex:
+                    message = (
+                        "Exception getting result of calc_stats_image_gdf for "
+                        f"{calc_stats_batch_info}"
+                    )
+                    logger.exception(message)
+                    nb_errors += 1
+                    if nb_errors > nb_errors_max:
+                        raise Exception(message) from ex
+
+            # For images with calc busy, check if they are ready
+
+            # Loop through busy image calculations
+            for image_path_str in _filter_on_status(image_dict, "IMAGE_CALC_BUSY"):
+                # If still batches busy for this image, continue to next image
+                batches_busy = False
+                for calc_stats_batch_id in _filter_on_status(
                     calc_stats_batch_dict, "BATCH_CALC_BUSY"
                 ):
-                    # If it is not done yet, continue
-                    calc_stats_batch_info = calc_stats_batch_dict[calc_stats_batch_id]
-                    if calc_stats_batch_info["calc_stats_future"].done() is False:
-                        continue
-
-                    try:
-                        # Get the result
-                        result = calc_stats_batch_info["calc_stats_future"].result()
-                        if not result:
-                            raise Exception("Returned False?")
-
-                        # Set the processed flag to True
-                        calc_stats_batch_info["status"] = "BATCH_CALC_DONE"
-                        logger.debug(
-                            "Ready processing batch of "
-                            f"{calc_stats_batch_info['nb_items_batch']} for features "
-                            f"image: {calc_stats_batch_info['image_path']}"
-                        )
-
-                    except Exception as ex:
-                        message = (
-                            "Exception getting result of calc_stats_image_gdf for "
-                            f"{calc_stats_batch_info}"
-                        )
-                        logger.exception(message)
-                        nb_errors += 1
-                        if nb_errors > nb_errors_max:
-                            raise Exception(message) from ex
-
-                # For images with calc busy, check if they are ready
-
-                # Loop through busy image calculations
-                for image_path_str in filter_on_status(image_dict, "IMAGE_CALC_BUSY"):
-                    # If still batches busy for this image, continue to next image
-                    batches_busy = False
-                    for calc_stats_batch_id in filter_on_status(
-                        calc_stats_batch_dict, "BATCH_CALC_BUSY"
-                    ):
-                        if (
-                            calc_stats_batch_dict[calc_stats_batch_id]["image_path"]
-                            == image_path_str
-                        ):
-                            batches_busy = True
-                            break
-                    if batches_busy:
-                        continue
-
-                    # All batches are done, so the image is done
-                    image = image_dict[image_path_str]
-                    image["status"] = "IMAGE_CALC_DONE"
-
-                    # If the preprocessing created a temp image file, clean it up
-                    image_prepared_path = image["image_prepared_path"]
-                    if image_prepared_path != image_path_str:
-                        logger.info(
-                            f"Remove local temp image copy: {image_prepared_path}"
-                        )
-                        if image_prepared_path.isdir():
-                            shutil.rmtree(image_prepared_path, ignore_errors=True)
-                        else:
-                            os.remove(image_prepared_path)
+                    path = calc_stats_batch_dict[calc_stats_batch_id]["image_path"]
+                    if str(path) == image_path_str:
+                        batches_busy = True
+                        break
+                if batches_busy:
+                    continue
 
-                    # If the preprocessing created temp pickle files with features,
-                    # clean them up
-                    shutil.rmtree(image["temp_features_dir"], ignore_errors=True)
+                # All batches are done, so the image is done
+                image = image_dict[image_path_str]
+                image["status"] = "IMAGE_CALC_DONE"
 
-                    # Rename the (completed) output files
-                    output_base_path_noext, output_base_ext = os.path.splitext(
-                        image["output_base_path"]
+                # If the preprocessing created a temp image file, clean it up
+                image_prepared_path = image["image_prepared_path"]
+                if str(image_prepared_path) != str(image_path_str):
+                    logger.info(f"Remove local temp image copy: {image_prepared_path}")
+                    if image_prepared_path.is_dir():
+                        shutil.rmtree(image_prepared_path, ignore_errors=True)
+                    else:
+                        image_prepared_path.unlink()
+
+                # If the preprocessing created temp pickle files with features,
+                # clean them up
+                # shutil.rmtree(image["temp_features_dir"], ignore_errors=True)
+
+                # Move the (completed) output files
+                output_base_path = Path(image["output_base_path"])
+                output_base_busy_path = Path(image["output_base_busy_path"])
+                for band in image["bands_todo"]:
+                    # TODO: creating the output paths should probably be cleaner
+                    # centralised
+                    output_band_path = (
+                        output_base_path.parent
+                        / f"{output_base_path.stem}_{band}{output_base_path.suffix}"
                     )
-                    (
-                        output_base_busy_noext,
-                        output_base_busy_ext,
-                    ) = os.path.splitext(image["output_base_busy_path"])
-                    for band in bands:
-                        # TODO: creating the output paths should probably be cleaner
-                        # centralised
-                        output_band_busy_path = Path(
-                            f"{output_base_busy_noext}_{band}{output_base_busy_ext}"
-                        )
-                        output_band_path = Path(
-                            f"{output_base_path_noext}_{band}{output_base_ext}"
-                        )
-
-                        # If BUSY output file exists, rename it
-                        if output_band_busy_path.exists():
-                            os.rename(output_band_busy_path, output_band_path)
-                        else:
-                            # If BUSY output file doesn't exist, create empty file
-                            logger.info(
-                                "No features found overlapping image after processing, "
-                                f"create done file: {output_band_path}"
-                            )
-                            io_util.create_file_atomic(output_band_path)
-
-                    # Log the progress and prediction speed
-                    logger.info(f"Ready processing image: {image_path_str}")
-                    nb_done_latestbatch = 1
-                    nb_done_total += nb_done_latestbatch
-                    progress_msg = format_progress_message(
-                        nb_todo,
-                        nb_done_total,
-                        nb_done_latestbatch,
-                        start_time,
-                        image["calc_starttime"],
+                    suffix = output_base_busy_path.suffix
+                    output_band_busy_path = (
+                        output_base_busy_path.parent
+                        / f"{output_base_busy_path.stem}_{band}{suffix}"
                     )
-                    logger.info(progress_msg)
-
-                # Check if we are completely ready
-
-                # This is the case if:
-                #     - no processing is needed (= empty image_dict)
-                #     - OR if all processing is started + is done
-                #       (= status IMAGE_CALC_DONE)
-                if len(image_dict) == 0 or (
-                    i == len(image_paths)
-                    and len(image_dict)
-                    == len(filter_on_status(image_dict, "IMAGE_CALC_DONE"))
-                ):
-                    if nb_errors == 0:
-                        # Ready, so jump out of unending while loop
-                        break
+
+                    # If BUSY output file exists, move it (rename cannot move to
+                    # different file volume)
+                    if output_band_busy_path.exists():
+                        shutil.move(output_band_busy_path, output_band_path)
                     else:
-                        raise Exception(
-                            f"Ready processing, but there were {nb_errors} errors!"
+                        # If BUSY output file doesn't exist, create empty file
+                        logger.info(
+                            "No features found overlapping image after processing, "
+                            f"create done file: {output_band_path}"
                         )
+                        io_util.create_file_atomic(output_band_path)
+
+                # Log the progress and prediction speed
+                logger.info(f"Ready processing image: {image_path_str}")
+                nb_done_latestbatch = 1
+                nb_done_total += nb_done_latestbatch
+                progress_msg = _format_progress_message(
+                    nb_todo,
+                    nb_done_total,
+                    nb_done_latestbatch,
+                    start_time,
+                    image["calc_starttime"],
+                )
+                logger.info(progress_msg)
+
+            # Check if we are completely ready
+
+            # This is the case if:
+            #     - no processing is needed (= empty image_dict)
+            #     - OR if all (possible) processing is started + everything is done
+            #       (= status IMAGE_CALC_DONE)
+            if len(image_dict) == 0 or (
+                image_idx == len(images_bands)
+                and len(image_dict)
+                == len(_filter_on_status(image_dict, "IMAGE_CALC_DONE"))
+            ):
+                if nb_errors == 0:
+                    # Ready, so jump out of unending while loop
+                    break
+                else:
+                    raise Exception(
+                        f"Ready processing, but there were {nb_errors} errors!"
+                    )
 
-                # Some sleep before starting next iteration...
-                time.sleep(0.1)
+            # Some sleep before starting next iteration...
+            time.sleep(0.1)
 
-        except KeyboardInterrupt:
-            # If CTRL+C is used, shut down pool and kill children
-            print("You pressed Ctrl+C")
-            print("Worker processes are being stopped, followed by exit!")
+    except KeyboardInterrupt:
+        # If CTRL+C is used, shut down pool and kill children
+        print("You pressed Ctrl+C")
+        print("Worker processes are being stopped, followed by exit!")
 
-            # Stop process pool + kill children + exit
-            try:
-                pool.shutdown(wait=False)
-                parent = psutil.Process(os.getpid())
-                children = parent.children(recursive=True)
-                for process_pid in children:
-                    print(f"Kill child with pid {process_pid}")
-                    process_pid.send_signal(signal.SIGTERM)
-            finally:
-                sys.exit(1)
-
-        logger.info(
-            f"Time taken to calculate data for {nb_todo} images: "
-            f"{(datetime.now()-start_time).total_seconds()} sec"
-        )
+        # Stop process pool + kill children + exit
+        try:
+            pool.shutdown(wait=False)
+            parent = psutil.Process(os.getpid())
+            children = parent.children(recursive=True)
+            for process_pid in children:
+                print(f"Kill child with pid {process_pid}")
+                process_pid.send_signal(signal.SIGTERM)
+        finally:
+            sys.exit(1)
+
+    finally:
+        pool.shutdown()
+
+    logger.info(
+        f"Time taken to calculate data for {nb_todo} images: "
+        f"{(datetime.now()-start_time).total_seconds()} sec"
+    )
 
 
-def filter_on_status(dict: dict, status_to_check: str) -> List[str]:
+def _filter_on_status(dict: dict, status_to_check: str) -> List[str]:
     """
     Function to check the number of images that are being prepared for processing
     """
@@ -451,7 +470,7 @@ def filter_on_status(dict: dict, status_to_check: str) -> List[str]:
     return keys_with_status
 
 
-def format_output_path(
+def _format_output_path(
     features_path: Path,
     image_path: Path,
     output_dir: Path,
@@ -461,14 +480,6 @@ def format_output_path(
     """
     Prepare the output path.
     """
-    # Filename of the features
-    features_filename = os.path.splitext(os.path.basename(features_path))[0]
-    # Filename of the image -> remove .zip if it is a zip file
-    image_filename = os.path.basename(image_path)
-    image_filename_noext, image_ext = os.path.splitext(image_filename)
-    if image_ext.lower() == ".zip":
-        image_filename = image_filename_noext
-
     # Interprete the orbit...
     if orbit_properties_pass is not None:
         if orbit_properties_pass == "ASCENDING":
@@ -482,21 +493,15 @@ def format_output_path(
     else:
         orbit = ""
 
-    # Prepare basic path
-    output_path_noext = os.path.join(
-        output_dir, f"{features_filename}__{image_filename}{orbit}"
-    )
-
-    # If a band was specified
+    # Format + return output path
+    output_stem = f"{features_path.stem}__{image_path.stem}{orbit}"
     if band is not None:
-        output_path_noext = f"{output_path_noext}_{band}"
-
-    # Add extension
-    output_path = Path(f"{output_path_noext}.sqlite")
+        output_stem = f"{output_stem}_{band}"
+    output_path = output_dir / f"{output_stem}.sqlite"
     return output_path
 
 
-def format_progress_message(
+def _format_progress_message(
     nb_todo: int,
     nb_done_total: int,
     nb_done_latestbatch: int,
@@ -543,11 +548,10 @@ def format_progress_message(
     return message
 
 
-def prepare_calc(
+def _prepare_calc(
     features_path: Path,
     id_column: str,
     image_path: Path,
-    output_path: Path,
     temp_dir: Path,
     log_dir: Path,
     log_level: Union[str, int],
@@ -562,48 +566,51 @@ def prepare_calc(
     """
     # When running in parallel processes, the logging needs to be write to seperate
     # files + no console logging
+    if len(logging.getLogger().handlers) == 0:
+        log_path = (
+            log_dir
+            / f"{datetime.now():%Y-%m-%d_%H-%M-%S}_prepare_calc_{os.getpid()}.log"
+        )
+        filehandler = logging.FileHandler(filename=log_path)
+        filehandler.setLevel(log_level)
+        log_format = "%(asctime)s|%(levelname)s|%(name)s|%(funcName)s|%(message)s"
+        logging.basicConfig(
+            level=logging.DEBUG,
+            format=log_format,
+            handlers=[filehandler],
+        )
+
     global logger
     logger = logging.getLogger("prepare_calc")
-    logger.propagate = False
-    log_path = (
-        log_dir / f"{datetime.now():%Y-%m-%d_%H-%M-%S}_prepare_calc_{os.getpid()}.log"
-    )
-    fh = logging.FileHandler(filename=log_path)
-    fh.setLevel(log_level)
-    fh.setFormatter(
-        logging.Formatter("%(asctime)s|%(levelname)s|%(name)s|%(funcName)s|%(message)s")
-    )
-    logger.addHandler(fh)
-
     ret_val = {}
 
     # Prepare the image
     logger.info(f"Start prepare_image for {image_path} to {temp_dir}")
-    image_prepared_path = raster_helper.prepare_image(image_path, temp_dir)
+    image_prepared_path = _image_helper.prepare_image(image_path, temp_dir)
     logger.debug(f"Preparing ready, result: {image_prepared_path}")
     ret_val["image_prepared_path"] = image_prepared_path
 
     # Get info about the image
     logger.info(f"Start get_image_info for {image_prepared_path}")
-    image_info = raster_helper.get_image_info(image_prepared_path)
+    image_info = _image_helper.get_image_info(image_prepared_path)
     logger.info(f"image_info: {image_info}")
 
     # Load the features that overlap with the image.
     # TODO: passing both bbox and poly is double, or not?
     # footprint epsg should be passed as well, or reproject here first?
     footprint_shape = None
-    if "footprint" in image_info:
-        logger.info(f"poly: {image_info['footprint']['shape']}")
-        footprint_shape = image_info["footprint"]["shape"]
+    if image_info.footprint is not None:
+        logger.info(f"poly: {image_info.footprint['shape']}")
+        footprint_shape = image_info.footprint["shape"]
 
-    if image_info["image_epsg"] == "NONE":
+    if image_info.image_epsg == "NONE":
         raise Exception(f"target_epsg == NONE: {image_info}")
 
-    features_gdf = load_features_file(
+    features_gdf = _load_features_file(
         features_path=features_path,
-        target_epsg=image_info["image_epsg"],
+        target_epsg=image_info.image_epsg,
         columns_to_retain=[id_column, "geometry"],
-        bbox=image_info["image_bounds"],
+        bbox=image_info.image_bounds,
         polygon=footprint_shape,
     )
 
@@ -627,25 +634,18 @@ def prepare_calc(
     # Pickle the batches to temporary files
     # Create temp dir to put the pickles in... and clean or create it.
     # TODO: change dir so it is always unique
-    temp_features_dirname = os.path.splitext(os.path.basename(image_path))[0]
-    temp_features_dir = os.path.join(temp_dir, temp_features_dirname)
+    temp_features_dir = temp_dir / image_path.stem
     ret_val["temp_features_dir"] = temp_features_dir
-    if os.path.exists(temp_features_dir):
-        logger.info(f"Remove dir {temp_features_dir + os.sep}")
-        shutil.rmtree(temp_features_dir + os.sep)
-    else:
-        os.makedirs(temp_features_dir, exist_ok=True)
-
-    """
-    if os.path.exists(temp_features_dir):
-        logger.info(f"Temp dir exists: {temp_features_dir}")
-    """
+    if temp_features_dir.exists():
+        logger.info(f"Remove dir {str(temp_features_dir)}{os.sep}")
+        shutil.rmtree(f"{str(temp_features_dir)}{os.sep}")
+    temp_features_dir.mkdir(parents=True, exist_ok=True)
 
     # Loop over the batches, pickle them and add the paths to the result...
     ret_val["feature_batches"] = []
-    for i, features_gdf_batch in enumerate(features_gdf_batches):
+    for batch_idx, features_gdf_batch in enumerate(features_gdf_batches):
         batch_info = {}
-        pickle_path = os.path.join(temp_features_dir, f"{i}.pkl")
+        pickle_path = temp_features_dir / f"{batch_idx}.pkl"
         logger.info(
             f"Write pkl of {len(features_gdf_batch.index)} features: {pickle_path}"
         )
@@ -655,16 +655,17 @@ def prepare_calc(
             raise Exception(f"Exception writing pickle: {pickle_path}: {ex}") from ex
         batch_info["path"] = pickle_path
         batch_info["nb_items"] = len(features_gdf_batch.index)
+        batch_info["id"] = f"{image_path.stem}_{batch_idx}"
         ret_val["feature_batches"].append(batch_info)
 
     # Ready... return
     return ret_val
 
 
-def load_features_file(
+def _load_features_file(
     features_path: Path,
     columns_to_retain: List[str],
-    target_epsg: str,
+    target_epsg: int,
     bbox=None,
     polygon=None,
 ) -> gpd.GeoDataFrame:
@@ -685,10 +686,21 @@ def load_features_file(
     # By convention, the features filename should end on the projection... so extract
     # epsg from filename
     start_time = datetime.now()
+    """
     features_epsg = None
     if features_path.stem.find("_") != -1:
         splitted = features_path.stem.split("_")
         features_epsg = splitted[len(splitted) - 1]
+        try:
+            features_epsg = int(features_epsg)
+        except Exception:
+            logger.info(
+                f"features_epsg {features_epsg } could not be cast to int for "
+                f"{features_path.name}"
+            )
+            features_epsg = None
+    """
+    features_epsg = gfo.get_layerinfo(features_path).crs.to_epsg()
 
     # Determine the correct filename for the input features in the correct projection.
     if features_epsg != target_epsg:
@@ -844,15 +856,16 @@ def load_features_file(
     return features_gdf
 
 
-def calc_stats_image_gdf(
-    features_gdf,
+def _zonal_stats_image_gdf(
+    features,
     id_column: str,
     image_path: Path,
-    bands: List[str],
+    bands: Union[List[str], str],
     output_base_path: Path,
     log_dir: Path,
     log_level: Union[str, int],
     future_start_time=None,
+    cloud_filter_band: Optional[str] = None,
 ) -> bool:
     """
     Calculate stats for an image.
@@ -866,19 +879,22 @@ def calc_stats_image_gdf(
     # this function should process only one band!
     # When running in parallel processes, the logging needs to be write to seperate
     # files + no console logging
+    if len(logging.getLogger().handlers) == 0:
+        log_path = (
+            log_dir
+            / f"{datetime.now():%Y-%m-%d_%H-%M-%S}_calc_stats_image_{os.getpid()}.log"
+        )
+        filehandler = logging.FileHandler(filename=log_path)
+        filehandler.setLevel(log_level)
+        log_format = "%(asctime)s|%(levelname)s|%(name)s|%(funcName)s|%(message)s"
+        logging.basicConfig(
+            level=logging.DEBUG,
+            format=log_format,
+            handlers=[filehandler],
+        )
+
     global logger
-    logger = logging.getLogger("calc_stats_image_gdf")
-    logger.propagate = False
-    log_path = (
-        log_dir
-        / f"{datetime.now():%Y-%m-%d_%H-%M-%S}_calc_stats_image_gdf_{os.getpid()}.log"
-    )
-    fh = logging.FileHandler(filename=str(log_path))
-    fh.setLevel(log_level)
-    fh.setFormatter(
-        logging.Formatter("%(asctime)s|%(levelname)s|%(name)s|%(funcName)s|%(message)s")
-    )
-    logger.addHandler(fh)
+    logger = logging.getLogger("calc_stats_image")
 
     # Log the time between scheduling the future and acually run...
     if future_start_time is not None:
@@ -888,20 +904,24 @@ def calc_stats_image_gdf(
         )
 
     # If the features_gdf is a string, use it as file path to unpickle geodataframe...
-    if isinstance(features_gdf, str):
-        features_gdf_pkl_path = features_gdf
+    if isinstance(features, str) or isinstance(features, Path):
+        features_gdf_pkl_path = features
         logger.info(f"Read pickle: {features_gdf_pkl_path}")
-        features_gdf = pd.read_pickle(features_gdf_pkl_path)
-        logger.info(f"Read pickle with {len(features_gdf.index)} features ready")
+        features = pd.read_pickle(features_gdf_pkl_path)
+        logger.info(f"Read pickle with {len(features.index)} features ready")
+
+    if isinstance(bands, str):
+        bands = [bands]
 
     # Init some variables
-    quality_band = "SCL-20m"
-    features_total_bounds = features_gdf.total_bounds
+    features_total_bounds = features.total_bounds
     output_base_path_noext, output_ext = os.path.splitext(output_base_path)
 
     # If the image has a quality band, check that one first so parcels with
     # bad pixels can be removed as the data can't be trusted anyway
-    if quality_band in bands:
+    if cloud_filter_band is not None:
+        if cloud_filter_band != "SCL-20m":
+            raise ValueError(f"cloud_filter_band not supported: {cloud_filter_band}")
         # Specific for Scene Classification (SCL) band, interprete already
         # Folowing values are considered "bad":
         #   -> 0 (no_data), 1 (saturated or defective), 3 (cloud shadows),
@@ -914,8 +934,8 @@ def calc_stats_image_gdf(
 
         # Get the image data and calculate
         logger.info(
-            f"Calculate categorical counts for band {quality_band} on "
-            f"{len(features_gdf.index)} features"
+            f"Calculate categorical counts for band {cloud_filter_band} on "
+            f"{len(features.index)} features"
         )
         category_map = {
             0.0: "nodata",
@@ -940,19 +960,19 @@ def calc_stats_image_gdf(
             "cloud_mediumproba",
             "cloud_highproba",
         ]
-        image_data = raster_helper.get_image_data(
+        image_data = _image_helper.get_image_data(
             image_path,
             bounds=features_total_bounds,
-            bands=[quality_band],
+            bands=[cloud_filter_band],
             pixel_buffer=1,
         )
         # Before zonal_stats, do reset_index so index is "clean", otherwise the
         # concat/insert/... later on gives wrong results
-        features_gdf.reset_index(drop=True, inplace=True)
-        features_stats = zonal_stats(
-            features_gdf,
-            image_data[quality_band]["data"],
-            affine=image_data[quality_band]["transform"],
+        features.reset_index(drop=True, inplace=True)
+        features_stats = rasterstats.zonal_stats(
+            features,
+            image_data[cloud_filter_band]["data"],
+            affine=image_data[cloud_filter_band]["transform"],
             prefix="",
             nodata=0,
             categorical=True,
@@ -979,8 +999,11 @@ def calc_stats_image_gdf(
         ].sum(axis=1)
 
         # Add index and write to file
-        features_stats_df.insert(loc=0, column=id_column, value=features_gdf[id_column])
-        output_band_path = Path(f"{output_base_path_noext}_{quality_band}{output_ext}")
+        features_stats_df.insert(loc=0, column=id_column, value=features[id_column])
+        output_band_path = output_base_path.with_stem(
+            f"{output_base_path.stem}_{cloud_filter_band}"
+        )
+
         logger.info(
             f"Write data for {len(features_stats_df.index)} parcels found to "
             f"{output_band_path}"
@@ -988,34 +1011,30 @@ def calc_stats_image_gdf(
         pdh.to_file(features_stats_df, output_band_path, index=False, append=True)
 
         # Use the nb_bad_pixels column to filter only parcels without bad pixels
-        features_gdf.insert(
+        features.insert(
             loc=0,
             column=nb_bad_pixels_column,
             value=features_stats_df[nb_bad_pixels_column],
         )
-        features_gdf = features_gdf.loc[features_gdf[nb_bad_pixels_column] == 0]
-        features_gdf.drop(columns=[nb_bad_pixels_column], inplace=True)
+        features = features.loc[features[nb_bad_pixels_column] == 0]
+        features.drop(columns=[nb_bad_pixels_column], inplace=True)
 
         # Check if there are still features to be calculated
-        if len(features_gdf.index) == 0:
+        if len(features.index) == 0:
             logger.info(
                 "After checking quality band, no more features to be calculated: stop"
             )
             return True
 
     # Loop over image bands
-    features_total_bounds = features_gdf.total_bounds
+    features_total_bounds = features.total_bounds
     for band in bands:
-        # The quality band is already treated, so skip it here
-        if band == "SCL-20m":
-            continue
-
         # Check if output file exists already...
-        output_band_path = Path(f"{output_base_path_noext}_{band}{output_ext}")
+        output_band_path = output_base_path.with_stem(f"{output_base_path.stem}_{band}")
 
         # Get the image data and calculate statistics
         logger.info(f"Read band {band} for bounds {features_total_bounds}")
-        image_data = raster_helper.get_image_data(
+        image_data = _image_helper.get_image_data(
             image_path, bounds=features_total_bounds, bands=[band], pixel_buffer=1
         )
 
@@ -1033,17 +1052,17 @@ def calc_stats_image_gdf(
 
         logger.info(
             f"Calculate zonal statistics for band {band} on "
-            f"{len(features_gdf.index)} features"
+            f"{len(features.index)} features"
         )
         # Before zonal_stats, do reset_index so index is "clean", otherwise the
         # concat/insert/... later on gives wrong results
-        features_gdf.reset_index(drop=True, inplace=True)
-        features_stats = zonal_stats(
-            features_gdf,
+        features.reset_index(drop=True, inplace=True)
+        features_stats = rasterstats.zonal_stats(
+            features,
             image_data_upsampled,
             affine=affine_upsampled,
             prefix="",
-            nodata=0,
+            nodata=image_data[band]["nodata"],
             all_touched=False,
             stats=["count", "mean", "median", "std", "min", "max"],
         )
@@ -1053,7 +1072,7 @@ def calc_stats_image_gdf(
         )
 
         # Add original id column to statistics dataframe
-        features_stats_df.insert(loc=0, column=id_column, value=features_gdf[id_column])
+        features_stats_df.insert(loc=0, column=id_column, value=features[id_column])
 
         # Remove rows with empty data
         features_stats_df.dropna(inplace=True)
diff --git a/cropclassification/util/zonal_stats_bulk/_image_helper.py b/cropclassification/util/zonal_stats_bulk/_image_helper.py
new file mode 100644
index 00000000..e12a5c32
--- /dev/null
+++ b/cropclassification/util/zonal_stats_bulk/_image_helper.py
@@ -0,0 +1,776 @@
+import logging
+import os
+from pathlib import Path
+import time
+from typing import Dict, List, Optional, Tuple
+import xml.etree.ElementTree as ET
+
+from osgeo import gdal
+
+# Suppress gdal warnings/errors
+gdal.PushErrorHandler("CPLQuietErrorHandler")
+
+import rasterio
+import shapely.geometry as sh_geom
+
+from cropclassification.util import io_util
+
+# General init
+logger = logging.getLogger(__name__)
+
+
+class BandInfo:
+    def __init__(
+        self,
+        path: str,
+        relative_path: str,
+        filename: str,
+        bandindex: int,
+        bounds: Optional[Tuple[float, float, float]] = None,
+        affine=None,
+        crs: Optional[str] = None,
+        epsg: Optional[int] = None,
+    ):
+        self.path = path
+        self.relative_path = relative_path
+        self.filename = filename
+        self.bandindex = bandindex
+        self.bounds = bounds
+        self.affine = affine
+        self.crs = crs
+        self.epsg = epsg
+
+
+class ImageInfo:
+    def __init__(
+        self,
+        imagetype: str,
+        filetype: str,
+        image_id: str,
+        filename: str,
+        footprint: dict,
+        image_epsg: int,
+        image_crs: str,
+        image_bounds: Tuple[float, float, float],
+        image_affine,
+        bands: Dict[str, BandInfo],
+        extra: dict,
+    ):
+        self.imagetype = imagetype
+        self.filetype = filetype
+        self.image_id = image_id
+        self.filename = filename
+        self.footprint = footprint
+        self.image_epsg = image_epsg
+        self.image_crs = image_crs
+        self.image_bounds = image_bounds
+        self.image_affine = image_affine
+        self.bands = bands
+        self.extra = extra
+
+
+def get_image_data(
+    image_path: Path,
+    bounds: Tuple[float, float, float, float],
+    bands: List[str],
+    pixel_buffer: int = 0,
+) -> dict:
+    """
+    Reads the data from the image.
+
+    Adds a small buffer around the bounds asked to evade possible rounding issues.
+
+    Returns a dict of the following structure:
+    imagedata[band]['data']: the data read as numpy array
+                   ['transform']: the Affine transform of the band read
+
+    Args
+        image_path: the path the image
+        bounds: the bounds to be read, in coordinates in the projection of the image
+        bands: list of bands to be read, eg. "VV", "VH",...
+        pixel_buffer: number to pixels to take as buffer around the bounds provided in
+            pixels
+    """
+    # Get info about the image
+    logger.debug("get_image_data start")
+    image_info = get_image_info(image_path)
+
+    # Now get the data
+    image_data = {}  # Dict for the transforms and the data per band
+    if image_info.filetype in ("CARD", "SAFE", "TIF"):
+        # Loop over bands and get data
+        for band in bands:
+            band_relative_path = image_info.bands[band].relative_path
+            if band_relative_path is not None:
+                image_band_path = image_path / band_relative_path
+            else:
+                image_band_path = image_path
+            band_index = image_info.bands[band].bandindex
+            logger.info(
+                f"Read image data from {image_band_path}, with band_index: {band_index}"
+            )
+            image_data[band] = {}
+            with rasterio.open(str(image_band_path)) as src:
+                # Determine the window we need to read from the image:
+                window_to_read = projected_bounds_to_window(
+                    bounds, src.transform, src.width, src.height, pixel_buffer
+                )
+                image_data[band][
+                    "transform"
+                ] = rasterio.windows.transform(  # type: ignore
+                    window_to_read, src.transform
+                )
+                image_data[band]["nodata"] = src.nodata
+                # Read!
+                # Remark: bandindex in rasterio is 1-based instead of 0-based -> +1
+                logger.debug(f"Read image data from {image_band_path}")
+                image_data[band]["data"] = src.read(
+                    band_index + 1, window=window_to_read
+                )
+                logger.debug("Image data read")
+    else:
+        message = f"Format currently not supported: {image_path}"
+        logger.error(message)
+        raise NotImplementedError(message)
+
+    logger.debug("get_image_data ready")
+
+    return image_data
+
+
+def get_image_info(image_path: Path) -> ImageInfo:
+    """
+    Returns basic information about an image.
+
+    Args:
+        image_path (Path): the path to the image.
+
+    Raises:
+        NotImplementedError: image is of an unsupported image type.
+
+    Returns:
+        ImageInfo: basic information about the image
+    """
+
+    image_info = {}
+
+    # Specific code per image type
+    suffix_upper = image_path.suffix.upper()
+    if suffix_upper == ".CARD":
+        image_info = _get_image_info_card(image_path)
+    elif suffix_upper == ".SAFE":
+        # This is a level 2 sentinel 2 file
+        image_info = _get_image_info_safe(image_path)
+    elif suffix_upper == ".TIF":
+        image_info = _get_image_info_tif(image_path)
+    else:
+        message = f"Not a supported image format: {image_path}"
+        logger.error(message)
+        raise NotImplementedError(message)
+
+    return image_info
+
+
+def projected_bounds_to_window(
+    projected_bounds,
+    image_transform,
+    image_pixel_width: int,
+    image_pixel_height: int,
+    pixel_buffer: int = 0,
+):
+    """
+    Returns a rasterio.windows.Window to be used in rasterio to read the part of the
+    image specified.
+
+    Args
+        projected_bounds: bounds to created the window from, in projected coordinates
+        image_transform: Affine transform of the image you want to create the pixel
+            window for
+        image_pixel_width: total width of the image you want to create the pixel window
+            for, in pixels
+        image_pixel_height: total height of the image you want to create the pixel
+            window for, in pixels
+        pixel_buffer: number to pixels to take as buffer around the bounds provided in
+            pixels
+    """
+    # Take bounds of the features + convert to image pixels
+    xmin, ymin, xmax, ymax = projected_bounds
+    window_to_read_raw = rasterio.windows.from_bounds(  # type: ignore
+        xmin, ymin, xmax, ymax, image_transform
+    )
+
+    # Round
+    window_to_read = window_to_read_raw.round_offsets().round_lengths()
+
+    # Now some general math on window properties, but as they are readonly properties,
+    # work on copy
+    col_off, row_off, width, height = window_to_read.flatten()
+    # Add buffer of 1 pixel extra around
+    col_off -= pixel_buffer
+    row_off -= pixel_buffer
+    width += 2 * pixel_buffer
+    height += 2 * pixel_buffer
+
+    # Make sure the offsets aren't negative, as the pixels that are 'read' there acually
+    # get some value instead of eg. nan...!
+    if col_off < 0:
+        width -= abs(col_off)
+        col_off = 0
+    if row_off < 0:
+        height -= abs(row_off)
+        row_off = 0
+
+    # Make sure there won't be extra pixels to the top and right that will be read
+    if (col_off + width) > image_pixel_width:
+        width = image_pixel_width - col_off
+    if (row_off + height) > image_pixel_height:
+        height = image_pixel_height - row_off
+
+    # Ready... prepare to return...
+    window_to_read = rasterio.windows.Window(  # type: ignore
+        col_off, row_off, width, height
+    )
+
+    # Debug info
+    """
+    bounds_to_read = rasterio.windows.bounds(window_to_read, image_transform)
+    logger.debug(f"projected_bounds: {projected_bounds}, "
+        f"window_to_read_raw: {window_to_read_raw}, window_to_read: {window_to_read}, "
+        f"image_pixel_width: {image_pixel_width}, "
+        f"image_pixel_height: {image_pixel_height}, "
+        f"file transform: {image_transform}, bounds_to_read: {bounds_to_read}"
+    )
+    """
+
+    return window_to_read
+
+
+def prepare_image(image_path: Path, temp_dir: Path) -> Path:
+    """
+    Prepares the input image for usages.
+
+    In case of a zip file, the file is unzipped to the temp dir specified.
+
+    Returns the path to the prepared file/directory.
+    """
+
+    # If the input path is not a zip file, don't make local copy + just return path
+    if image_path.suffix.lower() != ".zip":
+        return image_path
+    else:
+        # It is a zip file, so it needs to be unzipped first...
+        # Create destination file path
+        image_basename_withzipext = os.path.basename(image_path)
+        image_basename = os.path.splitext(image_basename_withzipext)[0]
+        image_unzipped_path = temp_dir / image_basename
+
+        image_unzipped_path_busy = Path(f"{image_unzipped_path}_busy")
+        # If the input is a zip file, unzip file to temp local location if it doesn't
+        # exist yet. If the file doesn't exist yet in right projection, read original
+        # input file to reproject/write to new file with correct epsg
+        if not (image_unzipped_path_busy.exists() or image_unzipped_path.exists()):
+            # Create temp dir if it doesn't exist yet
+            os.makedirs(temp_dir, exist_ok=True)
+
+            # Create lock file in an atomic way, so we are sure we are the only process
+            # working on it. If function returns true, there isn't any other
+            # thread/process already working on it
+            if io_util.create_file_atomic(image_unzipped_path_busy):
+                try:
+                    logger.info(
+                        f"Unzip image {image_path} to local location {temp_dir}"
+                    )
+
+                    # Create the dest dir where the file will be unzipped to + unzip!
+                    if not image_unzipped_path.exists():
+                        import zipfile
+
+                        with zipfile.ZipFile(image_path, "r") as zippedfile:
+                            zippedfile.extractall(temp_dir)
+                finally:
+                    # Remove lock file when we are ready
+                    os.remove(image_unzipped_path_busy)
+
+        # If a "busy file" still exists, the file isn't ready yet, but another process
+        # is working on it, so wait till it disappears
+        while image_unzipped_path_busy.exists():
+            time.sleep(1)
+
+        # Now we are ready to return the path...
+        return image_unzipped_path
+
+
+def _get_image_info_card(image_path: Path) -> ImageInfo:
+    # This is a sentinel 1 image (GRD or coherence)
+    # First extract and fill out some basic info
+    imagetype = image_path.stem.split("_")[0]
+    filetype = "CARD"
+    image_id = image_path.stem
+
+    # Read info from the metadata file
+    metadata_xml_path = image_path / "metadata.xml"
+    metadata = ET.parse(str(metadata_xml_path))
+    metadata_root = metadata.getroot()
+
+    logger.debug(f"Parse metadata info from {metadata_xml_path}")
+
+    try:
+        # Get the filename
+        filename = metadata_root.find("productFilename").text  # type: ignore
+
+        # Get the footprint
+        footprint = {}
+        footprint_elem = metadata_root.find("footprint")
+        assert footprint_elem is not None
+        poly = footprint_elem.find("{http://www.opengis.net/gml}Polygon")
+        assert poly is not None
+        footprint["srsname"] = poly.attrib.get("srsName")
+        linear_ring = []
+        # for coord in poly.findall(
+        #   "{http://www.opengis.net/gml}outerBoundaryIs/"
+        #   "{http://www.opengis.net/gml}LinearRing/"
+        #   "{http://www.opengis.net/gml}coord"
+        # ):
+        #    linear_ring.append(
+        #        (float(coord.findtext(
+        #            "{http://www.opengis.net/gml}X"
+        #        )),
+        #        float(coord.findtext("{http://www.opengis.net/gml}Y")))
+        #    )
+        coord_str = poly.find(
+            "{http://www.opengis.net/gml}outerBoundaryIs/"
+            "{http://www.opengis.net/gml}LinearRing/"
+            "{http://www.opengis.net/gml}coordinates"
+        )
+        assert coord_str is not None and coord_str.text is not None
+        logger.debug(f"coord_str: {coord_str}, coord_str.text: {coord_str.text}")
+        coord_list = coord_str.text.split(" ")
+        for coord in coord_list:
+            # Watch out, latitude (~y) is first, than longitude (~x)
+            # TODO: add check if the projection is in degrees (latlon) or
+            # coordinates (xy) instead of hardcoded latlon
+            y, x = coord.split(",")
+            linear_ring.append((float(x), float(y)))
+        footprint["shape"] = sh_geom.polygon.Polygon(linear_ring)
+
+        # get epsg
+        epsg = metadata_root.find("imageProjection/EPSG")
+        assert epsg is not None and epsg.text is not None
+        image_epsg = float(epsg.text)
+        image_crs = f"EPSG:{epsg.text}"
+    except Exception as ex:
+        raise Exception(f"Exception extracting info from {metadata_xml_path}") from ex
+
+    # Read info from the manifest.safe file
+    manifest_xml_searchstring = "*_manifest.safe"
+    manifest_xml_paths = list(image_path.glob(manifest_xml_searchstring))
+
+    # The number of .safe indicates whether it is a GRD or a Coherence image
+    extra = {}
+    nb_safefiles = len(manifest_xml_paths)
+    if nb_safefiles == 1:
+        # Now parse the .safe file
+        manifest_xml_path = manifest_xml_paths[0]
+
+        try:
+            manifest = ET.parse(str(manifest_xml_path))
+            manifest_root = manifest.getroot()
+
+            # Define namespaces...
+            ns = {
+                "safe": "http://www.esa.int/safe/sentinel-1.0",
+                "s1": "http://www.esa.int/safe/sentinel-1.0/sentinel-1",
+                "s1sarl1": (
+                    "http://www.esa.int/safe/sentinel-1.0/sentinel-1/sar/level-1"
+                ),
+            }
+
+            logger.debug(f"Parse manifest info from {metadata_xml_path}")
+            extra["transmitter_receiver_polarisation"] = []
+            for polarisation in manifest_root.findall(
+                "metadataSection/metadataObject/metadataWrap/xmlData/"
+                "s1sarl1:standAloneProductInformation/"
+                "s1sarl1:transmitterReceiverPolarisation",
+                ns,
+            ):
+                extra["transmitter_receiver_polarisation"].append(polarisation.text)
+            extra["productTimelinessCategory"] = manifest_root.find(
+                "metadataSection/metadataObject/metadataWrap/xmlData/"
+                "s1sarl1:standAloneProductInformation/"
+                "s1sarl1:productTimelinessCategory",
+                ns,
+            ).text  # type: ignore
+
+            extra["instrument_mode"] = manifest_root.find(
+                "metadataSection/metadataObject/metadataWrap/xmlData/"
+                "safe:platform/safe:instrument/safe:extension/"
+                "s1sarl1:instrumentMode/s1sarl1:mode",
+                ns,
+            ).text  # type: ignore
+            extra["orbit_properties_pass"] = manifest_root.find(
+                "metadataSection/metadataObject/metadataWrap/xmlData/"
+                "safe:orbitReference/safe:extension/s1:orbitProperties/s1:pass",
+                ns,
+            ).text  # type: ignore
+            extra["acquisition_date"] = manifest_root.find(
+                "metadataSection/metadataObject/metadataWrap/xmlData/"
+                "safe:acquisitionPeriod/safe:startTime",
+                ns,
+            ).text  # type: ignore
+
+            # Now have a look in the files themselves to get band info,...
+            # TODO: probably cleaner/easier to read from metadata files...
+            image_basename_noext_nodate = image_path.stem[
+                : image_path.stem.rfind("_", 0)
+            ]
+            image_datadirname = f"{image_basename_noext_nodate}.data"
+            image_datadir = image_path / image_datadirname
+            band_paths = list(image_datadir.glob("*.img"))
+
+            # If no files were found, check if the images are in .tif format!
+            if len(band_paths) == 0:
+                tif_path = image_path / "Gamma0_VH.tif"
+                if tif_path.exists():
+                    band_paths.append(tif_path)
+                tif_path = image_path / "Gamma0_VV.tif"
+                if tif_path.exists():
+                    band_paths.append(tif_path)
+
+            # If no files were found, error
+            if len(band_paths) == 0:
+                message = f"No image files found in {image_datadir} or {image_path}!"
+                logger.error(message)
+                raise Exception(message)
+
+            bands = {}
+            for i, band_path in enumerate(band_paths):
+                # Extract bound,... info from the first file only
+                # (they are all the same)
+                if i == 0:
+                    # logger.debug(f"Read image metadata from {band_path}")
+                    with rasterio.open(str(band_path)) as src:
+                        image_bounds = src.bounds
+                        image_affine = src.transform
+                        if src.crs is not None:
+                            image_crs = src.crs.to_string()
+                            image_epsg = src.crs.to_epsg()
+                    # logger.debug(f"Image metadata read: {image_info}")
+                if band_path.stem == "Gamma0_VH":
+                    band = "VH"
+                elif band_path.stem == "Gamma0_VV":
+                    band = "VV"
+                else:
+                    raise NotImplementedError(f"Filename not supported: {band_path}")
+
+                # Add specific info about the band
+                bands[band] = BandInfo(
+                    path=str(band_path),
+                    relative_path=str(band_path).replace(
+                        str(image_path) + os.pathsep, ""
+                    ),
+                    filename=band_path.name,
+                    bandindex=0,
+                )
+
+        except Exception as ex:
+            raise Exception(
+                f"Exception extracting info from {manifest_xml_path}"
+            ) from ex
+
+    elif nb_safefiles == 2 or nb_safefiles == 3:
+        # 2 safe files -> coherence
+        # Now parse the first .safe file
+        # TODO: maybe check if the info in all safe files  are the same or?.?
+        manifest_xml_path = manifest_xml_paths[0]
+
+        try:
+            manifest = ET.parse(str(manifest_xml_path))
+            manifest_root = manifest.getroot()
+
+            # Define namespaces...
+            ns = {
+                "safe": "http://www.esa.int/safe/sentinel-1.0",
+                "s1": "http://www.esa.int/safe/sentinel-1.0/sentinel-1",
+                "s1sarl1": (
+                    "http://www.esa.int/safe/sentinel-1.0/sentinel-1/sar/level-1"
+                ),
+            }
+
+            # logger.debug(f"Parse manifest info from {metadata_xml_path}")
+            extra["transmitter_receiver_polarisation"] = []
+            for polarisation in manifest_root.findall(
+                "metadataSection/metadataObject/metadataWrap/xmlData/"
+                "s1sarl1:standAloneProductInformation/"
+                "s1sarl1:transmitterReceiverPolarisation",
+                ns,
+            ):
+                extra["transmitter_receiver_polarisation"].append(polarisation.text)
+            extra["productTimelinessCategory"] = manifest_root.find(
+                "metadataSection/metadataObject/metadataWrap/xmlData/"
+                "s1sarl1:standAloneProductInformation/"
+                "s1sarl1:productTimelinessCategory",
+                ns,
+            ).text  # type: ignore
+            extra["instrument_mode"] = manifest_root.find(
+                "metadataSection/metadataObject/metadataWrap/xmlData/"
+                "safe:platform/safe:instrument/safe:extension/"
+                "s1sarl1:instrumentMode/s1sarl1:mode",
+                ns,
+            ).text  # type: ignore
+            extra["orbit_properties_pass"] = manifest_root.find(
+                "metadataSection/metadataObject/metadataWrap/xmlData/"
+                "safe:orbitReference/safe:extension/s1:orbitProperties/s1:pass",
+                ns,
+            ).text  # type: ignore
+            extra["acquisition_date"] = manifest_root.find(
+                "metadataSection/metadataObject/metadataWrap/xmlData/"
+                "safe:acquisitionPeriod/safe:startTime",
+                ns,
+            ).text  # type: ignore
+
+            # Now have a look in the files themselves to get band info,...
+            # TODO: probably cleaner/easier to read from metadata files...
+            # For coherence filename, remove extra .LCO1 extension + date
+            # image_datafilebasename, _ = os.path.splitext(image_basename_noext)
+            # image_datafilebasename = image_datafilebasename[
+            #     :image_datafilebasename.rindex('_')
+            # ]
+
+            metadata = ET.parse(str(image_path / "metadata.xml"))
+            metadata_root = metadata.getroot()
+
+            image_basename_noext_nodate = image_path.stem[
+                : image_path.stem.rfind("_", 0)
+            ]
+            image_datafilename = f"{image_basename_noext_nodate}_byte.tif"
+            image_datapath = image_path / image_datafilename
+
+            # logger.debug(f"Read image metadata from {image_datapath}")
+            with rasterio.open(str(image_datapath)) as src:
+                image_bounds = tuple(src.bounds)
+                image_affine = src.transform
+                if src.crs is not None:
+                    image_crs = src.crs.to_string()
+                    image_epsg = src.crs.to_epsg()
+
+            # Add specific info about the bands
+            bands = {}
+            bands["VH"] = BandInfo(
+                path=str(image_datapath),
+                relative_path=image_datafilename,
+                filename=image_datafilename,
+                bandindex=0,
+            )
+            bands["VV"] = BandInfo(
+                path=str(image_datapath),
+                relative_path=image_datafilename,
+                filename=image_datafilename,
+                bandindex=1,
+            )
+
+        except Exception as ex:
+            raise Exception(
+                f"Exception extracting info from {manifest_xml_path}"
+            ) from ex
+    else:
+        message = (
+            f"Error: found {nb_safefiles} .safe files doing glob in {image_path} "
+            f"with {manifest_xml_searchstring}"
+        )
+        logger.error(message)
+        raise Exception(message)
+
+    return ImageInfo(
+        imagetype=imagetype,
+        filetype=filetype,
+        image_id=image_id,
+        filename=filename,
+        footprint=footprint,
+        image_epsg=image_epsg,
+        image_crs=image_crs,
+        image_bounds=image_bounds,
+        image_affine=image_affine,
+        bands=bands,
+        extra=extra,
+    )
+
+
+def _get_image_info_safe(image_path: Path) -> ImageInfo:
+    # First extract and fill out some basic info
+    imagetype = image_path.stem.split("_")[0]
+    filetype = "SAFE"
+    image_id = image_path.stem
+    extra = {}
+
+    # Read info from the manifest.safe file
+    metadata_xml_path = image_path / "MTD_MSIL2A.xml"
+    try:
+        metadata = ET.parse(str(metadata_xml_path))
+        metadata_root = metadata.getroot()
+
+        # Define namespaces...
+        # xsi:schemaLocation=
+        #     "https://psd-14.sentinel2.eo.esa.int/PSD/User_Product_Level-2A.xsd"
+        ns = {
+            "n1": ("https://psd-14.sentinel2.eo.esa.int/PSD/User_Product_Level-2A.xsd"),
+            "xsi": "http://www.w3.org/2001/XMLSchema-instance",
+        }
+
+        # logger.debug(f"Parse metadata info from {metadata_xml_path}")
+        extra["Cloud_Coverage_Assessment"] = float(
+            metadata_root.find(
+                "n1:Quality_Indicators_Info/Cloud_Coverage_Assessment", ns
+            ).text  # type: ignore
+        )
+        extra["acquisition_date"] = metadata_root.find(
+            "n1:General_Info/Product_Info/PRODUCT_START_TIME", ns
+        ).text  # type: ignore
+
+    except Exception as ex:
+        raise Exception(f"Exception extracting info from {metadata_xml_path}") from ex
+
+    # Now have a look in the files themselves to get band info,...
+    image_datadir = image_path / "GRANULE"
+    band_paths = list(image_datadir.rglob("*.jp2"))
+
+    metadata_path = list(image_datadir.rglob("MTD_TL.xml"))[0]
+    try:
+        # read epsg
+        metadata = ET.parse(str(metadata_path))
+        metadata_root = metadata.getroot()
+
+        ns = {
+            "n1": (
+                "https://psd-14.sentinel2.eo.esa.int/PSD/"
+                "S2_PDI_Level-2A_Tile_Metadata.xsd"
+            ),
+            "xsi": "http://www.w3.org/2001/XMLSchema-instance",
+        }
+
+        image_crs = metadata_root.find(
+            "n1:Geometric_Info/Tile_Geocoding/HORIZONTAL_CS_CODE", ns
+        ).text  # type: ignore
+        image_epsg = int(image_crs.replace("EPSG:", ""))  # type: ignore
+    except Exception as ex:
+        raise Exception(f"Exception extracting info from {metadata_path}") from ex
+
+    # If no files were found, error!
+    if len(band_paths) == 0:
+        message = f"No image files found in {image_datadir}!"
+        logger.error(message)
+        raise Exception(message)
+
+    bands = {}
+    for i, band_path in enumerate(band_paths):
+        band_filename = os.path.basename(band_path)
+        band_filename_noext, _ = os.path.splitext(band_filename)
+        band_filename_noext_split = band_filename_noext.split("_")
+
+        if len(band_filename_noext_split) == 4:
+            # IMG_DATA files
+            band = f"{band_filename_noext_split[2]}-{band_filename_noext_split[3]}"
+        elif len(band_filename_noext_split) == 5:
+            # IMG_DATA files
+            band = f"{band_filename_noext_split[3]}-{band_filename_noext_split[4]}"
+        elif len(band_filename_noext_split) == 3:
+            # QI_DATA files
+            band = band_filename_noext
+        else:
+            message = f"Filename of band doesn't have supported format: {band_path}"
+            logger.error(message)
+            raise Exception(message)
+
+        # Extract bound,... info for band
+        logger.debug(f"Read image metadata from {band_path}")
+        with rasterio.open(str(band_path)) as src:
+            band_bounds = tuple(src.bounds)
+            band_affine = src.transform
+            if src.crs is None:
+                band_crs = image_crs
+                band_epsg = image_epsg
+            else:
+                band_crs = src.crs.to_string()
+                band_epsg = src.crs.to_epsg()
+
+        # Store the crs also on image level + check if all bands have the same crs
+        if i == 0:
+            image_bounds = band_bounds
+        else:
+            if image_crs != band_crs or image_epsg != band_epsg:
+                message = f"Not all bands have the same crs for {image_path}"
+                logger.error(message)
+                raise Exception(message)
+
+        bands[band] = BandInfo(
+            path=str(band_path),
+            relative_path=str(band_path).replace(str(image_path), "")[1:],
+            filename=band_filename,
+            bandindex=0,
+            bounds=band_bounds,
+            affine=band_affine,
+            crs=band_crs,
+            epsg=band_epsg,
+        )
+
+    return ImageInfo(
+        imagetype=imagetype,
+        filetype=filetype,
+        image_id=image_id,
+        filename=image_path.name,
+        footprint=None,
+        image_epsg=image_epsg,
+        image_crs=image_crs,
+        image_bounds=image_bounds,
+        image_affine=None,
+        bands=bands,
+        extra=extra,
+    )
+
+
+def _get_image_info_tif(image_path: Path) -> ImageInfo:
+    with rasterio.open(image_path) as src:
+        image_bounds = tuple(src.bounds)
+        image_affine = src.transform
+        if src.crs is not None:
+            image_crs = src.crs.to_string()
+            image_epsg = src.crs.to_epsg()
+
+        # Add all bands
+        bands = {}
+        for band_idx in src.indexes:
+            band_tags = src.tags(band_idx)
+            band_name = band_tags.get("DESCRIPTION", None)
+            band_name = band_name if band_name is not None else str(band_idx)
+            bands[band_name] = BandInfo(
+                path=str(image_path),
+                relative_path=None,
+                filename=image_path.name,
+                bandindex=band_idx - 1,
+            )
+
+    footprint = {}
+    footprint["shape"] = None
+    footprint["crs"] = None
+    footprint["epsg"] = None
+
+    extra = {}
+    imagetype = image_path.stem.split("_")[0].upper()
+    if imagetype in ("S1-ASC", "S1-DESC"):
+        orbit = imagetype.split("-")[1]
+        if orbit == "ASC":
+            extra["orbit_properties_pass"] = "ASCENDING"
+        elif orbit == "DESC":
+            extra["orbit_properties_pass"] = "DESCENDING"
+    return ImageInfo(
+        imagetype=imagetype,
+        filetype="TIF",
+        image_id=image_path.stem,
+        filename=image_path.name,
+        footprint=footprint,
+        image_epsg=image_epsg,
+        image_crs=image_crs,
+        image_bounds=image_bounds,
+        image_affine=image_affine,
+        bands=bands,
+        extra=extra,
+    )
diff --git a/environment-dev.yml b/environment-dev.yml
index cc67d5ea..6386d733 100644
--- a/environment-dev.yml
+++ b/environment-dev.yml
@@ -10,7 +10,9 @@ dependencies:
   - geopandas>=0.11
   - h5py <3
   - numpy
-  - openeo
+  # Openeo >=0.14 +- doesn't work on windows.
+  # More info here: https://github.com/Open-EO/openeo-python-client/issues/387
+  - openeo <0.14
   - pandas
   - psutil
   - pyproj
diff --git a/openeo-client-config.ini b/openeo-client-config.ini
new file mode 100644
index 00000000..d05fcdd0
--- /dev/null
+++ b/openeo-client-config.ini
@@ -0,0 +1,3 @@
+[Connection]
+default_backend = https://openeo.vito.be
+default_backend.auto_authenticate = oidc
\ No newline at end of file
diff --git a/setup.py b/setup.py
index 8e00116c..943b4340 100644
--- a/setup.py
+++ b/setup.py
@@ -20,6 +20,7 @@
     install_requires=[
         "geofileops",
         "geopandas",
+        "openeo",
         "psutil",
         "rasterio",
         "rasterstats",


Test patch:
diff --git a/tests/test_calc_timeseries.py b/tests/test_calc_timeseries.py
index e51a8b1c..689725dc 100644
--- a/tests/test_calc_timeseries.py
+++ b/tests/test_calc_timeseries.py
@@ -2,7 +2,7 @@
 import logging
 from pathlib import Path
 
-from cropclassification.preprocess import timeseries_calc_dias_onda_per_image as calc_ts
+from cropclassification.util import zonal_stats_bulk
 
 
 def get_testdata_dir() -> Path:
@@ -25,11 +25,11 @@ def test_calc_stats_per_image_s1_bs(tmpdir):
     tmp_dir = Path(tmpdir) / "raw"
     try:
         start_time = datetime.now()
-        calc_ts.calc_stats_per_image(
+        images_bands = [(path, ["VV", "VH"]) for path in input_image_paths]
+        zonal_stats_bulk.zonal_stats(
             features_path=input_features_path,
             id_column="UID",
-            image_paths=input_image_paths,
-            bands=["VV", "VH"],
+            images_bands=images_bands,
             output_dir=tmp_dir,
             temp_dir=tmp_dir / "tmp",
             log_dir=tmp_dir / "log",
@@ -54,11 +54,11 @@ def test_calc_stats_per_image_s1_bs(tmpdir):
     tmp_dir = Path(tmpdir) / "tif"
     try:
         start_time = datetime.now()
-        calc_ts.calc_stats_per_image(
+        images_bands = [(path, ["VV", "VH"]) for path in input_image_paths]
+        zonal_stats_bulk.zonal_stats(
             features_path=input_features_path,
             id_column="UID",
-            image_paths=input_image_paths,
-            bands=["VV", "VH"],
+            images_bands=images_bands,
             output_dir=tmp_dir,
             temp_dir=tmp_dir / "tmp",
             log_dir=tmp_dir / "log",