diff --git a/slurp/eomultiprocessing/slurp_executor.py b/slurp/eomultiprocessing/slurp_executor.py index 122a4d4..b40c435 100644 --- a/slurp/eomultiprocessing/slurp_executor.py +++ b/slurp/eomultiprocessing/slurp_executor.py @@ -352,19 +352,39 @@ def mp_n_to_m_images( ) -def _find_spatial_axes(arr: np.ndarray, h: int, w: int): +def _find_spatial_axes(arr: np.ndarray, h: int, w: int) -> tuple[int, int]: """ - Find which axes correspond to spatial dimensions (H, W) - regardless of axis order. - - Returns: - (h_axis, w_axis) + Find spatial axes (H, W) inside an arbitrary tensor layout. + + The function does NOT assume spatial axes are the last dimensions. + + Examples + -------- + (H, W) -> (0, 1) + (C, H, W) -> (1, 2) + (H, W, C) -> (0, 1) + (T, C, H, W) -> (2, 3) + + Parameters + ---------- + arr : np.ndarray + Input tensor. + h : int + Expected spatial height. + w : int + Expected spatial width. + + Returns + ------- + tuple[int, int] + Indices of spatial axes. """ matches = [] for i, size_i in enumerate(arr.shape): for j, size_j in enumerate(arr.shape): + if i == j: continue @@ -373,12 +393,13 @@ def _find_spatial_axes(arr: np.ndarray, h: int, w: int): if not matches: raise ValueError( - f"Cannot find spatial axes matching ({h},{w}) " + f"Cannot find spatial axes matching ({h}, {w}) " f"in array shape {arr.shape}" ) - # take first valid match - return matches[0] + # Prefer the last matching pair + # because spatial axes are usually deeper + return matches[-1] # ============================================================ diff --git a/slurp/masks/watermask.py b/slurp/masks/watermask.py index e58923f..95e42af 100644 --- a/slurp/masks/watermask.py +++ b/slurp/masks/watermask.py @@ -241,10 +241,7 @@ def build_samples( valid_stack: np.ndarray, mask_hand: np.ndarray, mask_pekel: np.ndarray, - vhr1: np.ndarray, - vhr2: np.ndarray, - vhr3: np.ndarray, - vhr4: np.ndarray, + vhr: np.ndarray, ndvi: np.ndarray, ndwi: np.ndarray, ndwi_threshold: float, @@ -299,7 +296,6 @@ def build_samples( """ if aux_inputs is None: aux_inputs = [] - # ---- validity mask ---- validity_mask = valid_stack == 0 @@ -359,7 +355,13 @@ def build_samples( cols = np.concatenate((cols_pekel, cols_hand)) # ---- stack features ---- - base_features = [mask_pekel, vhr1, vhr2, vhr3, vhr4, ndvi, ndwi] + vhr_features = [vhr[i] for i in range(vhr.shape[0])] + base_features = [ + mask_pekel, + *vhr_features, + ndvi, + ndwi, + ] features = base_features + list(aux_inputs) im_stack = np.stack(features, axis=0) @@ -373,10 +375,7 @@ def build_samples( def rf_prediction( valid_stack: np.ndarray, - vhr1: np.ndarray, - vhr2: np.ndarray, - vhr3: np.ndarray, - vhr4: np.ndarray, + vhr: np.ndarray, ndvi: np.ndarray, ndwi: np.ndarray, aux_inputs=None, @@ -384,29 +383,7 @@ def rf_prediction( debug: bool = False, ): """ - Random Forest prediction - - Parameters - ---------- - valid_stack : np.ndarray - Validity mask (1, H, W) - vhr : np.ndarray - VHR features - ndvi : np.ndarray - NDVI layer(s) - ndwi : np.ndarray - NDWI layer(s) - classifier : sklearn-like estimator - Trained classifier - debug : bool - Enable memory debug logs - aux_inputs : np.ndarray - Additional feature layers - - Returns - ------- - np.ndarray - Predicted mask (uint8) + Random Forest prediction on multiband VHR tensor. """ # ------------------------------------------------- # AUX INPUT NORMALIZATION @@ -416,33 +393,54 @@ def rf_prediction( if isinstance(aux_inputs, np.ndarray): aux_inputs = [aux_inputs] - # ---- build feature stack ---------------------------------------------- - feature_layers = [ - vhr1, - vhr2, - vhr3, - vhr4, - ndvi, - ndwi, - ] + aux_inputs - # concatenate exactly like legacy algo + # ------------------------------------------------- + # MULTIBAND FEATURE EXTRACTION + # ------------------------------------------------- + vhr_features = [ + vhr[i] + for i in range(vhr.shape[0]) + ] + + # ------------------------------------------------- + # BUILD FEATURE STACK + # ------------------------------------------------- + feature_layers = ( + vhr_features + + [ndvi, ndwi] + + list(aux_inputs) + ) + im_stack = np.concatenate( [layer[np.newaxis, ...] for layer in feature_layers], axis=0, ) - # valid_stack shape expected: (1, H, W) + + # ------------------------------------------------- + # VALID MASK + # ------------------------------------------------- valid_mask = np.logical_not(valid_stack) - # ---- reshape for sklearn ---------------------------------------------- - buffer_to_predict = np.transpose(im_stack[:, valid_mask]) + # ------------------------------------------------- + # SKLEARN BUFFER + # ------------------------------------------------- + buffer_to_predict = np.transpose( + im_stack[:, valid_mask] + ) - prediction = np.zeros(im_stack.shape[1:], dtype=np.uint8) + prediction = np.zeros( + im_stack.shape[1:], + dtype=np.uint8, + ) if buffer_to_predict.shape[0] > 0: - prediction[valid_mask] = classifier.predict(buffer_to_predict) + prediction[valid_mask] = classifier.predict( + buffer_to_predict + ) - # ---- debug ------------------------------------------------------------- + # ------------------------------------------------- + # DEBUG + # ------------------------------------------------- utils.display_mem_usage( debug, f"RF Prediction on buffer " @@ -905,10 +903,7 @@ def nominal_case_predict( valid_stack[0][0], mask_hand[0], local_mask_pekel, - vhr[0][0], - vhr[0][1], - vhr[0][2], - vhr[0][3], + vhr[0], ndvi[0][0], ndwi[0][0], ] @@ -963,10 +958,7 @@ def nominal_case_predict( # -- Predict full tile -- # input_for_prediction = [ valid_stack[0][0], - vhr[0][0], - vhr[0][1], - vhr[0][2], - vhr[0][3], + vhr[0], # <-- tensor multibande complet ndvi[0][0], ndwi[0][0], ] diff --git a/slurp/tools/random_forest_utils.py b/slurp/tools/random_forest_utils.py index b2b08e1..24aecc9 100644 --- a/slurp/tools/random_forest_utils.py +++ b/slurp/tools/random_forest_utils.py @@ -31,21 +31,58 @@ logger = logging.getLogger("slurp") +import numpy as np + + def print_feature_importance(classifier, layers): - """Compute feature importance.""" - feature_names = ["R", "G", "B", "NIR", "NDVI", "NDWI"] + layers + """ + Print feature importance ranking. + + Features are assumed to be ordered as: + B1..BN, NDVI, NDWI, auxiliary layers. + """ + + n_features = len(classifier.feature_importances_) + n_aux_layers = len(layers) + + # NDVI + NDWI + n_indices = 2 + + n_bands = n_features - n_aux_layers - n_indices + + if n_bands <= 0: + raise ValueError( + f"Unable to infer spectral band count: " + f"{n_features=} {n_aux_layers=}." + ) + + feature_names = ( + [f"B{i}" for i in range(1, n_bands + 1)] + + ["NDVI", "NDWI"] + + list(layers) + ) + + if len(feature_names) != n_features: + raise ValueError( + f"Feature count mismatch: " + f"{len(feature_names)} names generated for " + f"{n_features} model features." + ) importances = classifier.feature_importances_ indices = np.argsort(importances)[::-1] std = np.std( - [tree.feature_importances_ for tree in classifier.estimators_], axis=0 + [tree.feature_importances_ for tree in classifier.estimators_], + axis=0, ) logger.info("Feature ranking:") for idx in indices: logger.info( - f" {feature_names[idx]:4s} ({importances[idx]:f}) (std={std[idx]:f})" + f"{feature_names[idx]:<10s} " + f"({importances[idx]:.6f}) " + f"(std={std[idx]:.6f})" ) diff --git a/tests/config_test_n_bands.json b/tests/config_test_n_bands.json new file mode 100644 index 0000000..d123394 --- /dev/null +++ b/tests/config_test_n_bands.json @@ -0,0 +1,6 @@ +{ + "main_config": "tests/main_n_bands_config.json", + "features_test_img": "/work/CAMPUS/etudes/Masques_CO3D/Validation/Input-Data/uc_toulouse_with_n_bands/xt_Toulouse_with_n_bands.tif", + "ref_dir": "/work/CAMPUS/etudes/Masques_CO3D/Validation/Input-Data/uc_toulouse_with_n_bands/ref", + "valid_stack": "/work/CAMPUS/etudes/Masques_CO3D/Validation/Input-Data/uc_toulouse_with_n_bands/ref/valid_stack.tif" +} diff --git a/tests/main_n_bands_config.json b/tests/main_n_bands_config.json new file mode 100644 index 0000000..77c408b --- /dev/null +++ b/tests/main_n_bands_config.json @@ -0,0 +1,137 @@ +{ + "input": { + "file_vhr": "/work/CAMPUS/etudes/Masques_CO3D/Validation/Input-Data/uc_toulouse_with_n_bands/xt_Toulouse_with_n_bands.tif", + "sensor_mode": false + }, + "aux_layers": { + "valid_stack": "/work/CAMPUS/etudes/Masques_CO3D/Validation/Input-Data/uc_toulouse_with_n_bands/ref/valid_stack.tif", + "file_ndvi": "/work/CAMPUS/etudes/Masques_CO3D/Validation/Input-Data/uc_toulouse_with_n_bands/ref/ndvi.tif", + "file_ndwi": "/work/CAMPUS/etudes/Masques_CO3D/Validation/Input-Data/uc_toulouse_with_n_bands/ref/ndwi.tif", + "extracted_pekel": "/work/CAMPUS/etudes/Masques_CO3D/Validation/Input-Data/uc_toulouse_with_n_bands/ref/pekel.tif", + "extracted_hand": "/work/CAMPUS/etudes/Masques_CO3D/Validation/Input-Data/uc_toulouse_with_n_bands/ref/hand.tif", + "extracted_wsf": "/work/CAMPUS/etudes/Masques_CO3D/Validation/Input-Data/uc_toulouse_with_n_bands/ref/wsf.tif", + "extracted_wbm": "/work/CAMPUS/etudes/Masques_CO3D/Validation/Input-Data/uc_toulouse_with_n_bands/ref/wbm.tif", + "file_texture": "/work/CAMPUS/etudes/Masques_CO3D/Validation/Input-Data/uc_toulouse_with_n_bands/ref/texture.tif", + "mnh": null, + "file_cloud_gml": null + }, + "masks": { + "watermask": "/work/CAMPUS/etudes/Masques_CO3D/Validation/Input-Data/uc_toulouse_with_n_bands/ref/watermask.tif", + "urbanmask": "/work/CAMPUS/etudes/Masques_CO3D/Validation/Input-Data/uc_toulouse_with_n_bands/ref/urbanmask.tif", + "vegetationmask": "/work/CAMPUS/etudes/Masques_CO3D/Validation/Input-Data/uc_toulouse_with_n_bands/ref/vegetationmask.tif", + "shadowmask": "/work/CAMPUS/etudes/Masques_CO3D/Validation/Input-Data/uc_toulouse_with_n_bands/ref/shadowmask.tif", + "stackmask": "/work/CAMPUS/etudes/Masques_CO3D/Validation/Input-Data/uc_toulouse_with_n_bands/ref/stackmask.tif" + }, + "resources": { + "n_workers": 8, + "tile_max_size": 1000, + "multiproc_context": "spawn", + "n_jobs": 1, + "save_mode": "none" + }, + "prepare": { + "red": 1, + "green": 2, + "nir": 4, + "cloud_mask": null, + "pekel_method": "all", + "pekel": "/work/datalake/static_aux/MASQUES/PEKEL/data2021/occurrence/occurrence.vrt", + "pekel_monthly_occurrence": "/work/datalake/static_aux/MASQUES/PEKEL/data2021/MonthlyRecurrence", + "pekel_obs": null, + "hand": "/work/datalake/static_aux/MASQUES/HAND_MERIT/hnd.vrt", + "wsf": "/work/datalake/static_aux/MASQUES/WSF/WSF2019_v1/WSF2019_v1.vrt", + "texture_rad": 5, + "dtm": null, + "geoid_file": "/softs/projets/3d/geoids/egm96.grd", + "analyse_glcm": true, + "land_cover_map": "/work/CAMPUS/DATA/ESA_WORLDCOVER/ESA_WorldCover.vrt", + "cropped_land_cover_map": false, + "effective_used_config": "out/effective_used_config.json", + "wbm": "/work/CAMPUS/etudes/Masques_CO3D/Data/WBM/wbm.vrt" + }, + "post_process": { + "binary_opening": 2, + "binary_closing": 2, + "binary_dilation": 2, + "remove_small_objects": 100, + "remove_small_holes": 100, + "area_closing": null + }, + "shadows": { + "th_rgb": 0.2, + "th_nir": 0.2, + "percentile": 2, + "absolute_threshold": false + }, + "urban": { + "files_layers": [], + "vegmask_min_value": 21, + "veg_binary_dilation": 5, + "value_classif": 255, + "gt_binary_erosion": 5, + "nb_samples_other": 5000, + "nb_samples_urban": 1000, + "max_depth": 8, + "nb_estimators": 100 + }, + "vegetation": { + "texture_mode": "yes", + "filter_texture": 90, + "slic_seg_size": 100, + "slic_compactness": 0.1, + "nb_clusters_veg": 2, + "min_ndvi_veg": null, + "max_ndvi_noveg": null, + "non_veg_clusters": null, + "nb_clusters_low_veg": 2, + "max_texture_th": null, + "debug": false, + "pct_veg": 0.19173363949483352, + "pct_low_veg": 0.023727516264829697, + "pct_high_veg": 0.16800612323000383, + "pct_non_veg": 0.14772292384232683 + }, + "water": { + "files_layers": [], + "thresh_pekel": 50, + "thresh_hand": 25, + "hand_strict": false, + "strict_thresh": 50, + "simple_ndwi_threshold": false, + "ndwi_threshold": 0.1, + "samples_method": "grid", + "nb_samples_water": 2000, + "nb_samples_other": 10000, + "nb_samples_auto": false, + "auto_pct": 0.0002, + "smart_area_pct": 50, + "smart_minimum": 10, + "grid_spacing": 40, + "max_depth": 8, + "nb_estimators": 100, + "no_pekel_filter": false, + "hand_filter": false, + "value_classif": 1 + }, + "stack": { + "building_threshold": 70, + "building_erosion": 2, + "bonus_gt": 10, + "malus_shadow": 10, + "value_classif_low_veg": 1, + "value_classif_high_veg": 2, + "value_classif_water": 3, + "value_classif_buildings": 4, + "value_classif_bare_ground": 6, + "value_classif_sea": 7, + "value_classif_lake": 8, + "value_classif_river": 9, + "value_classif_false_positive_buildings": 10, + "value_classif_background": 11, + "vegmask_min_value": 21, + "minimal_size_water_area": 10000, + "binary_closing": 0, + "binary_opening": 0, + "categorized_watermask": true + } +} diff --git a/tests/test_features/test_features_watermask.py b/tests/test_features/test_features_watermask.py index f401d19..f518b87 100644 --- a/tests/test_features/test_features_watermask.py +++ b/tests/test_features/test_features_watermask.py @@ -7,6 +7,7 @@ # """Test water mask with differents features and different arguments values""" +import json import sys import pytest @@ -214,3 +215,32 @@ def test_hand_filter_ci( command = f"{cmd} -hand_filter".split() sys.argv = command slurp.masks.watermask.main() + +@pytest.mark.features +@pytest.mark.parametrize( + "params_file", + [ + "tests/config_test_n_bands.json", + ], +) +def test_n_bands(params_file, output_dir): + """Tests the water mask computation when recieving n bands > 4""" + + # Charger le JSON + with open(params_file, "r") as f: + params = json.load(f) + + cmd = ( + write_command_compute_watermask( + 1, + params["main_config"], + params["features_test_img"], + output_dir, + params["ref_dir"], + params["valid_stack"], + ) + ) + command = f"{cmd} -nb_samples_auto".split() + + sys.argv = command + slurp.masks.watermask.main() \ No newline at end of file