Performance_metrics_calculation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performance Metrics Calculation for NWM Data\n",
    "This script calculates several performance metrics for the NWM streamflow data comparing it to the USGS streamflow data. \n",
    "# Input:\n",
    "- File with list of COMIDs where the metrics are to be calculated\n",
    "- NWM data (streamflow time  series) for the COMIDs\n",
    "- USGS data (streamflow time series) for the station on/near the COMIDs. \n",
    "- File with the annual maximum streamflow value (for thresholds) at the locations\n",
    "\n",
    "# Output:\n",
    "- Excel file with performance metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from concurrent.futures import ThreadPoolExecutor, as_completed\n",
    "import cProfile\n",
    "import pstats\n",
    "\n",
    "# Function to read NWM and USGS data\n",
    "def read_data(comid, gauge_name, year, start_date, end_date):\n",
    "    nwm_discharge_path = f'{year}/discharge_{comid}_0{gauge_name}.csv'\n",
    "    usgs_data_path = f'usgs_data_1979_2023_parquet/discharge_0{gauge_name}.parquet'\n",
    "    \n",
    "    # Read NWM discharge data\n",
    "    nwm_discharge = pd.read_csv(nwm_discharge_path)\n",
    "    nwm_discharge['Time'] = nwm_discharge['time']\n",
    "    nwm_discharge = nwm_discharge[(nwm_discharge['Time'] >= start_date) & (nwm_discharge['Time'] <= end_date)]\n",
    "    nwm_discharge['Discharge_NWM'] = nwm_discharge['streamflow']\n",
    "\n",
    "    # Read USGS data\n",
    "    usgs_data = pd.read_parquet(usgs_data_path)\n",
    "    \n",
    "    return nwm_discharge, usgs_data\n",
    "\n",
    "\n",
    "# Function to filter USGS data based on the threshold value\n",
    "def filter_and_save_data(year, usgs_data, comid, gauge_name, threshold, start_date, end_date):\n",
    "    usgs_data = usgs_data.copy()\n",
    "\n",
    "    if usgs_data.empty:\n",
    "        return pd.DataFrame()\n",
    "\n",
    "    usgs_data['Time'] = pd.to_datetime(usgs_data['Time'], format='%Y-%m-%d %H:%M:%S', errors='coerce')\n",
    "    # IF USGS data is downloaded from the NWIS site no time conversion is required. Here we are using USGS data downloaded from website which is in CDT, so conversion is required\n",
    "    usgs_data['Time'] = usgs_data['Time'].dt.tz_localize('US/Central', ambiguous='NaT')\n",
    "    usgs_data['Time'] = usgs_data['Time'].dt.tz_convert('UTC')\n",
    "    usgs_data = usgs_data[usgs_data['Time'].dt.minute == 0]\n",
    "    usgs_data = usgs_data[(usgs_data['Time'] >= start_date) & (usgs_data['Time'] <= end_date)]\n",
    "    usgs_data['Discharge'] = pd.to_numeric(usgs_data['Discharge'], errors='coerce') / 35.3147  # Convert from CFS to CMS\n",
    "    usgs_filtered = usgs_data[usgs_data['Discharge'] >= threshold]\n",
    "\n",
    "    return usgs_filtered\n",
    "\n",
    "\n",
    "# Function to calculate Kling-Gupta Efficiency (KGE)\n",
    "def calculate_kge(observed, modeled):\n",
    "    mean_observed = np.mean(observed)\n",
    "    std_observed = np.std(observed)\n",
    "    mean_modeled = np.mean(modeled)\n",
    "    std_modeled = np.std(modeled)\n",
    "    if std_observed == 0 or std_modeled == 0:\n",
    "            return np.nan\n",
    "    r = np.corrcoef(observed, modeled)[0, 1]\n",
    "    alpha = np.std(modeled) / np.std(observed)\n",
    "    beta = np.mean(modeled) / np.mean(observed)\n",
    "    kge = 1 - np.sqrt((r - 1) ** 2 + (alpha - 1) ** 2 + (beta - 1) ** 2)\n",
    "    return kge\n",
    "\n",
    "# Function to calculate Nash-Sutcliffe Efficiency (NSE)\n",
    "def calculate_nse(observed, simulated):\n",
    "    mean_observed = np.mean(observed)\n",
    "    denominator = np.sum((observed - mean_observed)**2)\n",
    "    \n",
    "    if denominator == 0:\n",
    "        return np.nan  \n",
    "    \n",
    "    try:\n",
    "        nse = 1 - (np.sum((observed - simulated)**2) / denominator)\n",
    "    except:\n",
    "        return np.nan\n",
    "    return nse\n",
    "\n",
    "# Function to calculate Normalized Nash-Sutcliffe Efficiency (NNSE)\n",
    "def calculate_nnse(nse, observed):\n",
    "    mean_observed = np.mean(observed)\n",
    "    nnse = 1/(2-nse)\n",
    "    return nnse\n",
    "\n",
    "\n",
    "# Function to process each row\n",
    "def process_row(row, year, start_date, end_date, read_data, filter_and_save_data):\n",
    "    comid, gauge_name, threshold = row['COMID'], row['station_no'], row['Q_cms']\n",
    "    try:\n",
    "        # Read NWM and USGS data\n",
    "        nwm_discharge, usgs_data = read_data(comid, gauge_name, year, start_date, end_date)\n",
    "        usgs_filtered = filter_and_save_data(year, usgs_data, comid, gauge_name, threshold, start_date, end_date)\n",
    "\n",
    "        if usgs_filtered.empty:\n",
    "            return None\n",
    "\n",
    "        # Convert the time to UTC and remove the timezone information (convert to timezone-naive datetime)\n",
    "        nwm_discharge['Time'] = pd.to_datetime(nwm_discharge['Time'], utc=True).dt.tz_localize(None).astype('datetime64[ns]')\n",
    "        usgs_filtered['Time'] = pd.to_datetime(usgs_filtered['Time'], utc=True).dt.tz_localize(None).astype('datetime64[ns]')\n",
    "\n",
    "        # Using merge_asof to match nearest timestamps between NWM and USGS data\n",
    "        merged_df = pd.merge_asof(\n",
    "            usgs_filtered.sort_values('Time'),  # USGS data (sorted by time)\n",
    "            nwm_discharge.sort_values('Time'),  # NWM data (sorted by time)\n",
    "            on='Time',  # The time column to merge on\n",
    "            direction='nearest',  # Match the nearest timestamp\n",
    "            tolerance=pd.Timedelta('15min')  # Set a tolerance (considered 15 minutes here)\n",
    "        )\n",
    "\n",
    "        \n",
    "        merged_df = merged_df.dropna()\n",
    "\n",
    "        # Creating the final dataframe for comparison\n",
    "        new_df_final = pd.DataFrame({'Date': merged_df['Time']})\n",
    "        new_df_final['Observed'] = merged_df['Discharge']  # USGS discharge\n",
    "\n",
    "        new_df_final['NWM'] = merged_df['Discharge_NWM']  # NWM discharge\n",
    "        new_df_final['Observed - Model'] = new_df_final['Observed'] - new_df_final['NWM']\n",
    "\n",
    "        # Save final data to a CSV file\n",
    "        output_folder_path = f'result_{year}'\n",
    "        if not os.path.exists(output_folder_path):\n",
    "            os.makedirs(output_folder_path)\n",
    "        new_df_final.to_csv(f'{output_folder_path}/{gauge_name}.csv')\n",
    "\n",
    "        # Calculate performance metrics\n",
    "        observed = new_df_final['Observed']\n",
    "        modeled = new_df_final['NWM']\n",
    "\n",
    "        # RMSE\n",
    "        rmse = np.sqrt(np.mean((modeled - observed) ** 2))\n",
    "\n",
    "        # KGE\n",
    "        kge = calculate_kge(observed, modeled)\n",
    "\n",
    "        # PBIAS (Percent Bias)\n",
    "        pbias = 100 * (np.sum(modeled - observed) / np.sum(observed)) if np.sum(observed) != 0 else np.nan\n",
    "\n",
    "        # Normalized Error\n",
    "        # normalized_error = rmse / np.mean(observed) if np.mean(observed) != 0 else np.nan \n",
    "\n",
    "        nse = calculate_nse(observed, modeled)\n",
    "        # Calculate NNSE\n",
    "        nnse = calculate_nnse(nse, observed)\n",
    "\n",
    "        # Save results to the original row\n",
    "        result = row.copy()\n",
    "        result['RMSE'] = rmse\n",
    "        result['KGE'] = kge\n",
    "        result['PBIAS'] = pbias\n",
    "        result['Discharge'] = np.mean(observed)\n",
    "        # result['Normalized Error'] = normalized_error\n",
    "        result['NSE'] = nse\n",
    "        result['NNSE'] = nnse\n",
    "\n",
    "        return result\n",
    "\n",
    "    except Exception as e:\n",
    "        print(f\"An error occurred for COMID {comid}, Gauge {gauge_name}: {e}\")\n",
    "        return None\n",
    "\n",
    "\n",
    "# Main processing function with parallel execution using ThreadPoolExecutor\n",
    "def main_parallel_processing(comid_stn, year, start_date, end_date):\n",
    "    results = []\n",
    "    \n",
    "    # Use ThreadPoolExecutor\n",
    "    with ThreadPoolExecutor(max_workers=4) as executor:\n",
    "        future_to_row = {executor.submit(process_row, row, year, start_date, end_date, read_data, filter_and_save_data): row for idx, row in comid_stn.iterrows()}\n",
    "        \n",
    "        for future in as_completed(future_to_row):\n",
    "            result = future.result()\n",
    "            if result is not None:\n",
    "                results.append(result)\n",
    "\n",
    "    return pd.DataFrame(results) if results else pd.DataFrame()\n",
    "\n",
    "\n",
    "# Main processing function without parallelism (sequential)\n",
    "def main_sequential_processing(comid_stn, year, start_date, end_date):\n",
    "    results = []\n",
    "    \n",
    "    for idx, row in comid_stn.iterrows():\n",
    "        result = process_row(row, year, start_date, end_date, read_data, filter_and_save_data)\n",
    "        if result is not None:\n",
    "            results.append(result)\n",
    "\n",
    "    return pd.DataFrame(results) if results else pd.DataFrame()\n",
    "\n",
    "\n",
    "# Define the main function\n",
    "def main(parallel=True):\n",
    "    year = '1979'   # All data starting from 1979 till 2023 are stored in the folder 1979. This is for retrospective 3.0\n",
    "    start_date = \"1979-02-01\"\n",
    "    end_date = \"2023-01-31\"\n",
    "    comid_stn = pd.read_excel('usgs_with_comid.xlsx')  # This is the file that contains the information of the COMIDs that we want to perform the analysis for \n",
    "    rpq_filtered = pd.read_csv('consolidated_annual_min_retro_3.csv')  # This file contains the annual maximum discharge for all the COMIDs sorted in ascending order\n",
    "    comid_stn = pd.merge(comid_stn, rpq_filtered, on='COMID', how='inner')\n",
    "    # comid_stn = comid_stn[0:10]  # Limiting for testing\n",
    "\n",
    "    if parallel:\n",
    "        # Run parallel processing\n",
    "        final_results = main_parallel_processing(comid_stn, year, start_date, end_date)\n",
    "    else:\n",
    "        # Run sequential processing\n",
    "        final_results = main_sequential_processing(comid_stn, year, start_date, end_date)\n",
    "\n",
    "    # Save final results\n",
    "    if not final_results.empty:\n",
    "        output_folder_path = f'result_{year}'\n",
    "        if not os.path.exists(output_folder_path):\n",
    "            os.makedirs(output_folder_path)\n",
    "        final_results.to_csv(f'error_{year}_nwm_3_0_with_metrics.csv')\n",
    "\n",
    "\n",
    "# Profile the main function\n",
    "if __name__ == \"__main__\":\n",
    "    # Choose whether to run parallel or sequential processing\n",
    "    cProfile.run('main(parallel=False)', 'performance_profile')  # If you have a large set of dataset, use parallel or sequential processing as per your choice\n",
    "    # Analyze profiling stats\n",
    "    p = pstats.Stats('performance_profile')\n",
    "    p.strip_dirs().sort_stats('cumulative').print_stats(20)  \n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.18 ('assimilation')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.18"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "4bc09726d3eb0feb4d73d3636dad7388c6ae5ba645f208395a639ba59e5d5f59"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}