ENH: add feature to exclude tree crowns from agricultural parcels

.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

cropclassification/calc_cropclass.py

@@ -114,6 +114,7 @@ def run_cropclass(
     input_model_to_use_relativepath = conf.calc_marker_params.getpath(
         "input_model_to_use_relativepath"
     )
+    erase_layer_filename = conf.calc_marker_params.getpath("erase_layer_filename")
 
     # Prepare input paths
     if input_model_to_use_relativepath is not None:
@@ -138,6 +139,11 @@ def run_cropclass(
     refe_dir = conf.paths.getpath("refe_dir")
     classes_refe_path = refe_dir / classes_refe_filename
 
+    if erase_layer_filename is not None:
+        erase_layer_path = input_dir / erase_layer_filename
+    else:
+        erase_layer_path = None
+
     # Check if the necessary input files exist...
     for path in [classes_refe_path, input_parcel_path]:
         if path is not None and not path.exists():
@@ -157,8 +163,9 @@ def run_cropclass(
     # -------------------------------------------------------------
 
     # Prepare the input data for optimal image data extraction:
-    #    1) apply a negative buffer on the parcel to evade mixels
-    #    2) remove features that became null because of buffer
+    #    1) exlude erase_layer from agricultural parcels
+    #    2) apply a negative buffer on the parcel to evade mixels
+    #    3) remove features that became null because of buffer
     input_preprocessed_dir = conf.paths.getpath("input_preprocessed_dir")
     buffer = conf.timeseries.getfloat("buffer")
     input_parcel_nogeo_path = (
@@ -174,6 +181,8 @@ def run_cropclass(
         input_parcel_path=input_parcel_path,
         output_imagedata_parcel_input_path=imagedata_input_parcel_path,
         output_parcel_nogeo_path=input_parcel_nogeo_path,
+        erase_layer_path=erase_layer_path,
+        classes_refe_path=classes_refe_path,
     )
 
     # STEP 2: Get the timeseries data needed for the classification

cropclassification/general.ini

@@ -36,6 +36,8 @@ input_parcel_filename = MUST_OVERRIDE
 input_parcel_filetype = MUST_OVERRIDE
 # The lookup table (LUT) file where classes to classify to are specified
 classes_refe_filename = MUST_OVERRIDE
+# The filename of the erase_layer geofile
+erase_layer_filename
 # A groundtruth file for extra reporting, optional
 input_groundtruth_filename
 # The model to use if you want to reuse one instead of training one, optional

cropclassification/preprocess/_timeseries_helper.py

@@ -23,14 +23,18 @@ def prepare_input(
     input_parcel_path: Path,
     output_imagedata_parcel_input_path: Path,
     output_parcel_nogeo_path: Path | None = None,
+    erase_layer_path: Path | None = None,
+    classes_refe_path: Path | None = None,
     force: bool = False,
 ) -> bool:
     """Prepare a file so it is ready for timeseries extraction of sentinel images.
 
     Args:
         input_parcel_path (Path): input file
         output_imagedata_parcel_input_path (Path): prepared output file
-        output_parcel_nogeo_path (Path): output file with a copy of the non-geo data
+        output_parcel_nogeo_path (Path | None): output file with a copy of non-geo data
+        erase_layer_path (Path | None): erase layer path
+        classes_refe_path: (Path | None): classes reference path
         force: force creation, even if output file(s) exist already
     """
     # Check if parameters are OK and init some extra params
@@ -101,6 +105,17 @@ def prepare_input(
         )
         return False
 
+    # Exclude erase_layer from agricultural parcels
+    if erase_layer_path is not None and classes_refe_path is not None:
+        output_imagedata_parcel_input_path = exclude_erase_layer(
+            input_parcel_path=input_parcel_path,
+            output_imagedata_parcel_input_path=output_imagedata_parcel_input_path,
+            erase_layer_path=erase_layer_path,
+            classes_refe_path=classes_refe_path,
+        )
+    else:
+        logger.info("No erase layer file found, erase_layer will not be excluded")
+
     # Apply buffer
     parceldata_buf_gdf = parceldata_gdf.copy()
     # resolution = number of segments per circle
@@ -272,3 +287,56 @@ def get_fileinfo_timeseries_periods(path: Path) -> dict:
         }
 
     return get_fileinfo_timeseries(path)
+
+
+def exclude_erase_layer(
+    input_parcel_path: Path,
+    output_imagedata_parcel_input_path: Path,
+    classes_refe_path: Path,
+    erase_layer_path: Path,
+) -> Path:
+    """This function excludes the erase layer from the agricultural parcels."""
+    # Read reference data
+    # Select all parcels where ignore_erase_layer is True
+    classes_refe_df = pdh.read_file(classes_refe_path)
+
+    # Check if ignore_erase_layer_df has columns IGNORE_ERASE_LAYER
+    if "IGNORE_ERASE_LAYER" not in classes_refe_df.columns:
+        raise Exception("IGNORE_ERASE_LAYER column not found in classes reference file")
+
+    ignore_erase_layer_df = classes_refe_df[classes_refe_df["IGNORE_ERASE_LAYER"] == 1][
+        "CROPCODE"
+    ]
+    gewascodes = tuple(map(str, (ignore_erase_layer_df)))
+
+    logger.info("Erase layer will be excluded")
+
+    # Create temp dir to store temporary data for tracebility
+    temp_output_dir = output_imagedata_parcel_input_path.parent / "temp"
+    temp_output_dir.mkdir(parents=True, exist_ok=True)
+
+    # Select all parcels where gwscod_h is in gewascodes
+    sql_stmt = f"""
+        SELECT *
+        FROM "{{input_layer}}"
+        WHERE gwscod_h IN {gewascodes}
+    """
+    gfo.select(
+        input_path=input_parcel_path,
+        output_path=temp_output_dir / "ignore_erase_layer.gpkg",
+        sql_stmt=sql_stmt,
+    )
+    # Create erase_layer_prc file by erasing all
+    # agricultural parcels with "ignore_erase_layer" True from the erase_layer" file
+    gfo.erase(
+        input_path=erase_layer_path,
+        erase_path=temp_output_dir / "ignore_erase_layer.gpkg",
+        output_path=temp_output_dir / "erase_layer_prc.gpkg",
+    )
+    # Erase the erase_layer_prc from the input agricultural parcels
+    gfo.erase(
+        input_path=input_parcel_path,
+        erase_path=temp_output_dir / "erase_layer_prc.gpkg",
+        output_path=output_imagedata_parcel_input_path,
+    )
+    return output_imagedata_parcel_input_path

markers/_inputdata/Prc_BEFL_2023_2023_07_24_groundtruth.tsv

@@ -4,7 +4,7 @@ UID	HOOFDTEELT_CTRL_COD	HOOFDTEELT_CTRL_COD_ORIG	HOOFDTEELT_CTRL_NAM_ORIG	GESP_P
 000028105F1645A400000001	91	91	Suikerbieten		
 000028104AD7C78000000007	60	60	Grasland	MAA	BG
 000028104AD7CB7600000004	60	60	Grasland	MAA	BG
-000028104C44425B00000001	60	60	Grasland	MAA	BG
+000028104C44425B00000001	8915	8915	Bebossing		
 000028104C4440C700000001	60	60	Grasland	MAA	BG
 00002810694D925500000002	60	60	Grasland	MAA	BG
 000028104AD7C7E500000001	60	60	Grasland	MAA	BG