Repository: theroggy/cropclassification (language: python)

Issue:
ENH: Improve performance of zonal_stats_bulk
Calculating the zonal stats now takes a long time. It would be great if that could be made faster.

Steps to take:
  - [x] make list of existing relevant libraries: see below
  - [x] benchmark them: [benchmarks](
https://github.com/geofileops/geobenchmark/blob/main/results_zonalstats/GeoBenchmark.png): based on the benchmarks, the pyqgis zonal_stats was the (only) one with with significantly better performance.
  - [x] first rough implementation: #42
  - [ ] tests still need to be added 
  - [ ] decouple implementation further for rasterstats one + make implementation more generic

### Listing of relevant libraries

- rasterio: [manual](https://pysal.org/scipy2019-intermediate-gds/deterministic/gds2-rasters.html)
- pygeoprocessing (open?): [doc](https://pygeoprocessing.readthedocs.io/en/latest/api/pygeoprocessing.html#pygeoprocessing.zonal_statistics)
- pyjeo (GPL): [doc](https://pyjeo.readthedocs.io/en/latest/2_tutorial.html#tutorial-on-extract-calculating-regional-statistics)
- pyQGIS (GPL): [post](https://gis.stackexchange.com/questions/421556/performance-problem-with-getting-average-pixel-values-within-buffered-circles)
    - fast, but less clean on server.
    - [relevant src](https://github.com/qgis/QGIS/blob/d5626d92360efffb4b8085389c8d64072ef65833/src/analysis/vector/qgszonalstatistics.cpp#L266)
- ?: [doc](https://pygis.io/docs/f_rs_extraction.html)

These libraries also take partial pixels in account, so they are not fit for purpose.
- xagg: [repo](https://github.com/ks905383/xagg)
- exactextract (apache2): [repo](https://github.com/isciences/exactextract)



Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 27536908..fb309015 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,7 @@
 
 ### Improvements
 
+- Improve performance of zonal_stats_bulk (#38)
 - Use black to comply to pep8 + minor general improvements (#13)
 - Upgrade all dependencies (#12)
 - Add support for pandas 2.0 (#41)
diff --git a/README.md b/README.md
index 2c2f7cd6..9e89b468 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" openeo psutil rasterio rasterstats scikit-learn
+      conda install python=3.9 geopandas geofileops "h5py<3" openeo psutil qgis 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/cropclassification/calc_timeseries.py b/cropclassification/calc_timeseries.py
index 40376ebe..123afbcb 100644
--- a/cropclassification/calc_timeseries.py
+++ b/cropclassification/calc_timeseries.py
@@ -175,9 +175,9 @@ def calc_timeseries_task(config_paths: List[Path], default_basedir: Path):
     try:
         images_bands = [(path, ["VV", "VH"]) for path in input_image_paths]
         zonal_stats_bulk.zonal_stats(
-            features_path=input_features_path,
+            vector_path=input_features_path,
             id_column=conf.columns["id"],
-            images_bands=images_bands,
+            rasters_bands=images_bands,
             output_dir=output_dir,
             temp_dir=temp_dir,
             log_dir=log_dir,
@@ -248,9 +248,9 @@ def calc_timeseries_task(config_paths: List[Path], default_basedir: Path):
         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,
+            vector_path=input_features_path,
             id_column=conf.columns["id"],
-            images_bands=images_bands,
+            rasters_bands=images_bands,
             output_dir=output_dir,
             temp_dir=temp_dir,
             log_dir=log_dir,
@@ -292,9 +292,9 @@ def calc_timeseries_task(config_paths: List[Path], default_basedir: Path):
     try:
         images_bands = [(path, ["VV", "VH"]) for path in input_image_paths]
         zonal_stats_bulk.zonal_stats(
-            features_path=input_features_path,
+            vector_path=input_features_path,
             id_column=conf.columns["id"],
-            images_bands=images_bands,
+            rasters_bands=images_bands,
             output_dir=output_dir,
             temp_dir=temp_dir,
             log_dir=log_dir,
diff --git a/cropclassification/predict/classification_keras.py b/cropclassification/predict/classification_keras.py
index 7fb28203..ff13f34f 100644
--- a/cropclassification/predict/classification_keras.py
+++ b/cropclassification/predict/classification_keras.py
@@ -335,7 +335,9 @@ def predict_proba(
         logger.info(f"Resulting Columns for predicting data: {parcel_data_df.columns}")
 
     # Load the classifier and predict
-    model = tf.keras.models.load_model(classifier_path)
+    model = tf.keras.models.load_model(str(classifier_path))
+    if model is None:
+        raise RuntimeError(f"Error loading model {classifier_path}")
     logger.info(f"Predict classes with probabilities: {len(parcel_df.index)} rows")
     class_proba = model.predict(parcel_data_df)
     logger.info("Predict classes with probabilities ready")
diff --git a/cropclassification/preprocess/_timeseries_calc_openeo.py b/cropclassification/preprocess/_timeseries_calc_openeo.py
index b54d97c0..7eb5e172 100644
--- a/cropclassification/preprocess/_timeseries_calc_openeo.py
+++ b/cropclassification/preprocess/_timeseries_calc_openeo.py
@@ -55,13 +55,12 @@ def calculate_periodic_timeseries(
     if temp_dir == "None":
         temp_dir = Path(tempfile.gettempdir())
     zonal_stats_bulk.zonal_stats(
-        features_path=input_parcel_path,
+        vector_path=input_parcel_path,
         id_column=conf.columns["id"],
-        images_bands=images_bands,
+        rasters_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"),
+        stats=["count", "mean", "median", "std", "min", "max"],
+        engine="pyqgis",
     )
 
     """
diff --git a/cropclassification/util/zonal_stats_bulk/__init__.py b/cropclassification/util/zonal_stats_bulk/__init__.py
index 15089570..8e97d454 100644
--- a/cropclassification/util/zonal_stats_bulk/__init__.py
+++ b/cropclassification/util/zonal_stats_bulk/__init__.py
@@ -1,17 +1,15 @@
 from pathlib import Path
-from typing import List, Optional, Tuple, Union
+from typing import List, Literal, Optional, Tuple
 
 from ._raster_helper import *  # noqa: F401, F403
 
 
 def zonal_stats(
-    features_path: Path,
+    vector_path: Path,
     id_column: str,
-    images_bands: List[Tuple[Path, List[str]]],
+    rasters_bands: List[Tuple[Path, List[str]]],
     output_dir: Path,
-    temp_dir: Path,
-    log_dir: Path,
-    log_level: Union[str, int],
+    stats: Literal["count", "mean", "median", "std", "min", "max"],
     cloud_filter_band: Optional[str] = None,
     calc_bands_parallel: bool = True,
     engine: str = "rasterstats",
@@ -21,44 +19,39 @@ def zonal_stats(
     Calculate zonal statistics.
 
     Args:
-        features_path (Path): _description_
+        vector_path (Path): _description_
         id_column (str): _description_
-        images_bands (List[Tuple[Path, List[str]]]): _description_
+        rasters_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_
     """
     if engine == "pyqgis":
-        from _zonal_stats_bulk_pyqgis import zonal_stats
+        if cloud_filter_band is not None:
+            raise ValueError(
+                'cloud_filter_band parameter not supported for engine "pyqgis"'
+            )
+        from . import _zonal_stats_bulk_pyqgis
 
-        return zonal_stats(
-            features_path=features_path,
-            id_column=id_column,
-            images_bands=images_bands,
+        return _zonal_stats_bulk_pyqgis.zonal_stats(
+            vector_path=vector_path,
+            rasters_bands=rasters_bands,
             output_dir=output_dir,
-            temp_dir=temp_dir,
-            log_dir=log_dir,
-            log_level=log_level,
-            cloud_filter_band=cloud_filter_band,
-            calc_bands_parallel=calc_bands_parallel,
+            stats=stats,
+            columns=[id_column],
             force=force,
         )
     elif engine == "rasterstats":
-        from _zonal_stats_bulk_rs import zonal_stats
+        from . import _zonal_stats_bulk_rs
 
-        return zonal_stats(
-            features_path=features_path,
+        return _zonal_stats_bulk_rs.zonal_stats(
+            vector_path=vector_path,
             id_column=id_column,
-            images_bands=images_bands,
+            rasters_bands=rasters_bands,
             output_dir=output_dir,
-            temp_dir=temp_dir,
-            log_dir=log_dir,
-            log_level=log_level,
+            stats=stats,
             cloud_filter_band=cloud_filter_band,
             calc_bands_parallel=calc_bands_parallel,
             force=force,
diff --git a/cropclassification/util/zonal_stats_bulk/_general_helper.py b/cropclassification/util/zonal_stats_bulk/_general_helper.py
index 3cc960ee..cc49839c 100644
--- a/cropclassification/util/zonal_stats_bulk/_general_helper.py
+++ b/cropclassification/util/zonal_stats_bulk/_general_helper.py
@@ -1,7 +1,7 @@
 from datetime import datetime
 import logging
 from pathlib import Path
-from typing import Optional
+from typing import Optional, Union
 
 logger = logging.getLogger(__name__)
 
@@ -11,7 +11,7 @@ def _format_output_path(
     image_path: Path,
     output_dir: Path,
     orbit_properties_pass: Optional[str],
-    band: Optional[str],
+    band: Optional[Union[str, int]],
 ) -> Path:
     """
     Prepare the output path.
@@ -40,9 +40,9 @@ def _format_output_path(
 def _format_progress_message(
     nb_todo: int,
     nb_done_total: int,
-    nb_done_latestbatch: int,
     start_time: datetime,
-    start_time_latestbatch: datetime,
+    nb_done_latestbatch: Optional[int] = None,
+    start_time_latestbatch: Optional[datetime] = None,
 ) -> str:
     """
     Returns a progress message based on the input.
@@ -55,9 +55,11 @@ def _format_progress_message(
         start_time_latestbatch: datetime the latest batch started
     """
     time_passed_s = (datetime.now() - start_time).total_seconds()
-    time_passed_latestbatch_s = (
-        datetime.now() - start_time_latestbatch
-    ).total_seconds()
+    time_passed_latestbatch_s = None
+    if start_time_latestbatch is not None:
+        time_passed_latestbatch_s = (
+            datetime.now() - start_time_latestbatch
+        ).total_seconds()
 
     # Calculate the overall progress
     large_number = 9999999999
@@ -68,17 +70,20 @@ def _format_progress_message(
     hours_to_go = (int)((nb_todo - nb_done_total) / nb_per_hour)
     min_to_go = (int)((((nb_todo - nb_done_total) / nb_per_hour) % 1) * 60)
 
-    # Calculate the speed of the latest batch
-    if time_passed_latestbatch_s > 0:
+    # Format message
+    message = (
+        f"{hours_to_go}:{min_to_go} left for {nb_todo-nb_done_total} todo at "
+        f"{nb_per_hour:0.0f}/h"
+    )
+    # Add speed of the latest batch to message if appropriate
+    if (
+        time_passed_latestbatch_s is not None
+        and nb_done_latestbatch is not None
+        and time_passed_latestbatch_s > 0
+    ):
         nb_per_hour_latestbatch = (
             nb_done_latestbatch / time_passed_latestbatch_s
         ) * 3600
-    else:
-        nb_per_hour_latestbatch = large_number
+        message = f"{message} ({nb_per_hour_latestbatch:0.0f}/h last batch)"
 
-    # Return formatted message
-    message = (
-        f"{hours_to_go}:{min_to_go} left for {nb_todo-nb_done_total} todo at "
-        f"{nb_per_hour:0.0f}/h ({nb_per_hour_latestbatch:0.0f}/h last batch)"
-    )
     return message
diff --git a/cropclassification/util/zonal_stats_bulk/_vector_helper.py b/cropclassification/util/zonal_stats_bulk/_vector_helper.py
index ee2f9ce3..ce2c093d 100644
--- a/cropclassification/util/zonal_stats_bulk/_vector_helper.py
+++ b/cropclassification/util/zonal_stats_bulk/_vector_helper.py
@@ -2,7 +2,7 @@
 import logging
 from pathlib import Path
 import time
-from typing import List
+from typing import List, Optional
 
 import geofileops as gfo
 import geopandas as gpd
@@ -12,6 +12,121 @@
 logger = logging.getLogger(__name__)
 
 
+def reproject_synced(
+    path: Path,
+    columns: List[str],
+    target_epsg: int,
+    dst_dir: Optional[Path] = None,
+) -> Path:
+    """
+    Reproject the input file. It is process locked and if another process is already
+    busy, wait till it is ready.
+
+    Args:
+        features_path (Path): _description_
+        columns_to_retain (List[str]): _description_
+        target_epsg (int): _description_
+
+    Raises:
+        Exception: _description_
+        Exception: _description_
+
+    Returns:
+        Path: _description_
+    """
+    vector_info = gfo.get_layerinfo(path)
+    assert vector_info.crs is not None
+    vector_epsg = vector_info.crs.to_epsg()
+
+    if dst_dir is None:
+        dst_dir = path.parent
+
+    # Determine the correct filename for the input features in the correct projection.
+    if vector_epsg != target_epsg:
+        vector_prepr_path = dst_dir / f"{path.stem}_{target_epsg:.0f}.gpkg"
+    else:
+        vector_prepr_path = path
+
+    # Prepare filename for a "busy file" to ensure proper behaviour in a parallel
+    # processing context
+    vector_prepr_path_busy = Path(f"{str(vector_prepr_path)}_busy")
+
+    # If the file exists already or if a busy file exists, return
+    vector_gdf = None
+    if vector_prepr_path.exists():
+        return vector_prepr_path
+
+    # 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(vector_prepr_path_busy):
+        try:
+            # Read (all) original features + remove unnecessary columns...
+            logger.info(f"Read original file {path}")
+            start_time = datetime.now()
+            logging.getLogger("fiona.ogrext").setLevel(logging.INFO)
+            vector_gdf = gfo.read_file(path)
+            logger.info(
+                f"Read ready, found {len(vector_gdf.index)} features, "
+                f"crs: {vector_gdf.crs}, took "
+                f"{(datetime.now()-start_time).total_seconds()} s"
+            )
+            for column in vector_gdf.columns:
+                if column not in columns and column not in [
+                    "geometry",
+                    "x_ref",
+                ]:
+                    vector_gdf.drop(columns=column, inplace=True)
+
+            # Reproject them
+            logger.info(
+                f"Reproject features from {vector_gdf.crs} to epsg {target_epsg}"
+            )
+            vector_gdf = vector_gdf.to_crs(epsg=target_epsg)
+            logger.info("Reprojected, now sort on x_ref")
+
+            if vector_gdf is None:
+                raise Exception("features_gdf is None")
+            # Order features on x coordinate
+            if "x_ref" not in vector_gdf.columns:
+                vector_gdf["x_ref"] = vector_gdf.geometry.bounds.minx
+            vector_gdf.sort_values(by=["x_ref"], inplace=True)
+            vector_gdf.reset_index(inplace=True)
+
+            # Cache the file for future use
+            logger.info(
+                f"Write {len(vector_gdf.index)} reprojected features to "
+                f"{vector_prepr_path}"
+            )
+            gfo.to_file(vector_gdf, vector_prepr_path, index=False)  # type: ignore
+            logger.info("Reprojected features written")
+
+        except Exception as ex:
+            # If an exception occurs...
+            message = f"Delete possibly incomplete file: {vector_prepr_path}"
+            logger.exception(message)
+            gfo.remove(vector_prepr_path)
+            raise Exception(message) from ex
+        finally:
+            # Remove lock file as everything is ready for other processes to use it
+            vector_prepr_path_busy.unlink()
+
+    # If a "busy file" exists, the file isn't ready yet, but another process
+    # is working on it, so wait till it disappears
+    wait_secs_max = 600
+    wait_start_time = datetime.now()
+    while vector_prepr_path_busy.exists():
+        time.sleep(1)
+        wait_secs = (datetime.now() - wait_start_time).total_seconds()
+        if wait_secs > wait_secs_max:
+            raise Exception(
+                f"Waited {wait_secs} for busy file "
+                f"{vector_prepr_path_busy} and it is still there!"
+            )
+
+    return vector_prepr_path
+
+
 def _load_features_file(
     features_path: Path,
     columns_to_retain: List[str],
@@ -33,142 +148,31 @@ def _load_features_file(
         bbox: bounds of the area to be loaded, in the target_epsg
     """
     # Load parcel file and preprocess it: remove excess columns + reproject if needed.
-    # 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_info = gfo.get_layerinfo(features_path)
-    assert features_info.crs is not None
-    features_epsg = features_info.crs.to_epsg()
-
-    # Determine the correct filename for the input features in the correct projection.
-    if features_epsg != target_epsg:
-        features_prepr_path = (
-            features_path.parent / f"{features_path.stem}_{target_epsg:.0f}.gpkg"
-        )
-    else:
-        features_prepr_path = features_path
-
-    # Prepare filename for a "busy file" to ensure proper behaviour in a parallel
-    # processing context
-    features_prepr_path_busy = Path(f"{str(features_prepr_path)}_busy")
-
-    # If the file doesn't exist yet in right projection, read original input file to
-    # reproject/write to new file with correct epsg
-    features_gdf = None
-    if not (features_prepr_path_busy.exists() or features_prepr_path.exists()):
-        # 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(features_prepr_path_busy):
-            try:
-                # Read (all) original features + remove unnecessary columns...
-                logger.info(f"Read original file {features_path}")
-                start_time = datetime.now()
-                logging.getLogger("fiona.ogrext").setLevel(logging.INFO)
-                features_gdf = gfo.read_file(features_path)
-                logger.info(
-                    f"Read ready, found {len(features_gdf.index)} features, "
-                    f"crs: {features_gdf.crs}, took "
-                    f"{(datetime.now()-start_time).total_seconds()} s"
-                )
-                for column in features_gdf.columns:
-                    if column not in columns_to_retain and column not in [
-                        "geometry",
-                        "x_ref",
-                    ]:
-                        features_gdf.drop(columns=column, inplace=True)
-
-                # Reproject them
-                logger.info(
-                    f"Reproject features from {features_gdf.crs} to epsg {target_epsg}"
-                )
-                features_gdf = features_gdf.to_crs(epsg=target_epsg)
-                logger.info("Reprojected, now sort on x_ref")
-
-                if features_gdf is None:
-                    raise Exception("features_gdf is None")
-                # Order features on x coordinate
-                if "x_ref" not in features_gdf.columns:
-                    features_gdf["x_ref"] = features_gdf.geometry.bounds.minx
-                features_gdf.sort_values(by=["x_ref"], inplace=True)
-                features_gdf.reset_index(inplace=True)
-
-                # Cache the file for future use
-                logger.info(
-                    f"Write {len(features_gdf.index)} reprojected features to "
-                    f"{features_prepr_path}"
-                )
-                gfo.to_file(
-                    features_gdf, features_prepr_path, index=False  # type: ignore
-                )
-                logger.info("Reprojected features written")
-
-            except Exception as ex:
-                # If an exception occurs...
-                message = f"Delete possibly incomplete file: {features_prepr_path}"
-                logger.exception(message)
-                gfo.remove(features_prepr_path)
-                raise Exception(message) from ex
-            finally:
-                # Remove lock file as everything is ready for other processes to use it
-                features_prepr_path_busy.unlink()
-
-            # Now filter the parcels that are in bbox provided
-            if bbox is not None:
-                logger.info(f"bbox provided, so filter features in the bbox of {bbox}")
-                xmin, ymin, xmax, ymax = bbox
-                features_gdf = features_gdf.cx[xmin:xmax, ymin:ymax]
-                logger.info(f"Found {len(features_gdf.index)} features in bbox")
+    features_prepr_path = reproject_synced(
+        path=features_path, columns=columns_to_retain, target_epsg=target_epsg
+    )
 
     # If there exists already a file with the features in the right projection, we can
     # just read the data
-    if features_gdf is None:
-        # If a "busy file" still exists, the file isn't ready yet, but another process
-        # is working on it, so wait till it disappears
-        wait_secs_max = 600
-        wait_start_time = datetime.now()
-        while features_prepr_path_busy.exists():
-            time.sleep(1)
-            wait_secs = (datetime.now() - wait_start_time).total_seconds()
-            if wait_secs > wait_secs_max:
-                raise Exception(
-                    f"Waited {wait_secs} for busy file "
-                    f"{features_prepr_path_busy} and it is still there!"
-                )
-
-        logger.info(f"Read {features_prepr_path}")
-        start_time = datetime.now()
-        features_gdf = gfo.read_file(features_prepr_path, bbox=bbox)
-        logger.info(
-            f"Read ready, found {len(features_gdf.index)} features, crs: "
-            f"{features_gdf.crs}, took {(datetime.now()-start_time).total_seconds()} s"
-        )
+    logger.info(f"Read {features_prepr_path}")
+    start_time = datetime.now()
+    features_gdf = gfo.read_file(features_prepr_path, bbox=bbox)
+    logger.info(
+        f"Read ready, found {len(features_gdf.index)} features, crs: "
+        f"{features_gdf.crs}, took {(datetime.now()-start_time).total_seconds()} s"
+    )
 
-        # Order features on x_ref to (probably) have more clustering of features in
-        # further action...
-        if "x_ref" not in features_gdf.columns:
-            features_gdf["x_ref"] = features_gdf.geometry.bounds.minx
-        features_gdf.sort_values(by=["x_ref"], inplace=True)
-        features_gdf.reset_index(inplace=True)
-
-        # To be sure, remove the columns anyway...
-        for column in features_gdf.columns:
-            if column not in columns_to_retain and column not in ["geometry"]:
-                features_gdf.drop(columns=column, inplace=True)
+    # Order features on x_ref to (probably) have more clustering of features in
+    # further action...
+    if "x_ref" not in features_gdf.columns:
+        features_gdf["x_ref"] = features_gdf.geometry.bounds.minx
+    features_gdf.sort_values(by=["x_ref"], inplace=True)
+    features_gdf.reset_index(inplace=True)
+
+    # To be sure, remove the columns anyway...
+    for column in features_gdf.columns:
+        if column not in columns_to_retain and column not in ["geometry"]:
+            features_gdf.drop(columns=column, inplace=True)
 
     # If there is a polygon provided, filter on the polygon (as well)
     if polygon is not None:
diff --git a/cropclassification/util/zonal_stats_bulk/_zonal_stats_bulk_pyqgis.py b/cropclassification/util/zonal_stats_bulk/_zonal_stats_bulk_pyqgis.py
index 8b2e0e20..083299de 100644
--- a/cropclassification/util/zonal_stats_bulk/_zonal_stats_bulk_pyqgis.py
+++ b/cropclassification/util/zonal_stats_bulk/_zonal_stats_bulk_pyqgis.py
@@ -2,27 +2,53 @@
 from datetime import datetime
 import logging
 import multiprocessing
+import os
 from pathlib import Path
-from typing import List, Optional, Tuple, Union
+import sys
+import tempfile
+from typing import List, Literal, Tuple, Union
 
+# Set path for qgis
+qgis_path = Path(os.environ["CONDA_PREFIX"]) / "Library/python"
+sys.path.insert(0, str(qgis_path))
+
+import geofileops as gfo
+import geopandas as gpd
+import pandas as pd
 import qgis.core  # type: ignore
 import qgis.analysis  # type: ignore
 
 from cropclassification.util import io_util
+from cropclassification.helpers import pandas_helper as pdh
+from . import _general_helper as general_helper
+from . import _raster_helper as raster_helper
+from . import _vector_helper as vector_helper
 
 logger = logging.getLogger(__name__)
 
 
+Statistic = Literal[
+    "count",
+    "sum",
+    "mean",
+    "median",
+    "std",
+    "min",
+    "max",
+    "range",
+    "minority",
+    "majority",
+    "variance",
+]
+DEFAULT_STATS = ["count", "median"]
+
+
 def zonal_stats(
-    features_path: Path,
-    id_column: str,
-    images_bands: List[Tuple[Path, List[str]]],
+    vector_path: Path,
+    columns: List[str],
+    rasters_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,
+    stats: List[Statistic] = DEFAULT_STATS,
     force: bool = False,
 ):
     """
@@ -32,13 +58,257 @@ def zonal_stats(
         features_path (Path): _description_
         id_column (str): _description_
         images_bands (List[Tuple[Path, List[str]]]): _description_
+        stats (List[Statistic])
         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_
     """
-    raise Exception("Not implemented")
+    # Some checks on the input parameters
+    if len(rasters_bands) == 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 = io_util.create_tempdir("zonal_stats_pyqgis")
+
+    start_time = datetime.now()
+    nb_todo = 0
+    nb_done_total = 0
+    nb_parallel_max = multiprocessing.cpu_count()
+    nb_parallel = nb_parallel_max
+
+    # Loop over all images and bands to calculate zonal stats in parallel...
+    calc_queue = {}
+    pool = futures.ProcessPoolExecutor(nb_parallel)
+    try:
+        for raster_path, bands in rasters_bands:
+            raster_info = raster_helper.get_image_info(raster_path)
+            # Create base output filename
+            # TODO: hoort hier niet echt thuis
+            orbit = None
+            if raster_info.imagetype.lower() == "s1-grd-sigma0-asc":
+                orbit = "ASCENDING"
+            elif raster_info.imagetype.lower() == "s1-grd-sigma0-desc":
+                orbit = "DESCENDING"
+            output_base_path = general_helper._format_output_path(
+                vector_path,
+                raster_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 = raster_info.bands
+                assert bands is not None
+            bands_todo = {}
+
+            for band in bands:
+                # Prepare the output paths...
+                output_band_path = general_helper._format_output_path(
+                    vector_path,
+                    raster_path,
+                    output_dir,
+                    orbit,
+                    band,
+                )
+                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}"
+                    )
+                    continue
+
+                nb_todo += 1
+                future = pool.submit(
+                    zonal_stats_band_tofile,
+                    vector_path=vector_path,
+                    raster_path=raster_path,
+                    band=band,
+                    stats=stats,
+                    columns=columns,
+                    output_band_path=output_band_path,
+                )
+                calc_queue[future] = {
+                    "vector_path": vector_path,
+                    "raster_path": raster_path,
+                    "bands": bands,
+                    "bands_todo": bands_todo,
+                    "prepare_calc_future": future,
+                    "image_info": raster_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",
+                }
+
+        # Write progress for each completed calculation
+        for future in futures.as_completed(calc_queue):
+            try:
+                _ = future.result()
+            except Exception as ex:
+                raise Exception(f"Error calculating {calc_queue[future]}: {ex}") from ex
+            nb_done_total += 1
+            progress_msg = general_helper._format_progress_message(
+                nb_todo,
+                nb_done_total,
+                start_time,
+            )
+            logger.info(progress_msg)
+    finally:
+        pool.shutdown()
+
+
+def zonal_stats_band(
+    vector_path,
+    raster_path: Path,
+    band: str,
+    stats: List[Statistic],
+    columns: List[str] = [],
+) -> Union[pd.DataFrame, gpd.GeoDataFrame]:
+    # Init
+    stats_mask = 0
+    for stat in stats:
+        stats_mask |= stat_to_qgisstat(stat)
+
+    # Get the image info
+    image_info = raster_helper.get_image_info(raster_path)
+
+    # Reproject the vector data
+    tmp_dir = Path(tempfile.gettempdir()) / "zonal_stats_bulk"
+    tmp_dir.mkdir(exist_ok=True, parents=True)
+    vector_proj_path = vector_helper.reproject_synced(
+        path=vector_path,
+        columns=columns + ["geometry"],
+        target_epsg=image_info.image_epsg,
+        dst_dir=tmp_dir,
+    )
+    layer = gfo.get_only_layer(vector_proj_path)
+
+    # Init qgis
+    qgs = qgis.core.QgsApplication([], False)
+    qgs.initQgis()
+    # Read the vector file + copy to memory layer for:
+    #   - improved performance
+    #   - QgsZonalStatistics actually adds the data to the input file, so copy needed
+    vector = qgis.core.QgsVectorLayer(str(vector_proj_path), layer, "ogr")
+    vector_mem = vector.materialize(
+        qgis.core.QgsFeatureRequest().setFilterFids(vector.allFeatureIds())
+    )
+    del vector
+
+    # Calculates zonal stats with raster
+    raster = qgis.core.QgsRasterLayer(image_info.bands[band].path)
+
+    zoneStats = qgis.analysis.QgsZonalStatistics(
+        vector_mem,
+        raster,
+        stats=stats_mask,
+        rasterBand=image_info.bands[band].bandindex,
+    )
+    if zoneStats.calculateStatistics(None) != 0:
+        raise RuntimeError(
+            "Error: calculateStatistics didn't return 0 for zonal stats between "
+            f"{vector_proj_path.name} and {raster_path.name}"
+        )
+
+    # Convert result to (geo)dataframe
+    if "geometry" in columns:
+        columns = [f.name() for f in vector_mem.fields()] + ["geometry"]
+        data = [
+            dict(zip(columns, f.attributes() + [f.geometry().asWkt()]))
+            for f in vector_mem.getFeatures()
+        ]
+    else:
+        columns = [f.name() for f in vector_mem.fields()]
+        data = [dict(zip(columns, f.attributes())) for f in vector_mem.getFeatures()]
+    stats_df = pd.DataFrame(data, columns=columns)
+
+    if "geometry" in columns:
+        stats_df["geometry"] = gpd.GeoSeries.from_wkt(stats_df["geometry"])
+        stats_df = gpd.GeoDataFrame(
+            stats_df, geometry="geometry", crs=vector_mem.crs().toWkt()
+        )  # type: ignore
+
+    return stats_df
+
+
+def zonal_stats_band_tofile(
+    vector_path,
+    raster_path: Path,
+    output_band_path: Path,
+    band: str,
+    stats: List[Statistic],
+    columns: List[str],
+    force: bool = False,
+) -> Path:
+    if output_band_path.exists():
+        if force:
+            output_band_path.unlink()
+        else:
+            return output_band_path
+
+    stats_df = zonal_stats_band(
+        vector_path=vector_path,
+        raster_path=raster_path,
+        band=band,
+        stats=stats,
+        columns=columns,
+    )
+    # Remove rows with empty data
+    # stats_df.dropna(inplace=True)
+
+    logger.info(
+        f"Write data for {len(stats_df.index)} parcels found to " f"{output_band_path}"
+    )
+    pdh.to_file(stats_df, output_band_path, index=False)
+    return output_band_path
+
+
+def stat_to_qgisstat(stat: str):
+    if stat == "count":
+        return qgis.analysis.QgsZonalStatistics.Count
+    elif stat == "sum":
+        return qgis.analysis.QgsZonalStatistics.Sum
+    elif stat == "mean":
+        return qgis.analysis.QgsZonalStatistics.Mean
+    elif stat == "median":
+        return qgis.analysis.QgsZonalStatistics.Median
+    elif stat == "std":
+        return qgis.analysis.QgsZonalStatistics.StDev
+    elif stat == "min":
+        return qgis.analysis.QgsZonalStatistics.Min
+    elif stat == "max":
+        return qgis.analysis.QgsZonalStatistics.Max
+    elif stat == "range":
+        return qgis.analysis.QgsZonalStatistics.Range
+    elif stat == "minority":
+        return qgis.analysis.QgsZonalStatistics.Minority
+    elif stat == "majority":
+        return qgis.analysis.QgsZonalStatistics.Majority
+    elif stat == "variance":
+        return qgis.analysis.QgsZonalStatistics.Variance
+    else:
+        raise ValueError(f"unsupported value in stats: {stat}")
diff --git a/cropclassification/util/zonal_stats_bulk/_zonal_stats_bulk_rs.py b/cropclassification/util/zonal_stats_bulk/_zonal_stats_bulk_rs.py
index 55680eb6..25b730f8 100644
--- a/cropclassification/util/zonal_stats_bulk/_zonal_stats_bulk_rs.py
+++ b/cropclassification/util/zonal_stats_bulk/_zonal_stats_bulk_rs.py
@@ -14,7 +14,7 @@
 import signal  # To catch CTRL-C explicitly and kill children
 import sys
 import time
-from typing import List, Optional, Tuple, Union
+from typing import List, Literal, Optional, Tuple, Union
 
 from osgeo import gdal
 
@@ -37,13 +37,11 @@
 
 
 def zonal_stats(
-    features_path: Path,
+    vector_path: Path,
     id_column: str,
-    images_bands: List[Tuple[Path, List[str]]],
+    rasters_bands: List[Tuple[Path, List[str]]],
     output_dir: Path,
-    temp_dir: Path,
-    log_dir: Path,
-    log_level: Union[str, int],
+    stats: Literal["count", "mean", "median", "std", "min", "max"],
     cloud_filter_band: Optional[str] = None,
     calc_bands_parallel: bool = True,
     force: bool = False,
@@ -72,15 +70,16 @@ def zonal_stats(
     # after calculation together with the processing logic
 
     # Some checks on the input parameters
-    nb_todo = len(images_bands)
+    nb_todo = len(rasters_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)
-    log_dir.mkdir(parents=True, exist_ok=True)
-    temp_dir = io_util.create_tempdir("calc_stats_per_image", parent_dir=temp_dir)
+    temp_dir = io_util.create_tempdir("zonal_stats_rs")
+    log_dir = temp_dir
+    log_level = logger.level
 
     # Create process pool for parallelisation...
     nb_parallel_max = multiprocessing.cpu_count()
@@ -106,13 +105,13 @@ def zonal_stats(
             # 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)
+                image_idx < len(rasters_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, bands = rasters_bands[image_idx]
                 image_path_str = str(image_path)
                 image_idx += 1
 
@@ -125,7 +124,7 @@ def zonal_stats(
                 elif image_info.imagetype.lower() == "s1-grd-sigma0-desc":
                     orbit = "DESCENDING"
                 output_base_path = _general_helper._format_output_path(
-                    features_path,
+                    vector_path,
                     image_path,
                     output_dir,
                     orbit,
@@ -141,7 +140,7 @@ def zonal_stats(
                 for band in bands:
                     # Prepare the output paths...
                     output_band_path = _general_helper._format_output_path(
-                        features_path,
+                        vector_path,
                         image_path,
                         output_dir,
                         orbit,
@@ -179,16 +178,16 @@ def zonal_stats(
                 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,
+                    features_path=vector_path,
+                    id_column=id_column,
+                    image_path=image_path,
+                    temp_dir=temp_dir,
+                    log_dir=log_dir,
+                    log_level=log_level,
+                    nb_parallel_max=nb_parallel_max,
                 )
                 image_dict[image_path_str] = {
-                    "features_path": features_path,
+                    "features_path": vector_path,
                     "image_path": image_path,
                     "bands": bands,
                     "bands_todo": bands_todo,
@@ -264,15 +263,15 @@ def zonal_stats(
                             start_time_batch = datetime.now()
                             future = pool.submit(
                                 _zonal_stats_image_gdf,
-                                features_batch["path"],
-                                id_column,
-                                image["image_prepared_path"],
-                                band,
-                                image["output_base_busy_path"],
-                                log_dir,
-                                log_level,
-                                start_time_batch,
-                                cloud_filter_band,
+                                features=features_batch["path"],
+                                id_column=id_column,
+                                image_path=image["image_prepared_path"],
+                                band=band,
+                                output_base_path=image["output_base_busy_path"],
+                                log_dir=log_dir,
+                                log_level=log_level,
+                                future_start_time=start_time_batch,
+                                cloud_filter_band=cloud_filter_band,
                             )
                             band_id = band if isinstance(band, str) else "-".join(band)
                             batch_id = f"{features_batch['id']}_{band_id}"
@@ -393,9 +392,9 @@ def zonal_stats(
                 progress_msg = _general_helper._format_progress_message(
                     nb_todo,
                     nb_done_total,
-                    nb_done_latestbatch,
                     start_time,
-                    image["calc_starttime"],
+                    nb_done_latestbatch,
+                    start_time_latestbatch=image["calc_starttime"],
                 )
                 logger.info(progress_msg)
 
@@ -406,7 +405,7 @@ def zonal_stats(
             #     - 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)
+                image_idx == len(rasters_bands)
                 and len(image_dict)
                 == len(_filter_on_status(image_dict, "IMAGE_CALC_DONE"))
             ):
@@ -463,7 +462,7 @@ def _prepare_calc(
     image_path: Path,
     temp_dir: Path,
     log_dir: Path,
-    log_level: Union[str, int],
+    log_level: int,
     nb_parallel_max: int = 16,
 ) -> dict:
     """
@@ -476,10 +475,8 @@ 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"
-        )
+        log_name = f"{datetime.now():%Y-%m-%d_%H-%M-%S}_prepare_calc_{os.getpid()}.log"
+        log_path = log_dir / log_name
         filehandler = logging.FileHandler(filename=log_path)
         filehandler.setLevel(log_level)
         log_format = "%(asctime)s|%(levelname)s|%(name)s|%(funcName)s|%(message)s"
diff --git a/environment-dev.yml b/environment-dev.yml
index af4f70d4..c994b579 100644
--- a/environment-dev.yml
+++ b/environment-dev.yml
@@ -17,6 +17,7 @@ dependencies:
   - psutil
   - pyproj
   - scikit-learn
+  - qgis
   - rasterio
   - rasterstats
   # testing
@@ -27,4 +28,5 @@ dependencies:
   - flake8
   - pre-commit
   - pip:
-    - tensorflow-cpu
+    # Loading the classification model gives an error with 2.12?
+    - tensorflow <2.12


Test patch:
diff --git a/tests/test_calc_timeseries.py b/tests/test_calc_timeseries.py
index 689725dc..31ee7e86 100644
--- a/tests/test_calc_timeseries.py
+++ b/tests/test_calc_timeseries.py
@@ -27,9 +27,9 @@ def test_calc_stats_per_image_s1_bs(tmpdir):
         start_time = datetime.now()
         images_bands = [(path, ["VV", "VH"]) for path in input_image_paths]
         zonal_stats_bulk.zonal_stats(
-            features_path=input_features_path,
+            vector_path=input_features_path,
             id_column="UID",
-            images_bands=images_bands,
+            rasters_bands=images_bands,
             output_dir=tmp_dir,
             temp_dir=tmp_dir / "tmp",
             log_dir=tmp_dir / "log",
@@ -56,9 +56,9 @@ def test_calc_stats_per_image_s1_bs(tmpdir):
         start_time = datetime.now()
         images_bands = [(path, ["VV", "VH"]) for path in input_image_paths]
         zonal_stats_bulk.zonal_stats(
-            features_path=input_features_path,
+            vector_path=input_features_path,
             id_column="UID",
-            images_bands=images_bands,
+            rasters_bands=images_bands,
             output_dir=tmp_dir,
             temp_dir=tmp_dir / "tmp",
             log_dir=tmp_dir / "log",