Implements an alternative to random forest or absolute threshold

scripts/compute_all_masks.sh

@@ -36,6 +36,7 @@ otbcli_Superimpose -inr ${VHR_IM} -inm /work/datalake/static_aux/MASQUES/PEKEL/d
 otbcli_Superimpose -inr ${VHR_IM} -inm /work/datalake/static_aux/MASQUES/HAND_MERIT/hnd.vrt -out "out/hand.tif?&gdal:co:TILED=YES&gdal:co:COMPRESS=DEFLATE" 
 otbcli_Superimpose -inr ${VHR_IM} -inm /work/datalake/static_aux/MASQUES/WSF/WSF2019_v1/WSF2019_v1.vrt -out "out/wsf.tif?&gdal:co:TILED=YES&gdal:co:COMPRESS=DEFLATE" uint8 -interpolator nn
 otbcli_Superimpose -inr ${VHR_IM} -inm /work/CAMPUS/etudes/Masques_CO3D/Data/WBM/wbm.vrt -out "out/wbm.tif?&gdal:co:TILED=YES&gdal:co:COMPRESS=DEFLATE" uint8 -interpolator nn
+otbcli_Superimpose -inr ${VHR_IM} -inm /work/CAMPUS/DATA/ESA_WORLDCOVER/ESA_WorldCover.vrt -out "out/esa_lc.tif?&gdal:co:TILED=YES&gdal:co:COMPRESS=DEFLATE" uint8 -interpolator nn
 
 
 # Prepare

slurp/masks/stack_masks.py

@@ -603,7 +603,7 @@ def slurp_stackmask(
                 image_filter=post_process,
                 filter_parameters=vars(args),
                 generate_output_profiles=eo_utils.three_uint8_profile,
-                stable_margin=100,  # TODO : add a stability margin parameter ?
+                stable_margin=0,  # TODO : add a stability margin parameter ?
                 context_manager=eoscale_manager,
                 multiproc_context=args.multiproc_context,
                 filter_desc="Post processing...",

slurp/masks/watermask.py

@@ -68,7 +68,7 @@ def compute_pekel_mask(
         return [mask_pekel, mask_pekelxx]
 
     if not params["no_pekel_filter"]:
-        mask_pekel0 = utils.compute_mask(input_buffer[0], [0])
+        mask_pekel0 = utils.compute_mask(input_buffer[0], [30])
     else:
         mask_pekel0 = np.zeros(input_buffer[0].shape)
 
@@ -227,7 +227,7 @@ def build_samples(
     # But we propagate a simple boolean mask (heigth, width)
     validity_mask = (input_buffer[0] == 0)[0]
     valid_water_pixels = np.logical_and(
-        input_buffer[2], input_buffer[5] > params["ndwi_threshold"]
+        input_buffer[2], input_buffer[5] > 1000*params["ndwi_threshold"]
     )
     valid_water_pixels = np.logical_and(valid_water_pixels, validity_mask)
 
@@ -338,6 +338,8 @@ def mask_filter(im_in, mask_ref):
     return im_filtered
 
 
+"""
+TODO remove
 def apply_ndwi_thresh(args, eoscale_manager, key_ndwi, key_valid_stack):
     logger.info("Simple threshold mask NDWI > " + str(args.ndwi_threshold))
     key_predict = eoexe.n_images_to_m_images_filter(
@@ -353,6 +355,7 @@ def apply_ndwi_thresh(args, eoscale_manager, key_ndwi, key_valid_stack):
     time_samples = time_random_forest
     do_post_process = False
     return do_post_process, key_predict, time_random_forest, time_samples
+"""
 
 
 def post_process(
@@ -526,7 +529,7 @@ def display_computation_info(
     logger.info(
         "- Build_stack           :\t" + utils.convert_time(time_stack - t0)
     )
-    if not args.simple_ndwi_threshold and not not_enough_water_samples:
+    if not not_enough_water_samples:
         logger.info(
             "- Build_samples         :\t"
             + utils.convert_time(time_samples - time_stack)
@@ -585,12 +588,11 @@ def process_pekel(args, eoscale_manager, margin):
     return local_mask_pekel, mask_pekel, not_enough_water_samples
 
 
-def process_hand(args, eoscale_manager, margin):
+def process_hand(args, key_hand, eoscale_manager, margin):
     """
     Processes a HAND (Height Above Nearest Drainage) raster and applies the
     `compute_hand_mask` filter to create a usable mask for further processing.
     """
-    key_hand = eoscale_manager.open_raster(raster_path=args.extracted_hand)
     hand_profile = eoscale_manager.get_profile(key_hand)
     args.hand_nodata = hand_profile["nodata"]
     # Create HAND mask
@@ -606,16 +608,7 @@ def process_hand(args, eoscale_manager, margin):
         filter_desc="Hand valid mask processing...",
     )
 
-    not_enough_ground_samples = False
-    # Check Hand mask : if there are too few valid pixels, we propose to threshold NDWI
-    local_mask_hand = eoscale_manager.get_array(mask_hand[0])
-    if len(np.where(local_mask_hand)[0]) < args.nb_samples_other:
-        not_enough_ground_samples = True
-        logger.warning(
-            "** WARNING ** not enough ground samples are found in Hand : return a simple mask"
-        )
-
-    return mask_hand, not_enough_ground_samples
+    return mask_hand
 
 
 def nominal_case_predict(
@@ -712,7 +705,7 @@ def nominal_case_predict(
     return key_predict, time_random_forest, time_samples
 
 
-def launch_postprocess(
+def launch_post_process(
     args,
     eoscale_manager,
     key_predict,
@@ -759,7 +752,15 @@ def getarguments():
     parser.add_argument(
         "main_config", help="First JSON file, load basis arguments"
     )
-
+    parser.add_argument(
+        "-mode",
+        choices=["nominal", "no_HAND", "relative_th", "absolute_th"],
+        dest="mode",
+        default="nominal",
+        help="Watermask algorithm : nominal refers to Random Forest learning/prediction, "
+        "no HAND won't use HAND to select samples, relative_th computes a NDWI threshold based "
+        "on Pekel, and absolute_th simply filters NDWI values above a threshold",
+    )
     parser.add_argument(
         "-log_f",
         "--logs_to_file",
@@ -825,11 +826,6 @@ def getarguments():
         choices=["none", "debug"],
         help="Save all files (debug) or only output mask (none)",
     )
-    group2.add_argument(
-        "-simple_ndwi_threshold",
-        help="Compute water mask as a simple NDWI threshold - "
-        "useful in arid places where no water is known by Peckel",
-    )
     group2.add_argument(
         "-ndwi_threshold",
         type=float,
@@ -964,6 +960,7 @@ def getarguments():
 
 def slurp_watermask(
     main_config: str,
+    mode: str,
     logs_to_file: bool,
     debug: bool,
     user_config: str,
@@ -980,7 +977,6 @@ def slurp_watermask(
     thresh_hand: int,
     strict_thresh: float,
     save_mode: str,
-    simple_ndwi_threshold: bool,
     ndwi_threshold: float,
     samples_method: str,
     nb_samples_water: int,
@@ -1045,78 +1041,153 @@ def slurp_watermask(
             t0 = time.time()
 
             utils.display_mem_usage(args.debug, "Start computation")
-
-            # --Build stack with all layers-- #
-
-            key_ndvi, key_ndwi, key_phr, key_valid_stack, margin = (
+            """
+            Build stack with all layers (4 bands VHR image, NDVI, NDWI, valid stack)
+            TODO : do we separate read NDVI from read NDWI ?
+            """
+            key_ndvi, key_ndwi, key_vhr, key_valid_stack, margin = (
                 build_stack_water(args, eoscale_manager)
             )
-
             time_stack = time.time()
 
-            # --Build samples-- #
-
-            # Pekel
-            local_mask_pekel, mask_pekel, not_enough_water_samples = (
-                process_pekel(args, eoscale_manager, margin)
-            )
-
-            # HAND
-            mask_hand, not_enough_ground_samples = process_hand(
-                args, eoscale_manager, margin
-            )
-
-            # Flag to command post-process
-            do_post_process = True
-
-            if args.simple_ndwi_threshold or not_enough_ground_samples:
-                # Simple NDWI threshold,
-                # but taking account valid stack to take care of NO_DATA values
-                (
-                    do_post_process,
-                    key_predict,
-                    time_random_forest,
-                    time_samples,
-                ) = apply_ndwi_thresh(
-                    args, eoscale_manager, key_ndwi, key_valid_stack
+            if args.mode == "absolute_th":
+                logger.debug(
+                    f"Absolute threshold mode : filter NDWI values above "
+                    f"{args.ndwi_threshold=}"
                 )
 
-            elif not_enough_water_samples:
-                # We compute a void mask (0 everywhere, except for NO DATA values)
-                # Tips : we threshold NDWI > 1000 : no pixel should be detected.
                 key_predict = eoexe.n_images_to_m_images_filter(
                     inputs=[key_ndwi, key_valid_stack],
                     image_filter=utils.compute_mask_threshold,
-                    filter_parameters={"threshold": 1000},
+                    filter_parameters={"threshold": 1000 * args.ndwi_threshold},
                     context_manager=eoscale_manager,
                     generate_output_profiles=eo_utils.single_uint8_profile,
                     multiproc_context=args.multiproc_context,
-                    filter_desc="Void mask",
+                    filter_desc="Simple NDWI threshold",
                 )
-
+                # TODO : allow a post-processing even with NDWI thresholding?
                 do_post_process = False
-
+                time_random_forest = time.time()
+                time_samples = time_random_forest  # TODO : remove timers
+                end_time = time.time()
             else:
-                # Nominal case : select samples, train, predict
-                #
-                # Taking optional layers into account
-                key_predict, time_random_forest, time_samples = (
-                    nominal_case_predict(
-                        args,
-                        eoscale_manager,
-                        key_ndvi,
-                        key_ndwi,
-                        key_phr,
-                        key_valid_stack,
-                        local_mask_pekel,
-                        margin,
-                        mask_hand,
-                        mask_pekel,
-                    )
+                """
+                Read Pekel and check if there are enough water samples to learn/predict
+                a watermask
+                """
+                # TODO : rename local_mask_pekel and mask_pekel
+                local_mask_pekel, mask_pekel, not_enough_water_samples = (
+                    process_pekel(args, eoscale_manager, margin)
                 )
+
+                if args.mode == "relative_th":
+                    """
+                    Compute relative threshold
+                    """
+                    logger.debug(
+                        f"Relative threshold mode : compute a relative "
+                        "threshold based on Pekel"
+                    )
+
+                    args.ndwi_threshold = 0.1  # compute_relative_NDWI_threshold(key_ndwi, percentile, mask_pekel)
+
+                    key_predict = eoexe.n_images_to_m_images_filter(
+                        inputs=[key_ndwi, key_valid_stack],
+                        image_filter=utils.compute_mask_threshold,
+                        filter_parameters={
+                            "threshold": 1000 * args.ndwi_threshold
+                        },
+                        context_manager=eoscale_manager,
+                        generate_output_profiles=eo_utils.single_uint8_profile,
+                        multiproc_context=args.multiproc_context,
+                        filter_desc="Simple NDWI threshold",
+                    )
+                    do_post_process = True
+                    time_random_forest = time.time()
+                    time_samples = time_random_forest  # TODO : remove timers
+                    end_time = time.time()
+                else:
+                    # args.mode = nominal or "no HAND"
+                    if args.mode == "no_HAND":
+                        """
+                        Create a void HAND layer to select samples
+                        """
+                        logger.debug(f"no HAND mode : won't read HAND file")
+                        profile = eoscale_manager.get_profile(key_vhr)
+                        profile["count"] = 1
+                        profile["dtype"] = int  # np.float64
+                        key_hand = eoscale_manager.create_image(profile)
+                        hand = eoscale_manager.get_array(key=key_hand)
+                        hand.fill(0.0)
+                    else:
+                        key_hand = eoscale_manager.open_raster(
+                            raster_path=args.extracted_hand
+                        )
+                        logger.debug(f"nominal mode : read HAND file")
+
+                    mask_hand = process_hand(
+                        args, key_hand, eoscale_manager, margin
+                    )
+                    local_mask_hand = eoscale_manager.get_array(mask_hand[0])
+
+                    not_enough_ground_samples = False
+                    # Check Hand mask : if there are too few valid pixels, we propose to threshold NDWI
+
+                    if (
+                        len(np.where(local_mask_hand)[0])
+                        < args.nb_samples_other
+                    ):
+                        not_enough_ground_samples = True
+                        logger.warning(
+                            "** WARNING ** not enough ground samples are found in Hand : return a simple mask"
+                        )
+
+                    """
+                    select samples
+                    learn / predict (or filter NDWI, corner case)
+                    """
+                    logger.debug("Learn/Predict watermask")
+
+                    if not_enough_ground_samples or not_enough_water_samples:
+                        # We compute a simple mask based on NDWI threshold
+                        #
+                        key_predict = eoexe.n_images_to_m_images_filter(
+                            inputs=[key_ndwi, key_valid_stack],
+                            image_filter=utils.compute_mask_threshold,
+                            filter_parameters={
+                                "threshold": args.ndwi_threshold
+                            },
+                            context_manager=eoscale_manager,
+                            generate_output_profiles=eo_utils.single_uint8_profile,
+                            multiproc_context=args.multiproc_context,
+                            filter_desc="Simple NDWI threshold",
+                        )
+                        # Allow post-processing to remove false positive
+                        do_post_process = True
+                        time_random_forest, time_samples = 0, 0
+                    else:
+                        # Nominal case : select samples, train, predict
+                        key_predict, time_random_forest, time_samples = (
+                            nominal_case_predict(
+                                args,
+                                eoscale_manager,
+                                key_ndvi,
+                                key_ndwi,
+                                key_vhr,
+                                key_valid_stack,
+                                local_mask_pekel,
+                                margin,
+                                mask_hand,
+                                mask_pekel,
+                            )
+                        )
+                        do_post_process = True
+
             if do_post_process:
+                """Post-process"""
+                logger.debug("Post-process watermask")
                 # --Post_processing-- #
-                launch_postprocess(
+                launch_post_process(
                     args,
                     eoscale_manager,
                     key_predict,
@@ -1125,8 +1196,8 @@ def slurp_watermask(
                     mask_hand,
                     mask_pekel,
                 )
-
             else:
+                # Write raw prediction
                 eoscale_manager.write(
                     key=key_predict[0], img_path=args.watermask
                 )  # classif
@@ -1135,7 +1206,7 @@ def slurp_watermask(
             end_time = time.time()
 
             display_global_infos(args, end_time, t0, time_stack)
-            if not args.simple_ndwi_threshold and not not_enough_water_samples:
+            if args.mode == "nominal" or args.mode == "no_HAND":
                 display_rf_infos(
                     end_time, time_random_forest, time_samples, time_stack
                 )

slurp/tools/pydantic_class.py

@@ -250,11 +250,6 @@ class Water(BaseModel):
     strict_thresh: Optional[int] = Field(
         default=50, description="Pekel Threshold float if hand_strict"
     )
-    simple_ndwi_threshold: Optional[bool] = Field(
-        default=False,
-        description="Compute water mask as a simple NDWI threshold - useful in arid places "
-        "where no water is known by Peckel",
-    )
     ndwi_threshold: Optional[float] = Field(
         default=100,
         description="Threshold used when Pekel is empty in the area",