Merge pull request #15 from supathdhitalGEO/main

added USGS, GeoGLOWS module, update statistics, update the vizualization and code usage note

Author: supathdhitalGEO

README.md

@@ -38,6 +38,8 @@ FIMserv/
 │   └── fimserve/         
 │       ├── streamflowdata/ # Handles streamflow data
 │       │   ├── nwmretrospectivedata.py   # Processes NWM retrospective data
+│       │   ├── geoglows.py   # Module to retrieve geoglows streamflow data 
+│       │   ├── usgsdata.py   # Retrieve USGS gauge station data
 │       │   └── forecasteddata.py        # Processes all range forecasted streamflow data
 │       ├── plots/          # Vizualization functionalities
 │       ├── FIMsubset/      # Subsetting functionalities for FIM
@@ -61,6 +63,8 @@ FIMserv/
 
 Although not mandatory, 
 **we strongly recommend users create a virtual environment and install this package on that virtual environment to avoid the conflict between system dependencies and package dependencies.**
+
+**‼️ If your system doesnot have git, install it first. Download link of git for windows or MacOS: https://git-scm.com/downloads**
 ```bash
 #creating a virtual environment using conda
 conda create --name fimserve python==3.10

dist/fimserve-0.1.77-py3-none-any.whl

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "fimserve"
-version = "0.1.76"
+version = "0.1.77"
 description = "Framework which is developed with the purpose of quickly generating Flood Inundation Maps (FIM) for emergency response and risk assessment. It is developed under Surface Dynamics Modeling Lab (SDML)."
 authors = ["Supath Dhital <sdhital@crimson.ua.edu>"]
 license = "GPL-3.0"

dist/fimserve-0.1.77.tar.gz

@@ -7,8 +7,7 @@
 
 from .datadownload import setup_directories
 
-
-def runfim(code_dir, output_dir, HUC_code, data_dir):
+def runfim(code_dir, output_dir, HUC_code, data_dir, depth=False):
     original_dir = os.getcwd()
     try:
         tools_path = os.path.join(code_dir, "tools")
@@ -18,12 +17,11 @@ def runfim(code_dir, output_dir, HUC_code, data_dir):
         load_dotenv(dotenv_path)
         sys.path.append(src_path)
         sys.path.append(code_dir)
+
         HUC_code = str(HUC_code)
         HUC_dir = os.path.join(output_dir, f"flood_{HUC_code}")
-
         csv_path = data_dir
 
-        # Get the inundation file path
         discharge_basename = os.path.basename(data_dir).split(".")[0]
         inundation_dir = os.path.join(HUC_dir, f"{HUC_code}_inundation")
         temp_dir = os.path.join(inundation_dir, "temp")
@@ -32,24 +30,21 @@ def runfim(code_dir, output_dir, HUC_code, data_dir):
             os.makedirs(temp_dir)
 
         inundation_file = os.path.join(temp_dir, f"{discharge_basename}_inundation.tif")
-        # depth_file = os.path.join(
-        #     temp_dir, f"{discharge_basename}_depth.tif"
-        # )
-
         Command = [
             sys.executable,
             "inundate_mosaic_wrapper.py",
-            "-y",
-            HUC_dir,
-            "-u",
-            HUC_code,
-            "-f",
-            csv_path,
-            "-i",
-            inundation_file,
-            # "-d",
-            # depth_file,
+            "-y", HUC_dir,
+            "-u", HUC_code,
+            "-f", csv_path,
+            "-i", inundation_file,
         ]
+
+        if depth:
+            depth_file = os.path.join(temp_dir, f"{discharge_basename}_depth.tif")
+            Command += ["-d", depth_file]
+        else:
+            depth_file = None
+
         env = os.environ.copy()
         env["PYTHONPATH"] = f"{src_path}{os.pathsep}{code_dir}"
 
@@ -61,34 +56,29 @@ def runfim(code_dir, output_dir, HUC_code, data_dir):
             stderr=subprocess.PIPE,
         )
 
-        # Print the output and error (if any)
         print(result.stdout.decode())
         if result.stderr:
             print(result.stderr.decode())
 
-        # Check if the command was successful
         if result.returncode == 0:
             print(f"Inundation mapping for {HUC_code} completed successfully.")
 
-            # Move the generated files to inundation_dir
-            # if os.path.exists(inundation_file) and os.path.exists(depth_file):
             if os.path.exists(inundation_file):
                 shutil.move(inundation_file, inundation_dir)
-                # shutil.move(depth_file, inundation_dir)
 
-            # Delete the temp directory
+            if depth and depth_file and os.path.exists(depth_file):
+                shutil.move(depth_file, inundation_dir)
+
             if os.path.exists(temp_dir):
                 shutil.rmtree(temp_dir)
-
         else:
             print(f"Failed to complete inundation mapping for {HUC_code}.")
 
     finally:
         os.chdir(original_dir)
 
-
-def runOWPHANDFIM(huc):
+def runOWPHANDFIM(huc, depth=False):
     code_dir, data_dir, output_dir = setup_directories()
     discharge = glob.glob(os.path.join(data_dir, f"*{huc}*.csv"))
     for file in discharge:
-        runfim(code_dir, output_dir, huc, file)
+        runfim(code_dir, output_dir, huc, file, depth=depth)

docs/code_usage.ipynb

@@ -21,23 +21,21 @@ def calculate_metrics(nwm_data, usgs_data):
         (usgs_data - np.mean(usgs_data)) ** 2
     )
     apb = np.sum(np.abs(nwm_data - usgs_data)) / np.sum(usgs_data) * 100
-    r2 = r2_score(usgs_data, nwm_data)
 
-    return {"KGE": kge, "PBias (%)": apb, "R²": r2, "NSE": nse}
+    return {"KGE": kge, "PBias (%)": apb, "NSE": nse}
 
 
 # Dual-axis visualization
 def visualize_comparison(metrics, output_dir, usgs_site, huc):
-    metric_names = ["KGE", "R²", "NSE", "PBias (%)"]
+    metric_names = ["KGE", "NSE", "PBias (%)"]
     metric_values = [
         metrics["KGE"],
-        metrics["R²"],
         metrics["NSE"],
         metrics["PBias (%)"],
     ]
 
     fig, ax1 = plt.subplots(figsize=(6, 4))
-    bar_colors = ["tab:blue", "tab:green", "tab:purple", "tab:orange"]
+    bar_colors = ["tab:blue", "tab:purple", "tab:orange"]
 
     # Bar plot for the first three metrics
     ax1.bar(
@@ -54,7 +52,7 @@ def visualize_comparison(metrics, output_dir, usgs_site, huc):
     ax1.grid(axis="y", linestyle="--", alpha=1, zorder=0)
 
     # Vertical line between the third and fourth bars
-    plt.axvline(x=2.5, color="black", linestyle="--", linewidth=1.5)
+    plt.axvline(x=1.5, color="black", linestyle="--", linewidth=1.5)
 
     # Bar plot for PBias on secondary y-axis
     ax2 = ax1.twinx()
@@ -70,22 +68,22 @@ def visualize_comparison(metrics, output_dir, usgs_site, huc):
     ax2.tick_params(axis="y", labelcolor="red", labelsize=14)
     ax2.set_ylim(0, metric_values[-1] * 1.2)
 
-    # Add text on bars
     padding = 0.01
     for i, val in enumerate(metric_values):
-        ax = ax1 if i < 3 else ax2
-        text_color = "red" if i == 3 else "black"  # Red text for PBias
+        ax = ax1 if i < 2 else ax2
+        text_color = "red" if i == 2 else "black"
+        # Adjust padding depending on axis
+        pad = padding if i < 2 else val * 0.05  # 5% of PBias value
         ax.text(
             i,
-            val + padding,
+            val + pad,
             f"{val:.2f}",
             ha="center",
             va="bottom",
             fontsize=12,
             fontweight="bold",
             color=text_color,
         )
-
     # Title and labels
     ax1.set_xlabel("Statistical Metrics", fontsize=16)
     plt.xticks(ticks=range(len(metric_names)), labels=metric_names, fontsize=14)

poetry.lock

@@ -18,16 +18,16 @@ def getdischargeforspecifiedtime(
         df = pd.read_parquet(file)
         all_data = pd.concat([all_data, df], ignore_index=True)
 
-    df["value_time"] = pd.to_datetime(df["value_time"])
+    all_data["value_time"] = pd.to_datetime(all_data["value_time"])
 
     locationID_df = pd.read_csv(location_ids)
     location_ids = [f"nwm30-{int(fid)}" for fid in locationID_df["feature_id"]]
 
     specific_date = pd.to_datetime(specific_date)
     if date_type == "date":
-        filtered_df = df[
-            (df["location_id"].isin(location_ids))
-            & (df["value_time"].dt.date == specific_date.date())
+        filtered_df = all_data[
+            (all_data["location_id"].isin(location_ids))
+            & (all_data["value_time"].dt.date == specific_date.date())
         ].copy()
         filtered_df["feature_id"] = filtered_df["location_id"].str.replace("nwm30-", "")
         discharge_data = (
@@ -38,8 +38,8 @@ def getdischargeforspecifiedtime(
         )
         formatted_datetime = specific_date.strftime("%Y%m%d")
     else:
-        filtered_df = df[
-            (df["location_id"].isin(location_ids)) & (df["value_time"] == specific_date)
+        filtered_df = all_data[
+            (all_data["location_id"].isin(location_ids)) & (all_data["value_time"] == specific_date)
         ].copy()
         filtered_df.loc[:, "feature_id"] = filtered_df["location_id"].str.replace(
             "nwm30-", ""

pyproject.toml

@@ -28,7 +28,6 @@ def getusgs_discharge(
         output_parquet_dir=output_dir,
         overwrite_output=True,
     )
-    print(f"USGS discharge data saved to {output_dir}.")
 
 #If value_times is mentioned and user need the discharge for specific time
 def getdischargeforspecifiedtime(
@@ -73,7 +72,9 @@ def getdischargeforspecifiedtime(
         .reset_index()
         .rename(columns={"value": "discharge"})
     )
-
+    #Make sure the columns are in int and float
+    discharge_data = discharge_data.astype({"feature_id": int, "discharge": float})
+    
     formatted_datetime = specific_date.strftime("%Y%m%d") if date_type == "date" else specific_date.strftime("%Y%m%d%H%M%S")
     finalHANDdischarge_dir = os.path.join(data_dir, f"USGS_{formatted_datetime}_{huc}.csv")
     discharge_data.to_csv(finalHANDdischarge_dir, index=False)
@@ -87,6 +88,9 @@ def getUSGSsitedata(huc=None, start_date = None, end_date= None, usgs_sites=None
     assign a feature_id corresponding to the usgs_sites and save in the data/inputs as required by the FIMserv.
     """
     code_dir, data_dir, output_dir = setup_directories()
+    HUC_dir = os.path.join(output_dir, f"flood_{huc}")
+    featureID_dir = os.path.join(HUC_dir, f"feature_IDs.csv")
+    
     def process_value_times(huc_key, value_times_list, allow_cleanup=False):
         site_data = GetUSGSIDandCorrFID(huc_key)
         usgs_ids = site_data["USGS gauge station ID"].tolist()