improve the robustness of the FIM evaluation with building footprint

Author: supathdhitalGEO

dist/fimeval-0.1.46-py3-none-any.whl

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "fimeval"
-version = "0.1.45"
+version = "0.1.46"
 description = "A Framework for Automatic Evaluation of Flood Inundation Mapping Predictions Evaluation"
 authors = [
     "Surface Dynamics Modeling Lab",

dist/fimeval-0.1.46.tar.gz

@@ -21,18 +21,6 @@ def Changeintogpkg(input_path, output_dir, layer_name):
         gdf.to_file(output_gpkg, driver="GPKG")
         return output_gpkg
 
-
-def Changeintogpkg(input_path, output_dir, layer_name):
-    input_path = str(input_path)
-    if input_path.endswith(".gpkg"):
-        return input_path
-    else:
-        gdf = gpd.read_file(input_path)
-        output_gpkg = os.path.join(output_dir, f"{layer_name}.gpkg")
-        gdf.to_file(output_gpkg, driver="GPKG")
-        return output_gpkg
-
-
 def GetFloodedBuildingCountInfo(
     building_fp_path,
     study_area_path,
@@ -51,8 +39,17 @@ def GetFloodedBuildingCountInfo(
     building_gdf = gpd.read_file(building_fp_gpkg)
     study_area_gdf = gpd.read_file(study_area_path)
 
-    if building_gdf.crs != study_area_gdf.crs:
-        building_gdf = building_gdf.to_crs(study_area_gdf.crs)
+    with rasterio.open(raster1_path) as src:
+        target_crs = str(src.crs)
+
+    # Reproject all GeoDataFrames to the target CRS
+    if building_gdf.crs != target_crs:
+        building_gdf = building_gdf.to_crs(target_crs)
+        print("reproject building_gdf")
+
+    if study_area_gdf.crs != target_crs:
+        study_area_gdf = study_area_gdf.to_crs(target_crs)
+        print("reproject study_area_gdf")
 
     clipped_buildings = gpd.overlay(building_gdf, study_area_gdf, how="intersection")
     clipped_buildings["centroid"] = clipped_buildings.geometry.centroid
@@ -85,7 +82,7 @@ def count_centroids_in_raster(raster_path, label):
                         elif pixel_value == 4:
                             centroid_counts["True Positive"] += 1
 
-    if "benchmark" in str(raster1_path).lower():
+    if "bm" in str(raster1_path).lower():
         count_centroids_in_raster(raster1_path, "Benchmark")
         count_centroids_in_raster(raster2_path, "Candidate")
     elif "candidate" in str(raster2_path).lower():
@@ -220,46 +217,32 @@ def count_centroids_in_raster(raster_path, label):
     print(f"Performance metrics chart is saved as PNG at {output_path}")
     fig.show()
 
-
 def process_TIFF(
     tif_files, contingency_files, building_footprint, boundary, method_path
 ):
     benchmark_path = None
-    candidate_path = []
-
-    if len(tif_files) == 2:
-        for tif_file in tif_files:
-            if "benchmark" in tif_file.name.lower():
-                benchmark_path = tif_file
-            else:
-                candidate_path.append(tif_file)
+    candidate_paths = []
 
-    elif len(tif_files) > 2:
-        for tif_file in tif_files:
-            if "benchmark" in tif_file.name.lower():
+    for tif_file in tif_files:
+        if "bm" in tif_file.name.lower() or "benchmark" in tif_file.name.lower():
+            if benchmark_path is None:
                 benchmark_path = tif_file
-                print(f"---Benchmark: {tif_file.name}---")
             else:
-                candidate_path.append(tif_file)
-
-    if benchmark_path and candidate_path:
-        for candidate in candidate_path:
+                candidate_paths.append(tif_file)
+        else:
+            candidate_paths.append(tif_file)
 
+    if benchmark_path and candidate_paths:
+        for candidate in candidate_paths:
             matching_contingency_map = None
             candidate_base_name = candidate.stem.replace("_clipped", "")
 
             for contingency_file in contingency_files:
                 if candidate_base_name in contingency_file.name:
                     matching_contingency_map = contingency_file
-                    print(
-                        f"Found matching contingency map for candidate {candidate.name}: {contingency_file.name}"
-                    )
                     break
 
             if matching_contingency_map:
-                print(
-                    f"---FIM evaluation with Building Footprint starts for {candidate.name}---"
-                )
                 GetFloodedBuildingCountInfo(
                     building_footprint,
                     boundary,
@@ -273,8 +256,11 @@ def process_TIFF(
                 print(
                     f"No matching contingency map found for candidate {candidate.name}. Skipping..."
                 )
-
-
+    elif not benchmark_path:
+        print("Warning: No benchmark file found.")
+    elif not candidate_paths:
+        print("Warning: No candidate files found.")
+        
 def find_existing_footprint(out_dir):
     gpkg_files = list(Path(out_dir).glob("*.gpkg"))
     return gpkg_files[0] if gpkg_files else None
@@ -337,7 +323,6 @@ def EvaluationWithBuildingFootprint(
                         )
                     else:
                         building_footprintMS = EX_building_footprint
-
                 process_TIFF(
                     tif_files,
                     contingency_files,

poetry.lock

@@ -18,12 +18,12 @@ def getContingencyMap(raster_path, method_path):
     combined_flood = np.full_like(band1, fill_value=1, dtype=int)
 
     # Map pixel values to colors
-    combined_flood[band1 == 5] = 0
-    combined_flood[band1 == 0] = 1
-    combined_flood[band1 == 1] = 2
-    combined_flood[band1 == 2] = 3
-    combined_flood[band1 == 3] = 4
-    combined_flood[band1 == 4] = 5
+    combined_flood[band1 == 5] = 5
+    combined_flood[band1 == 0] = 0
+    combined_flood[band1 == 1] = 1
+    combined_flood[band1 == 2] = 2
+    combined_flood[band1 == 3] = 3
+    combined_flood[band1 == 4] = 4
 
     # Handle NoData explicitly, mapping it to "No Data" class (1)
     if nodata_value is not None:
@@ -42,7 +42,7 @@ def getContingencyMap(raster_path, method_path):
     ys_dd = np.array(latitudes).reshape(ys.shape)
 
     # Define the color map and normalization
-    flood_colors = ["black", "white", "grey", "green", "blue", "red"]  # 6 classes
+    flood_colors = ["white", "grey", "green", "blue", "red", "black"]  # 6 classes
     flood_cmap = mcolors.ListedColormap(flood_colors)
     flood_norm = mcolors.BoundaryNorm(
         boundaries=np.arange(-0.5, 6.5, 1), ncolors=len(flood_colors)
@@ -60,6 +60,7 @@ def getContingencyMap(raster_path, method_path):
 
     # Create legend patches
     value_labels = {
+        0: "No data",
         1: "True negative",
         2: "False positive",
         3: "False negative",

pyproject.toml

@@ -132,6 +132,7 @@ def MakeFIMsUniform(fim_dir, target_crs=None, target_resolution=None):
             for src_path in tif_files:
                 dst_path = processing_folder / src_path.name
                 reprojectFIMs(str(src_path), str(dst_path), target_crs)
+                compress_tif_lzw(dst_path)
         else:
             all_within_conus = all(is_within_conus(bounds_list[i], crs_list[i]) for i in range(len(bounds_list)))
 
@@ -148,7 +149,6 @@ def MakeFIMsUniform(fim_dir, target_crs=None, target_resolution=None):
         for src_path in tif_files:
             dst_path = processing_folder / src_path.name
             shutil.copy(src_path, dst_path)
-            compress_tif_lzw(dst_path)
 
     # Resolution check and resampling
     processed_tifs = list(processing_folder.glob('*.tif'))
@@ -170,7 +170,6 @@ def MakeFIMsUniform(fim_dir, target_crs=None, target_resolution=None):
                 for src_path in processed_tifs:
                     resample_to_resolution(str(src_path), target_resolution, target_resolution)
             else:
-                print("FIMs are in different resolution after projection. \n")
                 coarser_x = max(res[0] for res in resolutions)
                 coarser_y = max(res[1] for res in resolutions)
                 print(f"Using coarser resolution: X={coarser_x}, Y={coarser_y}. Resampling all FIMS to this resolution.")

src/fimeval/BuildingFootprint/evaluationwithBF.py

@@ -1,53 +1,51 @@
-import fimpef as fp
+import fimeval as fe
 from pathlib import Path
 
 Main_dir = (
     # "../docs/sampledata"
-    '/Users/Supath/Downloads/SDML/FIMeval/FIMS'
 )
 PWD_dir = "./path/to/PWB"
-# output_dir = "../docs/Output"
-output_dir = '/Users/Supath/Downloads/SDML/FIMeval/output'
+output_dir = "../docs/Output"
 target_crs = "EPSG:5070" # Target CRS for reprojecting the FIMs, need to be in EPSG code of Projected CRS
-target_resolution = 30 #This will be in meters, if it passes the FIMS will be resampled to this resolution else, it will find the coarser resolution among all FIMS for this case and use that to resample!
+target_resolution = 10 #This will be in meters, if it passes the FIMS will be resampled to this resolution else, it will find the coarser resolution among all FIMS for this case and use that to resample!
 
 
-building_footprint = "./path/to/building_footprint"
+building_footprint = "../path/to/building_footprint.shp" # If user is working with user defined building footprint shapefile
 
 # If user is working with user defined shapefile
-# AOI = "./path/to/shapefile"
+AOI = "path/to/shapefile.shp"  # This shapefile should be in projected CRS, if not, it will be reprojected to the target CRS
 
 # Three methods are available to evaluate the FIM,
 # 1. Smallest extent
 # 2. Convex Hull
 # 3. AOI (User defined shapefile)
-method_name = "smallest_extent"
+method_name = "AOI"
 
 # 3 letter country ISO code
 countryISO = "USA"
 
 # Run the evaluation
 # It has the Permanent Water Bodies (PWB) dataset as default for United States
-fp.EvaluateFIM(Main_dir, method_name, output_dir)
+# fe.EvaluateFIM(Main_dir, method_name, output_dir)
 
 # OR, If the Evaluation Study Area is outside the US or, user has their own PWB dataset
-# fp.EvaluateFIM(Main_dir, method_name, output_dir, PWB_dir=PWD_dir)
+# fe.EvaluateFIM(Main_dir, method_name, output_dir, PWB_dir=PWD_dir)
 
 
 #If the FIMS are not in projected crs or are in different spatial resolution
-fp.EvaluateFIM(Main_dir, method_name, output_dir, PWB_dir=PWD_dir, target_crs=target_crs, target_resolution=target_resolution)
+fe.EvaluateFIM(Main_dir, method_name, output_dir, target_crs=target_crs, shapefile_dir = AOI, target_resolution=target_resolution)
 
 # Once the FIM evaluation is done, print the contingency map
-fp.PrintContingencyMap(Main_dir, method_name, output_dir)
+fe.PrintContingencyMap(Main_dir, method_name, output_dir)
 
 
 # Plot rhe evaluation metrics after the FIM evaluation
-fp.PlotEvaluationMetrics(Main_dir, method_name, output_dir)
+fe.PlotEvaluationMetrics(Main_dir, method_name, output_dir)
 
-# FIM Evaluation with Building Footprint (by default, it uses the Microsoft Building Footprint dataset)
-fp.EvaluationWithBuildingFootprint(
+# # FIM Evaluation with Building Footprint (by default, it uses the Microsoft Building Footprint dataset)
+fe.EvaluationWithBuildingFootprint(
     Main_dir, method_name, output_dir, country=countryISO
 )
 
 # If user have their own building footprint dataset, they can use it as well
-# fp.EvaluationWithBuildingFootprint(Main_dir, method_name, output_dir, building_footprint = building_footprint, shapefile_dir=AOI)
+fe.EvaluationWithBuildingFootprint(Main_dir, method_name, output_dir, building_footprint = building_footprint, shapefile_dir=AOI)