ENH: Use exactextract by default for zonalstats calculation

.github/workflows/tests.yml

@@ -28,7 +28,7 @@ jobs:
       matrix:
         os: [ubuntu-latest]
         dev: [false]
-        python: ["3.10", "3.11", "3.12"]
+        python: ["3.10", "3.11", "3.12", "3.13"]
         env: ["latest"]
         # Use openblas instead of mkl saves 600 MB. Linux OK, 50% slower on Windows and OSX!
         extra: ["nomkl"]
@@ -44,7 +44,7 @@ jobs:
           - env: latest
             os: windows-latest
             dev: false
-            python: "3.12"
+            python: "3.13"
 
     steps:
       - uses: actions/checkout@v6

CHANGELOG.md

@@ -52,6 +52,7 @@
 - Add first version of a cover/bare soil marker (#168, #200)
 - Avoid high committed memory for zonal stats calculation (#197)
 - General small improvements, e.g. save randomforest models compressed,.. (#144)
+- Use Exactextract as default engine for zonalstats calculation (#181)
 
 ### Bugs fixed
 

cropclassification/__init__.py

@@ -1 +1,7 @@
 """Package with functionalities to support agricultural parcel monitoring."""
+
+import logging
+
+# Disable info logging pf botocore.credentials
+logger_botocore_credentials = logging.getLogger("botocore.credentials")
+logger_botocore_credentials.setLevel(logging.WARNING)

cropclassification/calc_cropclass.py

@@ -189,6 +189,7 @@ def run_cropclass(
     parceldata_aggregations_to_use = conf.marker.getlist(
         "parceldata_aggregations_to_use"
     )
+
     ts.calc_timeseries_data(
         input_parcel_path=imagedata_input_parcel_path,
         roi_bounds=tuple(conf.roi.getlistfloat("roi_bounds")),

cropclassification/general.ini

@@ -138,14 +138,30 @@ on_missing_image = calculate_raise
 # Configuration on how/which periodic images/timeseries statistics should be generated.
 [timeseries]
 
+# Engine to use for the timeseries calculation.
+# Possible values are pyqgis, rasterstats and exactextract.
+engine = exactextract
+
+# Stats to calculate for the timeseries. The following stats are available:
+#     - "rasterstats": documentation: https://pythonhosted.org/rasterstats/manual.html#statistics
+#     - "pyqgis": "count", "sum", "mean", "median", "std", "min", "max", "range", "minority", "majority" and "variance".
+#     - "exactextract": documentation: https://isciences.github.io/exactextract/operations.html
+stats = [
+        "count(min_coverage_frac=1,coverage_weight=none)",
+        "mean(min_coverage_frac=1,coverage_weight=none)",
+        "median(min_coverage_frac=1,coverage_weight=none)",
+        "stdev(min_coverage_frac=1,coverage_weight=none)",
+        "min(min_coverage_frac=1,coverage_weight=none)",
+        "max(min_coverage_frac=1,coverage_weight=none)"
+    ]
 # Negative buffer to apply to input parcels to account for mixels.
-buffer = 5
+buffer = 0
 
 # 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 = 80
+min_parcels_with_data_pct = 75
 
 # Configuration specific to the marker being calculated.
 [marker]

cropclassification/helpers/config_helper.py

@@ -146,6 +146,7 @@ def read_config(
             "list": lambda x: [i.strip() for i in x.split(",")],
             "listint": lambda x: [int(i.strip()) for i in x.split(",")],
             "listfloat": lambda x: [float(i.strip()) for i in x.split(",")],
+            "jsonlist": lambda x: None if x is None else json.loads(x),
             "dict": lambda x: None if x is None else json.loads(x),
             "path": lambda x: None if x is None else Path(x),
         },

cropclassification/preprocess/_timeseries_calc_openeo.py

@@ -25,6 +25,8 @@ def calculate_periodic_timeseries(
     imageprofiles: dict[str, ImageProfile],
     images_periodic_dir: Path,
     timeseries_periodic_dir: Path,
+    engine: str,
+    stats: list[str],
     nb_parallel: int,
     on_missing_image: str,
     force: bool = False,
@@ -44,6 +46,16 @@ def calculate_periodic_timeseries(
         images_periodic_dir (Path): directory where the images are stored.
         timeseries_periodic_dir (Path): directory where the timeseries data will be
             saved.
+        engine (str): the engine to use for the calculation. Options are
+            "exactextract", "rasterstats" and "pyqgis".
+        stats (list[str]): statistics to calculate. Available statistics and
+            special options are dependent on the `engine` specified:
+
+                - "rasterstats": `rasterstats documentation <https://pythonhosted.org/rasterstats/manual.html#statistics>`_
+                - "pyqgis": "count", "sum", "mean", "median", "std", "min", "max",
+                        "range", "minority", "majority" and "variance".
+                - "exactextract": `exactextract documentation <https://isciences.github.io/exactextract/operations.html>`_
+
         nb_parallel (int): number of parallel processes to use.
         on_missing_image (str): what to do when an image is missing. Options are:
 
@@ -88,14 +100,13 @@ def calculate_periodic_timeseries(
     if temp_dir == "None":
         temp_dir = Path(tempfile.gettempdir())
 
-    logger.info(f"Calculating timeseries for {len(images_bands)} images")
     zonal_stats_bulk.zonal_stats(
         vector_path=input_parcel_path,
         id_column=conf.columns["id"],
         rasters_bands=images_bands,
         output_dir=timeseries_periodic_dir,
-        stats=["count", "mean", "median", "std", "min", "max"],
-        engine="pyqgis",
+        stats=stats,
+        engine=engine,
         nb_parallel=nb_parallel,
         force=force,
     )

cropclassification/preprocess/timeseries.py

@@ -64,6 +64,8 @@ def calc_timeseries_data(
             imageprofiles=conf.image_profiles,
             images_periodic_dir=conf.paths.getpath("images_periodic_dir"),
             timeseries_periodic_dir=timeseries_periodic_dir,
+            engine=conf.timeseries.get("engine"),
+            stats=conf.timeseries.getjsonlist("stats"),
             nb_parallel=conf.general.getint("nb_parallel", -1),
             on_missing_image=conf.images.get("on_missing_image", "calculate_raise"),
             force=force,

cropclassification/util/zonal_stats_bulk/__init__.py

@@ -1,5 +1,6 @@
 """Calculate zonal statistics for a vector file with many raster files."""
 
+import logging
 from pathlib import Path
 
 from . import (
@@ -9,16 +10,18 @@
 )
 from ._raster_helper import *  # noqa: F403
 
+logger = logging.getLogger(__name__)
+
 
 def zonal_stats(
     vector_path: Path,
     id_column: str,
     rasters_bands: list[tuple[Path, list[str]]],
     output_dir: Path,
+    engine: str,
     stats: list[str] | str | None = None,
     cloud_filter_band: str | None = None,
     calc_bands_parallel: bool = True,
-    engine: str = "rasterstats",
     nb_parallel: int = -1,
     force: bool = False,
 ) -> None:
@@ -31,21 +34,19 @@ def zonal_stats(
         rasters_bands (List[Tuple[Path, List[str]]]): List of tuples with the path to
             the raster files and the bands to calculate the zonal statistics on.
         output_dir (Path): directory to write the results to.
+        engine (str): the engine to use for the calculation. Options are
+            "exactextract", "rasterstats" and "pyqgis".
         stats (List[str]): statistics to calculate. Default to ["count", "median"].
             Available statistics and special options are dependent on the `engine`
             specified:
-
                 - "rasterstats": `rasterstats documentation <https://pythonhosted.org/rasterstats/manual.html#statistics>`_
                 - "pyqgis": "count", "sum", "mean", "median", "std", "min", "max",
                         "range", "minority", "majority" and "variance".
                 - "exactextract": `exactextract documentation <https://isciences.github.io/exactextract/operations.html>`_
-
         cloud_filter_band (str, optional): the band to use as a cloud filter. Only
             supported for engine "rasterstats". Defaults to None.
         calc_bands_parallel (bool, optional): True to calculate the bands in parallel.
             Only supported for engine "rasterstats". Defaults to True.
-        engine (str, optional): the engine to use for the calculation. Options are
-            "exactextract", "rasterstats" and "pyqgis". Defaults to "rasterstats".
         nb_parallel (int, optional): the number of parallel processes to use.
             Defaults to -1: use all available processors.
         force (bool, optional): False to skip calculating existing output files. True to

cropclassification/util/zonal_stats_bulk/_zonal_stats_bulk_exactextract.py

@@ -137,6 +137,7 @@ def zonal_stats(
                 tmp_dir=tmp_dir,
                 include_cols=columns,
                 output_paths=output_paths,
+                force=force,
             )
             calc_queue[future] = {
                 "vector_path": vector_path,
@@ -236,9 +237,10 @@ def zonal_stats_band(
             output="pandas",
             include_cols=include_cols,
         )
-
-    except Exception:
-        raise
+    except Exception as ex:
+        message = f"Error calculating zonal stats {stats}: {ex}"
+        logger.error(message)
+        raise ValueError(message) from ex
 
     return stats_df
 
@@ -253,13 +255,6 @@ def zonal_stats_band_tofile(
     include_cols: list[str],
     force: bool = False,
 ) -> dict[str, Path]:
-    # Init
-    if all(output_path.exists() for output_path in output_paths.values()):
-        if force:
-            for output_path in output_paths.values():
-                output_path.unlink(missing_ok=True)
-        return output_paths
-
     stats_df = zonal_stats_band(
         vector_path=vector_path,
         raster_path=raster_path,
@@ -273,17 +268,26 @@ def zonal_stats_band_tofile(
     for band in bands:
         index = raster_info.bands[band].band_index
         band_columns = include_cols.copy()
-        band_columns.extend(
-            [f"band_{index}_{stat}" for stat in [stat.split("(")[0] for stat in stats]]
-        )
-        band_stats_df = stats_df[band_columns].copy()
-        band_stats_df.rename(
-            columns={
-                f"band_{index}_{stat}": stat
-                for stat in [stat.split("(")[0] for stat in stats]
-            },
-            inplace=True,
-        )
+        if len(bands) == 1:
+            band_columns.extend(
+                [f"{stat}" for stat in [stat.split("(")[0] for stat in stats]]
+            )
+            band_stats_df = stats_df[band_columns].copy()
+        else:
+            band_columns.extend(
+                [
+                    f"band_{index}_{stat}"
+                    for stat in [stat.split("(")[0] for stat in stats]
+                ]
+            )
+            band_stats_df = stats_df[band_columns].copy()
+            band_stats_df.rename(
+                columns={
+                    f"band_{index}_{stat}": stat
+                    for stat in [stat.split("(")[0] for stat in stats]
+                },
+                inplace=True,
+            )
         # Add fid column to the beginning of the dataframe
         band_stats_df.insert(0, "fid", range(len(band_stats_df)))
 

markers/_images_periodic/roi_test/s1-grd-sigma0-vvdvh-asc-weekly/s1-grd-sigma0-vvdvh-asc-weekly_2024-03-04_2024-03-10_vvdvh_last.tif.json

@@ -0,0 +1,31 @@
+{
+    "imageprofile": "s1-grd-sigma0-vvdvh-asc-weekly",
+    "collection": null,
+    "index_type": "vvdvh",
+    "max_cloud_cover": null,
+    "image_source": "local",
+    "base_imageprofile": "s1-grd-sigma0-asc-weekly",
+    "pixel_type": "FLOAT32",
+    "satellite": "s1",
+    "roi_bounds": [
+        161400.0,
+        188000.0,
+        161900.0,
+        188500.0
+    ],
+    "roi_crs": "31370",
+    "start_date": "2024-03-04",
+    "end_date_incl": "2024-03-10",
+    "end_date": "2024-03-11",
+    "period_name": "weekly",
+    "weeks": [
+        10
+    ],
+    "bands": [
+        "vvdvh"
+    ],
+    "time_reducer": "last",
+    "path": "C:/Users/local_KRIWAY/Temp/1/pytest-of-KRIWAY/pytest-5/test_task_calc_periodic_mosaic0/markers/_images_periodic/roi_test/s1-grd-sigma0-vvdvh-asc-weekly/s1-grd-sigma0-vvdvh-asc-weekly_2024-03-04_2024-03-10_vvdvh_last.tif",
+    "job_options": null,
+    "process_options": null
+}

markers/_images_periodic/roi_test/s1-grd-sigma0-vvdvh-asc-weekly/s1-grd-sigma0-vvdvh-asc-weekly_2024-03-11_2024-03-17_vvdvh_last.tif.json

@@ -0,0 +1,31 @@
+{
+    "imageprofile": "s1-grd-sigma0-vvdvh-asc-weekly",
+    "collection": null,
+    "index_type": "vvdvh",
+    "max_cloud_cover": null,
+    "image_source": "local",
+    "base_imageprofile": "s1-grd-sigma0-asc-weekly",
+    "pixel_type": "FLOAT32",
+    "satellite": "s1",
+    "roi_bounds": [
+        161400.0,
+        188000.0,
+        161900.0,
+        188500.0
+    ],
+    "roi_crs": "31370",
+    "start_date": "2024-03-11",
+    "end_date_incl": "2024-03-17",
+    "end_date": "2024-03-18",
+    "period_name": "weekly",
+    "weeks": [
+        11
+    ],
+    "bands": [
+        "vvdvh"
+    ],
+    "time_reducer": "last",
+    "path": "C:/Users/local_KRIWAY/Temp/1/pytest-of-KRIWAY/pytest-5/test_task_calc_periodic_mosaic0/markers/_images_periodic/roi_test/s1-grd-sigma0-vvdvh-asc-weekly/s1-grd-sigma0-vvdvh-asc-weekly_2024-03-11_2024-03-17_vvdvh_last.tif",
+    "job_options": null,
+    "process_options": null
+}

markers/_images_periodic/roi_test/s1-grd-sigma0-vvdvh-desc-weekly/s1-grd-sigma0-vvdvh-desc-weekly_2024-03-04_2024-03-10_vvdvh_last.tif.json

@@ -0,0 +1,31 @@
+{
+    "imageprofile": "s1-grd-sigma0-vvdvh-desc-weekly",
+    "collection": null,
+    "index_type": "vvvh",
+    "max_cloud_cover": null,
+    "image_source": "local",
+    "base_imageprofile": "s1-grd-sigma0-desc-weekly",
+    "pixel_type": "FLOAT32",
+    "satellite": "s1",
+    "roi_bounds": [
+        161400.0,
+        188000.0,
+        161900.0,
+        188500.0
+    ],
+    "roi_crs": "31370",
+    "start_date": "2024-03-04",
+    "end_date_incl": "2024-03-10",
+    "end_date": "2024-03-11",
+    "period_name": "weekly",
+    "weeks": [
+        10
+    ],
+    "bands": [
+        "vvdvh"
+    ],
+    "time_reducer": "last",
+    "path": "C:/Users/local_KRIWAY/Temp/1/pytest-of-KRIWAY/pytest-5/test_task_calc_periodic_mosaic0/markers/_images_periodic/roi_test/s1-grd-sigma0-vvdvh-desc-weekly/s1-grd-sigma0-vvdvh-desc-weekly_2024-03-04_2024-03-10_vvdvh_last.tif",
+    "job_options": null,
+    "process_options": null
+}

markers/_images_periodic/roi_test/s1-grd-sigma0-vvdvh-desc-weekly/s1-grd-sigma0-vvdvh-desc-weekly_2024-03-11_2024-03-17_vvdvh_last.tif.json

@@ -0,0 +1,31 @@
+{
+    "imageprofile": "s1-grd-sigma0-vvdvh-desc-weekly",
+    "collection": null,
+    "index_type": "vvvh",
+    "max_cloud_cover": null,
+    "image_source": "local",
+    "base_imageprofile": "s1-grd-sigma0-desc-weekly",
+    "pixel_type": "FLOAT32",
+    "satellite": "s1",
+    "roi_bounds": [
+        161400.0,
+        188000.0,
+        161900.0,
+        188500.0
+    ],
+    "roi_crs": "31370",
+    "start_date": "2024-03-11",
+    "end_date_incl": "2024-03-17",
+    "end_date": "2024-03-18",
+    "period_name": "weekly",
+    "weeks": [
+        11
+    ],
+    "bands": [
+        "vvdvh"
+    ],
+    "time_reducer": "last",
+    "path": "C:/Users/local_KRIWAY/Temp/1/pytest-of-KRIWAY/pytest-5/test_task_calc_periodic_mosaic0/markers/_images_periodic/roi_test/s1-grd-sigma0-vvdvh-desc-weekly/s1-grd-sigma0-vvdvh-desc-weekly_2024-03-11_2024-03-17_vvdvh_last.tif",
+    "job_options": null,
+    "process_options": null
+}