diff --git a/examples/parflow/parflow_streamflow_point_scale_evaluation.ipynb b/examples/parflow/parflow_streamflow_point_scale_evaluation.ipynb new file mode 100644 index 0000000..bd7e8c2 --- /dev/null +++ b/examples/parflow/parflow_streamflow_point_scale_evaluation.ipynb @@ -0,0 +1,722 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Use HydroData to Retrieve Modeled and Observed Streamflow Data for a Watershed of Interest with ParFlow-CONUS Outputs vs USGS Observed Streamflow - Full Evaluation Workflow\n", + "Authors: Irene Garousi-Nejad (igarousi@cuahsi.org), Danielle Tijerina-Kreuzer (dtijerina@cuahsi.org) \n", + "Last updated: March 2026" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Introduction: \n", + "This notebook demonstrates how to perform a point-scale analysis comparing modeled and observed streamflow at selected USGS stream gages. We focus on analyzing model performance for **site-by-site streamflow comparisons**. " + ] + }, + { + "cell_type": "markdown", + "id": "e86aae63", + "metadata": {}, + "source": [ + "## 1. Setup" + ] + }, + { + "cell_type": "markdown", + "id": "e88ed1ef", + "metadata": {}, + "source": [ + "### 1a. Python Environment \n", + "\n", + "Ensure that the `cssi_evaluation` conda environment is selected as your Jupyter kernel. This environment should already be created if you followed the instructions under section \"Set up the recommended local environment\" in the [GettingStarted.md](https://github.com/hydroframe/cssi_evaluation/blob/main/GettingStarted.md) file." + ] + }, + { + "cell_type": "markdown", + "id": "c0f30927", + "metadata": {}, + "source": [ + "Import the libraries needed to run this notebook:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0dfc3fb3", + "metadata": {}, + "outputs": [], + "source": [ + "import datetime\n", + "from pathlib import Path\n", + "import pandas as pd\n", + "import numpy as np\n", + "import hf_hydrodata as hf\n", + "import parflow as pf\n", + "import plotly.express as px\n", + "import plotly.io as pio\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.dates as mdates\n", + "\n", + "from cssi_evaluation.utils import plot_utils\n", + "from cssi_evaluation.utils import evaluation_utils\n", + "from cssi_evaluation.variables import streamflow_utils\n", + "\n", + "%load_ext autoreload\n", + "%autoreload 2\n" + ] + }, + { + "cell_type": "markdown", + "id": "8e058228", + "metadata": {}, + "source": [ + "### 1b. Register Pin and Access HydroData\n", + "\n", + "To access the HydroData catalog you will need to sign up for a [HydroFrame account](https://hydrogen.princeton.edu/signup) (do this only once), [create a 4-digit PIN](https://hydrogen.princeton.edu/pin), and register your PIN in order to have access to the HydroData datasets (you will do this in the next code cell below). To note, you PIN will expire after 7 days and will need to recreate it after that time. The [Getting Started page](https://hf-hydrodata.readthedocs.io/en/latest/getting_started.html) has full details on how to sign up for an account and set up a PIN." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a365996d", + "metadata": {}, + "outputs": [], + "source": [ + "# Once you create a PIN, you'll need to register your PIN.\n", + "hf.register_api_pin(\"your_email\", \"your_api_pin\")" + ] + }, + { + "cell_type": "markdown", + "id": "b8620cfc", + "metadata": {}, + "source": [ + "## 2. Define geographic domain and date range\n", + "\n", + "### 2a. Define date range and (optional) file directories" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e50dc99d", + "metadata": {}, + "outputs": [], + "source": [ + "# Start and end times of a water year (to note, these dates were chosen to align with the PFCONUS1 early 2000s runs)\n", + "StartDate = '2002-10-01'\n", + "EndDate = '2003-09-30'\n", + "\n", + "domain_data_path = 'examples/parflow/domain_data/' # path to the model domain data\n", + "\n", + "# Output folders for modeled and observed data\n", + "OBS_OutputFolder = Path(\"./obs_outputs\")\n", + "MOD_OutputFolder = Path(\"./mod_outputs\")\n", + "\n", + "OBS_OutputFolder.mkdir(parents=True, exist_ok=True)\n", + "MOD_OutputFolder.mkdir(parents=True, exist_ok=True)" + ] + }, + { + "cell_type": "markdown", + "id": "31d78679", + "metadata": {}, + "source": [ + "### 2b. Define the watershed of interest\n", + "\n", + "One of the simplest ways to gather data and model output from HydroData is by specifying a [Hydrologic Unit Code](https://www.usgs.gov/national-hydrography/watershed-boundary-dataset). Before we retrieve any hydrologic information, we need to indicate a HUC8 code and use it to gather streamflow observations from USGS stream gages. \n", + "\n", + "✏️ If you have a specific HUC in mind, you can change the variable `huc_code` below." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "53c003cb", + "metadata": {}, + "outputs": [], + "source": [ + "# ✏️ Specify HUC8 ID and Name for watershed of interest\n", + "# (The Bayou De Chien HUC '08010201' has only 1 gage, which will run faster than the following Patoka-White HUC '051202')\n", + "huc_code = '051202'\n", + "print(f'HUC-8 ID: {huc_code}')\n", + "\n", + "huc_name = 'Patoka-White'\n", + "print(f'HUC-8 name: {huc_name}')" + ] + }, + { + "cell_type": "markdown", + "id": "feb58871", + "metadata": {}, + "source": [ + "## 3. Retrieve Observed USGS Streamflow Data via `hf_hydrodata`" + ] + }, + { + "cell_type": "markdown", + "id": "0f6cdb0a", + "metadata": {}, + "source": [ + "Here we are providing the code needed to access streamflow observations data from the `hf_hydrodata` Python package for our HUC. If you'd like more detail on this code and the `hf_hydrodata` package, please see the notebook `examples/collect_observations/hydrodata_collect_observations.ipynb`. \n", + "\n", + "One of the dataframes we collect is the metadata_df. This is an important addition to the observations and it is recommended to always gather and save this for the observations you are using (particularly to support reproducibility within an open-science workflow). The saved file has useful attributes like site names, first and last date of available data, lat/lon, and the query URL. \n", + "\n", + "Additionally, the metadata contains **ParFlow-CONUS1 and ParFlow-CONUS2 `i,j` indices, which indicate the exact model domain grid cell the observation aligns with**. This is a useful HydroData feature that removes the need for users to manually match station latitude/longitude coordinates to the appropriate model grid cell, as this spatial mapping is handled directly within HydroData. We will use these indices below to extract PF-CONUS1 modeled streamflow for each USGS gage in the section below. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10ae5690", + "metadata": {}, + "outputs": [], + "source": [ + "# Create dictionary of options to pass to the HydroData API for point-level streamflow data\n", + "point_options = {\n", + " \"dataset\": \"usgs_nwis\",\n", + " \"variable\": \"streamflow\",\n", + " \"temporal_resolution\": \"daily\",\n", + " \"aggregation\": \"mean\",\n", + " \"huc_id\": [huc_code], \n", + " \"grid\": \"conus1\",\n", + " \"date_start\": StartDate,\n", + " \"date_end\": EndDate\n", + "}\n", + "\n", + "# Let's request site-level metadata first. \n", + "streamflow_metadata_df = hf.get_point_metadata(point_options)\n", + "\n", + "print(f\"Number of streamflow sites: {len(streamflow_metadata_df)}\")\n", + "streamflow_metadata_df.head(5)" + ] + }, + { + "cell_type": "markdown", + "id": "068c31ab", + "metadata": {}, + "source": [ + "Plot the locations of the streamflow gages onto a map:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a16038bd", + "metadata": {}, + "outputs": [], + "source": [ + "# Let's define a helper plotting function to use throughout the notebook\n", + "def plot_sites_on_map(metadata_df, color_by=None):\n", + " \"\"\"\n", + " Plots site locations on an interactive map using Plotly.\n", + "\n", + " Parameters:\n", + " - metadata_df: DataFrame with 'latitude', 'longitude', 'site_name', 'site_id'\n", + " - color_by: column name to color by (e.g., 'nse', 'kge', 'site_type')\n", + " \"\"\"\n", + " pio.renderers.default = \"notebook\"\n", + "\n", + " center = {\n", + " \"lat\": metadata_df[\"latitude\"].mean(),\n", + " \"lon\": metadata_df[\"longitude\"].mean()\n", + " }\n", + "\n", + " fig = px.scatter_map(\n", + " metadata_df,\n", + " lat=\"latitude\",\n", + " lon=\"longitude\",\n", + " color=color_by, \n", + " hover_name=\"site_name\",\n", + " hover_data=[\"site_id\"],\n", + " center=center,\n", + " zoom=8,\n", + " height=600,\n", + " color_continuous_scale=\"RdYlGn\", # good default for performance metrics\n", + " )\n", + "\n", + " fig.update_layout(mapbox_style=\"carto-voyager\")\n", + " fig.update_traces(marker=dict(size=12))\n", + "\n", + " return fig\n", + "\n", + "plot_sites_on_map(streamflow_metadata_df)" + ] + }, + { + "cell_type": "markdown", + "id": "6e2b894b", + "metadata": {}, + "source": [ + "##### Request USGS Streamflow data through the HydroData API at each of the gages located in the given HUC and time period using the `hf.get_point_data()` function. This is based on our `point_options` dictionary above. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fccadf5e", + "metadata": {}, + "outputs": [], + "source": [ + "# Now let's request data.\n", + "streamflow_data_df = hf.get_point_data(point_options)\n", + "print(f\"Streamflow data shape: {streamflow_data_df.shape}\")" + ] + }, + { + "cell_type": "markdown", + "id": "c89b9f6e", + "metadata": {}, + "source": [ + "Clean up any sites that have no data for both `streamflow_data_df` and `streamflow_metadata_df`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "32cee48a", + "metadata": {}, + "outputs": [], + "source": [ + "# Identify columns with NO NaNs\n", + "valid_sites = streamflow_data_df.columns[streamflow_data_df.notna().all()]\n", + "\n", + "# Subset the data dataframe\n", + "streamflow_data_df = streamflow_data_df[valid_sites]\n", + "\n", + "# Subset the metadata dataframe\n", + "streamflow_metadata_df = streamflow_metadata_df[\n", + " streamflow_metadata_df[\"site_id\"].isin(valid_sites)\n", + "]\n", + "\n", + "print(f\"Number of valid streamflow sites: {len(valid_sites)}\")\n", + "streamflow_data_df.head(5)" + ] + }, + { + "cell_type": "markdown", + "id": "c7cea09a", + "metadata": {}, + "source": [ + "## 4. Retrieve ParFlow-CONUS1 Modeled Streamflow Data" + ] + }, + { + "cell_type": "markdown", + "id": "4dd99610", + "metadata": {}, + "source": [ + "
\n", + "

📖 Note

\n", + "

\n", + "As a project, we have pre-computed the ParFlow grid cell mappings for every stream gage, groundwater well, SNOTEL station, SCAN station, and AmeriFlux site. Those are stored in the variables conus1_i and conus1_j for the ParFlow-CONUS1 grid and conus2_i and conus2_j for the ParFlow-CONUS2 grid. In this example, we will use those grid indices directly. Please see the parflow_grid_mapping.ipynb notebook in this same directory for more details on how we constructed those mappings.\n", + "

\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "998994b5", + "metadata": {}, + "source": [ + "Now that we know where our sites are located, we're going to extract only those grid cells. The ParFlow-CONUS1 simulation results are a gridded dataset, so we'll use the `get_gridded_data` function. This function takes inputs similar to the `get_point_data` function; sometimes fewer parameters are needed. We'll request from the [ParFlow-CONUS1 modern simulations](https://hf-hydrodata.readthedocs.io/en/latest/gen_conus1_baseline_mod.html) dataset.\n", + "- `dataset=\"conus1_baseline_mod\"`\n", + "- `variable=\"streamflow\"`\n", + "- `temporal_resolution=\"daily\"`\n", + "\n", + "**Note**: The `get_gridded_data` function follows NumPy indexing convention. The end date will be exclusive of the request. We are working on updating the point module to be consistent with this convention. For now, set the start and end dates separately within each query.\n", + "\n", + "Because there are only specific locations we want to extract, we'll use the `grid_point` parameter to extract the grid cell associated with each site. As a reminder, the [Python API reference](https://hf-hydrodata.readthedocs.io/en/latest/api_reference.html) has complete information on available parameter options." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b0df123a", + "metadata": {}, + "outputs": [], + "source": [ + "date_start = datetime.datetime.strptime(StartDate, '%Y-%m-%d')\n", + "date_end = datetime.datetime.strptime(EndDate, '%Y-%m-%d')\n", + "\n", + "# Create DataFrame for storing ParFlow results\n", + "timesteps = np.arange(\n", + " date_start,\n", + " date_end + datetime.timedelta(days=1),\n", + " datetime.timedelta(days=1)).astype(\"datetime64[D]\")\n", + "\n", + "model_df = pd.DataFrame(columns=[\"date\"], data=timesteps)\n", + "model_df.head(5)" + ] + }, + { + "cell_type": "markdown", + "id": "c9c9542b", + "metadata": {}, + "source": [ + "### Compute simulated streamflow from PF-CONUS1 pressure outputs " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "310c8fc1", + "metadata": {}, + "outputs": [], + "source": [ + "# Gather static fields needed for streamflow calculation\n", + "mannings = 5.0e-5\n", + "dx, dy = 1000, 1000\n", + "\n", + "slope_x = hf.get_gridded_data({\n", + " \"dataset\": \"conus1_domain\",\n", + " \"variable\": \"slope_x\",\n", + " })\n", + "slope_y = hf.get_gridded_data({\n", + " \"dataset\": \"conus1_domain\",\n", + " \"variable\": \"slope_y\",\n", + " })" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9ebe3c48", + "metadata": {}, + "outputs": [], + "source": [ + "# Loop over each station in obs_metadata_df\n", + "for idx, row in streamflow_metadata_df.iterrows():\n", + " site_id = row[\"site_id\"] # original USGS site_id\n", + " conus_i = int(row[\"conus1_i\"])\n", + " conus_j = int(row[\"conus1_j\"])\n", + " print(f\"Processing site {site_id} at PF-CONUS1 grid point ({conus_i}, {conus_j})\")\n", + " \n", + " # Build options dict for this station\n", + " parflow_options = {\n", + " \"dataset\": \"conus1_baseline_mod\",\n", + " \"variable\": \"pressure_head\",\n", + " \"temporal_resolution\": \"daily\",\n", + " \"start_time\": date_start,\n", + " \"end_time\": date_end + datetime.timedelta(days=1), ### NOTE: the gridded function has exclusive end date\n", + " \"grid_point\": [conus_i, conus_j]\n", + " }\n", + " print(f\"Requesting pressure data for site {site_id}\")\n", + " pressure = hf.get_gridded_data(parflow_options)\n", + " \n", + " # Loop through timesteps to calculate streamflow at each timestep\n", + " print(f\"Calculating streamflow for site {site_id}\")\n", + " streamflow = np.empty(len(pressure))\n", + " for t in range(len(pressure)):\n", + " streamflow[t] = pf.tools.hydrology.calculate_overland_flow(pressure[t].reshape(5, 1, 1), \n", + " slope_x[conus_j, conus_i].reshape(1, 1),\n", + " slope_y[conus_j, conus_i].reshape(1, 1),\n", + " mannings, dx, dy, mask=np.ones((5, 1, 1)), \n", + " flow_method=\"OverlandFlow\") / 3600 # convert from m^3/h to cms\n", + " \n", + " # Fill column in model_df\n", + " model_df[site_id] = np.squeeze(np.array(streamflow))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7551a523", + "metadata": {}, + "outputs": [], + "source": [ + "model_df.head(5)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "263a9ca5", + "metadata": {}, + "outputs": [], + "source": [ + "streamflow_data_df.head(5)" + ] + }, + { + "cell_type": "markdown", + "id": "c8af8c51", + "metadata": {}, + "source": [ + "Great! Now we have our ParFlow results in the same format as our observations data.\n", + "\n", + "*Optional*: Save observations data and matching ParFlow outputs to .csv files" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a6d957d2", + "metadata": {}, + "outputs": [], + "source": [ + "# ensure 'date' columns are datetime \n", + "streamflow_data_df['date'] = pd.to_datetime(streamflow_data_df['date'])\n", + "model_df['date'] = pd.to_datetime(model_df['date'])\n", + "\n", + "# Save CSVs to the output folders\n", + "streamflow_data_df.to_csv(f'{OBS_OutputFolder}/streamflow_obs_data.csv', index=False)\n", + "streamflow_metadata_df.to_csv(f'{OBS_OutputFolder}/streamflow_obs_metadata.csv', index=False)\n", + "\n", + "model_df.to_csv(f'{MOD_OutputFolder}/streamflow_parflow_model_data.csv', index=False)" + ] + }, + { + "cell_type": "markdown", + "id": "a59c41c2", + "metadata": {}, + "source": [ + "## 5. Evaluate ParFlow-CONUS1 outputs against observations" + ] + }, + { + "cell_type": "markdown", + "id": "933fd0a9", + "metadata": {}, + "source": [ + "### 5a. Create hydrograph plots and compute metrics for each stream gage \n", + "#### Plot the observed and modeled values together over time. This helps identify:\n", + " - Periods of systematic bias\n", + " - Timing differences in peaks and lows\n", + " - General agreement in trends" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1929bd27", + "metadata": {}, + "outputs": [], + "source": [ + "# The resulting plots will be saved to the current working directory.\n", + "plot_utils.plot_time_series(\n", + " streamflow_data_df,\n", + " model_df,\n", + " streamflow_metadata_df,\n", + " variable=\"streamflow\")" + ] + }, + { + "cell_type": "markdown", + "id": "a68b9667", + "metadata": {}, + "source": [ + ">Optional: Copy the following code to a Python code block to view saved plots in the notebook. \n", + "To note, this will display all of the saved plots, so be mindful if you have many sites. \n", + "```from IPython.display import display, Image\n", + "\n", + "for site_id in list(streamflow_metadata_df[\"site_id\"]):\n", + " img = f\"streamflow_{site_id}.png\"\n", + " display(Image(filename=img))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1f44982d", + "metadata": {}, + "outputs": [], + "source": [ + "# look at one of the saved hydrographs in the notebook\n", + "from IPython.display import display, Image\n", + "\n", + "img = \"streamflow_03354000.png\"\n", + "display(Image(filename=img))" + ] + }, + { + "cell_type": "markdown", + "id": "ea4971e5", + "metadata": {}, + "source": [ + "#### Calculate site-level statistical metrics of comparison" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cfac2e30", + "metadata": {}, + "outputs": [], + "source": [ + "metrics_df = evaluation_utils.calculate_metrics(streamflow_data_df, model_df, streamflow_metadata_df)\n", + "metrics_df.head(5)" + ] + }, + { + "cell_type": "markdown", + "id": "fa36a124", + "metadata": {}, + "source": [ + "### 5b. Multi-Site Evaluation \n", + "We begin by examining model performance across all gages in the watershed." + ] + }, + { + "cell_type": "markdown", + "id": "dcb98bc1", + "metadata": {}, + "source": [ + "#### Scatter Plot for all gages with 1:1 Line:\n", + "Plots each modeled value against its corresponding observed value. This highlights:\n", + " - Accuracy across the full range of values\n", + " - Over- or under-prediction patterns\n", + " - Outliers or extreme events" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "48a4c1bc", + "metadata": {}, + "outputs": [], + "source": [ + "# Scatter plot for all sites comparing values\n", + "plot_utils.plot_compare_scatter(streamflow_data_df, model_df, variable=\"streamflow\", log_scale=False)\n" + ] + }, + { + "cell_type": "markdown", + "id": "15cab20f", + "metadata": {}, + "source": [ + "#### Visualize Metric Distribution \n", + "Histograms of individual performance metrics provide a high-level evaluation of model performance across the domain." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db73f69a", + "metadata": {}, + "outputs": [], + "source": [ + "metrics_df[[\"nse\", \"kge\", \"spearman_rho\", \"bias\"]].hist(bins=20)\n", + "plt.suptitle(f'Distribution of Selected Performance Metrics, n={len(metrics_df)}')" + ] + }, + { + "cell_type": "markdown", + "id": "f94591f5", + "metadata": {}, + "source": [ + "#### Spatial mapping of gage with colored by performance metric " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8b1cb373", + "metadata": {}, + "outputs": [], + "source": [ + "plot_sites_on_map(metrics_df, color_by=\"bias\")" + ] + }, + { + "cell_type": "markdown", + "id": "764e1fd3", + "metadata": {}, + "source": [ + "#### Combination Metrics \n", + "The Condon diagram combines Spearman's Rho and Absolute Relative Bias to represent the performance of streamflow magnitude and timing. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "674a3a28", + "metadata": {}, + "outputs": [], + "source": [ + "# Map Condon diagram\n", + "plot_utils.plot_condon_diagram(metrics_df, variable=\"streamflow\")" + ] + }, + { + "cell_type": "markdown", + "id": "23d1553a", + "metadata": {}, + "source": [ + "### 5c. Single-Site Evaluation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c54992a5", + "metadata": {}, + "outputs": [], + "source": [ + "# pick a gage \n", + "single_gage = \"03354000\"\n", + "\n", + "gage_metrics = metrics_df[metrics_df[\"site_id\"] == single_gage].iloc[0]\n", + "display(gage_metrics.to_frame(name=\"value\"))" + ] + }, + { + "cell_type": "markdown", + "id": "79c4a3ea", + "metadata": {}, + "source": [ + "#### Plot streamflow diagnostic plots using the function `plot_streamflow_diagnostics()` in `streamflow_utils.py` to show modeled vs. observed hydrograph, flow duration curve, and Q-Q plot: " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2540851c", + "metadata": {}, + "outputs": [], + "source": [ + "streamflow_utils.plot_streamflow_diagnostics(\n", + " streamflow_data_df,\n", + " model_df,\n", + " single_gage,\n", + " gage_metrics\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "30f069c4", + "metadata": {}, + "source": [ + "FDC: plots streamflow against exceedance probability; compares the distribution of flows between observed and modeled data \n", + "Q-Q: compares modeled and observed flows at the same percentiles (quantiles)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "cssi_evaluation", + "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.10.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pyproject.toml b/pyproject.toml index 98b4310..40aaa1b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -40,6 +40,7 @@ notebooks = [ "plotly", "pyproj", "subsettools>=2.0.0", + "nbformat>=4.2.0", ] [dependency-groups] diff --git a/src/cssi_evaluation/utils/plot_utils.py b/src/cssi_evaluation/utils/plot_utils.py index 23a34f6..87caff4 100644 --- a/src/cssi_evaluation/utils/plot_utils.py +++ b/src/cssi_evaluation/utils/plot_utils.py @@ -408,11 +408,10 @@ def plot_condon_diagram(metrics_df, variable, output_dir=".", adaptive_xlim=True ax.scatter( df_plot["abs_rel_bias"], df_plot["spearman_rho"], - c=[CONDON_COLORS[c] for c in df_plot["condon"]], - s=8, - zorder=2, - alpha=0.5, - edgecolors="none", + c=df_plot["condon"].map(CONDON_COLORS), + s=25, + zorder=1, + alpha=0.4, ) custom = [ @@ -458,16 +457,65 @@ def plot_condon_diagram(metrics_df, variable, output_dir=".", adaptive_xlim=True # Add percentages in quadrants total_obs = df_plot.shape[0] - if total_obs > 0: - def pct(label): - return round(df_plot[df_plot["condon"] == label].shape[0] / total_obs * 100) - - ax.text(0.02 * xmax, 0.90, f"{pct('Low bias, good shape')}%", weight="bold", fontsize=10) - ax.text(0.90 * xmax, 0.90, f"{pct('High bias, good shape')}%", weight="bold", fontsize=10) - ax.text(0.02 * xmax, -0.95, f"{pct('Low bias, poor shape')}%", weight="bold", fontsize=10) - ax.text(0.90 * xmax, -0.95, f"{pct('High bias, poor shape')}%", weight="bold", fontsize=10) + ax.text( + 0.1, + 0.9, + str( + round( + df_plot[df_plot["condon"] == "Low bias, good shape"].shape[0] + / total_obs + * 100 + ) + ) + + "%", + weight="bold", + fontsize=12, + ) + ax.text( + 9.3, + 0.9, + str( + round( + df_plot[df_plot["condon"] == "High bias, good shape"].shape[0] + / total_obs + * 100 + ) + ) + + "%", + weight="bold", + fontsize=12, + ) + ax.text( + 0.1, + -0.99, + str( + round( + df_plot[df_plot["condon"] == "Low bias, poor shape"].shape[0] + / total_obs + * 100 + ) + ) + + "%", + weight="bold", + fontsize=12, + ) + ax.text( + 9.3, + -0.99, + str( + round( + df_plot[df_plot["condon"] == "High bias, poor shape"].shape[0] + / total_obs + * 100 + ) + ) + + "%", + weight="bold", + fontsize=12, + ) - ax.set_title(f"{variable.capitalize()} Performance Category", fontsize=14, pad=28) + plt.title(f"{variable.capitalize()} Performance Category") + plt.savefig(f"{output_dir}/{variable}_condon_diagram.png", bbox_inches="tight", dpi=300) # Leave room at top for outside legend fig.subplots_adjust(top=0.82) diff --git a/src/cssi_evaluation/variables/streamflow_utils.py b/src/cssi_evaluation/variables/streamflow_utils.py index fece850..7b68da8 100644 --- a/src/cssi_evaluation/variables/streamflow_utils.py +++ b/src/cssi_evaluation/variables/streamflow_utils.py @@ -1,3 +1,123 @@ """ Streamflow evaluation utility functions """ + +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt + +def plot_streamflow_diagnostics( + streamflow_data_df, + model_df, + site_id, + metrics_row=None, + output_dir="." +): + """ + Create a 3-panel diagnostic plot: + 1) Hydrograph + 2) Flow Duration Curve (FDC) + 3) Q-Q plot + + Parameters + ---------- + streamflow_data_df : DataFrame + Observed data with 'date' column + model_df : DataFrame + Modeled data with 'date' column + site_id : str + Gage/site ID column name + metrics_row : Series (optional) + Row containing metrics (rho, bias, nse, condon) + output_dir : str; default="." + String path to where plots should be saved. Default is current working directory. + """ + + # ensure datetime + dates = pd.to_datetime(streamflow_data_df['date']) + + obs = streamflow_data_df[site_id].values + mod = model_df[site_id].values + + # remove any NaNs in observed and modeled dataframes + mask = ~np.isnan(obs) & ~np.isnan(mod) + obs = obs[mask] + mod = mod[mask] + dates = dates[mask] + + # create figure + fig, axes = plt.subplots(1, 3, figsize=(18, 4)) + + # ========================================================= + # 1. HYDROGRAPH + # ========================================================= + ax = axes[0] + ax.plot(dates, obs, label="Observed", linewidth=2) + ax.plot(dates, mod, label="Modeled", linewidth=2, alpha=0.8) + + ax.set_title(f"Hydrograph at gage {site_id}") + ax.set_ylabel("Streamflow") + ax.legend() + + plt.setp(ax.get_xticklabels(), rotation=45, ha="right") + + # add metrics annotation + if metrics_row is not None: + textstr = ( + f"Srho: {metrics_row['spearman_rho']:.2f}\n" + f"Bias: {metrics_row['bias']:.2f}\n" + f"NSE: {metrics_row['nse']:.2f}\n" + f"Condon: {metrics_row['condon']}" + ) + ax.text( + 0.02, 0.98, textstr, + transform=ax.transAxes, + fontsize=9, + verticalalignment="top", + bbox=dict(boxstyle="round", alpha=0.2) + ) + + # ========================================================= + # 2. FLOW DURATION CURVE (FDC) + # ========================================================= + ax = axes[1] + + # sort descending + obs_sorted = np.sort(obs)[::-1] + mod_sorted = np.sort(mod)[::-1] + + # exceedance probability + p = np.arange(1, len(obs_sorted)+1) / (len(obs_sorted)+1) + + ax.plot(p, obs_sorted, label="Observed", linewidth=2) + ax.plot(p, mod_sorted, label="Modeled", linewidth=2, alpha=0.8) + + ax.set_yscale("log") + ax.set_xlabel("Exceedance Probability") + ax.set_ylabel("Streamflow") + ax.set_title("Flow Duration Curve") + ax.legend() + + # ========================================================= + # 3. Q–Q PLOT + # ========================================================= + ax = axes[2] + + q = np.linspace(0, 1, 100) + obs_q = np.quantile(obs, q) + mod_q = np.quantile(mod, q) + + ax.scatter(obs_q, mod_q, s=20, alpha=0.7) + + # 1:1 line + min_val = min(obs_q.min(), mod_q.min()) + max_val = max(obs_q.max(), mod_q.max()) + ax.plot([min_val, max_val], [min_val, max_val], 'k--') + + ax.set_xlabel("Observed Quantiles") + ax.set_ylabel("Modeled Quantiles") + ax.set_title("Q-Q Plot") + + # ========================================================= + plt.tight_layout() + plt.savefig(f"{output_dir}/streamflow_diagnostics_3panel_{site_id}.png", bbox_inches="tight", dpi=300) \ No newline at end of file