diff --git a/scripts/compute_all_masks.sh b/scripts/compute_all_masks.sh index 225d6b7..34e8a90 100755 --- a/scripts/compute_all_masks.sh +++ b/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 diff --git a/slurp/masks/stack_masks.py b/slurp/masks/stack_masks.py index d372b63..24f50ed 100644 --- a/slurp/masks/stack_masks.py +++ b/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...", diff --git a/slurp/masks/watermask.py b/slurp/masks/watermask.py index fcae5ec..dd60069 100644 --- a/slurp/masks/watermask.py +++ b/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 ) diff --git a/slurp/tools/pydantic_class.py b/slurp/tools/pydantic_class.py index 8fae874..85da7ef 100644 --- a/slurp/tools/pydantic_class.py +++ b/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", diff --git a/tests/test_features/test_features_watermask.py b/tests/test_features/test_features_watermask.py index f401d19..73c4c36 100644 --- a/tests/test_features/test_features_watermask.py +++ b/tests/test_features/test_features_watermask.py @@ -63,14 +63,16 @@ def test_hand_strict_ci( sys.argv = command slurp.masks.watermask.main() +""" +TODO : replace these tests by new tests on different modes @pytest.mark.features def test_simple_ndwi_threshold( main_config, features_test_img, output_dir, ref_dir ): - """Tests the water mask computation with a simple NDWI threshold enabled. + Tests the water mask computation with a simple NDWI threshold enabled. simple_ndwi_threshold: Compute water mask as a simple NDWI threshold, - useful in arid places where no water is known by Peckel""" + useful in arid places where no water is known by Peckel cmd = write_command_compute_watermask( 1, main_config, features_test_img, output_dir, ref_dir ) @@ -83,14 +85,14 @@ def test_simple_ndwi_threshold( def test_simple_ndwi_threshold_ci( main_config, features_test_img, output_dir, ref_dir, valid_stack ): - """Run test_simple_nwdi_threshold with specified valid_stack (for GithubCI).""" +Run test_simple_nwdi_threshold with specified valid_stack (for GithubCI). cmd = write_command_compute_watermask( 1, main_config, features_test_img, output_dir, ref_dir, valid_stack ) command = f"{cmd} -simple_ndwi_threshold True".split() sys.argv = command slurp.masks.watermask.main() - +""" @pytest.mark.features @pytest.mark.parametrize("samples_method", ["random", "smart", "grid"])