Cam 968 snow notebook 03 - Spatial Snow Analysis

examples/notebooks/snow/02_point_scale_comparison_ParFlow.ipynb

@@ -112,14 +112,14 @@
    "outputs": [],
    "source": [
     "# Observations \n",
-    "obs_df = pd.read_csv(os.path.join('obs_outputs', 'df_East-Taylor_14020001_SNOTEL.csv'))\n",
+    "obs_df = pd.read_csv(os.path.join('obs_outputs_PF', 'df_East-Taylor_14020001_SNOTEL.csv'))\n",
     "obs_df = obs_df.set_index(\"date\")\n",
     "\n",
     "# Observations metadata \n",
-    "metadata_df = pd.read_csv(os.path.join('obs_outputs', 'df_East-Taylor_14020001_SNOTEL_metadata.csv'))\n",
+    "metadata_df = pd.read_csv(os.path.join('obs_outputs_PF', 'df_East-Taylor_14020001_SNOTEL_metadata.csv'))\n",
     "\n",
     "# Model outputs\n",
-    "mod_df = pd.read_csv(os.path.join('mod_outputs', 'df_East-Taylor_14020001_PFCONUS1.csv'))\n",
+    "mod_df = pd.read_csv(os.path.join('mod_outputs_PF', 'df_East-Taylor_14020001_PFCONUS1.csv'))\n",
     "mod_df = mod_df.set_index(\"date\")"
    ]
   },

examples/notebooks/snow/nwm_env.yml

@@ -37,4 +37,6 @@ dependencies:
       - planetary-computer
       - python-dotenv
       - requests-cache
-      - hf_hydrodata
+      - hf_hydrodata
+      - subsettools
+      - jupyter_bokeh

examples/notebooks/snow/utils/nwm_utils.py

@@ -12,6 +12,7 @@
 import holoviews as hv
 import hvplot.xarray
 from holoviews import opts
+import xarray as xr
 import geopandas as gpd
 import matplotlib.pyplot as plt
 import matplotlib.dates as mdates
@@ -444,10 +445,8 @@ def compute_stats(df, ts1, ts2):
 
     return stats_table
 
-def compute_stats_period(
-    df, ts_obs, ts_mod, months
-):
-    df_sub = df[df.index.month.isin(months)]
+def compute_stats_period(df, ts_obs, ts_mod, months):
+
     """
     Create a wrapper for compute_stats to filter dataframe by months before computing stats.
     For example, to compute stats for melt season (April to July), use months=[4,5,6,7].
@@ -461,26 +460,27 @@ def compute_stats_period(
     Returns:
     - DataFrame with computed statistics.
     """
-
+    df_sub = df[df.index.month.isin(months)]
+
     return compute_stats(df_sub, ts_obs, ts_mod)
 
 
-def comparison_plots(df, ts1, ts2):
+def comparison_plots(df, ts_obs, ts_mod):
     '''
     Create a set of comparison plots (timeseries overlay and scatter plot with 1:1 line)
 
     Parameters:
     df: dataframe with combined observed and modeled timeseries for each site  
-    ts1 (str): column heading for observed timeseries in df
-    ts2 (str): column heading for modeled timeseries in df
+    ts_obs (str): column heading for observed timeseries in df
+    ts_mod (str): column heading for modeled timeseries in df
     '''
 
     df = df.copy()
     df.index.name = "date" # change the index name to "Date" for better hover tooltip display
 
     # Timeseries plot (Overlay)
     observed_plot = df.hvplot.line(
-        y=ts1,
+        y=ts_obs,
         ylabel='Snow Water Equivalent (mm)',
         xlabel='',
         label='Observed SWE',
@@ -491,7 +491,7 @@ def comparison_plots(df, ts1, ts2):
     )
 
     modeled_plot = df.hvplot.line(
-    y=ts2,
+    y=ts_mod,
     ylabel='Snow Water Equivalent (mm)',
     xlabel='',
     label='Modeled SWE',
@@ -509,8 +509,8 @@ def comparison_plots(df, ts1, ts2):
 
     # Scatter plot
     scatter_plot = df.hvplot.scatter(
-        x=ts1,
-        y=ts2,
+        x=ts_obs,
+        y=ts_mod,
         xlabel='Observed SWE (mm)',
         ylabel='Modeled SWE (mm)',
         color='black',
@@ -521,7 +521,7 @@ def comparison_plots(df, ts1, ts2):
     )
 
     # Add 1:1 line (perfect match line)
-    swe_max = max(df[ts1].max(), df[ts2].max())
+    swe_max = max(df[ts_obs].max(), df[ts_mod].max())
     one_to_one_line = hv.Curve(([0, swe_max], [0, swe_max])).opts(
         color='gray',
         line_dash='solid',

src/cssi_evaluation/snow/utils.py

@@ -3,6 +3,7 @@
 """
 
 import pandas as pd
+import numpy as np
 from typing import Any, Union
 
 
@@ -338,3 +339,58 @@ def compute_melt_period_statistics(
 
     return pd.DataFrame(result)
 
+def compute_snow_timing_grid(grid_swe_da, water_year, threshold=1.0, smooth_window=5):
+    """
+    Compute peak SWE timing, melt-out timing, and snow duration
+    for a single water year.
+
+    Parameters
+    ----------
+    grid_swe_da : xarray.DataArray
+        SWE with dims (time, ny, nx) and coordinate 'water_year'
+    water_year : int
+        Water year to analyze (e.g., 2004)
+    threshold : float
+        SWE threshold for defining snow-free (default = 1.0 mm)
+    smooth_window : int
+        Rolling window size for smoothing (default = 5 days)
+
+    Returns
+    -------
+    peak_wy_day : DataArray (ny, nx)
+    melt_wy_day : DataArray (ny, nx)
+    snow_duration : DataArray (ny, nx)
+    """
+
+    # Select only the specified water year
+    swe_wy = grid_swe_da.sel(time=grid_swe_da.water_year == water_year)
+
+    if swe_wy.time.size == 0:
+        raise ValueError(f"No data found for water year {water_year}")
+
+    # Peak SWE date
+    peak_date = swe_wy.idxmax(dim="time")
+    wy_start = swe_wy.time.values[0]
+    print(f"Water year {water_year} starts on {pd.to_datetime(wy_start)}")
+    # Store the peak date as a water-year day number for easier comparison across grid cells (The number of days after October 1 when peak SWE occurs)
+    peak_wy_day = (peak_date - wy_start) / np.timedelta64(1, "D")
+
+    # Smooth SWE 
+    swe_smooth = swe_wy.rolling(time=smooth_window, center=True).mean()
+
+    # Define snow-free condition; where the 5-day rolling mean SWE falls below the threshold to identify snow-free conditions
+    snow_free = swe_smooth <= threshold
+
+    # Mask snow free and make sure it's the after peak (i.e., not considering the snow-free period in the fall)
+    after_peak = swe_wy.time >= peak_date
+    snow_free_after_peak = snow_free.where(after_peak)
+
+    # Melt-out date
+    melt_date = snow_free_after_peak.idxmax(dim="time")
+    melt_wy_day = (melt_date - wy_start) / np.timedelta64(1, "D")
+
+    # Snow duration 
+    snow_duration = melt_wy_day - peak_wy_day
+
+    return peak_wy_day, melt_wy_day, snow_duration
+