From 39111f17a58f79cc8a9fe00ef6ce2007c5c664b9 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Wed, 18 Feb 2026 15:17:50 -0500
Subject: [PATCH 01/37] move some vector plotting functionality into
 vector_utils.py (#2321)

---
 .../invest/coastal_vulnerability/reporter.py  | 77 +++----------------
 src/natcap/invest/reports/vector_utils.py     | 57 ++++++++++++++
 2 files changed, 69 insertions(+), 65 deletions(-)
 create mode 100644 src/natcap/invest/reports/vector_utils.py

diff --git a/src/natcap/invest/coastal_vulnerability/reporter.py b/src/natcap/invest/coastal_vulnerability/reporter.py
index 2c36f8622a..3a98cbb140 100644
--- a/src/natcap/invest/coastal_vulnerability/reporter.py
+++ b/src/natcap/invest/coastal_vulnerability/reporter.py
@@ -9,7 +9,7 @@
 from natcap.invest import __version__
 from natcap.invest import gettext
 import natcap.invest.spec
-from natcap.invest.reports import jinja_env
+from natcap.invest.reports import jinja_env, vector_utils
 
 
 LOGGER = logging.getLogger(__name__)
@@ -28,61 +28,6 @@
 POINT_SIZE = 20
 MAP_WIDTH = 450  # pixels
 
-LEGEND_CONFIG = {
-    'labelFontSize': 14,
-    'titleFontSize': 14,
-    'orient': 'left',
-    'gradientLength': 120
-}
-AXIS_CONFIG = {
-    'labelFontSize': 12,
-    'titleFontSize': 12,
-}
-
-
-def _get_geojson_bbox(geodataframe):
-    """Get the bounding box of a GeoDataFrame as a GeoJSON feature.
-
-    Also calculate its aspect ratio. These are useful for cropping
-    other layers in altair plots.
-
-    Args:
-        geodataframe (geopandas.GeoDataFrame):
-    Returns:
-        tuple: A 2-tuple containing:
-            - extent_feature (dict): A GeoJSON feature representing the bounding
-            box of the input GeoDataFrame.
-            - xy_ratio (float): The aspect ratio of the bounding box
-            (width/height).
-
-    """
-    xmin, ymin, xmax, ymax = geodataframe.total_bounds
-    xy_ratio = (xmax - xmin) / (ymax - ymin)
-    extent_feature = {
-        "type": "Feature",
-        "geometry": {"type": "Polygon",
-                     "coordinates": [[
-                         [xmax, ymax],
-                         [xmax, ymin],
-                         [xmin, ymin],
-                         [xmin, ymax],
-                         [xmax, ymax]]]},
-        "properties": {}
-    }
-    return extent_feature, xy_ratio
-
-
-def _chart_landmass(geodataframe, clip=False, extent_feature=None):
-    landmass = altair.Chart(geodataframe).mark_geoshape(
-        clip=clip,
-        fill='lightgrey'
-    ).project(
-        type='identity',
-        reflectY=True,  # Canvas and SVG treats positive y as down
-        fit=extent_feature
-    )
-    return landmass
-
 
 def _chart_base_points(geodataframe):
     # Plot points using mark_circle instead of mark_geoshape
@@ -134,13 +79,15 @@ def concat_habitats(row):
         ]
     )
 
-    _, xy_ratio = _get_geojson_bbox(exposure_geodf)
+    _, xy_ratio = vector_utils.get_geojson_bbox(exposure_geodf)
     habitat_map = landmass_chart + habitat_points
     habitat_map = habitat_map.properties(
         width=MAP_WIDTH,
         height=MAP_WIDTH / xy_ratio,
         title=gettext('The role of habitat in reducing coastal exposure')
-    ).configure_legend(**LEGEND_CONFIG).configure_axis(**AXIS_CONFIG)
+    ).configure_legend(
+        **vector_utils.LEGEND_CONFIG
+    ).configure_axis(**vector_utils.AXIS_CONFIG)
     return habitat_map
 
 
@@ -172,8 +119,8 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
 
     landmass_geo = geopandas.read_file(
         file_registry['clipped_projected_landmass'])
-    extent_feature, xy_ratio = _get_geojson_bbox(exposure_geo)
-    landmass_chart = _chart_landmass(
+    extent_feature, xy_ratio = vector_utils.get_geojson_bbox(exposure_geo)
+    landmass_chart = vector_utils.chart_landmass(
         landmass_geo, clip=True, extent_feature=extent_feature)
     base_points = _chart_base_points(exposure_geo)
 
@@ -247,7 +194,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
         width=MAP_WIDTH,
         height=MAP_WIDTH / xy_ratio,
         title='coastal exposure'
-    ).configure_legend(**LEGEND_CONFIG)
+    ).configure_legend(**vector_utils.LEGEND_CONFIG)
     exposure_map_json = exposure_map.to_json()
     exposure_map_caption = [model_spec.get_output(
         'coastal_exposure').get_field('exposure').about]
@@ -284,7 +231,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
     ).properties(
         width=MAP_WIDTH,
         height=200
-    ).configure_axis(**AXIS_CONFIG)
+    ).configure_axis(**vector_utils.AXIS_CONFIG)
     exposure_histogram_json = exposure_histogram.to_json()
 
     base_rank_vars_chart = base_points.mark_circle(
@@ -308,7 +255,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
     rank_vars_figure = altair.vconcat(
         altair.hconcat(*rank_vars_chart_list[:n_cols]),
         altair.hconcat(*rank_vars_chart_list[n_cols:])
-    ).configure_axis(**AXIS_CONFIG)
+    ).configure_axis(**vector_utils.AXIS_CONFIG)
     rank_vars_figure_json = rank_vars_figure.to_json()
     rank_vars_figure_caption = gettext(
         """
@@ -345,7 +292,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
         )
         histograms.append(hist)
     facetted_histograms = altair.hconcat(
-        *histograms).configure_axis(**AXIS_CONFIG)
+        *histograms).configure_axis(**vector_utils.AXIS_CONFIG)
     facetted_histograms_json = facetted_histograms.to_json()
     facetted_histograms_caption = model_spec.get_output(
         'intermediate_exposure').about
@@ -380,7 +327,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
         width=MAP_WIDTH + 30,  # extra space for legend
         height=MAP_WIDTH / xy_ratio,
         title=gettext('local wind-driven waves vs. open ocean waves')
-    ).configure_legend(**LEGEND_CONFIG)
+    ).configure_legend(**vector_utils.LEGEND_CONFIG)
     wave_energy_map_json = wave_energy_map.to_json()
 
     wave_energy_map_caption = [model_spec.get_output(
diff --git a/src/natcap/invest/reports/vector_utils.py b/src/natcap/invest/reports/vector_utils.py
new file mode 100644
index 0000000000..06d32042e1
--- /dev/null
+++ b/src/natcap/invest/reports/vector_utils.py
@@ -0,0 +1,57 @@
+import altair
+
+
+LEGEND_CONFIG = {
+    'labelFontSize': 14,
+    'titleFontSize': 14,
+    'orient': 'left',
+    'gradientLength': 120
+}
+AXIS_CONFIG = {
+    'labelFontSize': 12,
+    'titleFontSize': 12,
+}
+
+
+def get_geojson_bbox(geodataframe):
+    """Get the bounding box of a GeoDataFrame as a GeoJSON feature.
+
+    Also calculate its aspect ratio. These are useful for cropping
+    other layers in altair plots.
+
+    Args:
+        geodataframe (geopandas.GeoDataFrame):
+    Returns:
+        tuple: A 2-tuple containing:
+            - extent_feature (dict): A GeoJSON feature representing the bounding
+            box of the input GeoDataFrame.
+            - xy_ratio (float): The aspect ratio of the bounding box
+            (width/height).
+
+    """
+    xmin, ymin, xmax, ymax = geodataframe.total_bounds
+    xy_ratio = (xmax - xmin) / (ymax - ymin)
+    extent_feature = {
+        "type": "Feature",
+        "geometry": {"type": "Polygon",
+                     "coordinates": [[
+                         [xmax, ymax],
+                         [xmax, ymin],
+                         [xmin, ymin],
+                         [xmin, ymax],
+                         [xmax, ymax]]]},
+        "properties": {}
+    }
+    return extent_feature, xy_ratio
+
+
+def chart_landmass(geodataframe, clip=False, extent_feature=None):
+    landmass = altair.Chart(geodataframe).mark_geoshape(
+        clip=clip,
+        fill='lightgrey'
+    ).project(
+        type='identity',
+        reflectY=True,  # Canvas and SVG treats positive y as down
+        fit=extent_feature
+    )
+    return landmass

From a92d06958dbbd988154c1045caf932b529e78ad6 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Wed, 18 Feb 2026 15:18:48 -0500
Subject: [PATCH 02/37] groundwork for a swy report (#2321)

---
 .../models/seasonal_water_yield.html          |  91 ++++++++++++
 .../invest/reports/templates/styles.html      |   1 +
 .../invest/seasonal_water_yield/reporter.py   | 138 ++++++++++++++++++
 .../seasonal_water_yield.py                   |  16 ++
 4 files changed, 246 insertions(+)
 create mode 100644 src/natcap/invest/reports/templates/models/seasonal_water_yield.html
 create mode 100644 src/natcap/invest/seasonal_water_yield/reporter.py

diff --git a/src/natcap/invest/reports/templates/models/seasonal_water_yield.html b/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
new file mode 100644
index 0000000000..43215af53d
--- /dev/null
+++ b/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
@@ -0,0 +1,91 @@
+{% extends 'base.html' %}
+
+{% block content %}
+  {{ super() }}
+  {% from 'args-table.html' import args_table %}
+  {% from 'caption.html' import caption %}
+  {% from 'content-grid.html' import content_grid %}
+  {% from 'metadata.html' import list_metadata %}
+  {% from 'raster-plot-img.html' import raster_plot_img %}
+
+  <h2 class="section-header">Results</h2>
+  {{ accordion_section(
+    'Primary Outputs',
+    content_grid([
+      (caption(raster_group_caption, pre_caption=True), 100),
+      (raster_plot_img(outputs_img_src, 'Primary Outputs'), 100),
+      (caption(outputs_caption, definition_list=True), 100)
+    ])
+  ) }}
+
+  {{ accordion_section(
+    'Results aggregated by watershed',
+    content_grid([
+      (content_grid([
+        ('<div id="qb_map"></div>', 100),
+        (caption(qb_map_caption, aggregate_map_source_list), 100)
+      ]), 50),
+      (content_grid([
+        ('<div id="vri_sum_map"></div>', 100),
+        (caption(vri_sum_map_caption, aggregate_map_source_list), 100)
+      ]), 50),
+    ])
+  ) }}
+
+  {#
+  <h2 class="section-header">Intermediate Outputs</h2>
+  {{ accordion_section(
+    'Ranked exposure variables',
+    content_grid([
+      ('<div id="rank_vars_figure"></div>', 100),
+      (caption(rank_vars_figure_caption, rank_vars_figure_source_list), 100)
+    ])
+  ) }}
+
+  {{ accordion_section(
+    'Pre-ranked variables',
+    content_grid([
+      ('<div id="facetted_histograms_figure"></div>', 100),
+      (caption(facetted_histograms_caption, facetted_histograms_source_list), 100)
+    ])
+  ) }}
+
+  {{ accordion_section(
+    'Wave Exposure',
+    content_grid([
+      ('<div id="wave_energy_map"></div>', 65),
+      (caption(wave_energy_map_caption, wave_energy_map_source_list), 35)
+    ])
+  ) }}
+  #}
+
+  <h2 class="section-header">Inputs</h2>
+  {{ accordion_section(
+    'Arguments',
+    args_table(args_dict),
+  )}}
+
+  <h2 class="section-header">Metadata</h2>
+  {{
+    accordion_section(
+      'Output Filenames and Descriptions',
+      list_metadata(model_spec_outputs),
+      expanded=False
+    )
+  }}
+
+{% endblock content %}
+
+{% from 'vegalite-plot.html' import embed_vega %}
+{% block scripts %}
+  {{ super() }}
+  <script src="https://cdn.jsdelivr.net/npm/vega@5"></script>
+  <script src="https://cdn.jsdelivr.net/npm/vega-lite@5.20.1"></script>
+  <script src="https://cdn.jsdelivr.net/npm/vega-embed@6"></script>
+  {% include 'vega-embed-js.html' %}
+  {% set chart_spec_id_list = [
+    (qb_map_json, 'qb_map'),
+    (vri_sum_map_json, 'vri_sum_map'),
+  ] %}
+  {{ embed_vega(chart_spec_id_list) }}
+{% endblock scripts %}
diff --git a/src/natcap/invest/reports/templates/styles.html b/src/natcap/invest/reports/templates/styles.html
index ae649de10b..def11c7204 100644
--- a/src/natcap/invest/reports/templates/styles.html
+++ b/src/natcap/invest/reports/templates/styles.html
@@ -11,6 +11,7 @@
     --accent-color-carbon: #bed4f0;
     --accent-color-coastal_vulnerability: #fff6d8;
     --accent-color-ndr: #c2da9d;
+    --accent-color-seasonal_water_yield: #d1d1f3;
     --accent-color-sdr: #e4b5db;
 
     --accent-color: var(--accent-color-{{ accent_color_suffix }});
diff --git a/src/natcap/invest/seasonal_water_yield/reporter.py b/src/natcap/invest/seasonal_water_yield/reporter.py
new file mode 100644
index 0000000000..cd525f628b
--- /dev/null
+++ b/src/natcap/invest/seasonal_water_yield/reporter.py
@@ -0,0 +1,138 @@
+import logging
+import os
+import time
+
+import altair
+import geopandas
+import pandas
+
+from natcap.invest import __version__
+from natcap.invest import gettext
+import natcap.invest.spec
+from natcap.invest.reports import jinja_env, raster_utils, report_constants, vector_utils
+from natcap.invest.unit_registry import u
+
+
+LOGGER = logging.getLogger(__name__)
+
+TEMPLATE = jinja_env.get_template('models/seasonal_water_yield.html')
+
+MAP_WIDTH = 450 # pixels
+
+OUTPUT_RASTER_PLOT_TUPLES = [
+    ('l_aligned', 'continuous'),
+    ('l_avail', 'continuous'),
+    ('b', 'continuous'),
+    ('b_sum', 'continuous'),
+    ('qf', 'continuous')
+]
+
+
+def _create_aggregate_map(geodataframe, extent_feature, xy_ratio, attribute,
+                          title):
+    attr_map = altair.Chart(geodataframe).mark_geoshape(
+        clip=True,
+        stroke="white",
+        strokeWidth=0.5
+    ).project(
+        type='identity',
+        reflectY=True,
+        fit=extent_feature
+    ).encode(
+        color=attribute,
+        tooltip=[altair.Tooltip(attribute, title=attribute)]
+    ).properties(
+        width=MAP_WIDTH,
+        height=MAP_WIDTH / xy_ratio,
+        title=title
+    ).configure_legend(**vector_utils.LEGEND_CONFIG)
+
+    return attr_map.to_json()
+#    qb_map = altair.Chart(aggregated_results).mark_geoshape(
+#        clip=True,
+#        stroke="white",
+#        strokeWidth=0.5
+#    ).project(
+#        type='identity',
+#        reflectY=True,
+#        fit=extent_feature
+#    ).encode(
+#        color='qb',
+#        tooltip=[altair.Tooltip("qb", title="qb")]
+#    ).properties(
+#        width=MAP_WIDTH,
+#        height=MAP_WIDTH / xy_ratio,
+#        title='Mean local recharge value within the watershed'
+#    ).configure_legend(**vector_utils.LEGEND_CONFIG)
+#    qb_map_caption = [model_spec.get_output(
+#        'aggregate_vector').get_field('qb').about]
+
+
+def report(file_registry, args_dict, model_spec, target_html_filepath):
+    """Generate an html summary of Coastal Vulnerability results.
+
+    Args:
+        file_registry (dict): The ``natcap.invest.FileRegistry.registry``
+            that was returned by ``natcap.invest.seasonal_water_yield.execute``.
+        args_dict (dict): The arguments that were passed to
+            ``natcap.invest.seasonal_water_yield.execute``.
+        model_spec (natcap.invest.spec.ModelSpec):
+            ``natcap.invest.seasonal_water_yield.MODEL_SPEC``
+        target_html_filepath (str): path to an html file generated by this
+            function.
+
+    Returns:
+        None
+    """
+
+    aggregated_results = geopandas.read_file(file_registry['aggregate_vector'])
+    extent_feature, xy_ratio = vector_utils.get_geojson_bbox(aggregated_results)
+
+    qb_map_json = _create_aggregate_map(
+        aggregated_results, extent_feature, xy_ratio, 'qb',
+        gettext('Mean local recharge value within the watershed'))
+    qb_map_caption = [
+        model_spec.get_output('aggregate_vector').get_field('qb').about,
+        gettext('Values are in millimeters, but should be interpreted as '
+                'relative values, not absolute values.')]
+
+    vri_sum_map_json = _create_aggregate_map(
+        aggregated_results, extent_feature, xy_ratio, 'vri_sum',
+        gettext('Total recharge contribution of the watershed'))
+    vri_sum_map_caption = [
+        model_spec.get_output('aggregate_vector').get_field('vri_sum').about,
+        gettext('The sum of ``Vri_[suffix].tif`` pixel values within the watershed.')]
+
+    vector_map_source_list = [model_spec.get_output('aggregate_vector').path]
+
+    output_raster_plot_configs = raster_utils.build_raster_plot_configs(
+            file_registry, OUTPUT_RASTER_PLOT_TUPLES)
+    outputs_img_src = raster_utils.plot_and_base64_encode_rasters(
+            output_raster_plot_configs)
+    output_raster_caption = raster_utils.generate_caption_from_raster_list(
+            [(id, 'output') for (id, _) in OUTPUT_RASTER_PLOT_TUPLES],
+            args_dict, file_registry, model_spec)
+
+    with open(target_html_filepath, 'w', encoding='utf-8') as target_file:
+        target_file.write(TEMPLATE.render(
+            report_script=model_spec.reporter,
+            invest_version=__version__,
+            report_filepath=target_html_filepath,
+            model_id=model_spec.model_id,
+            model_name=model_spec.model_title,
+            model_description=model_spec.about,
+            userguide_page=model_spec.userguide,
+            timestamp=time.strftime('%Y-%m-%d %H:%M'),
+            args_dict=args_dict,
+            raster_group_caption=report_constants.RASTER_GROUP_CAPTION,
+            outputs_img_src=outputs_img_src,
+            outputs_caption=output_raster_caption,
+            qb_map_json=qb_map_json,
+            qb_map_caption=qb_map_caption,
+            vri_sum_map_json=vri_sum_map_json,
+            vri_sum_map_caption=vri_sum_map_caption,
+            aggregate_map_source_list=vector_map_source_list,
+            model_spec_outputs=model_spec.outputs,
+        ))
+
+    LOGGER.info(f'Created {target_html_filepath}')
diff --git a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
index d0f1ba984f..c96452b211 100644
--- a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
+++ b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
@@ -25,10 +25,26 @@
     'jan', 'feb', 'mar', 'apr', 'may', 'jun', 'jul', 'aug', 'sep', 'oct',
     'nov', 'dec']
 
+_model_description = gettext(
+    """
+    The Seasonal Water Yield (SWY) model estimates the amount of water produced
+    by a watershed, arriving in streams over the course of a year. The primary
+    outputs of the model are quickflow, local recharge, and baseflow. Quickflow
+    represents the amount of precipitation that runs off of the land directly,
+    during and soon after a rain event, and local recharge represents the amount
+    of rainfall that infiltrates into soil, minus what is evaporated or used by
+    vegetation. Baseflow is the amount of precipitation that enters streams more
+    gradually through sub-surface flow, including during the dry season. The model
+    is based on inputs of topography (DEM), soils, land cover and management,
+    rainfall, and vegetation water demand.
+    """)
+
 MODEL_SPEC = spec.ModelSpec(
     model_id="seasonal_water_yield",
     model_title=gettext("Seasonal Water Yield"),
     userguide="seasonal_water_yield.html",
+    reporter="natcap.invest.seasonal_water_yield.reporter",
+    about=_model_description,
     validate_spatial_overlap=True,
     different_projections_ok=True,
     aliases=("swy",),

From 03ef51ffa1def3c4b7eb9f9b1b83c74d3c9512c6 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Tue, 3 Mar 2026 17:14:20 -0500
Subject: [PATCH 03/37] second iteration of swy report (#2321)

---
 .../models/seasonal_water_yield.html          |  72 ++++--
 .../invest/seasonal_water_yield/reporter.py   | 218 ++++++++++++++++--
 2 files changed, 260 insertions(+), 30 deletions(-)

diff --git a/src/natcap/invest/reports/templates/models/seasonal_water_yield.html b/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
index 43215af53d..6f1a9ac8e1 100644
--- a/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
+++ b/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
@@ -9,6 +9,15 @@
   {% from 'raster-plot-img.html' import raster_plot_img %}
 
   <h2 class="section-header">Results</h2>
+  {{ accordion_section(
+    stream_outputs_heading,
+    content_grid([
+      (caption(raster_group_caption, pre_caption=True), 100),
+      (raster_plot_img(stream_img_src, stream_outputs_heading), 100),
+      (caption(stream_caption, definition_list=True), 100)
+    ])
+  )}}
+
   {{ accordion_section(
     'Primary Outputs',
     content_grid([
@@ -19,7 +28,33 @@ <h2 class="section-header">Results</h2>
   ) }}
 
   {{ accordion_section(
-    'Results aggregated by watershed',
+    'Annual Quickflow (QF)',
+    content_grid([
+      (caption(raster_group_caption, pre_caption=True), 100),
+      (raster_plot_img(qf_img_src, 'Quickflow values'), 100),
+      (caption(qf_caption, definition_list=True), 100)
+    ])
+  )}}
+
+  {% set qf_section_contents %}
+  <div class="content-grid">
+    <div class="content-grid-item" style="width: 100%">
+      {{ caption(qf_plots_caption) }}
+    </div>
+    {% for item in qf_plots_json_tuples %}
+    <div class="content-grid-item" style="width: 33%">
+      <div id={{ item[1] }}></div>
+    </div>
+    {% endfor %}
+  </div>
+  {% endset %}
+  {{ accordion_section(
+    'Quickflow by Month per Polygon',
+    qf_section_contents
+  ) }}
+
+  {{ accordion_section(
+    'Results Aggregated by Polygon',
     content_grid([
       (content_grid([
         ('<div id="qb_map"></div>', 100),
@@ -32,34 +67,36 @@ <h2 class="section-header">Results</h2>
     ])
   ) }}
 
-  {#
-  <h2 class="section-header">Intermediate Outputs</h2>
   {{ accordion_section(
-    'Ranked exposure variables',
+    'Intermediate Outputs',
     content_grid([
-      ('<div id="rank_vars_figure"></div>', 100),
-      (caption(rank_vars_figure_caption, rank_vars_figure_source_list), 100)
+      (caption(raster_group_caption, pre_caption=True), 100),
+      (raster_plot_img(intermediate_img_src, 'Intermediate Outputs'), 100),
+      (caption(intermediate_caption, definition_list=True), 100)
     ])
-  ) }}
+  )}}
 
   {{ accordion_section(
-    'Pre-ranked variables',
+    'Output Raster Stats',
     content_grid([
-      ('<div id="facetted_histograms_figure"></div>', 100),
-      (caption(facetted_histograms_caption, facetted_histograms_source_list), 100)
+      (stats_table_note, 100),
+      (wide_table(
+        output_raster_stats_table | safe,
+        font_size_px=16
+      ), 100)
     ])
-  ) }}
+  )}}
 
+  <h2 class="section-header">Inputs</h2>
   {{ accordion_section(
-    'Wave Exposure',
+    'Raster Inputs',
     content_grid([
-      ('<div id="wave_energy_map"></div>', 65),
-      (caption(wave_energy_map_caption, wave_energy_map_source_list), 35)
+      (caption(raster_group_caption, pre_caption=True), 100),
+      (raster_plot_img(inputs_img_src, 'Raster Inputs'), 100),
+      (caption(inputs_caption, definition_list=True), 100)
     ])
   ) }}
-  #}
 
-  <h2 class="section-header">Inputs</h2>
   {{ accordion_section(
     'Arguments',
     args_table(args_dict),
@@ -85,7 +122,8 @@ <h2 class="section-header">Metadata</h2>
   {% include 'vega-embed-js.html' %}
   {% set chart_spec_id_list = [
     (qb_map_json, 'qb_map'),
-    (vri_sum_map_json, 'vri_sum_map'),
+    (vri_sum_map_json, 'vri_sum_map')
   ] %}
   {{ embed_vega(chart_spec_id_list) }}
+  {{ embed_vega(qf_plots_json_tuples) }}
 {% endblock scripts %}
diff --git a/src/natcap/invest/seasonal_water_yield/reporter.py b/src/natcap/invest/seasonal_water_yield/reporter.py
index cd525f628b..6f449a4ed1 100644
--- a/src/natcap/invest/seasonal_water_yield/reporter.py
+++ b/src/natcap/invest/seasonal_water_yield/reporter.py
@@ -4,12 +4,17 @@
 
 import altair
 import geopandas
+import numpy
 import pandas
+import pygeoprocessing
+from osgeo import gdal
+from osgeo import ogr
 
 from natcap.invest import __version__
 from natcap.invest import gettext
 import natcap.invest.spec
 from natcap.invest.reports import jinja_env, raster_utils, report_constants, vector_utils
+from natcap.invest.reports.raster_utils import RasterDatatype, RasterPlotConfig
 from natcap.invest.unit_registry import u
 
 
@@ -19,13 +24,38 @@
 
 MAP_WIDTH = 450 # pixels
 
-OUTPUT_RASTER_PLOT_TUPLES = [
-    ('l_aligned', 'continuous'),
-    ('l_avail', 'continuous'),
-    ('b', 'continuous'),
-    ('b_sum', 'continuous'),
-    ('qf', 'continuous')
-]
+qf_label_month_map = {
+    f"qf_{month_index+1}": str(month_index+1) for month_index in range(12)
+}
+
+
+def _label_to_month(row):
+    return qf_label_month_map[row['MonthLabel']]
+
+
+def _make_qf_plots(input_vector_path, id_var):
+    df = geopandas.read_file(input_vector_path)
+    melted_df = df.melt(
+        id_vars=[id_var],
+        var_name="MonthLabel",
+        value_vars=[f'qf_{month_index+1}' for month_index in range(12)],
+        value_name="Quickflow")
+
+    melted_df['Month'] = melted_df.apply(_label_to_month, axis=1)
+
+    max_qf = round(melted_df.max()['Quickflow']) + 1
+    plots = []
+
+    for poly_id in set(melted_df[id_var]):
+        plt = altair.Chart(melted_df[melted_df[id_var] == poly_id]).mark_bar().encode(
+                altair.X("month(Month):T").title("Month"),
+                altair.Y("Quickflow").title("Quickflow (m^3/s)").scale(
+                    domain=(0, max_qf), clamp=True)
+        ).properties(title=f"Mean Quickflow for Polygon, FID {poly_id}")
+
+        plots.append(plt.to_json())
+
+    return plots
 
 
 def _create_aggregate_map(geodataframe, extent_feature, xy_ratio, attribute,
@@ -68,6 +98,59 @@ def _create_aggregate_map(geodataframe, extent_feature, xy_ratio, attribute,
 #        'aggregate_vector').get_field('qb').about]
 
 
+def _aggregate_monthly_qf_cubic_meters(aoi_path, file_registry,
+                                       aggregate_vector_path):
+    if os.path.exists(aggregate_vector_path):
+        print(
+            '%s exists, deleting and writing new output',
+            aggregate_vector_path)
+        os.remove(aggregate_vector_path)
+
+    original_aoi_vector = gdal.OpenEx(aoi_path, gdal.OF_VECTOR)
+
+    driver = gdal.GetDriverByName('ESRI Shapefile')
+    driver.CreateCopy(aggregate_vector_path, original_aoi_vector)
+    gdal.Dataset.__swig_destroy__(original_aoi_vector)
+    original_aoi_vector = None
+    aggregate_vector = gdal.OpenEx(aggregate_vector_path, 1)
+    aggregate_layer = aggregate_vector.GetLayer()
+
+    path_field_tuples = [(file_registry['qf_[MONTH]'][str(month_index +1)],
+                         f"qf_{month_index+1}") for month_index in range(12)]
+
+    raster_info = pygeoprocessing.get_raster_info(path_field_tuples[0][0])
+    pixel_area_m2 = numpy.prod([abs(x) for x in raster_info['pixel_size']])
+
+    for raster_path, aggregate_field_id in path_field_tuples:
+        aggregate_stats = pygeoprocessing.zonal_statistics(
+            (raster_path, 1), aggregate_vector_path)
+
+        aggregate_field = ogr.FieldDefn(aggregate_field_id, ogr.OFTReal)
+        aggregate_field.SetWidth(24)
+        aggregate_field.SetPrecision(11)
+        aggregate_layer.CreateField(aggregate_field)
+
+        aggregate_layer.ResetReading()
+        for poly_index, poly_feat in enumerate(aggregate_layer):
+            qf_sum_meters = aggregate_stats[poly_index]['sum'] * 0.001
+            total_qf = qf_sum_meters * pixel_area_m2
+
+            poly_feat.SetField(aggregate_field_id, float(total_qf))
+            aggregate_layer.SetFeature(poly_feat)
+
+    fid_field = ogr.FieldDefn('geom_fid', ogr.OFTInteger)
+    aggregate_layer.CreateField(fid_field)
+    for feature in aggregate_layer:
+        feature_id = feature.GetFID()
+        feature.SetField('geom_fid', feature_id)
+        aggregate_layer.SetFeature(feature)
+
+    aggregate_layer.SyncToDisk()
+    aggregate_layer = None
+    gdal.Dataset.__swig_destroy__(aggregate_vector)
+    aggregate_vector = None
+
+
 def report(file_registry, args_dict, model_spec, target_html_filepath):
     """Generate an html summary of Coastal Vulnerability results.
 
@@ -88,6 +171,19 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
     aggregated_results = geopandas.read_file(file_registry['aggregate_vector'])
     extent_feature, xy_ratio = vector_utils.get_geojson_bbox(aggregated_results)
 
+    aggregated_monthly_qf_path = os.path.join(args_dict['workspace_dir'],
+                                              'monthly_qf_mean_aggregate.shp')
+    _aggregate_monthly_qf_cubic_meters(args_dict['aoi_path'], file_registry,
+                                       aggregated_monthly_qf_path)
+    qf_plots_json = _make_qf_plots(aggregated_monthly_qf_path, 'geom_fid')
+    qf_plots_json_tuples = [(qf_json, f'qf_plot_{poly_id}') for poly_id, qf_json in
+                            enumerate(qf_plots_json)]
+    qf_plots_caption = gettext(
+        """
+
+        """
+    )
+
     qb_map_json = _create_aggregate_map(
         aggregated_results, extent_feature, xy_ratio, 'qb',
         gettext('Mean local recharge value within the watershed'))
@@ -105,13 +201,97 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
 
     vector_map_source_list = [model_spec.get_output('aggregate_vector').path]
 
-    output_raster_plot_configs = raster_utils.build_raster_plot_configs(
-            file_registry, OUTPUT_RASTER_PLOT_TUPLES)
+    # Raster config lists
+    stream_raster_config_list = [
+        RasterPlotConfig(
+            raster_path=file_registry['pit_filled_dem'],
+            datatype=RasterDatatype.continuous,
+            spec=model_spec.get_output('pit_filled_dem')),
+        RasterPlotConfig(
+            raster_path=file_registry['stream'],
+            datatype=RasterDatatype.continuous,
+            spec=model_spec.get_output('stream'))]
+
+    output_raster_config_list = [
+        RasterPlotConfig(
+            raster_path=file_registry['b'],
+            datatype=RasterDatatype.continuous,
+            spec=model_spec.get_output('b')),
+        RasterPlotConfig(
+            raster_path=file_registry['annual_precip'],
+            datatype=RasterDatatype.continuous,
+            spec=model_spec.get_output('annual_precip')),
+        RasterPlotConfig(
+            raster_path=file_registry['aet'],
+            datatype=RasterDatatype.continuous,
+            spec=model_spec.get_output('aet'))]
+
+    qf_raster_config_list = [
+        RasterPlotConfig(
+            raster_path=file_registry['qf'],
+            datatype=RasterDatatype.continuous,
+            spec=model_spec.get_output('qf'))]
+
+    intermediate_raster_config_list = [
+        RasterPlotConfig(
+            raster_path=file_registry['cn'],
+            datatype=RasterDatatype.continuous,
+            spec=model_spec.get_output('cn'))]
+
+    input_raster_config_list = [
+        RasterPlotConfig(
+            raster_path=args_dict['dem_raster_path'],
+            datatype=RasterDatatype.continuous,
+            spec=model_spec.get_input('dem_raster_path')),
+        RasterPlotConfig(
+            raster_path=args_dict['lulc_raster_path'],
+            datatype=RasterDatatype.nominal,
+            spec=model_spec.get_input('lulc_raster_path'))
+    ]
+
+    local_recharge_config = RasterPlotConfig(
+        raster_path=file_registry['l_aligned'],
+        datatype=RasterDatatype.continuous,
+        spec=model_spec.get_output('l_aligned'))
+
+    if args['user_defined_local_recharge']:
+        input_raster_config_list.append(local_recharge_config)
+    else:
+        output_raster_config_list.append(local_recharge_config)
+
+    # Create raster image sources and captions:
+    stream_img_src = raster_utils.plot_and_base64_encode_rasters(
+            stream_raster_config_list)
+    stream_raster_caption = raster_utils.caption_raster_list(
+            stream_raster_config_list)
+    stream_outputs_heading = gettext(
+            'Stream Network Maps (Flow Algorithm: '
+            f'{args_dict["flow_dir_algorithm"]}, '
+            'Threshold Flow Accumulation value: '
+            f'{args_dict["thresholdf_flow_accumulation"]})')
+
     outputs_img_src = raster_utils.plot_and_base64_encode_rasters(
-            output_raster_plot_configs)
-    output_raster_caption = raster_utils.generate_caption_from_raster_list(
-            [(id, 'output') for (id, _) in OUTPUT_RASTER_PLOT_TUPLES],
-            args_dict, file_registry, model_spec)
+            output_raster_config_list)
+    output_raster_caption = raster_utils.caption_raster_list(
+            output_raster_config_list)
+
+    qf_img_src = raster_utils.plot_and_base64_encode_rasters(
+            qf_raster_config_list)
+    qf_raster_caption = raster_utils.caption_raster_list(
+            qf_raster_config_list)
+
+    intermediate_img_src = raster_utils.plot_and_base64_encode_rasters(
+            intermediate_raster_config_list)
+    intermediate_raster_caption = raster_utils.caption_raster_list(
+            intermediate_raster_config_list)
+
+    output_raster_stats_table = raster_utils.raster_workspace_summary(
+            file_registry).to_html(na_rep='')
+
+    inputs_img_src = raster_utils.plot_and_base64_encode_rasters(
+            input_raster_config_list)
+    inputs_raster_caption = raster_utils.caption_raster_list(
+            input_raster_config_list)
 
     with open(target_html_filepath, 'w', encoding='utf-8') as target_file:
         target_file.write(TEMPLATE.render(
@@ -125,8 +305,20 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
             timestamp=time.strftime('%Y-%m-%d %H:%M'),
             args_dict=args_dict,
             raster_group_caption=report_constants.RASTER_GROUP_CAPTION,
+            stream_img_src=stream_img_src,
+            stream_caption=stream_raster_caption,
+            stream_outputs_heading=stream_outputs_heading,
             outputs_img_src=outputs_img_src,
             outputs_caption=output_raster_caption,
+            qf_img_src=qf_img_src,
+            qf_caption=qf_raster_caption,
+            intermediate_img_src=intermediate_img_src,
+            intermediate_caption=intermediate_raster_caption,
+            output_raster_stats_table=output_raster_stats_table,
+            inputs_img_src=inputs_img_src,
+            inputs_caption=inputs_raster_caption,
+            qf_plots_json_tuples=qf_plots_json_tuples,
+            qf_plots_caption=qf_plots_caption,
             qb_map_json=qb_map_json,
             qb_map_caption=qb_map_caption,
             vri_sum_map_json=vri_sum_map_json,

From eb61c0d8bd0b1dfd98631916a8a58b8322de0714 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Thu, 5 Mar 2026 17:41:18 -0500
Subject: [PATCH 04/37] interactive baseflow and quickflow plots, now in m3/s
 (#2321)

---
 .../models/seasonal_water_yield.html          |  37 ++-
 .../invest/seasonal_water_yield/reporter.py   | 287 +++++++++++-------
 .../seasonal_water_yield.py                   |  24 +-
 3 files changed, 217 insertions(+), 131 deletions(-)

diff --git a/src/natcap/invest/reports/templates/models/seasonal_water_yield.html b/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
index 6f1a9ac8e1..3983493ab8 100644
--- a/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
+++ b/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
@@ -7,17 +7,9 @@
   {% from 'content-grid.html' import content_grid %}
   {% from 'metadata.html' import list_metadata %}
   {% from 'raster-plot-img.html' import raster_plot_img %}
+  {% from 'wide-table.html' import wide_table %}
 
   <h2 class="section-header">Results</h2>
-  {{ accordion_section(
-    stream_outputs_heading,
-    content_grid([
-      (caption(raster_group_caption, pre_caption=True), 100),
-      (raster_plot_img(stream_img_src, stream_outputs_heading), 100),
-      (caption(stream_caption, definition_list=True), 100)
-    ])
-  )}}
-
   {{ accordion_section(
     'Primary Outputs',
     content_grid([
@@ -36,6 +28,17 @@ <h2 class="section-header">Results</h2>
     ])
   )}}
 
+  {{ accordion_section(
+    'Quickflow and Baseflow by Month',
+    content_grid([
+      (content_grid([
+        ('<div id="qf_b_charts"></div>', 100),
+        (caption(qf_b_charts_caption, qf_b_charts_source_list), 100)
+      ]), 100)
+    ])
+  ) }}
+
+  {#
   {% set qf_section_contents %}
   <div class="content-grid">
     <div class="content-grid-item" style="width: 100%">
@@ -52,9 +55,10 @@ <h2 class="section-header">Results</h2>
     'Quickflow by Month per Polygon',
     qf_section_contents
   ) }}
+  #}
 
   {{ accordion_section(
-    'Results Aggregated by Polygon',
+    'Results Aggregated by AOI Feature',
     content_grid([
       (content_grid([
         ('<div id="qb_map"></div>', 100),
@@ -67,6 +71,15 @@ <h2 class="section-header">Results</h2>
     ])
   ) }}
 
+  {{ accordion_section(
+    stream_outputs_heading,
+    content_grid([
+      (caption(raster_group_caption, pre_caption=True), 100),
+      (raster_plot_img(stream_img_src, stream_outputs_heading), 100),
+      (caption(stream_caption, definition_list=True), 100)
+    ])
+  )}}
+
   {{ accordion_section(
     'Intermediate Outputs',
     content_grid([
@@ -122,8 +135,8 @@ <h2 class="section-header">Metadata</h2>
   {% include 'vega-embed-js.html' %}
   {% set chart_spec_id_list = [
     (qb_map_json, 'qb_map'),
-    (vri_sum_map_json, 'vri_sum_map')
+    (vri_sum_map_json, 'vri_sum_map'),
+    (qf_b_charts_json, 'qf_b_charts')
   ] %}
   {{ embed_vega(chart_spec_id_list) }}
-  {{ embed_vega(qf_plots_json_tuples) }}
 {% endblock scripts %}
diff --git a/src/natcap/invest/seasonal_water_yield/reporter.py b/src/natcap/invest/seasonal_water_yield/reporter.py
index 6f449a4ed1..6a7c611aec 100644
--- a/src/natcap/invest/seasonal_water_yield/reporter.py
+++ b/src/natcap/invest/seasonal_water_yield/reporter.py
@@ -1,3 +1,4 @@
+import csv
 import logging
 import os
 import time
@@ -33,31 +34,6 @@ def _label_to_month(row):
     return qf_label_month_map[row['MonthLabel']]
 
 
-def _make_qf_plots(input_vector_path, id_var):
-    df = geopandas.read_file(input_vector_path)
-    melted_df = df.melt(
-        id_vars=[id_var],
-        var_name="MonthLabel",
-        value_vars=[f'qf_{month_index+1}' for month_index in range(12)],
-        value_name="Quickflow")
-
-    melted_df['Month'] = melted_df.apply(_label_to_month, axis=1)
-
-    max_qf = round(melted_df.max()['Quickflow']) + 1
-    plots = []
-
-    for poly_id in set(melted_df[id_var]):
-        plt = altair.Chart(melted_df[melted_df[id_var] == poly_id]).mark_bar().encode(
-                altair.X("month(Month):T").title("Month"),
-                altair.Y("Quickflow").title("Quickflow (m^3/s)").scale(
-                    domain=(0, max_qf), clamp=True)
-        ).properties(title=f"Mean Quickflow for Polygon, FID {poly_id}")
-
-        plots.append(plt.to_json())
-
-    return plots
-
-
 def _create_aggregate_map(geodataframe, extent_feature, xy_ratio, attribute,
                           title):
     attr_map = altair.Chart(geodataframe).mark_geoshape(
@@ -78,81 +54,169 @@ def _create_aggregate_map(geodataframe, extent_feature, xy_ratio, attribute,
     ).configure_legend(**vector_utils.LEGEND_CONFIG)
 
     return attr_map.to_json()
-#    qb_map = altair.Chart(aggregated_results).mark_geoshape(
-#        clip=True,
-#        stroke="white",
-#        strokeWidth=0.5
-#    ).project(
-#        type='identity',
-#        reflectY=True,
-#        fit=extent_feature
-#    ).encode(
-#        color='qb',
-#        tooltip=[altair.Tooltip("qb", title="qb")]
-#    ).properties(
-#        width=MAP_WIDTH,
-#        height=MAP_WIDTH / xy_ratio,
-#        title='Mean local recharge value within the watershed'
-#    ).configure_legend(**vector_utils.LEGEND_CONFIG)
-#    qb_map_caption = [model_spec.get_output(
-#        'aggregate_vector').get_field('qb').about]
-
-
-def _aggregate_monthly_qf_cubic_meters(aoi_path, file_registry,
-                                       aggregate_vector_path):
-    if os.path.exists(aggregate_vector_path):
-        print(
-            '%s exists, deleting and writing new output',
-            aggregate_vector_path)
-        os.remove(aggregate_vector_path)
+
+
+def create_aoi_copy_with_fid(aoi_path, output_vector_path):
+    if os.path.exists(output_vector_path):
+        print(f'{output_vector_path} exists, deleting and writing new output')
+        os.remove(output_vector_path)
 
     original_aoi_vector = gdal.OpenEx(aoi_path, gdal.OF_VECTOR)
 
     driver = gdal.GetDriverByName('ESRI Shapefile')
-    driver.CreateCopy(aggregate_vector_path, original_aoi_vector)
+    driver.CreateCopy(output_vector_path, original_aoi_vector)
     gdal.Dataset.__swig_destroy__(original_aoi_vector)
     original_aoi_vector = None
-    aggregate_vector = gdal.OpenEx(aggregate_vector_path, 1)
-    aggregate_layer = aggregate_vector.GetLayer()
 
-    path_field_tuples = [(file_registry['qf_[MONTH]'][str(month_index +1)],
-                         f"qf_{month_index+1}") for month_index in range(12)]
+    aoi_vector = gdal.OpenEx(output_vector_path, 1)
+    aoi_layer = aoi_vector.GetLayer()
 
-    raster_info = pygeoprocessing.get_raster_info(path_field_tuples[0][0])
+    fid_field = ogr.FieldDefn('fid', ogr.OFTInteger)
+    aoi_layer.CreateField(fid_field)
+    for feature in aoi_layer:
+        feature_id = feature.GetFID()
+        feature.SetField('fid', feature_id)
+        aoi_layer.SetFeature(feature)
+
+    aoi_layer.SyncToDisk()
+    aoi_layer = None
+    gdal.Dataset.__swig_destroy__(aoi_vector)
+    aoi_vector = None
+
+
+def create_monthly_stats_table(aoi_path, file_registry, output_table_path):
+    if os.path.exists(output_table_path):
+        print(f'{output_table_path} exists, deleting and writing new output')
+        os.remove(output_table_path)
+
+    seconds_per_month = {
+        1: 2678400,
+        2: 2440152,
+        3: 2678400,
+        4: 2592000,
+        5: 2678400,
+        6: 2592000,
+        7: 2678400,
+        8: 2678400,
+        9: 2592000,
+        10: 2678400,
+        11: 2592000,
+        12: 2678400}
+
+    annual_b_path = file_registry['b']
+    monthly_qf_path_list_tuples = [
+        (file_registry['qf_[MONTH]'][str(month_index +1)], month_index+1, "quickflow")
+        for month_index in range(12)]
+    monthly_precip_path_list_tuples = [
+        (file_registry['prcp_a[MONTH]'][str(month_index)], month_index+1, "precipitation")
+        for month_index in range(12)]
+
+    # Use the baseflow raster to get the pixel_size;
+    # all rasters should be aligned + the same size
+    raster_info = pygeoprocessing.get_raster_info(annual_b_path)
     pixel_area_m2 = numpy.prod([abs(x) for x in raster_info['pixel_size']])
+    pixel_area_m2
+
+    zonal_stats_b = pygeoprocessing.zonal_statistics((annual_b_path, 1), aoi_path)
+    b_avg_per_feat_per_month = {k: v['sum'] * 0.001 * pixel_area_m2 / 12
+                                for k, v in zonal_stats_b.items()}
+
+    values_dict = {fid: {month + 1: {'baseflow': b_val / seconds_per_month[month+1]}
+                         for month in range(12)}
+                   for fid, b_val in b_avg_per_feat_per_month.items()}
+
+    for raster_path, month_index, value_name in (
+            monthly_qf_path_list_tuples + monthly_precip_path_list_tuples):
+        zonal_stats = pygeoprocessing.zonal_statistics((raster_path, 1), aoi_path)
+        avg_per_feat_per_month = {k: v['sum'] * 0.001 * pixel_area_m2
+                                  for k, v in zonal_stats.items()}
+
+        for fid, value in avg_per_feat_per_month.items():
+            values_dict[fid][month_index][value_name] = value / seconds_per_month[month_index]
+
+    with open(output_table_path, 'w', newline='') as csvfile:
+        writer = csv.writer(csvfile, delimiter=',')
+        writer.writerow(['fid', 'month', 'quickflow', 'baseflow', 'precipitation'])
+        for fid, month_dicts in values_dict.items():
+            for month, val_dicts in month_dicts.items():
+                writer.writerow([fid,
+                                 month,
+                                 val_dicts['quickflow'],
+                                 val_dicts['baseflow'],
+                                 val_dicts['precipitation']])
+
+
+def create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
+    map_df = geopandas.read_file(aoi_vector_path)
+    values_df = pandas.read_csv(aggregate_csv_path)
+    values_df.month = values_df.month.astype(str)
+
+    extent_feature, xy_ratio = vector_utils.get_geojson_bbox(map_df)
+
+    feat_select = altair.selection_point(fields=["fid"], name="feat_select", value=0)
+
+    attr_map = altair.Chart(map_df).mark_geoshape(
+        clip=True, stroke="white", strokeWidth=0.5
+        ).project(
+            type='identity',
+            reflectY=True,
+            fit=extent_feature
+    ).encode(
+        color=altair.condition(
+            feat_select,
+            altair.value("seagreen"),
+            altair.value("lightgray")
+        ),
+        tooltip=[altair.Tooltip("fid", title="fid")]
+    ).properties(
+        width=MAP_WIDTH,
+        height=MAP_WIDTH / xy_ratio,
+        title="AOI"
+    ).add_params(
+        feat_select
+    )
 
-    for raster_path, aggregate_field_id in path_field_tuples:
-        aggregate_stats = pygeoprocessing.zonal_statistics(
-            (raster_path, 1), aggregate_vector_path)
-
-        aggregate_field = ogr.FieldDefn(aggregate_field_id, ogr.OFTReal)
-        aggregate_field.SetWidth(24)
-        aggregate_field.SetPrecision(11)
-        aggregate_layer.CreateField(aggregate_field)
+    base_chart = altair.Chart(values_df)
 
-        aggregate_layer.ResetReading()
-        for poly_index, poly_feat in enumerate(aggregate_layer):
-            qf_sum_meters = aggregate_stats[poly_index]['sum'] * 0.001
-            total_qf = qf_sum_meters * pixel_area_m2
+    bar_chart = base_chart.mark_bar().transform_fold(
+        ['baseflow', 'quickflow']
+    ).encode(
+        altair.X("month(month):T").title("Month"),
+        altair.Y("value:Q").title("Quickflow + Baseflow (m^3/s)"),
+        altair.Order(field='key', sort='ascending'),
+        color=altair.Color('key:N').scale(
+            domain=['quickflow', "baseflow"],
+            range=['#fdae6b', '#9ecae1']
+        ),
+        tooltip=[altair.Tooltip(val, type="quantitative", format='.5f')
+                 for val in ["quickflow", "baseflow"]]
+    )
 
-            poly_feat.SetField(aggregate_field_id, float(total_qf))
-            aggregate_layer.SetFeature(poly_feat)
+    precip_chart = base_chart.mark_line().encode(
+        altair.X("month(month):T").title("Month"),
+        altair.Y(
+            "precipitation",
+            axis=altair.Axis(orient="right")
+        ).title("Precipitation (m^3/s)"),
+        color=altair.value('#0500a3')
+    )
 
-    fid_field = ogr.FieldDefn('geom_fid', ogr.OFTInteger)
-    aggregate_layer.CreateField(fid_field)
-    for feature in aggregate_layer:
-        feature_id = feature.GetFID()
-        feature.SetField('geom_fid', feature_id)
-        aggregate_layer.SetFeature(feature)
+    combined_chart = altair.layer(bar_chart, precip_chart).resolve_scale(
+        y='independent'
+    ).transform_filter(
+        feat_select
+    ).properties(
+        title=altair.Title(altair.expr(
+            f'"Mean Quickflow + Baseflow for Feature, FID " + {feat_select.name}.fid')
+        )
+    )
 
-    aggregate_layer.SyncToDisk()
-    aggregate_layer = None
-    gdal.Dataset.__swig_destroy__(aggregate_vector)
-    aggregate_vector = None
+    chart = attr_map | combined_chart
+    return chart.to_json()
 
 
 def report(file_registry, args_dict, model_spec, target_html_filepath):
-    """Generate an html summary of Coastal Vulnerability results.
+    """Generate an html summary of Seasonal Water Yield results.
 
     Args:
         file_registry (dict): The ``natcap.invest.FileRegistry.registry``
@@ -168,22 +232,10 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
         None
     """
 
+    # qb and vri_sum plots from the output aggregate vector
     aggregated_results = geopandas.read_file(file_registry['aggregate_vector'])
     extent_feature, xy_ratio = vector_utils.get_geojson_bbox(aggregated_results)
 
-    aggregated_monthly_qf_path = os.path.join(args_dict['workspace_dir'],
-                                              'monthly_qf_mean_aggregate.shp')
-    _aggregate_monthly_qf_cubic_meters(args_dict['aoi_path'], file_registry,
-                                       aggregated_monthly_qf_path)
-    qf_plots_json = _make_qf_plots(aggregated_monthly_qf_path, 'geom_fid')
-    qf_plots_json_tuples = [(qf_json, f'qf_plot_{poly_id}') for poly_id, qf_json in
-                            enumerate(qf_plots_json)]
-    qf_plots_caption = gettext(
-        """
-
-        """
-    )
-
     qb_map_json = _create_aggregate_map(
         aggregated_results, extent_feature, xy_ratio, 'qb',
         gettext('Mean local recharge value within the watershed'))
@@ -201,6 +253,22 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
 
     vector_map_source_list = [model_spec.get_output('aggregate_vector').path]
 
+    # Monthly quickflow + baseflow plots and map
+    aoi_copy_path = os.path.join(args_dict['workspace_dir'], 'aoi_copy.shp')
+    create_aoi_copy_with_fid(args_dict['aoi_path'], aoi_copy_path)
+    qf_b_csv_path = os.path.join(args_dict['workspace_dir'], 'monthly_average_qf_b.shp')
+    create_monthly_stats_table(aoi_copy_path, file_registry, qf_b_csv_path)
+
+    qf_b_charts_json = create_linked_monthly_plots(aoi_copy_path, qf_b_csv_path)
+    qf_b_charts_caption = gettext(
+        """
+        This chart displays the monthly combined average baseflow + quickflow for
+        each feature within the AOI. Select a feature on the AOI map to see the
+        values for that feature.
+        """
+    )
+    qf_b_charts_source_list = [qf_b_csv_path, aoi_copy_path]
+
     # Raster config lists
     stream_raster_config_list = [
         RasterPlotConfig(
@@ -249,15 +317,18 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
             spec=model_spec.get_input('lulc_raster_path'))
     ]
 
-    local_recharge_config = RasterPlotConfig(
-        raster_path=file_registry['l_aligned'],
-        datatype=RasterDatatype.continuous,
-        spec=model_spec.get_output('l_aligned'))
-
-    if args['user_defined_local_recharge']:
-        input_raster_config_list.append(local_recharge_config)
+    if args_dict['user_defined_local_recharge']:
+        input_raster_config_list.append(
+            RasterPlotConfig(
+                raster_path=args_dict['l_path'],
+                datatype=RasterDatatype.continuous,
+                spec=model_spec.get_input('l_path')))
     else:
-        output_raster_config_list.append(local_recharge_config)
+        output_raster_config_list.append(
+            RasterPlotConfig(
+                raster_path=file_registry['l'],
+                datatype=RasterDatatype.continuous,
+                spec=model_spec.get_output('l')))
 
     # Create raster image sources and captions:
     stream_img_src = raster_utils.plot_and_base64_encode_rasters(
@@ -268,7 +339,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
             'Stream Network Maps (Flow Algorithm: '
             f'{args_dict["flow_dir_algorithm"]}, '
             'Threshold Flow Accumulation value: '
-            f'{args_dict["thresholdf_flow_accumulation"]})')
+            f'{args_dict["threshold_flow_accumulation"]})')
 
     outputs_img_src = raster_utils.plot_and_base64_encode_rasters(
             output_raster_config_list)
@@ -305,6 +376,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
             timestamp=time.strftime('%Y-%m-%d %H:%M'),
             args_dict=args_dict,
             raster_group_caption=report_constants.RASTER_GROUP_CAPTION,
+            stats_table_note=report_constants.STATS_TABLE_NOTE,
             stream_img_src=stream_img_src,
             stream_caption=stream_raster_caption,
             stream_outputs_heading=stream_outputs_heading,
@@ -317,8 +389,9 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
             output_raster_stats_table=output_raster_stats_table,
             inputs_img_src=inputs_img_src,
             inputs_caption=inputs_raster_caption,
-            qf_plots_json_tuples=qf_plots_json_tuples,
-            qf_plots_caption=qf_plots_caption,
+            qf_b_charts_json=qf_b_charts_json,
+            qf_b_charts_caption=qf_b_charts_caption,
+            qf_b_charts_source_list=qf_b_charts_source_list,
             qb_map_json=qb_map_json,
             qb_map_caption=qb_map_caption,
             vri_sum_map_json=vri_sum_map_json,
diff --git a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
index 1ea07cfc76..932d5ec161 100644
--- a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
+++ b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
@@ -368,6 +368,18 @@
             data_type=float,
             units=u.none
         ),
+        spec.SingleBandRasterOutput(
+            id="l",
+            path="L.tif",
+            about=gettext(
+                "Map of local recharge. If a user-defined local recharge input"
+                " is provided, this is a copy of that layer, aligned and clipped"
+                " to match the other spatial inputs. Otherwise, this is the"
+                " local recharge as calculated by the model."
+            ),
+            data_type=float,
+            units=u.millimeter
+        ),
         spec.SingleBandRasterOutput(
             id="l_avail",
             path="L_avail.tif",
@@ -549,18 +561,6 @@
             data_type=float,
             units=u.none
         ),
-        spec.SingleBandRasterOutput(
-            id="l",
-            path="L.tif",
-            about=gettext(
-                "Map of local recharge. If a user-defined local recharge input"
-                " is provided, this is a copy of that layer, aligned and clipped"
-                " to match the other spatial inputs. Otherwise, this is the"
-                " local recharge as calculated by the model."
-            ),
-            data_type=float,
-            units=u.millimeter
-        ),
         spec.SingleBandRasterOutput(
             id="cz_aligned",
             path="intermediate_outputs/cz_aligned.tif",

From 6e1bd1d1bb5c1f3d544a024f288d3f64ad0265b1 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Fri, 6 Mar 2026 15:27:48 -0500
Subject: [PATCH 05/37] adjust raster_utils.raster_inputs_summary to check
 input csvs for raster paths (#2321)

---
 src/natcap/invest/reports/raster_utils.py | 54 +++++++++++++++++------
 1 file changed, 41 insertions(+), 13 deletions(-)

diff --git a/src/natcap/invest/reports/raster_utils.py b/src/natcap/invest/reports/raster_utils.py
index b96c973eff..54e0937080 100644
--- a/src/natcap/invest/reports/raster_utils.py
+++ b/src/natcap/invest/reports/raster_utils.py
@@ -22,7 +22,8 @@
 from pydantic.dataclasses import dataclass
 
 from natcap.invest import gettext
-from natcap.invest.spec import ModelSpec, Input, Output
+from natcap.invest.spec import ModelSpec, Input, Output, \
+    CSVInput, SingleBandRasterInput
 
 LOGGER = logging.getLogger(__name__)
 
@@ -647,19 +648,46 @@ def raster_workspace_summary(file_registry):
     return pandas.DataFrame(raster_summary).T
 
 
-def raster_inputs_summary(args_dict):
+def raster_inputs_summary(args_dict, model_spec=None, check_csv_for_rasters=False):
     """Create a table of stats for all rasters in an args_dict."""
     raster_summary = {}
-    for v in args_dict.values():
-        if isinstance(v, str) and os.path.isfile(v):
-            resource = geometamaker.describe(v, compute_stats=True)
-            if isinstance(resource, geometamaker.models.RasterResource):
-                filename = os.path.basename(resource.path)
-                band = resource.get_band_description(1)
-                raster_summary[filename] = _build_stats_table_row(
-                    resource, band)
-                # Remove 'Units' column if all units are blank
-                if not any(raster_summary[filename][UNITS_COL_NAME]):
-                    del raster_summary[filename][UNITS_COL_NAME]
+
+    paths_to_check = [v for v in args_dict.values()
+                      if isinstance(v, str) and os.path.isfile(v)]
+
+    if model_spec and check_csv_for_rasters:
+        paths_to_check.extend(_parse_csv_paths_from_spec(args_dict, model_spec))
+
+    for v in paths_to_check:
+        resource = geometamaker.describe(v, compute_stats=True)
+        if isinstance(resource, geometamaker.models.RasterResource):
+            filename = os.path.basename(resource.path)
+            band = resource.get_band_description(1)
+            raster_summary[filename] = _build_stats_table_row(
+                resource, band)
+            # Remove 'Units' column if all units are blank
+            if not any(raster_summary[filename][UNITS_COL_NAME]):
+                del raster_summary[filename][UNITS_COL_NAME]
 
     return pandas.DataFrame(raster_summary).T
+
+
+def _parse_csv_paths_from_spec(args_dict, spec):
+    table_map_inputs = []
+    for input_ in spec.inputs:
+        if isinstance(input_, CSVInput):
+            table_map_inputs.extend([
+                (input_.id, col.id) for col in spec.get_input(input_.id).columns
+                if isinstance(col, SingleBandRasterInput)])
+
+    paths_to_check = []
+    for input_id, col_name in table_map_inputs:
+        if args_dict.get(input_id):
+            df = CSVInput.get_validated_dataframe(
+                    spec.get_input(input_id),
+                    csv_path=args_dict.get(input_id))
+            paths_to_check.extend([
+                v for v in df[col_name]
+                if isinstance(v, str) and os.path.isfile(v)])
+
+    return paths_to_check

From 4afb5e16a2716ba022752ef8f2ce1f010f931c5d Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Fri, 6 Mar 2026 15:29:37 -0500
Subject: [PATCH 06/37] reorganize report sections; include monthly qf rasters
 alongside annual (#2321)

---
 .../models/seasonal_water_yield.html          | 60 +++++--------
 .../invest/seasonal_water_yield/reporter.py   | 89 ++++++++++++-------
 2 files changed, 81 insertions(+), 68 deletions(-)

diff --git a/src/natcap/invest/reports/templates/models/seasonal_water_yield.html b/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
index 3983493ab8..1de32c8860 100644
--- a/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
+++ b/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
@@ -11,19 +11,11 @@
 
   <h2 class="section-header">Results</h2>
   {{ accordion_section(
-    'Primary Outputs',
+    'Annual and Monthly Quickflow (QF)',
     content_grid([
       (caption(raster_group_caption, pre_caption=True), 100),
-      (raster_plot_img(outputs_img_src, 'Primary Outputs'), 100),
-      (caption(outputs_caption, definition_list=True), 100)
-    ])
-  ) }}
-
-  {{ accordion_section(
-    'Annual Quickflow (QF)',
-    content_grid([
-      (caption(raster_group_caption, pre_caption=True), 100),
-      (raster_plot_img(qf_img_src, 'Quickflow values'), 100),
+      (raster_plot_img(annual_qf_img_src, 'Annual Quickflow values'), 100),
+      (raster_plot_img(monthly_qf_img_src, 'Monthly Quickflow values'), 100),
       (caption(qf_caption, definition_list=True), 100)
     ])
   )}}
@@ -38,24 +30,14 @@ <h2 class="section-header">Results</h2>
     ])
   ) }}
 
-  {#
-  {% set qf_section_contents %}
-  <div class="content-grid">
-    <div class="content-grid-item" style="width: 100%">
-      {{ caption(qf_plots_caption) }}
-    </div>
-    {% for item in qf_plots_json_tuples %}
-    <div class="content-grid-item" style="width: 33%">
-      <div id={{ item[1] }}></div>
-    </div>
-    {% endfor %}
-  </div>
-  {% endset %}
   {{ accordion_section(
-    'Quickflow by Month per Polygon',
-    qf_section_contents
+    'Other Raster Outputs',
+    content_grid([
+      (caption(raster_group_caption, pre_caption=True), 100),
+      (raster_plot_img(outputs_img_src, 'Primary Outputs'), 100),
+      (caption(outputs_caption, definition_list=True), 100)
+    ])
   ) }}
-  #}
 
   {{ accordion_section(
     'Results Aggregated by AOI Feature',
@@ -80,15 +62,6 @@ <h2 class="section-header">Results</h2>
     ])
   )}}
 
-  {{ accordion_section(
-    'Intermediate Outputs',
-    content_grid([
-      (caption(raster_group_caption, pre_caption=True), 100),
-      (raster_plot_img(intermediate_img_src, 'Intermediate Outputs'), 100),
-      (caption(intermediate_caption, definition_list=True), 100)
-    ])
-  )}}
-
   {{ accordion_section(
     'Output Raster Stats',
     content_grid([
@@ -101,6 +74,11 @@ <h2 class="section-header">Results</h2>
   )}}
 
   <h2 class="section-header">Inputs</h2>
+  {{ accordion_section(
+    'Arguments',
+    args_table(args_dict),
+  )}}
+
   {{ accordion_section(
     'Raster Inputs',
     content_grid([
@@ -111,8 +89,14 @@ <h2 class="section-header">Inputs</h2>
   ) }}
 
   {{ accordion_section(
-    'Arguments',
-    args_table(args_dict),
+    'Input Raster Stats',
+    content_grid([
+      (stats_table_note, 100),
+      (wide_table(
+        input_raster_stats_table | safe,
+        font_size_px=16
+      ), 100)
+    ])
   )}}
 
   <h2 class="section-header">Metadata</h2>
diff --git a/src/natcap/invest/seasonal_water_yield/reporter.py b/src/natcap/invest/seasonal_water_yield/reporter.py
index 6a7c611aec..27c5dab172 100644
--- a/src/natcap/invest/seasonal_water_yield/reporter.py
+++ b/src/natcap/invest/seasonal_water_yield/reporter.py
@@ -58,7 +58,7 @@ def _create_aggregate_map(geodataframe, extent_feature, xy_ratio, attribute,
 
 def create_aoi_copy_with_fid(aoi_path, output_vector_path):
     if os.path.exists(output_vector_path):
-        print(f'{output_vector_path} exists, deleting and writing new output')
+        LOGGER.info(f'{output_vector_path} exists, deleting and writing new output')
         os.remove(output_vector_path)
 
     original_aoi_vector = gdal.OpenEx(aoi_path, gdal.OF_VECTOR)
@@ -71,11 +71,11 @@ def create_aoi_copy_with_fid(aoi_path, output_vector_path):
     aoi_vector = gdal.OpenEx(output_vector_path, 1)
     aoi_layer = aoi_vector.GetLayer()
 
-    fid_field = ogr.FieldDefn('fid', ogr.OFTInteger)
+    fid_field = ogr.FieldDefn('geom_fid', ogr.OFTInteger)
     aoi_layer.CreateField(fid_field)
     for feature in aoi_layer:
         feature_id = feature.GetFID()
-        feature.SetField('fid', feature_id)
+        feature.SetField('geom_fid', feature_id)
         aoi_layer.SetFeature(feature)
 
     aoi_layer.SyncToDisk()
@@ -86,7 +86,7 @@ def create_aoi_copy_with_fid(aoi_path, output_vector_path):
 
 def create_monthly_stats_table(aoi_path, file_registry, output_table_path):
     if os.path.exists(output_table_path):
-        print(f'{output_table_path} exists, deleting and writing new output')
+        LOGGER.info(f'{output_table_path} exists, deleting and writing new output')
         os.remove(output_table_path)
 
     seconds_per_month = {
@@ -136,7 +136,7 @@ def create_monthly_stats_table(aoi_path, file_registry, output_table_path):
 
     with open(output_table_path, 'w', newline='') as csvfile:
         writer = csv.writer(csvfile, delimiter=',')
-        writer.writerow(['fid', 'month', 'quickflow', 'baseflow', 'precipitation'])
+        writer.writerow(['geom_fid', 'month', 'quickflow', 'baseflow', 'precipitation'])
         for fid, month_dicts in values_dict.items():
             for month, val_dicts in month_dicts.items():
                 writer.writerow([fid,
@@ -153,7 +153,7 @@ def create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
 
     extent_feature, xy_ratio = vector_utils.get_geojson_bbox(map_df)
 
-    feat_select = altair.selection_point(fields=["fid"], name="feat_select", value=0)
+    feat_select = altair.selection_point(fields=["geom_fid"], name="feat_select", value=0)
 
     attr_map = altair.Chart(map_df).mark_geoshape(
         clip=True, stroke="white", strokeWidth=0.5
@@ -167,7 +167,7 @@ def create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
             altair.value("seagreen"),
             altair.value("lightgray")
         ),
-        tooltip=[altair.Tooltip("fid", title="fid")]
+        tooltip=[altair.Tooltip("geom_fid", title="geom_fid")]
     ).properties(
         width=MAP_WIDTH,
         height=MAP_WIDTH / xy_ratio,
@@ -207,7 +207,7 @@ def create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
         feat_select
     ).properties(
         title=altair.Title(altair.expr(
-            f'"Mean Quickflow + Baseflow for Feature, FID " + {feat_select.name}.fid')
+            f'"Mean Quickflow + Baseflow for Feature, FID " + {feat_select.name}.geom_fid')
         )
     )
 
@@ -238,7 +238,8 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
 
     qb_map_json = _create_aggregate_map(
         aggregated_results, extent_feature, xy_ratio, 'qb',
-        gettext('Mean local recharge value within the watershed'))
+        gettext("Mean local recharge value within the watershed "
+                f"({model_spec.get_output('aggregate_vector').get_field('qb').units})"))
     qb_map_caption = [
         model_spec.get_output('aggregate_vector').get_field('qb').about,
         gettext('Values are in millimeters, but should be interpreted as '
@@ -246,7 +247,8 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
 
     vri_sum_map_json = _create_aggregate_map(
         aggregated_results, extent_feature, xy_ratio, 'vri_sum',
-        gettext('Total recharge contribution of the watershed'))
+        gettext("Total recharge contribution of the watershed "
+                f"({model_spec.get_output('aggregate_vector').get_field('vri_sum').units})"))
     vri_sum_map_caption = [
         model_spec.get_output('aggregate_vector').get_field('vri_sum').about,
         gettext('The sum of ``Vri_[suffix].tif`` pixel values within the watershed.')]
@@ -277,7 +279,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
             spec=model_spec.get_output('pit_filled_dem')),
         RasterPlotConfig(
             raster_path=file_registry['stream'],
-            datatype=RasterDatatype.continuous,
+            datatype=RasterDatatype.binary_high_contrast,
             spec=model_spec.get_output('stream'))]
 
     output_raster_config_list = [
@@ -292,19 +294,30 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
         RasterPlotConfig(
             raster_path=file_registry['aet'],
             datatype=RasterDatatype.continuous,
-            spec=model_spec.get_output('aet'))]
-
-    qf_raster_config_list = [
+            spec=model_spec.get_output('aet')),
         RasterPlotConfig(
+            raster_path=file_registry['cn'],
+            datatype=RasterDatatype.continuous,
+            spec=model_spec.get_output('cn'))]
+
+    annual_qf_raster_config = RasterPlotConfig(
             raster_path=file_registry['qf'],
             datatype=RasterDatatype.continuous,
-            spec=model_spec.get_output('qf'))]
+            spec=model_spec.get_output('qf'),
+            title=gettext(
+                f"Annual Quickflow ({os.path.basename(file_registry['qf'])})"
+            ))
 
-    intermediate_raster_config_list = [
+    monthly_qf_raster_config_list = [
         RasterPlotConfig(
-            raster_path=file_registry['cn'],
+            raster_path=file_registry['qf_[MONTH]'][str(month_index + 1)],
             datatype=RasterDatatype.continuous,
-            spec=model_spec.get_output('cn'))]
+            spec=model_spec.get_output('qf_[MONTH]'),
+            title=gettext(
+                f"Quickflow for month {month_index + 1} "
+                f"({os.path.basename(file_registry['qf_[MONTH]'][str(month_index + 1)])})"
+            )
+        ) for month_index in range(12)]
 
     input_raster_config_list = [
         RasterPlotConfig(
@@ -324,6 +337,11 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
                 datatype=RasterDatatype.continuous,
                 spec=model_spec.get_input('l_path')))
     else:
+        input_raster_config_list.append(
+            RasterPlotConfig(
+                raster_path=args_dict['soil_group_path'],
+                datatype=RasterDatatype.nominal,
+                spec=model_spec.get_input('soil_group_path')))
         output_raster_config_list.append(
             RasterPlotConfig(
                 raster_path=file_registry['l'],
@@ -346,19 +364,30 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
     output_raster_caption = raster_utils.caption_raster_list(
             output_raster_config_list)
 
-    qf_img_src = raster_utils.plot_and_base64_encode_rasters(
-            qf_raster_config_list)
-    qf_raster_caption = raster_utils.caption_raster_list(
-            qf_raster_config_list)
-
-    intermediate_img_src = raster_utils.plot_and_base64_encode_rasters(
-            intermediate_raster_config_list)
-    intermediate_raster_caption = raster_utils.caption_raster_list(
-            intermediate_raster_config_list)
+    annual_qf_img_src = raster_utils.plot_and_base64_encode_rasters(
+            [annual_qf_raster_config])
+    monthly_qf_plots = raster_utils.plot_raster_facets(
+            [raster_config.raster_path for raster_config
+             in monthly_qf_raster_config_list],
+            'continuous',
+            title_list=[raster_config.title for raster_config
+                        in monthly_qf_raster_config_list])
+    monthly_qf_img_src = raster_utils.base64_encode(monthly_qf_plots)
+    monthly_qf_displayname = os.path.basename(
+            monthly_qf_raster_config_list[0].raster_path).replace('1', '[MONTH]')
+    qf_raster_caption = [
+        (f'{annual_qf_raster_config.title}:'
+         f'{annual_qf_raster_config.spec.about}'),
+        (f'{monthly_qf_displayname}:'
+         f'{monthly_qf_raster_config_list[0].spec.about}')
+    ]
 
     output_raster_stats_table = raster_utils.raster_workspace_summary(
             file_registry).to_html(na_rep='')
 
+    input_raster_stats_table = raster_utils.raster_inputs_summary(
+            args_dict, model_spec, check_csv_for_rasters=True).to_html(na_rep='')
+
     inputs_img_src = raster_utils.plot_and_base64_encode_rasters(
             input_raster_config_list)
     inputs_raster_caption = raster_utils.caption_raster_list(
@@ -382,11 +411,11 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
             stream_outputs_heading=stream_outputs_heading,
             outputs_img_src=outputs_img_src,
             outputs_caption=output_raster_caption,
-            qf_img_src=qf_img_src,
+            annual_qf_img_src=annual_qf_img_src,
+            monthly_qf_img_src=monthly_qf_img_src,
             qf_caption=qf_raster_caption,
-            intermediate_img_src=intermediate_img_src,
-            intermediate_caption=intermediate_raster_caption,
             output_raster_stats_table=output_raster_stats_table,
+            input_raster_stats_table=input_raster_stats_table,
             inputs_img_src=inputs_img_src,
             inputs_caption=inputs_raster_caption,
             qf_b_charts_json=qf_b_charts_json,

From fc6b23b45f9095772216ec77f1ebb324cf128e1b Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Mon, 9 Mar 2026 12:09:38 -0400
Subject: [PATCH 07/37] fix qf+b graph to show sums when multiple features are
 selected (#2321)

---
 .../templates/models/seasonal_water_yield.html     |  2 +-
 src/natcap/invest/seasonal_water_yield/reporter.py | 14 ++++++++------
 2 files changed, 9 insertions(+), 7 deletions(-)

diff --git a/src/natcap/invest/reports/templates/models/seasonal_water_yield.html b/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
index 1de32c8860..4d163f4db4 100644
--- a/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
+++ b/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
@@ -21,7 +21,7 @@ <h2 class="section-header">Results</h2>
   )}}
 
   {{ accordion_section(
-    'Quickflow and Baseflow by Month',
+    'Average Quickflow and Baseflow by Month',
     content_grid([
       (content_grid([
         ('<div id="qf_b_charts"></div>', 100),
diff --git a/src/natcap/invest/seasonal_water_yield/reporter.py b/src/natcap/invest/seasonal_water_yield/reporter.py
index 27c5dab172..33f9072d6a 100644
--- a/src/natcap/invest/seasonal_water_yield/reporter.py
+++ b/src/natcap/invest/seasonal_water_yield/reporter.py
@@ -182,22 +182,23 @@ def create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
         ['baseflow', 'quickflow']
     ).encode(
         altair.X("month(month):T").title("Month"),
-        altair.Y("value:Q").title("Quickflow + Baseflow (m^3/s)"),
+        altair.Y("sum(value):Q").title("Quickflow + Baseflow (m3/s)"),
         altair.Order(field='key', sort='ascending'),
         color=altair.Color('key:N').scale(
             domain=['quickflow', "baseflow"],
             range=['#fdae6b', '#9ecae1']
         ),
-        tooltip=[altair.Tooltip(val, type="quantitative", format='.5f')
-                 for val in ["quickflow", "baseflow"]]
+        tooltip=[altair.Tooltip(val, aggregate="sum", type="quantitative",
+                                format='.5f', title=val)
+                 for val in ["quickflow", "baseflow", "precipitation"]]
     )
 
     precip_chart = base_chart.mark_line().encode(
         altair.X("month(month):T").title("Month"),
         altair.Y(
-            "precipitation",
+            "sum(precipitation)",
             axis=altair.Axis(orient="right")
-        ).title("Precipitation (m^3/s)"),
+        ).title("Precipitation (m3/s)"),
         color=altair.value('#0500a3')
     )
 
@@ -266,7 +267,8 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
         """
         This chart displays the monthly combined average baseflow + quickflow for
         each feature within the AOI. Select a feature on the AOI map to see the
-        values for that feature.
+        values for that feature. Selecting multiple features will display the sum
+        of their values.
         """
     )
     qf_b_charts_source_list = [qf_b_csv_path, aoi_copy_path]

From a56e945151b1277bbe237a9d022d7b3294bec99e Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Wed, 11 Mar 2026 17:08:36 -0400
Subject: [PATCH 08/37] move chart_landmass back into CV reporter (not used
 elsewhere) (#2321)

---
 .../invest/coastal_vulnerability/reporter.py       | 14 +++++++++++++-
 src/natcap/invest/reports/vector_utils.py          | 12 ------------
 2 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/src/natcap/invest/coastal_vulnerability/reporter.py b/src/natcap/invest/coastal_vulnerability/reporter.py
index 3a98cbb140..bf9349dbbd 100644
--- a/src/natcap/invest/coastal_vulnerability/reporter.py
+++ b/src/natcap/invest/coastal_vulnerability/reporter.py
@@ -29,6 +29,18 @@
 MAP_WIDTH = 450  # pixels
 
 
+def _chart_landmass(geodataframe, clip=False, extent_feature=None):
+    landmass = altair.Chart(geodataframe).mark_geoshape(
+        clip=clip,
+        fill='lightgrey'
+    ).project(
+        type='identity',
+        reflectY=True,  # Canvas and SVG treats positive y as down
+        fit=extent_feature
+    )
+    return landmass
+
+
 def _chart_base_points(geodataframe):
     # Plot points using mark_circle instead of mark_geoshape
     # so that they can get a size encoding later if needed.
@@ -120,7 +132,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
     landmass_geo = geopandas.read_file(
         file_registry['clipped_projected_landmass'])
     extent_feature, xy_ratio = vector_utils.get_geojson_bbox(exposure_geo)
-    landmass_chart = vector_utils.chart_landmass(
+    landmass_chart = _chart_landmass(
         landmass_geo, clip=True, extent_feature=extent_feature)
     base_points = _chart_base_points(exposure_geo)
 
diff --git a/src/natcap/invest/reports/vector_utils.py b/src/natcap/invest/reports/vector_utils.py
index 06d32042e1..5dca731962 100644
--- a/src/natcap/invest/reports/vector_utils.py
+++ b/src/natcap/invest/reports/vector_utils.py
@@ -43,15 +43,3 @@ def get_geojson_bbox(geodataframe):
         "properties": {}
     }
     return extent_feature, xy_ratio
-
-
-def chart_landmass(geodataframe, clip=False, extent_feature=None):
-    landmass = altair.Chart(geodataframe).mark_geoshape(
-        clip=clip,
-        fill='lightgrey'
-    ).project(
-        type='identity',
-        reflectY=True,  # Canvas and SVG treats positive y as down
-        fit=extent_feature
-    )
-    return landmass

From d77940595c839b023598e761f28c771234ce9abd Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Wed, 11 Mar 2026 17:10:02 -0400
Subject: [PATCH 09/37] fix alignment of monthly tickmarks; add FID to
 aggregate_vector.shp rather than creating another copy of the AOI (#2321)

---
 .../invest/seasonal_water_yield/reporter.py   | 50 ++++---------------
 .../seasonal_water_yield.py                   |  7 +++
 2 files changed, 18 insertions(+), 39 deletions(-)

diff --git a/src/natcap/invest/seasonal_water_yield/reporter.py b/src/natcap/invest/seasonal_water_yield/reporter.py
index 33f9072d6a..bc56315493 100644
--- a/src/natcap/invest/seasonal_water_yield/reporter.py
+++ b/src/natcap/invest/seasonal_water_yield/reporter.py
@@ -56,34 +56,6 @@ def _create_aggregate_map(geodataframe, extent_feature, xy_ratio, attribute,
     return attr_map.to_json()
 
 
-def create_aoi_copy_with_fid(aoi_path, output_vector_path):
-    if os.path.exists(output_vector_path):
-        LOGGER.info(f'{output_vector_path} exists, deleting and writing new output')
-        os.remove(output_vector_path)
-
-    original_aoi_vector = gdal.OpenEx(aoi_path, gdal.OF_VECTOR)
-
-    driver = gdal.GetDriverByName('ESRI Shapefile')
-    driver.CreateCopy(output_vector_path, original_aoi_vector)
-    gdal.Dataset.__swig_destroy__(original_aoi_vector)
-    original_aoi_vector = None
-
-    aoi_vector = gdal.OpenEx(output_vector_path, 1)
-    aoi_layer = aoi_vector.GetLayer()
-
-    fid_field = ogr.FieldDefn('geom_fid', ogr.OFTInteger)
-    aoi_layer.CreateField(fid_field)
-    for feature in aoi_layer:
-        feature_id = feature.GetFID()
-        feature.SetField('geom_fid', feature_id)
-        aoi_layer.SetFeature(feature)
-
-    aoi_layer.SyncToDisk()
-    aoi_layer = None
-    gdal.Dataset.__swig_destroy__(aoi_vector)
-    aoi_vector = None
-
-
 def create_monthly_stats_table(aoi_path, file_registry, output_table_path):
     if os.path.exists(output_table_path):
         LOGGER.info(f'{output_table_path} exists, deleting and writing new output')
@@ -167,7 +139,7 @@ def create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
             altair.value("seagreen"),
             altair.value("lightgray")
         ),
-        tooltip=[altair.Tooltip("geom_fid", title="geom_fid")]
+        tooltip=[altair.Tooltip("geom_fid", title="FID")]
     ).properties(
         width=MAP_WIDTH,
         height=MAP_WIDTH / xy_ratio,
@@ -181,12 +153,12 @@ def create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
     bar_chart = base_chart.mark_bar().transform_fold(
         ['baseflow', 'quickflow']
     ).encode(
-        altair.X("month(month):T").title("Month"),
+        altair.X("month(month):O").title("Month"),
         altair.Y("sum(value):Q").title("Quickflow + Baseflow (m3/s)"),
         altair.Order(field='key', sort='ascending'),
         color=altair.Color('key:N').scale(
-            domain=['quickflow', "baseflow"],
-            range=['#fdae6b', '#9ecae1']
+            domain=['quickflow', "baseflow", "precipitation"],
+            range=['#fdae6b', '#9ecae1', "#0500a3"]
         ),
         tooltip=[altair.Tooltip(val, aggregate="sum", type="quantitative",
                                 format='.5f', title=val)
@@ -194,7 +166,7 @@ def create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
     )
 
     precip_chart = base_chart.mark_line().encode(
-        altair.X("month(month):T").title("Month"),
+            altair.X("month(month):O").title("Month"),
         altair.Y(
             "sum(precipitation)",
             axis=altair.Axis(orient="right")
@@ -257,12 +229,12 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
     vector_map_source_list = [model_spec.get_output('aggregate_vector').path]
 
     # Monthly quickflow + baseflow plots and map
-    aoi_copy_path = os.path.join(args_dict['workspace_dir'], 'aoi_copy.shp')
-    create_aoi_copy_with_fid(args_dict['aoi_path'], aoi_copy_path)
-    qf_b_csv_path = os.path.join(args_dict['workspace_dir'], 'monthly_average_qf_b.shp')
-    create_monthly_stats_table(aoi_copy_path, file_registry, qf_b_csv_path)
+    qf_b_csv_path = os.path.join(args_dict['workspace_dir'], 'monthly_average_qf_b.csv')
+    create_monthly_stats_table(file_registry['aggregate_vector'],
+                               file_registry, qf_b_csv_path)
 
-    qf_b_charts_json = create_linked_monthly_plots(aoi_copy_path, qf_b_csv_path)
+    qf_b_charts_json = create_linked_monthly_plots(file_registry['aggregate_vector'],
+                                                   qf_b_csv_path)
     qf_b_charts_caption = gettext(
         """
         This chart displays the monthly combined average baseflow + quickflow for
@@ -271,7 +243,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
         of their values.
         """
     )
-    qf_b_charts_source_list = [qf_b_csv_path, aoi_copy_path]
+    qf_b_charts_source_list = [qf_b_csv_path, file_registry['aggregate_vector']]
 
     # Raster config lists
     stream_raster_config_list = [
diff --git a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
index 932d5ec161..7a8e116959 100644
--- a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
+++ b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
@@ -1456,6 +1456,13 @@ def _aggregate_recharge(
             poly_feat.SetField(aggregate_field_id, float(value))
             aggregate_layer.SetFeature(poly_feat)
 
+    fid_field = ogr.FieldDefn('geom_fid', ogr.OFTInteger)
+    aggregate_layer.CreateField(fid_field)
+    for feature in aggregate_layer:
+        feature_id = feature.GetFID()
+        feature.SetField('geom_fid', feature_id)
+        aggregate_layer.SetFeature(feature)
+
     aggregate_layer.SyncToDisk()
     aggregate_layer = None
     gdal.Dataset.__swig_destroy__(aggregate_vector)

From aa4dfb40f540ce5604af2c6323d7f2fe52356ec6 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Thu, 12 Mar 2026 10:17:42 -0400
Subject: [PATCH 10/37] move construction of monthly qf & b csv into execute
 (#2321)

---
 .../invest/seasonal_water_yield/reporter.py   |  77 +-------
 .../seasonal_water_yield.py                   | 166 +++++++++++++++++-
 2 files changed, 170 insertions(+), 73 deletions(-)

diff --git a/src/natcap/invest/seasonal_water_yield/reporter.py b/src/natcap/invest/seasonal_water_yield/reporter.py
index bc56315493..6bd84ed981 100644
--- a/src/natcap/invest/seasonal_water_yield/reporter.py
+++ b/src/natcap/invest/seasonal_water_yield/reporter.py
@@ -56,68 +56,6 @@ def _create_aggregate_map(geodataframe, extent_feature, xy_ratio, attribute,
     return attr_map.to_json()
 
 
-def create_monthly_stats_table(aoi_path, file_registry, output_table_path):
-    if os.path.exists(output_table_path):
-        LOGGER.info(f'{output_table_path} exists, deleting and writing new output')
-        os.remove(output_table_path)
-
-    seconds_per_month = {
-        1: 2678400,
-        2: 2440152,
-        3: 2678400,
-        4: 2592000,
-        5: 2678400,
-        6: 2592000,
-        7: 2678400,
-        8: 2678400,
-        9: 2592000,
-        10: 2678400,
-        11: 2592000,
-        12: 2678400}
-
-    annual_b_path = file_registry['b']
-    monthly_qf_path_list_tuples = [
-        (file_registry['qf_[MONTH]'][str(month_index +1)], month_index+1, "quickflow")
-        for month_index in range(12)]
-    monthly_precip_path_list_tuples = [
-        (file_registry['prcp_a[MONTH]'][str(month_index)], month_index+1, "precipitation")
-        for month_index in range(12)]
-
-    # Use the baseflow raster to get the pixel_size;
-    # all rasters should be aligned + the same size
-    raster_info = pygeoprocessing.get_raster_info(annual_b_path)
-    pixel_area_m2 = numpy.prod([abs(x) for x in raster_info['pixel_size']])
-    pixel_area_m2
-
-    zonal_stats_b = pygeoprocessing.zonal_statistics((annual_b_path, 1), aoi_path)
-    b_avg_per_feat_per_month = {k: v['sum'] * 0.001 * pixel_area_m2 / 12
-                                for k, v in zonal_stats_b.items()}
-
-    values_dict = {fid: {month + 1: {'baseflow': b_val / seconds_per_month[month+1]}
-                         for month in range(12)}
-                   for fid, b_val in b_avg_per_feat_per_month.items()}
-
-    for raster_path, month_index, value_name in (
-            monthly_qf_path_list_tuples + monthly_precip_path_list_tuples):
-        zonal_stats = pygeoprocessing.zonal_statistics((raster_path, 1), aoi_path)
-        avg_per_feat_per_month = {k: v['sum'] * 0.001 * pixel_area_m2
-                                  for k, v in zonal_stats.items()}
-
-        for fid, value in avg_per_feat_per_month.items():
-            values_dict[fid][month_index][value_name] = value / seconds_per_month[month_index]
-
-    with open(output_table_path, 'w', newline='') as csvfile:
-        writer = csv.writer(csvfile, delimiter=',')
-        writer.writerow(['geom_fid', 'month', 'quickflow', 'baseflow', 'precipitation'])
-        for fid, month_dicts in values_dict.items():
-            for month, val_dicts in month_dicts.items():
-                writer.writerow([fid,
-                                 month,
-                                 val_dicts['quickflow'],
-                                 val_dicts['baseflow'],
-                                 val_dicts['precipitation']])
-
-
 def create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
     map_df = geopandas.read_file(aoi_vector_path)
     values_df = pandas.read_csv(aggregate_csv_path)
@@ -125,7 +63,7 @@ def create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
 
     extent_feature, xy_ratio = vector_utils.get_geojson_bbox(map_df)
 
-    feat_select = altair.selection_point(fields=["geom_fid"], name="feat_select", value=0)
+    feat_select = altair.selection_point(fields=["geom_id"], name="feat_select", value=0)
 
     attr_map = altair.Chart(map_df).mark_geoshape(
         clip=True, stroke="white", strokeWidth=0.5
@@ -139,7 +77,7 @@ def create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
             altair.value("seagreen"),
             altair.value("lightgray")
         ),
-        tooltip=[altair.Tooltip("geom_fid", title="FID")]
+        tooltip=[altair.Tooltip("geom_id", title="FID")]
     ).properties(
         width=MAP_WIDTH,
         height=MAP_WIDTH / xy_ratio,
@@ -180,7 +118,7 @@ def create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
         feat_select
     ).properties(
         title=altair.Title(altair.expr(
-            f'"Mean Quickflow + Baseflow for Feature, FID " + {feat_select.name}.geom_fid')
+            f'"Mean Quickflow + Baseflow for Feature, FID " + {feat_select.name}.geom_id')
         )
     )
 
@@ -229,12 +167,8 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
     vector_map_source_list = [model_spec.get_output('aggregate_vector').path]
 
     # Monthly quickflow + baseflow plots and map
-    qf_b_csv_path = os.path.join(args_dict['workspace_dir'], 'monthly_average_qf_b.csv')
-    create_monthly_stats_table(file_registry['aggregate_vector'],
-                               file_registry, qf_b_csv_path)
-
     qf_b_charts_json = create_linked_monthly_plots(file_registry['aggregate_vector'],
-                                                   qf_b_csv_path)
+                                                   file_registry['monthly_qf_table'])
     qf_b_charts_caption = gettext(
         """
         This chart displays the monthly combined average baseflow + quickflow for
@@ -243,7 +177,8 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
         of their values.
         """
     )
-    qf_b_charts_source_list = [qf_b_csv_path, file_registry['aggregate_vector']]
+    qf_b_charts_source_list = [file_registry['monthly_qf_table'],
+                               file_registry['aggregate_vector']]
 
     # Raster config lists
     stream_raster_config_list = [
diff --git a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
index 7a8e116959..6a775c6f12 100644
--- a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
+++ b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
@@ -1,4 +1,5 @@
 """InVEST Seasonal Water Yield Model."""
+import csv
 import fractions
 import logging
 import os
@@ -454,6 +455,61 @@
                 )
             ]
         ),
+        spec.CSVOutput(
+            id="monthly_qf_table",
+            path="monthly_quickflow_baseflow.csv",
+            about=gettext(
+                "Table of average monthly baseflow, quickflow, and precipitation"
+                " values for each watershed (or feature) within the AOI."
+            ),
+            columns=[
+                spec.IntegerOutput(
+                    id="geom_id",
+                    about=gettext(
+                        "The Feature ID for the watershed. This will correspond to"
+                        " the FIDs in the Aggregate Results shapefile."
+                    ),
+                    units=u.none
+                ),
+                spec.NumberOutput(
+                    id="month",
+                    about=gettext(
+                        "Values are the numbers 1-12 corresponding to each month,"
+                        " January (1) through December (12)."
+                    ),
+                    units=u.none
+                ),
+                spec.NumberOutput(
+                    id="quickflow",
+                    about=gettext(
+                        "The average quickflow value for the month in the watershed,"
+                        " expressed in cubic meters per second."
+                    ),
+                    units=u.meter ** 3 / u.second
+                ),
+                spec.NumberOutput(
+                    id="baseflow",
+                    about=gettext(
+                        "The average baseflow value for the month in the watershed,"
+                        " expressed in cubic meters per second. Since baseflow is"
+                        " calculated on an annual scale, the values for each watershed"
+                        " have been distributed evenly across the year (annual average"
+                        " divided by 12)."
+                    ),
+                    units=u.meter ** 3 / u.second
+                ),
+                spec.NumberOutput(
+                    id="precipitation",
+                    about=gettext(
+                        "The average precipitation value for the month in the watershed,"
+                        " expressed in cubic meters per second. Values are based on"
+                        " the aligned input monthly precipitation rasters."
+                    ),
+                    units=u.meter ** 3 / u.second
+                )
+            ],
+            index_col="geom_id"
+        ),
         spec.SingleBandRasterOutput(
             id="aet",
             path="intermediate_outputs/aet.tif",
@@ -1043,6 +1099,20 @@ def execute(args):
         dependent_task_list=b_sum_dependent_task_list + [l_sum_task],
         task_name='calculate B_sum')
 
+    monthly_csv_task = task_graph.add_task(
+        func=_generate_monthly_qf_b_p_csv,
+        args=(
+            file_registry['aggregate_vector'],
+            file_registry['b'],
+            [file_registry['qf_[MONTH]', month] for month in range(1, 13)],
+            [file_registry['prcp_a[MONTH]', month] for month in range(12)],
+            file_registry['monthly_qf_table']
+        ),
+        target_path_list=[file_registry['monthly_qf_table']],
+        dependent_task_list=[
+            aggregate_recharge_task, align_task, b_sum_task] + quick_flow_task_list,
+        task_name='create monthly qf csv')
+
     task_graph.close()
     task_graph.join()
 
@@ -1456,11 +1526,11 @@ def _aggregate_recharge(
             poly_feat.SetField(aggregate_field_id, float(value))
             aggregate_layer.SetFeature(poly_feat)
 
-    fid_field = ogr.FieldDefn('geom_fid', ogr.OFTInteger)
+    fid_field = ogr.FieldDefn('geom_id', ogr.OFTInteger)
     aggregate_layer.CreateField(fid_field)
     for feature in aggregate_layer:
         feature_id = feature.GetFID()
-        feature.SetField('geom_fid', feature_id)
+        feature.SetField('geom_id', feature_id)
         aggregate_layer.SetFeature(feature)
 
     aggregate_layer.SyncToDisk()
@@ -1469,6 +1539,98 @@ def _aggregate_recharge(
     aggregate_vector = None
 
 
+def _generate_monthly_qf_b_p_csv(
+        aoi_path, annual_baseflow_path, monthly_quickflow_path_list,
+        monthly_precip_path_list, target_csv_path):
+    """Generate a CSV of average monthly Qf, B, and P values for the watersheds/AOI.
+
+    Args:
+        aoi_path (string): path to shapefile that will be used to
+            aggregate rasters
+        annual_baseflow_path (string): path to the annual baseflow raster
+        monthly_quickflow_path_list (list): list of paths to monthly quickflow
+            rasters
+        monthly_precip_path_list (list): list of paths to aligned monthly
+            precipitation rasters
+        target_csv_path (string): path to the output CSV. If this file exists
+            on disk prior to the call, it is overwritten with the result of
+            this call.
+
+    Returns:
+        None
+    """
+    if os.path.exists(target_csv_path):
+        LOGGER.warning(f'{target_csv_path} exists, deleting and writing new output')
+        os.remove(target_csv_path)
+
+    seconds_per_month = {
+        1: 2678400,
+        2: 2440152,
+        3: 2678400,
+        4: 2592000,
+        5: 2678400,
+        6: 2592000,
+        7: 2678400,
+        8: 2678400,
+        9: 2592000,
+        10: 2678400,
+        11: 2592000,
+        12: 2678400}
+
+    # Use the baseflow raster to get the pixel_size; all rasters should be aligned
+    raster_info = pygeoprocessing.get_raster_info(annual_baseflow_path)
+    pixel_area_m2 = numpy.prod([abs(x) for x in raster_info['pixel_size']])
+
+    # The baseflow raster is annual, so there will only be 1
+    raster_path_tuples = [(annual_baseflow_path, 0, "baseflow")]
+    raster_path_tuples.extend([
+            (qf_path, month, "quickflow") for qf_path, month
+            in zip(monthly_quickflow_path_list, range(1, 13))])
+    raster_path_tuples.extend([
+            (precip_path, month, "precipitation") for precip_path, month
+            in zip(monthly_precip_path_list, range(1, 13))])
+
+    raster_stats = pygeoprocessing.zonal_statistics(
+            [(path_tuple[0], 1) for path_tuple in raster_path_tuples],
+            aoi_path)
+
+    stats_by_field = {
+        "baseflow": {},
+        "quickflow": {},
+        "precipitation": {}
+    }
+    for (_, month, field), stats in zip(raster_path_tuples, raster_stats):
+        stats_by_field[field][month] = stats
+
+    # Baseflow is computed annually; distribute evenly over the year
+    b_avg_per_feat_per_month = {k: v['sum'] * 0.001 * pixel_area_m2 / 12
+                                for k, v in stats_by_field['baseflow'][0].items()}
+
+    values_dict = {fid: {month: {'baseflow': b_val / seconds_per_month[month]}
+                         for month in range(1, 13)}
+                   for fid, b_val in b_avg_per_feat_per_month.items()}
+
+    for value_name in ['quickflow', 'precipitation']:
+        for month in range(1, 13):
+            avg_per_feat_per_month = {
+                    k: v['sum'] * 0.001 * pixel_area_m2
+                    for k, v in stats_by_field[value_name][month].items()}
+            for fid, value in avg_per_feat_per_month.items():
+                value_per_second = value / seconds_per_month[month]
+                values_dict[fid][month][value_name] = value_per_second
+
+    with open(target_csv_path, 'w', newline='') as csvfile:
+        writer = csv.writer(csvfile, delimiter=',')
+        writer.writerow(['geom_id', 'month', 'quickflow', 'baseflow', 'precipitation'])
+        for fid, month_dicts in values_dict.items():
+            for month, val_dicts in month_dicts.items():
+                writer.writerow([fid,
+                                 month,
+                                 val_dicts['quickflow'],
+                                 val_dicts['baseflow'],
+                                 val_dicts['precipitation']])
+
+
 @validation.invest_validator
 def validate(args, limit_to=None):
     """Validate args to ensure they conform to `execute`'s contract.

From 70b252e36d66a52ae1728133dd62c061f1f83795 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Thu, 12 Mar 2026 11:44:19 -0400
Subject: [PATCH 11/37] update SWY created_if conditions; update reporter
 conditional logic to reflect these conditions (#2321)

---
 .../models/seasonal_water_yield.html          |  46 ++---
 .../invest/seasonal_water_yield/reporter.py   | 164 ++++++++++--------
 .../seasonal_water_yield.py                   |  66 ++++---
 3 files changed, 158 insertions(+), 118 deletions(-)

diff --git a/src/natcap/invest/reports/templates/models/seasonal_water_yield.html b/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
index 4d163f4db4..0850e75a9b 100644
--- a/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
+++ b/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
@@ -10,28 +10,32 @@
   {% from 'wide-table.html' import wide_table %}
 
   <h2 class="section-header">Results</h2>
-  {{ accordion_section(
-    'Annual and Monthly Quickflow (QF)',
-    content_grid([
-      (caption(raster_group_caption, pre_caption=True), 100),
-      (raster_plot_img(annual_qf_img_src, 'Annual Quickflow values'), 100),
-      (raster_plot_img(monthly_qf_img_src, 'Monthly Quickflow values'), 100),
-      (caption(qf_caption, definition_list=True), 100)
-    ])
-  )}}
+  {% if qf_rasters %}
+    {{ accordion_section(
+      'Annual and Monthly Quickflow (QF)',
+      content_grid([
+        (caption(raster_group_caption, pre_caption=True), 100),
+        (raster_plot_img(qf_rasters['annual_qf_img_src'], 'Annual Quickflow values'), 100),
+        (raster_plot_img(qf_rasters['monthly_qf_img_src'], 'Monthly Quickflow values'), 100),
+        (caption(qf_rasters['qf_caption'], definition_list=True), 100)
+      ])
+    )}}
+  {% endif %}
 
-  {{ accordion_section(
-    'Average Quickflow and Baseflow by Month',
-    content_grid([
-      (content_grid([
-        ('<div id="qf_b_charts"></div>', 100),
-        (caption(qf_b_charts_caption, qf_b_charts_source_list), 100)
-      ]), 100)
-    ])
-  ) }}
+  {% if qf_b_charts %}
+    {{ accordion_section(
+      'Average Quickflow and Baseflow by Month',
+      content_grid([
+        (content_grid([
+          ('<div id="qf_b_charts"></div>', 100),
+          (caption(qf_b_charts['caption'], qf_b_charts['sources']), 100)
+        ]), 100)
+      ])
+    ) }}
+  {% endif %}
 
   {{ accordion_section(
-    'Other Raster Outputs',
+    outputs_title,
     content_grid([
       (caption(raster_group_caption, pre_caption=True), 100),
       (raster_plot_img(outputs_img_src, 'Primary Outputs'), 100),
@@ -120,7 +124,9 @@ <h2 class="section-header">Metadata</h2>
   {% set chart_spec_id_list = [
     (qb_map_json, 'qb_map'),
     (vri_sum_map_json, 'vri_sum_map'),
-    (qf_b_charts_json, 'qf_b_charts')
   ] %}
   {{ embed_vega(chart_spec_id_list) }}
+  {% if qf_b_charts %}
+    {{ embed_vega([(qf_b_charts['json'], 'qf_b_charts')]) }}
+  {% endif %}
 {% endblock scripts %}
diff --git a/src/natcap/invest/seasonal_water_yield/reporter.py b/src/natcap/invest/seasonal_water_yield/reporter.py
index 6bd84ed981..c011164d4d 100644
--- a/src/natcap/invest/seasonal_water_yield/reporter.py
+++ b/src/natcap/invest/seasonal_water_yield/reporter.py
@@ -166,19 +166,30 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
 
     vector_map_source_list = [model_spec.get_output('aggregate_vector').path]
 
-    # Monthly quickflow + baseflow plots and map
-    qf_b_charts_json = create_linked_monthly_plots(file_registry['aggregate_vector'],
-                                                   file_registry['monthly_qf_table'])
-    qf_b_charts_caption = gettext(
-        """
-        This chart displays the monthly combined average baseflow + quickflow for
-        each feature within the AOI. Select a feature on the AOI map to see the
-        values for that feature. Selecting multiple features will display the sum
-        of their values.
-        """
-    )
-    qf_b_charts_source_list = [file_registry['monthly_qf_table'],
-                               file_registry['aggregate_vector']]
+    if args_dict['user_defined_local_recharge']:
+        # Quickflow isn't calculated if `user_defined_local_recharge`
+        # so we cannot construct monthly average qf + b charts
+        qf_b_charts = None
+    else:
+        # Monthly quickflow + baseflow plots and map
+        qf_b_charts_json = create_linked_monthly_plots(
+                file_registry['aggregate_vector'],
+                file_registry['monthly_qf_table'])
+        qf_b_charts_caption = gettext(
+            """
+            This chart displays the monthly combined average baseflow + quickflow for
+            each feature within the AOI. Select a feature on the AOI map to see the
+            values for that feature. Selecting multiple features will display the sum
+            of their values.
+            """
+        )
+        qf_b_charts_source_list = [file_registry['monthly_qf_table'],
+                                   file_registry['aggregate_vector']]
+        qf_b_charts = {
+            'json': qf_b_charts_json,
+            'caption': qf_b_charts_caption,
+            'sources': qf_b_charts_source_list
+        }
 
     # Raster config lists
     stream_raster_config_list = [
@@ -195,38 +206,8 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
         RasterPlotConfig(
             raster_path=file_registry['b'],
             datatype=RasterDatatype.continuous,
-            spec=model_spec.get_output('b')),
-        RasterPlotConfig(
-            raster_path=file_registry['annual_precip'],
-            datatype=RasterDatatype.continuous,
-            spec=model_spec.get_output('annual_precip')),
-        RasterPlotConfig(
-            raster_path=file_registry['aet'],
-            datatype=RasterDatatype.continuous,
-            spec=model_spec.get_output('aet')),
-        RasterPlotConfig(
-            raster_path=file_registry['cn'],
-            datatype=RasterDatatype.continuous,
-            spec=model_spec.get_output('cn'))]
-
-    annual_qf_raster_config = RasterPlotConfig(
-            raster_path=file_registry['qf'],
-            datatype=RasterDatatype.continuous,
-            spec=model_spec.get_output('qf'),
-            title=gettext(
-                f"Annual Quickflow ({os.path.basename(file_registry['qf'])})"
-            ))
-
-    monthly_qf_raster_config_list = [
-        RasterPlotConfig(
-            raster_path=file_registry['qf_[MONTH]'][str(month_index + 1)],
-            datatype=RasterDatatype.continuous,
-            spec=model_spec.get_output('qf_[MONTH]'),
-            title=gettext(
-                f"Quickflow for month {month_index + 1} "
-                f"({os.path.basename(file_registry['qf_[MONTH]'][str(month_index + 1)])})"
-            )
-        ) for month_index in range(12)]
+            spec=model_spec.get_output('b'))
+    ]
 
     input_raster_config_list = [
         RasterPlotConfig(
@@ -245,17 +226,77 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
                 raster_path=args_dict['l_path'],
                 datatype=RasterDatatype.continuous,
                 spec=model_spec.get_input('l_path')))
+        qf_rasters = None
+        raster_outputs_title = 'Annual Baseflow'
+
     else:
+        output_raster_config_list.extend([
+            RasterPlotConfig(
+                raster_path=file_registry['annual_precip'],
+                datatype=RasterDatatype.continuous,
+                spec=model_spec.get_output('annual_precip')),
+            RasterPlotConfig(
+                raster_path=file_registry['aet'],
+                datatype=RasterDatatype.continuous,
+                spec=model_spec.get_output('aet')),
+            RasterPlotConfig(
+                raster_path=file_registry['cn'],
+                datatype=RasterDatatype.continuous,
+                spec=model_spec.get_output('cn')),
+            RasterPlotConfig(
+                raster_path=file_registry['l'],
+                datatype=RasterDatatype.continuous,
+                spec=model_spec.get_output('l'))])
+        raster_outputs_title = 'Additional Raster Outputs'
+
         input_raster_config_list.append(
             RasterPlotConfig(
                 raster_path=args_dict['soil_group_path'],
                 datatype=RasterDatatype.nominal,
                 spec=model_spec.get_input('soil_group_path')))
-        output_raster_config_list.append(
+
+        # Quickflow outputs are only created if not `user_defined_local_recharge`
+        annual_qf_raster_config = RasterPlotConfig(
+                raster_path=file_registry['qf'],
+                datatype=RasterDatatype.continuous,
+                spec=model_spec.get_output('qf'),
+                title=gettext(
+                    f"Annual Quickflow ({os.path.basename(file_registry['qf'])})"
+                ))
+
+        monthly_qf_raster_config_list = [
             RasterPlotConfig(
-                raster_path=file_registry['l'],
+                raster_path=file_registry['qf_[MONTH]'][str(month)],
                 datatype=RasterDatatype.continuous,
-                spec=model_spec.get_output('l')))
+                spec=model_spec.get_output('qf_[MONTH]'),
+                title=gettext(
+                    f"Quickflow for month {month} "
+                    f"({os.path.basename(file_registry['qf_[MONTH]'][str(month)])})"
+                )
+            ) for month in range(1, 13)]
+
+        annual_qf_img_src = raster_utils.plot_and_base64_encode_rasters(
+                [annual_qf_raster_config])
+        monthly_qf_plots = raster_utils.plot_raster_facets(
+                [raster_config.raster_path for raster_config
+                 in monthly_qf_raster_config_list],
+                'continuous',
+                title_list=[raster_config.title for raster_config
+                            in monthly_qf_raster_config_list])
+        monthly_qf_img_src = raster_utils.base64_encode(monthly_qf_plots)
+        monthly_qf_displayname = os.path.basename(
+                monthly_qf_raster_config_list[0].raster_path).replace('1', '[MONTH]')
+        qf_raster_caption = [
+            (f'{annual_qf_raster_config.title}:'
+             f'{annual_qf_raster_config.spec.about}'),
+            (f'{monthly_qf_displayname}:'
+             f'{monthly_qf_raster_config_list[0].spec.about}')
+        ]
+
+        qf_rasters = {
+            'annual_qf_img_src': annual_qf_img_src,
+            'monthly_qf_img_src': monthly_qf_img_src,
+            'qf_caption': qf_raster_caption}
 
     # Create raster image sources and captions:
     stream_img_src = raster_utils.plot_and_base64_encode_rasters(
@@ -273,24 +314,6 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
     output_raster_caption = raster_utils.caption_raster_list(
             output_raster_config_list)
 
-    annual_qf_img_src = raster_utils.plot_and_base64_encode_rasters(
-            [annual_qf_raster_config])
-    monthly_qf_plots = raster_utils.plot_raster_facets(
-            [raster_config.raster_path for raster_config
-             in monthly_qf_raster_config_list],
-            'continuous',
-            title_list=[raster_config.title for raster_config
-                        in monthly_qf_raster_config_list])
-    monthly_qf_img_src = raster_utils.base64_encode(monthly_qf_plots)
-    monthly_qf_displayname = os.path.basename(
-            monthly_qf_raster_config_list[0].raster_path).replace('1', '[MONTH]')
-    qf_raster_caption = [
-        (f'{annual_qf_raster_config.title}:'
-         f'{annual_qf_raster_config.spec.about}'),
-        (f'{monthly_qf_displayname}:'
-         f'{monthly_qf_raster_config_list[0].spec.about}')
-    ]
-
     output_raster_stats_table = raster_utils.raster_workspace_summary(
             file_registry).to_html(na_rep='')
 
@@ -318,18 +341,15 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
             stream_img_src=stream_img_src,
             stream_caption=stream_raster_caption,
             stream_outputs_heading=stream_outputs_heading,
+            outputs_title=raster_outputs_title,
             outputs_img_src=outputs_img_src,
             outputs_caption=output_raster_caption,
-            annual_qf_img_src=annual_qf_img_src,
-            monthly_qf_img_src=monthly_qf_img_src,
-            qf_caption=qf_raster_caption,
+            qf_rasters=qf_rasters,
             output_raster_stats_table=output_raster_stats_table,
             input_raster_stats_table=input_raster_stats_table,
             inputs_img_src=inputs_img_src,
             inputs_caption=inputs_raster_caption,
-            qf_b_charts_json=qf_b_charts_json,
-            qf_b_charts_caption=qf_b_charts_caption,
-            qf_b_charts_source_list=qf_b_charts_source_list,
+            qf_b_charts=qf_b_charts,
             qb_map_json=qb_map_json,
             qb_map_caption=qb_map_caption,
             vri_sum_map_json=vri_sum_map_json,
diff --git a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
index 6a775c6f12..2076185a3f 100644
--- a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
+++ b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
@@ -367,7 +367,8 @@
             path="CN.tif",
             about=gettext("Map of curve number values."),
             data_type=float,
-            units=u.none
+            units=u.none,
+            created_if="not user_defined_local_recharge"
         ),
         spec.SingleBandRasterOutput(
             id="l",
@@ -396,7 +397,8 @@
                 " available for evapotranspiration by this pixel."
             ),
             data_type=float,
-            units=u.millimeter
+            units=u.millimeter,
+            created_if="not user_defined_local_recharge"
         ),
         spec.SingleBandRasterOutput(
             id="l_sum",
@@ -414,7 +416,8 @@
             path="QF.tif",
             about=gettext("Map of quickflow"),
             data_type=float,
-            units=u.millimeter / u.year
+            units=u.millimeter / u.year,
+            created_if="not user_defined_local_recharge"
         ),
         spec.STREAM,
         spec.SingleBandRasterOutput(
@@ -422,7 +425,8 @@
             path="P.tif",
             about=gettext("The total precipitation across all months on this pixel."),
             data_type=float,
-            units=u.millimeter / u.year
+            units=u.millimeter / u.year,
+            created_if="not user_defined_local_recharge"
         ),
         spec.SingleBandRasterOutput(
             id="vri",
@@ -508,14 +512,16 @@
                     units=u.meter ** 3 / u.second
                 )
             ],
-            index_col="geom_id"
+            index_col="geom_id",
+            created_if="not user_defined_local_recharge"
         ),
         spec.SingleBandRasterOutput(
             id="aet",
             path="intermediate_outputs/aet.tif",
             about=gettext("Map of actual evapotranspiration"),
             data_type=float,
-            units=u.millimeter
+            units=u.millimeter,
+            created_if="not user_defined_local_recharge"
         ),
         spec.SingleBandRasterOutput(
             id="flow_dir",
@@ -534,14 +540,16 @@
                 "Maps of monthly quickflow (1 = January… 12 = December)"
             ),
             data_type=float,
-            units=u.millimeter
+            units=u.millimeter,
+            created_if="not user_defined_local_recharge"
         ),
         spec.SingleBandRasterOutput(
             id="si",
             path="intermediate_outputs/Si.tif",
             about=gettext("Map of the S_i factor derived from CN"),
             data_type=float,
-            units=u.inch
+            units=u.inch,
+            created_if="not user_defined_local_recharge"
         ),
         spec.SingleBandRasterOutput(
             id="lulc_aligned",
@@ -578,7 +586,8 @@
                 " other spatial inputs"
             ),
             data_type=int,
-            units=None
+            units=None,
+            created_if="not user_defined_local_recharge"
         ),
         spec.FLOW_ACCUMULATION.model_copy(update=dict(
             id="flow_accum",
@@ -591,14 +600,16 @@
                 " other spatial inputs"
             ),
             data_type=float,
-            units=u.millimeter / u.year
+            units=u.millimeter / u.year,
+            created_if="not user_defined_local_recharge"
         ),
         spec.SingleBandRasterOutput(
             id="n_events[MONTH]",
             path="intermediate_outputs/n_events[MONTH].tif",
             about=gettext("Map of monthly rain events"),
             data_type=int,
-            units=None
+            units=None,
+            created_if="not user_defined_local_recharge"
         ),
         spec.SingleBandRasterOutput(
             id="et0_a[MONTH]",
@@ -608,14 +619,16 @@
                 " spatial inputs"
             ),
             data_type=float,
-            units=u.millimeter
+            units=u.millimeter,
+            created_if="not user_defined_local_recharge"
         ),
         spec.SingleBandRasterOutput(
             id="kc_[MONTH]",
             path="intermediate_outputs/kc_[MONTH].tif",
             about=gettext("Map of monthly KC values"),
             data_type=float,
-            units=u.none
+            units=u.none,
+            created_if="not user_defined_local_recharge"
         ),
         spec.SingleBandRasterOutput(
             id="cz_aligned",
@@ -1099,19 +1112,20 @@ def execute(args):
         dependent_task_list=b_sum_dependent_task_list + [l_sum_task],
         task_name='calculate B_sum')
 
-    monthly_csv_task = task_graph.add_task(
-        func=_generate_monthly_qf_b_p_csv,
-        args=(
-            file_registry['aggregate_vector'],
-            file_registry['b'],
-            [file_registry['qf_[MONTH]', month] for month in range(1, 13)],
-            [file_registry['prcp_a[MONTH]', month] for month in range(12)],
-            file_registry['monthly_qf_table']
-        ),
-        target_path_list=[file_registry['monthly_qf_table']],
-        dependent_task_list=[
-            aggregate_recharge_task, align_task, b_sum_task] + quick_flow_task_list,
-        task_name='create monthly qf csv')
+    if not args['user_defined_local_recharge']:
+        monthly_csv_task = task_graph.add_task(
+            func=_generate_monthly_qf_b_p_csv,
+            args=(
+                file_registry['aggregate_vector'],
+                file_registry['b'],
+                [file_registry['qf_[MONTH]', month] for month in range(1, 13)],
+                [file_registry['prcp_a[MONTH]', month] for month in range(12)],
+                file_registry['monthly_qf_table']
+            ),
+            target_path_list=[file_registry['monthly_qf_table']],
+            dependent_task_list=[
+                aggregate_recharge_task, align_task, b_sum_task] + quick_flow_task_list,
+            task_name='create monthly qf csv')
 
     task_graph.close()
     task_graph.join()

From f90e7d54bb4d597a13640e3f1f467f97100548c0 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Fri, 13 Mar 2026 12:43:37 -0400
Subject: [PATCH 12/37] mask qf rasters by stream network for improved contrast
 (#2321)

---
 .../models/seasonal_water_yield.html          |  3 +-
 .../invest/seasonal_water_yield/reporter.py   | 86 ++++++++++++++++---
 2 files changed, 78 insertions(+), 11 deletions(-)

diff --git a/src/natcap/invest/reports/templates/models/seasonal_water_yield.html b/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
index 0850e75a9b..7c8646b600 100644
--- a/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
+++ b/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
@@ -12,11 +12,12 @@
   <h2 class="section-header">Results</h2>
   {% if qf_rasters %}
     {{ accordion_section(
-      'Annual and Monthly Quickflow (QF)',
+      'Annual and Monthly Quickflow (QF), masked by stream network',
       content_grid([
         (caption(raster_group_caption, pre_caption=True), 100),
         (raster_plot_img(qf_rasters['annual_qf_img_src'], 'Annual Quickflow values'), 100),
         (raster_plot_img(qf_rasters['monthly_qf_img_src'], 'Monthly Quickflow values'), 100),
+        (caption(qf_rasters['qf_caption_note']), 100),
         (caption(qf_rasters['qf_caption'], definition_list=True), 100)
       ])
     )}}
diff --git a/src/natcap/invest/seasonal_water_yield/reporter.py b/src/natcap/invest/seasonal_water_yield/reporter.py
index c011164d4d..a854535f51 100644
--- a/src/natcap/invest/seasonal_water_yield/reporter.py
+++ b/src/natcap/invest/seasonal_water_yield/reporter.py
@@ -9,6 +9,7 @@
 import pandas
 import pygeoprocessing
 from osgeo import gdal
+from osgeo import gdal_array
 from osgeo import ogr
 
 from natcap.invest import __version__
@@ -24,6 +25,8 @@
 TEMPLATE = jinja_env.get_template('models/seasonal_water_yield.html')
 
 MAP_WIDTH = 450 # pixels
+TARGET_NODATA = -1
+TARGET_DTYPE = gdal.GDT_Float32
 
 qf_label_month_map = {
     f"qf_{month_index+1}": str(month_index+1) for month_index in range(12)
@@ -56,7 +59,7 @@ def _create_aggregate_map(geodataframe, extent_feature, xy_ratio, attribute,
     return attr_map.to_json()
 
 
-def create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
+def _create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
     map_df = geopandas.read_file(aoi_vector_path)
     values_df = pandas.read_csv(aggregate_csv_path)
     values_df.month = values_df.month.astype(str)
@@ -126,6 +129,42 @@ def create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
     return chart.to_json()
 
 
+def _mask_raster_list(source_raster_path_list, mask_stream_raster_path,
+                      target_masked_vrt_path_list):
+    """Using a raster of 1s and 0s, determine which pixels remain in output.
+
+    Args:
+        source_raster_path_list (list): List of paths to source rasters containing
+            pixel values to mask by the binary raster provided as a mask.
+        mask_stream_raster_path (str): The path to a stream network raster of
+            1s and 0s indicating whether a pixel should (0) or should not (1)
+            be copied to the target raster.
+        target_masked_vrt_path_list (list): List of paths to VRTs where the target
+            rasters should be written.
+
+    Returns:
+        None
+    """
+    source_raster_info = pygeoprocessing.get_raster_info(source_raster_path_list[0])
+    source_nodata = source_raster_info['nodata'][0]
+    target_numpy_dtype = gdal_array.GDALTypeCodeToNumericTypeCode(TARGET_DTYPE)
+
+    def _mask_op(mask, raster):
+        result = numpy.full(mask.shape, TARGET_NODATA,
+                            dtype=target_numpy_dtype)
+        valid_pixels = (
+            ~pygeoprocessing.array_equals_nodata(raster, source_nodata) &
+            (mask == 0))
+        result[valid_pixels] = raster[valid_pixels]
+        return result
+
+    for source_raster_path, target_vrt_path in zip(
+            source_raster_path_list, target_masked_vrt_path_list):
+        pygeoprocessing.raster_calculator(
+            [(mask_stream_raster_path, 1), (source_raster_path, 1)], _mask_op,
+            target_vrt_path, TARGET_DTYPE, TARGET_NODATA)
+
+
 def report(file_registry, args_dict, model_spec, target_html_filepath):
     """Generate an html summary of Seasonal Water Yield results.
 
@@ -172,7 +211,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
         qf_b_charts = None
     else:
         # Monthly quickflow + baseflow plots and map
-        qf_b_charts_json = create_linked_monthly_plots(
+        qf_b_charts_json = _create_linked_monthly_plots(
                 file_registry['aggregate_vector'],
                 file_registry['monthly_qf_table'])
         qf_b_charts_caption = gettext(
@@ -183,8 +222,9 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
             of their values.
             """
         )
-        qf_b_charts_source_list = [file_registry['monthly_qf_table'],
-                                   file_registry['aggregate_vector']]
+        qf_b_charts_source_list = [
+            model_spec.get_output('monthly_qf_table').path,
+            model_spec.get_output('aggregate_vector').path]
         qf_b_charts = {
             'json': qf_b_charts_json,
             'caption': qf_b_charts_caption,
@@ -256,8 +296,21 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
                 spec=model_spec.get_input('soil_group_path')))
 
         # Quickflow outputs are only created if not `user_defined_local_recharge`
+        masked_annual_qf_path = os.path.join(
+                args_dict['workspace_dir'], 'intermediate_outputs',
+                'annual_qf_masked' + (args_dict.get('results_suffix') or '') + '.vrt')
+        masked_monthly_qf_paths = [os.path.join(
+                args_dict['workspace_dir'], 'intermediate_outputs',
+                f'qf_{month}_masked' + (args_dict.get('results_suffix') or '') + '.vrt')
+                for month in range(1, 13)]
+        _mask_raster_list(
+                [file_registry['qf']] + [file_registry['qf_[MONTH]'][str(month)]
+                                         for month in range(1, 13)],
+                file_registry['stream'],
+                [masked_annual_qf_path] + masked_monthly_qf_paths)
+
         annual_qf_raster_config = RasterPlotConfig(
-                raster_path=file_registry['qf'],
+                raster_path=masked_annual_qf_path,
                 datatype=RasterDatatype.continuous,
                 spec=model_spec.get_output('qf'),
                 title=gettext(
@@ -266,14 +319,14 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
 
         monthly_qf_raster_config_list = [
             RasterPlotConfig(
-                raster_path=file_registry['qf_[MONTH]'][str(month)],
+                raster_path=masked_qf_path,
                 datatype=RasterDatatype.continuous,
                 spec=model_spec.get_output('qf_[MONTH]'),
                 title=gettext(
                     f"Quickflow for month {month} "
                     f"({os.path.basename(file_registry['qf_[MONTH]'][str(month)])})"
                 )
-            ) for month in range(1, 13)]
+            ) for month, masked_qf_path in zip(range(1, 13), masked_monthly_qf_paths)]
 
         annual_qf_img_src = raster_utils.plot_and_base64_encode_rasters(
                 [annual_qf_raster_config])
@@ -286,16 +339,29 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
         monthly_qf_img_src = raster_utils.base64_encode(monthly_qf_plots)
         monthly_qf_displayname = os.path.basename(
                 monthly_qf_raster_config_list[0].raster_path).replace('1', '[MONTH]')
+        monthly_qf_base_raster_name = os.path.basename(
+                file_registry['qf_[MONTH]']['1']).replace('1', '[MONTH]')
+        stream_path = os.path.basename(file_registry['stream'])
+        qf_caption_note = gettext(
+                f"The stream network in `{stream_path}`"
+                " has been used as a mask to help improve contrast for the"
+                " quickflow values of non-stream pixels."),
         qf_raster_caption = [
-            (f'{annual_qf_raster_config.title}:'
-             f'{annual_qf_raster_config.spec.about}'),
+            (f'{os.path.basename(annual_qf_raster_config.raster_path)}:'
+             f'{annual_qf_raster_config.spec.about}. This VRT is derived from'
+             f' `{os.path.basename(file_registry['qf'])}`, masked by `{stream_path}`.'
+             f' (Units: {annual_qf_raster_config.spec.units})'),
             (f'{monthly_qf_displayname}:'
-             f'{monthly_qf_raster_config_list[0].spec.about}')
+             f'{monthly_qf_raster_config_list[0].spec.about}. This VRT is derived from'
+             f' the monthly quickflow rasters, `{monthly_qf_base_raster_name}`,'
+             f' masked by `{stream_path}`.'
+             f' (Units: {monthly_qf_raster_config_list[0].spec.units})')
         ]
 
         qf_rasters = {
             'annual_qf_img_src': annual_qf_img_src,
             'monthly_qf_img_src': monthly_qf_img_src,
+            'qf_caption_note': qf_caption_note,
             'qf_caption': qf_raster_caption}
 
     # Create raster image sources and captions:

From faa6ed08da5bb4dfcaf7b38a8cc09976849d4129 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Fri, 13 Mar 2026 13:06:46 -0400
Subject: [PATCH 13/37] fix quotes in f-string (#2321)

---
 src/natcap/invest/seasonal_water_yield/reporter.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/natcap/invest/seasonal_water_yield/reporter.py b/src/natcap/invest/seasonal_water_yield/reporter.py
index a854535f51..b21dc38f61 100644
--- a/src/natcap/invest/seasonal_water_yield/reporter.py
+++ b/src/natcap/invest/seasonal_water_yield/reporter.py
@@ -349,7 +349,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
         qf_raster_caption = [
             (f'{os.path.basename(annual_qf_raster_config.raster_path)}:'
              f'{annual_qf_raster_config.spec.about}. This VRT is derived from'
-             f' `{os.path.basename(file_registry['qf'])}`, masked by `{stream_path}`.'
+             f' `{os.path.basename(file_registry["qf"])}`, masked by `{stream_path}`.'
              f' (Units: {annual_qf_raster_config.spec.units})'),
             (f'{monthly_qf_displayname}:'
              f'{monthly_qf_raster_config_list[0].spec.about}. This VRT is derived from'

From b7e061c89e2a223c10b0485dc676d274f8b08ace Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Fri, 13 Mar 2026 14:24:11 -0400
Subject: [PATCH 14/37] Revert "mask qf rasters by stream network for improved
 contrast (#2321)"

This reverts commit f90e7d54bb4d597a13640e3f1f467f97100548c0.
---
 .../models/seasonal_water_yield.html          |  3 +-
 .../invest/seasonal_water_yield/reporter.py   | 86 +++----------------
 2 files changed, 11 insertions(+), 78 deletions(-)

diff --git a/src/natcap/invest/reports/templates/models/seasonal_water_yield.html b/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
index 7c8646b600..0850e75a9b 100644
--- a/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
+++ b/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
@@ -12,12 +12,11 @@
   <h2 class="section-header">Results</h2>
   {% if qf_rasters %}
     {{ accordion_section(
-      'Annual and Monthly Quickflow (QF), masked by stream network',
+      'Annual and Monthly Quickflow (QF)',
       content_grid([
         (caption(raster_group_caption, pre_caption=True), 100),
         (raster_plot_img(qf_rasters['annual_qf_img_src'], 'Annual Quickflow values'), 100),
         (raster_plot_img(qf_rasters['monthly_qf_img_src'], 'Monthly Quickflow values'), 100),
-        (caption(qf_rasters['qf_caption_note']), 100),
         (caption(qf_rasters['qf_caption'], definition_list=True), 100)
       ])
     )}}
diff --git a/src/natcap/invest/seasonal_water_yield/reporter.py b/src/natcap/invest/seasonal_water_yield/reporter.py
index b21dc38f61..c011164d4d 100644
--- a/src/natcap/invest/seasonal_water_yield/reporter.py
+++ b/src/natcap/invest/seasonal_water_yield/reporter.py
@@ -9,7 +9,6 @@
 import pandas
 import pygeoprocessing
 from osgeo import gdal
-from osgeo import gdal_array
 from osgeo import ogr
 
 from natcap.invest import __version__
@@ -25,8 +24,6 @@
 TEMPLATE = jinja_env.get_template('models/seasonal_water_yield.html')
 
 MAP_WIDTH = 450 # pixels
-TARGET_NODATA = -1
-TARGET_DTYPE = gdal.GDT_Float32
 
 qf_label_month_map = {
     f"qf_{month_index+1}": str(month_index+1) for month_index in range(12)
@@ -59,7 +56,7 @@ def _create_aggregate_map(geodataframe, extent_feature, xy_ratio, attribute,
     return attr_map.to_json()
 
 
-def _create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
+def create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
     map_df = geopandas.read_file(aoi_vector_path)
     values_df = pandas.read_csv(aggregate_csv_path)
     values_df.month = values_df.month.astype(str)
@@ -129,42 +126,6 @@ def _create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
     return chart.to_json()
 
 
-def _mask_raster_list(source_raster_path_list, mask_stream_raster_path,
-                      target_masked_vrt_path_list):
-    """Using a raster of 1s and 0s, determine which pixels remain in output.
-
-    Args:
-        source_raster_path_list (list): List of paths to source rasters containing
-            pixel values to mask by the binary raster provided as a mask.
-        mask_stream_raster_path (str): The path to a stream network raster of
-            1s and 0s indicating whether a pixel should (0) or should not (1)
-            be copied to the target raster.
-        target_masked_vrt_path_list (list): List of paths to VRTs where the target
-            rasters should be written.
-
-    Returns:
-        None
-    """
-    source_raster_info = pygeoprocessing.get_raster_info(source_raster_path_list[0])
-    source_nodata = source_raster_info['nodata'][0]
-    target_numpy_dtype = gdal_array.GDALTypeCodeToNumericTypeCode(TARGET_DTYPE)
-
-    def _mask_op(mask, raster):
-        result = numpy.full(mask.shape, TARGET_NODATA,
-                            dtype=target_numpy_dtype)
-        valid_pixels = (
-            ~pygeoprocessing.array_equals_nodata(raster, source_nodata) &
-            (mask == 0))
-        result[valid_pixels] = raster[valid_pixels]
-        return result
-
-    for source_raster_path, target_vrt_path in zip(
-            source_raster_path_list, target_masked_vrt_path_list):
-        pygeoprocessing.raster_calculator(
-            [(mask_stream_raster_path, 1), (source_raster_path, 1)], _mask_op,
-            target_vrt_path, TARGET_DTYPE, TARGET_NODATA)
-
-
 def report(file_registry, args_dict, model_spec, target_html_filepath):
     """Generate an html summary of Seasonal Water Yield results.
 
@@ -211,7 +172,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
         qf_b_charts = None
     else:
         # Monthly quickflow + baseflow plots and map
-        qf_b_charts_json = _create_linked_monthly_plots(
+        qf_b_charts_json = create_linked_monthly_plots(
                 file_registry['aggregate_vector'],
                 file_registry['monthly_qf_table'])
         qf_b_charts_caption = gettext(
@@ -222,9 +183,8 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
             of their values.
             """
         )
-        qf_b_charts_source_list = [
-            model_spec.get_output('monthly_qf_table').path,
-            model_spec.get_output('aggregate_vector').path]
+        qf_b_charts_source_list = [file_registry['monthly_qf_table'],
+                                   file_registry['aggregate_vector']]
         qf_b_charts = {
             'json': qf_b_charts_json,
             'caption': qf_b_charts_caption,
@@ -296,21 +256,8 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
                 spec=model_spec.get_input('soil_group_path')))
 
         # Quickflow outputs are only created if not `user_defined_local_recharge`
-        masked_annual_qf_path = os.path.join(
-                args_dict['workspace_dir'], 'intermediate_outputs',
-                'annual_qf_masked' + (args_dict.get('results_suffix') or '') + '.vrt')
-        masked_monthly_qf_paths = [os.path.join(
-                args_dict['workspace_dir'], 'intermediate_outputs',
-                f'qf_{month}_masked' + (args_dict.get('results_suffix') or '') + '.vrt')
-                for month in range(1, 13)]
-        _mask_raster_list(
-                [file_registry['qf']] + [file_registry['qf_[MONTH]'][str(month)]
-                                         for month in range(1, 13)],
-                file_registry['stream'],
-                [masked_annual_qf_path] + masked_monthly_qf_paths)
-
         annual_qf_raster_config = RasterPlotConfig(
-                raster_path=masked_annual_qf_path,
+                raster_path=file_registry['qf'],
                 datatype=RasterDatatype.continuous,
                 spec=model_spec.get_output('qf'),
                 title=gettext(
@@ -319,14 +266,14 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
 
         monthly_qf_raster_config_list = [
             RasterPlotConfig(
-                raster_path=masked_qf_path,
+                raster_path=file_registry['qf_[MONTH]'][str(month)],
                 datatype=RasterDatatype.continuous,
                 spec=model_spec.get_output('qf_[MONTH]'),
                 title=gettext(
                     f"Quickflow for month {month} "
                     f"({os.path.basename(file_registry['qf_[MONTH]'][str(month)])})"
                 )
-            ) for month, masked_qf_path in zip(range(1, 13), masked_monthly_qf_paths)]
+            ) for month in range(1, 13)]
 
         annual_qf_img_src = raster_utils.plot_and_base64_encode_rasters(
                 [annual_qf_raster_config])
@@ -339,29 +286,16 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
         monthly_qf_img_src = raster_utils.base64_encode(monthly_qf_plots)
         monthly_qf_displayname = os.path.basename(
                 monthly_qf_raster_config_list[0].raster_path).replace('1', '[MONTH]')
-        monthly_qf_base_raster_name = os.path.basename(
-                file_registry['qf_[MONTH]']['1']).replace('1', '[MONTH]')
-        stream_path = os.path.basename(file_registry['stream'])
-        qf_caption_note = gettext(
-                f"The stream network in `{stream_path}`"
-                " has been used as a mask to help improve contrast for the"
-                " quickflow values of non-stream pixels."),
         qf_raster_caption = [
-            (f'{os.path.basename(annual_qf_raster_config.raster_path)}:'
-             f'{annual_qf_raster_config.spec.about}. This VRT is derived from'
-             f' `{os.path.basename(file_registry["qf"])}`, masked by `{stream_path}`.'
-             f' (Units: {annual_qf_raster_config.spec.units})'),
+            (f'{annual_qf_raster_config.title}:'
+             f'{annual_qf_raster_config.spec.about}'),
             (f'{monthly_qf_displayname}:'
-             f'{monthly_qf_raster_config_list[0].spec.about}. This VRT is derived from'
-             f' the monthly quickflow rasters, `{monthly_qf_base_raster_name}`,'
-             f' masked by `{stream_path}`.'
-             f' (Units: {monthly_qf_raster_config_list[0].spec.units})')
+             f'{monthly_qf_raster_config_list[0].spec.about}')
         ]
 
         qf_rasters = {
             'annual_qf_img_src': annual_qf_img_src,
             'monthly_qf_img_src': monthly_qf_img_src,
-            'qf_caption_note': qf_caption_note,
             'qf_caption': qf_raster_caption}
 
     # Create raster image sources and captions:

From cbfe9330b5633e564c97b0b6f0ba5293a9bcd702 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Mon, 16 Mar 2026 16:52:19 -0400
Subject: [PATCH 15/37] log-transform qf plots for increased contrast; add
 small plots and custom colormap functionality (#2321)

---
 src/natcap/invest/reports/raster_utils.py     | 47 +++++++++++++------
 .../invest/seasonal_water_yield/reporter.py   | 37 +++++++++------
 2 files changed, 55 insertions(+), 29 deletions(-)

diff --git a/src/natcap/invest/reports/raster_utils.py b/src/natcap/invest/reports/raster_utils.py
index 54e0937080..d24ea8c64b 100644
--- a/src/natcap/invest/reports/raster_utils.py
+++ b/src/natcap/invest/reports/raster_utils.py
@@ -158,6 +158,8 @@ class RasterPlotConfig:
     """For highly skewed data, a transformation can help reveal variation."""
     title: str | None = None
     """An optional plot title. If ``None``, the filename is used."""
+    custom_colormap: str | None = None
+    """The name of a custom matplotlib colormap."""
 
     def __post_init__(self):
         if self.title is None:
@@ -257,7 +259,7 @@ def _extra_wide_aoi(xy_ratio):
     return xy_ratio > EX_WIDE_AOI_THRESHOLD
 
 
-def _choose_n_rows_n_cols(xy_ratio, n_plots):
+def _choose_n_rows_n_cols(xy_ratio, n_plots, small_plots):
     if _extra_wide_aoi(xy_ratio):
         n_cols = 1
     elif _wide_aoi(xy_ratio):
@@ -265,14 +267,17 @@ def _choose_n_rows_n_cols(xy_ratio, n_plots):
     else:
         n_cols = 3
 
+    if small_plots:
+        n_cols += 1
+
     if n_cols > n_plots:
         n_cols = n_plots
     n_rows = int(math.ceil(n_plots / n_cols))
     return n_rows, n_cols
 
 
-def _figure_subplots(xy_ratio, n_plots):
-    n_rows, n_cols = _choose_n_rows_n_cols(xy_ratio, n_plots)
+def _figure_subplots(xy_ratio, n_plots, small_plots=False):
+    n_rows, n_cols = _choose_n_rows_n_cols(xy_ratio, n_plots, small_plots)
 
     figure_width = MAX_FIGURE_WIDTH_DEFAULT
     if (n_cols == 2) and (_wide_aoi(xy_ratio)):
@@ -391,7 +396,7 @@ def plot_raster_list(raster_list: list[RasterPlotConfig]):
         colorbar_kwargs = {}
         imshow_kwargs['norm'] = transform
         imshow_kwargs['interpolation'] = 'none'
-        cmap = COLORMAPS[dtype]
+        cmap = config.custom_colormap if config.custom_colormap else COLORMAPS[dtype]
         if dtype == 'divergent':
             if transform == 'log':
                 transform = matplotlib.colors.SymLogNorm(linthresh=0.03)
@@ -502,7 +507,8 @@ def plot_and_base64_encode_rasters(raster_list: list[RasterPlotConfig]) -> str:
     return base64_encode(figure)
 
 
-def plot_raster_facets(tif_list, datatype, transform=None, title_list=None):
+def plot_raster_facets(tif_list, datatype, transform=None, title_list=None,
+                       small_plots=False, custom_colormap=None):
     """Plot a list of rasters that will all share a fixed colorscale.
 
     When all the rasters have the same shape and represent the same variable,
@@ -518,15 +524,20 @@ def plot_raster_facets(tif_list, datatype, transform=None, title_list=None):
             to the colormap. Either 'linear' or 'log'.
         title_list (list): Optional list of strings to use as subplot titles.
             If ``None``, the raster filename is used as the title.
+        small_plots (bool): Defaults to False. If True, the typical number of
+            columns calculated for plotting facets will be increased by 1,
+            making the plots smaller so more can be viewed side-by-side.
+        custom_colormap (str): Optional name of a matplotlib colormap to use
+            in place of the default derived from the raster datatype.
 
     """
     raster_info = pygeoprocessing.get_raster_info(tif_list[0])
     bbox = raster_info['bounding_box']
     n_plots = len(tif_list)
     xy_ratio = _get_aspect_ratio(bbox)
-    fig, axes = _figure_subplots(xy_ratio, n_plots)
+    fig, axes = _figure_subplots(xy_ratio, n_plots, small_plots=small_plots)
 
-    cmap_str = COLORMAPS[datatype]
+    cmap_str = custom_colormap if custom_colormap else COLORMAPS[datatype]
     if transform is None:
         transform = 'linear'
     if title_list is None:
@@ -560,7 +571,6 @@ def plot_raster_facets(tif_list, datatype, transform=None, title_list=None):
         if numpy.isclose(vmin, 0.0):
             vmin = 1e-6
         normalizer = matplotlib.colors.LogNorm(vmin=vmin, vmax=vmax)
-        cmap.set_under(cmap.colors[0])  # values below vmin (0s) get this color
     else:
         normalizer = plt.Normalize(vmin=vmin, vmax=vmax)
     for arr, ax, raster_path, title in zip(
@@ -569,13 +579,20 @@ def plot_raster_facets(tif_list, datatype, transform=None, title_list=None):
         # all rasters are identical size; `resampled` will be the same for all
         title_line_width = _get_title_line_width(n_plots, xy_ratio)
         ax.set_title(**_get_title_kwargs(title, resampled, title_line_width))
-        units = _get_raster_units(raster_path)
-        if units:
-            (ylim_kwargs,
-             text_kwargs) = _get_units_text_kwargs(units, len(arr))
-            ax.set_ylim(**ylim_kwargs)
-            ax.text(**text_kwargs)
-        fig.colorbar(mappable, ax=ax)
+        if not small_plots:
+            units = _get_raster_units(raster_path)
+            if units:
+                (ylim_kwargs,
+                 text_kwargs) = _get_units_text_kwargs(units, len(arr))
+                ax.set_ylim(**ylim_kwargs)
+                ax.text(**text_kwargs)
+            fig.colorbar(mappable, ax=ax)
+
+    if small_plots:
+        units = _get_raster_units(tif_list[0])
+        legend_label = f"{UNITS_TEXT}: {units}" if units else None
+        fig.colorbar(mappable, ax=axes.ravel().tolist(), label=legend_label)
+
     [ax.set_axis_off() for ax in axes.flatten()]
     return fig
 
diff --git a/src/natcap/invest/seasonal_water_yield/reporter.py b/src/natcap/invest/seasonal_water_yield/reporter.py
index c011164d4d..e5be288deb 100644
--- a/src/natcap/invest/seasonal_water_yield/reporter.py
+++ b/src/natcap/invest/seasonal_water_yield/reporter.py
@@ -1,3 +1,4 @@
+import calendar
 import csv
 import logging
 import os
@@ -14,8 +15,13 @@
 from natcap.invest import __version__
 from natcap.invest import gettext
 import natcap.invest.spec
-from natcap.invest.reports import jinja_env, raster_utils, report_constants, vector_utils
-from natcap.invest.reports.raster_utils import RasterDatatype, RasterPlotConfig
+from natcap.invest.reports import jinja_env
+from natcap.invest.reports import raster_utils
+from natcap.invest.reports import report_constants
+from natcap.invest.reports import vector_utils
+from natcap.invest.reports.raster_utils import RasterDatatype
+from natcap.invest.reports.raster_utils import RasterPlotConfig
+from natcap.invest.reports.raster_utils import RasterTransform
 from natcap.invest.unit_registry import u
 
 
@@ -56,7 +62,7 @@ def _create_aggregate_map(geodataframe, extent_feature, xy_ratio, attribute,
     return attr_map.to_json()
 
 
-def create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
+def _create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
     map_df = geopandas.read_file(aoi_vector_path)
     values_df = pandas.read_csv(aggregate_csv_path)
     values_df.month = values_df.month.astype(str)
@@ -172,7 +178,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
         qf_b_charts = None
     else:
         # Monthly quickflow + baseflow plots and map
-        qf_b_charts_json = create_linked_monthly_plots(
+        qf_b_charts_json = _create_linked_monthly_plots(
                 file_registry['aggregate_vector'],
                 file_registry['monthly_qf_table'])
         qf_b_charts_caption = gettext(
@@ -183,8 +189,9 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
             of their values.
             """
         )
-        qf_b_charts_source_list = [file_registry['monthly_qf_table'],
-                                   file_registry['aggregate_vector']]
+        qf_b_charts_source_list = [
+            model_spec.get_output('monthly_qf_table').path,
+            model_spec.get_output('aggregate_vector').path]
         qf_b_charts = {
             'json': qf_b_charts_json,
             'caption': qf_b_charts_caption,
@@ -260,19 +267,18 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
                 raster_path=file_registry['qf'],
                 datatype=RasterDatatype.continuous,
                 spec=model_spec.get_output('qf'),
-                title=gettext(
-                    f"Annual Quickflow ({os.path.basename(file_registry['qf'])})"
-                ))
+                title=gettext("Annual Quickflow"),
+                transform=RasterTransform.log,
+                custom_colormap='GnBu')
 
         monthly_qf_raster_config_list = [
             RasterPlotConfig(
                 raster_path=file_registry['qf_[MONTH]'][str(month)],
                 datatype=RasterDatatype.continuous,
                 spec=model_spec.get_output('qf_[MONTH]'),
-                title=gettext(
-                    f"Quickflow for month {month} "
-                    f"({os.path.basename(file_registry['qf_[MONTH]'][str(month)])})"
-                )
+                title=gettext(f"{calendar.month_name[month]}"),
+                transform=RasterTransform.log,
+                custom_colormap='GnBu'
             ) for month in range(1, 13)]
 
         annual_qf_img_src = raster_utils.plot_and_base64_encode_rasters(
@@ -281,8 +287,11 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
                 [raster_config.raster_path for raster_config
                  in monthly_qf_raster_config_list],
                 'continuous',
+                transform=RasterTransform.log,
                 title_list=[raster_config.title for raster_config
-                            in monthly_qf_raster_config_list])
+                            in monthly_qf_raster_config_list],
+                small_plots=True,
+                custom_colormap='GnBu')
         monthly_qf_img_src = raster_utils.base64_encode(monthly_qf_plots)
         monthly_qf_displayname = os.path.basename(
                 monthly_qf_raster_config_list[0].raster_path).replace('1', '[MONTH]')

From 169be1e7f0d15452104474d8b5aeaec8fe14797c Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Wed, 18 Mar 2026 15:13:42 -0400
Subject: [PATCH 16/37] summarize output raster stats for nested paths; expand
 utils.fake_execute to provide a test case (#2321)

---
 src/natcap/invest/reports/raster_utils.py | 28 +++++++++++++----------
 tests/reports/test_raster_utils.py        |  4 ++--
 tests/test_spec.py                        |  2 +-
 tests/utils.py                            | 20 +++++++++++-----
 4 files changed, 33 insertions(+), 21 deletions(-)

diff --git a/src/natcap/invest/reports/raster_utils.py b/src/natcap/invest/reports/raster_utils.py
index d24ea8c64b..d56ffb0eeb 100644
--- a/src/natcap/invest/reports/raster_utils.py
+++ b/src/natcap/invest/reports/raster_utils.py
@@ -633,17 +633,13 @@ def _build_stats_table_row(resource, band):
 
 
 def _get_raster_metadata(filepath):
-    if isinstance(filepath, collections.abc.Mapping):
-        for path in filepath.values():
-            return _get_raster_metadata(path)
-    else:
-        try:
-            resource = geometamaker_load(f'{filepath}.yml')
-        except (FileNotFoundError, ValueError) as err:
-            LOGGER.debug(err)
-            return None
-        if isinstance(resource, geometamaker.models.RasterResource):
-            return resource
+    try:
+        resource = geometamaker_load(f'{filepath}.yml')
+    except (FileNotFoundError, ValueError) as err:
+        LOGGER.debug(err)
+        return None
+    if isinstance(resource, geometamaker.models.RasterResource):
+        return resource
 
 
 def _get_raster_units(filepath):
@@ -654,7 +650,8 @@ def _get_raster_units(filepath):
 def raster_workspace_summary(file_registry):
     """Create a table of stats for all rasters in a file_registry."""
     raster_summary = {}
-    for path in file_registry.values():
+
+    def _summarize_output(path):
         resource = _get_raster_metadata(path)
         band = resource.get_band_description(1) if resource else None
         if band:
@@ -662,6 +659,13 @@ def raster_workspace_summary(file_registry):
             raster_summary[filename] = _build_stats_table_row(
                 resource, band)
 
+    for item in file_registry.values():
+        if isinstance(item, collections.abc.Mapping):
+            for path in item.values():
+                _summarize_output(path)
+        else:
+            _summarize_output(item)
+
     return pandas.DataFrame(raster_summary).T
 
 
diff --git a/tests/reports/test_raster_utils.py b/tests/reports/test_raster_utils.py
index f4f04a2312..66469782e6 100644
--- a/tests/reports/test_raster_utils.py
+++ b/tests/reports/test_raster_utils.py
@@ -542,5 +542,5 @@ def test_raster_workspace_summary(self):
             file_registry, args_dict)
         dataframe = raster_utils.raster_workspace_summary(file_registry)
 
-        # There are 2 rasters in the sample output spec
-        self.assertEqual(dataframe.shape, (2, 7))
+        # There are 3 rasters in the sample output spec
+        self.assertEqual(dataframe.shape, (3, 7))
diff --git a/tests/test_spec.py b/tests/test_spec.py
index 654587c16e..a7c9ec5b6a 100644
--- a/tests/test_spec.py
+++ b/tests/test_spec.py
@@ -288,7 +288,7 @@ def test_write_metadata_for_outputs(self):
         SAMPLE_MODEL_SPEC.generate_metadata_for_outputs(file_registry, args_dict)
 
         files, messages = geometamaker.validate_dir(self.workspace_dir)
-        self.assertEqual(len(files), 4)
+        self.assertEqual(len(files), 5)
         self.assertFalse(any(messages))
 
         # Test some specific content of the metadata
diff --git a/tests/utils.py b/tests/utils.py
index 0a10cbd6ff..e74279410a 100644
--- a/tests/utils.py
+++ b/tests/utils.py
@@ -107,12 +107,7 @@ def fake_execute(output_spec, workspace):
         dict (FileRegistry.registry)
 
     """
-    file_registry = FileRegistry(output_spec, workspace)
-    for spec_data in output_spec:
-        reg_key = spec_data.id
-        if '[' in spec_data.id:
-            reg_key = (spec_data.id, 'A')
-        filepath = file_registry[reg_key]
+    def _create_file(spec_data, filepath):
         if isinstance(spec_data, spec.SingleBandRasterOutput):
             driver = gdal.GetDriverByName('GTIFF')
             raster = driver.Create(filepath, 2, 2, 1, gdal.GDT_Byte)
@@ -134,4 +129,17 @@ def fake_execute(output_spec, workspace):
             # Such as taskgraph.db, just create the file.
             with open(filepath, 'w') as file:
                 pass
+
+    file_registry = FileRegistry(output_spec, workspace)
+    for spec_data in output_spec:
+        reg_key = spec_data.id
+        if '[' in spec_data.id:
+            reg_keys = [(spec_data.id, 'A'), (spec_data.id, 'B')]
+            for reg_key in reg_keys:
+                filepath = file_registry[reg_key]
+                _create_file(spec_data, filepath)
+        else:
+            filepath = file_registry[reg_key]
+            _create_file(spec_data, filepath)
+
     return file_registry.registry

From 18047bca0174fff81cb3d374506a148db2b59bdc Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Wed, 18 Mar 2026 15:16:57 -0400
Subject: [PATCH 17/37] minor edits to swy reporter, including supertitle for
 raster facet group (#2321)

---
 src/natcap/invest/reports/raster_utils.py           |  7 ++++++-
 .../templates/models/seasonal_water_yield.html      |  2 +-
 src/natcap/invest/seasonal_water_yield/reporter.py  | 13 +++++++------
 3 files changed, 14 insertions(+), 8 deletions(-)

diff --git a/src/natcap/invest/reports/raster_utils.py b/src/natcap/invest/reports/raster_utils.py
index d56ffb0eeb..5986cbe1fa 100644
--- a/src/natcap/invest/reports/raster_utils.py
+++ b/src/natcap/invest/reports/raster_utils.py
@@ -508,7 +508,7 @@ def plot_and_base64_encode_rasters(raster_list: list[RasterPlotConfig]) -> str:
 
 
 def plot_raster_facets(tif_list, datatype, transform=None, title_list=None,
-                       small_plots=False, custom_colormap=None):
+                       small_plots=False, custom_colormap=None, supertitle=None):
     """Plot a list of rasters that will all share a fixed colorscale.
 
     When all the rasters have the same shape and represent the same variable,
@@ -529,6 +529,8 @@ def plot_raster_facets(tif_list, datatype, transform=None, title_list=None,
             making the plots smaller so more can be viewed side-by-side.
         custom_colormap (str): Optional name of a matplotlib colormap to use
             in place of the default derived from the raster datatype.
+        supertitle (str): Optional title to use for the entire group of
+            raster facets.
 
     """
     raster_info = pygeoprocessing.get_raster_info(tif_list[0])
@@ -593,6 +595,9 @@ def plot_raster_facets(tif_list, datatype, transform=None, title_list=None,
         legend_label = f"{UNITS_TEXT}: {units}" if units else None
         fig.colorbar(mappable, ax=axes.ravel().tolist(), label=legend_label)
 
+    if supertitle:
+        fig.suptitle(f"{supertitle}", fontsize=TITLE_FONT_SIZE)
+
     [ax.set_axis_off() for ax in axes.flatten()]
     return fig
 
diff --git a/src/natcap/invest/reports/templates/models/seasonal_water_yield.html b/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
index 0850e75a9b..4367d097ba 100644
--- a/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
+++ b/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
@@ -35,7 +35,7 @@ <h2 class="section-header">Results</h2>
   {% endif %}
 
   {{ accordion_section(
-    outputs_title,
+    outputs_heading,
     content_grid([
       (caption(raster_group_caption, pre_caption=True), 100),
       (raster_plot_img(outputs_img_src, 'Primary Outputs'), 100),
diff --git a/src/natcap/invest/seasonal_water_yield/reporter.py b/src/natcap/invest/seasonal_water_yield/reporter.py
index e5be288deb..c1c2506ca1 100644
--- a/src/natcap/invest/seasonal_water_yield/reporter.py
+++ b/src/natcap/invest/seasonal_water_yield/reporter.py
@@ -185,8 +185,8 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
             """
             This chart displays the monthly combined average baseflow + quickflow for
             each feature within the AOI. Select a feature on the AOI map to see the
-            values for that feature. Selecting multiple features will display the sum
-            of their values.
+            values for that feature. Shift+Click to select multiple features; the chart
+            will display the sum of their values.
             """
         )
         qf_b_charts_source_list = [
@@ -234,7 +234,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
                 datatype=RasterDatatype.continuous,
                 spec=model_spec.get_input('l_path')))
         qf_rasters = None
-        raster_outputs_title = 'Annual Baseflow'
+        raster_outputs_heading = 'Annual Baseflow'
 
     else:
         output_raster_config_list.extend([
@@ -254,7 +254,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
                 raster_path=file_registry['l'],
                 datatype=RasterDatatype.continuous,
                 spec=model_spec.get_output('l'))])
-        raster_outputs_title = 'Additional Raster Outputs'
+        raster_outputs_heading = 'Additional Raster Outputs'
 
         input_raster_config_list.append(
             RasterPlotConfig(
@@ -291,7 +291,8 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
                 title_list=[raster_config.title for raster_config
                             in monthly_qf_raster_config_list],
                 small_plots=True,
-                custom_colormap='GnBu')
+                custom_colormap='GnBu',
+                supertitle=gettext("Monthly Quickflow"))
         monthly_qf_img_src = raster_utils.base64_encode(monthly_qf_plots)
         monthly_qf_displayname = os.path.basename(
                 monthly_qf_raster_config_list[0].raster_path).replace('1', '[MONTH]')
@@ -350,7 +351,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
             stream_img_src=stream_img_src,
             stream_caption=stream_raster_caption,
             stream_outputs_heading=stream_outputs_heading,
-            outputs_title=raster_outputs_title,
+            outputs_heading=raster_outputs_heading,
             outputs_img_src=outputs_img_src,
             outputs_caption=output_raster_caption,
             qf_rasters=qf_rasters,

From 3b3b0f50fa7216f298083c15943a512843020172 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Wed, 18 Mar 2026 15:17:54 -0400
Subject: [PATCH 18/37] update SWY outputs units in model_spec (#2321)

---
 .../seasonal_water_yield.py                   | 34 +++++++++----------
 1 file changed, 17 insertions(+), 17 deletions(-)

diff --git a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
index 2076185a3f..406ac684c1 100644
--- a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
+++ b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
@@ -349,7 +349,7 @@
                 " (which is not evapotranspired before it reaches the stream)."
             ),
             data_type=float,
-            units=u.millimeter
+            units=u.millimeter / u.year
         ),
         spec.SingleBandRasterOutput(
             id="b_sum",
@@ -360,7 +360,7 @@
                 " stream."
             ),
             data_type=float,
-            units=u.millimeter
+            units=u.millimeter / u.year
         ),
         spec.SingleBandRasterOutput(
             id="cn",
@@ -380,14 +380,14 @@
                 " local recharge as calculated by the model."
             ),
             data_type=float,
-            units=u.millimeter
+            units=u.millimeter / u.year
         ),
         spec.SingleBandRasterOutput(
             id="l_avail",
             path="L_avail.tif",
             about=gettext("Map of available local recharge"),
             data_type=float,
-            units=u.millimeter
+            units=u.millimeter / u.year
         ),
         spec.SingleBandRasterOutput(
             id="l_sum_avail",
@@ -397,7 +397,7 @@
                 " available for evapotranspiration by this pixel."
             ),
             data_type=float,
-            units=u.millimeter,
+            units=u.millimeter / u.year,
             created_if="not user_defined_local_recharge"
         ),
         spec.SingleBandRasterOutput(
@@ -409,7 +409,7 @@
                 " evapotranspiration to downslope pixels."
             ),
             data_type=float,
-            units=u.millimeter
+            units=u.millimeter / u.year
         ),
         spec.SingleBandRasterOutput(
             id="qf",
@@ -436,7 +436,7 @@
                 " the total recharge."
             ),
             data_type=float,
-            units=u.millimeter
+            units=u.millimeter / u.year
         ),
         spec.VectorOutput(
             id="aggregate_vector",
@@ -447,7 +447,7 @@
                 spec.NumberOutput(
                     id="qb",
                     about=gettext("Mean local recharge value within the watershed"),
-                    units=u.millimeter
+                    units=u.millimeter / u.year
                 ),
                 spec.NumberOutput(
                     id="vri_sum",
@@ -455,7 +455,7 @@
                         "Total recharge contribution, (positive or negative) within the"
                         " watershed."
                     ),
-                    units=u.millimeter
+                    units=u.millimeter / u.year
                 )
             ]
         ),
@@ -520,7 +520,7 @@
             path="intermediate_outputs/aet.tif",
             about=gettext("Map of actual evapotranspiration"),
             data_type=float,
-            units=u.millimeter,
+            units=u.millimeter / u.year,
             created_if="not user_defined_local_recharge"
         ),
         spec.SingleBandRasterOutput(
@@ -531,7 +531,7 @@
                 " the option selected."
             ),
             data_type=int,
-            units=None
+            units=u.none
         ),
         spec.SingleBandRasterOutput(
             id="qf_[MONTH]",
@@ -540,7 +540,7 @@
                 "Maps of monthly quickflow (1 = January… 12 = December)"
             ),
             data_type=float,
-            units=u.millimeter,
+            units=u.millimeter / u.month,
             created_if="not user_defined_local_recharge"
         ),
         spec.SingleBandRasterOutput(
@@ -559,7 +559,7 @@
                 " spatial inputs"
             ),
             data_type=int,
-            units=None
+            units=u.none
         ),
         spec.SingleBandRasterOutput(
             id="dem_aligned",
@@ -586,7 +586,7 @@
                 " other spatial inputs"
             ),
             data_type=int,
-            units=None,
+            units=u.none,
             created_if="not user_defined_local_recharge"
         ),
         spec.FLOW_ACCUMULATION.model_copy(update=dict(
@@ -608,7 +608,7 @@
             path="intermediate_outputs/n_events[MONTH].tif",
             about=gettext("Map of monthly rain events"),
             data_type=int,
-            units=None,
+            units=u.none,
             created_if="not user_defined_local_recharge"
         ),
         spec.SingleBandRasterOutput(
@@ -619,7 +619,7 @@
                 " spatial inputs"
             ),
             data_type=float,
-            units=u.millimeter,
+            units=u.millimeter / u.month,
             created_if="not user_defined_local_recharge"
         ),
         spec.SingleBandRasterOutput(
@@ -639,7 +639,7 @@
             ),
             created_if="user_defined_climate_zones",
             data_type=int,
-            units=None
+            units=u.none
         ),
         spec.TASKGRAPH_CACHE
     ]

From d09660938cfa90a30cd6a7dbbd1ed6c8f4fc9fdb Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Wed, 18 Mar 2026 15:19:09 -0400
Subject: [PATCH 19/37] basic testing for swy reporter template (#2321)

---
 tests/test_seasonal_water_yield_regression.py | 85 +++++++++++++++++--
 1 file changed, 80 insertions(+), 5 deletions(-)

diff --git a/tests/test_seasonal_water_yield_regression.py b/tests/test_seasonal_water_yield_regression.py
index 777064409f..c75b7f2bc5 100644
--- a/tests/test_seasonal_water_yield_regression.py
+++ b/tests/test_seasonal_water_yield_regression.py
@@ -14,9 +14,6 @@
 from .utils import assert_complete_execute
 
 gdal.UseExceptions()
-REGRESSION_DATA = os.path.join(
-    os.path.dirname(__file__), '..', 'data', 'invest-test-data',
-    'seasonal_water_yield')
 
 
 def make_simple_shp(base_shp_path, origin):
@@ -508,10 +505,12 @@ class SeasonalWaterYieldRegressionTests(unittest.TestCase):
     def setUp(self):
         """Create temporary workspace."""
         self.workspace_dir = tempfile.mkdtemp()
+#        self.workspace_dir = "/Users/megannissel/dev/debugging/swy_reports/"
 
     def tearDown(self):
         """Remove temporary workspace."""
         shutil.rmtree(self.workspace_dir)
+#        pass
 
     @staticmethod
     def generate_base_args(workspace_dir):
@@ -601,7 +600,8 @@ def test_base_regression(self):
 
         execute_kwargs = {
             'generate_report': bool(seasonal_water_yield.MODEL_SPEC.reporter),
-            'save_file_registry': True
+            'save_file_registry': True,
+            'check_outputs': True
         }
         seasonal_water_yield.MODEL_SPEC.execute(args, **execute_kwargs)
         assert_complete_execute(
@@ -616,6 +616,23 @@ def test_base_regression(self):
             os.path.join(args['workspace_dir'], 'aggregated_results_swy.shp'),
             agg_results_csv_path)
 
+        # check the values in the avg monthly quickflow baseflow precip csv
+        actual_result_df = pandas.read_csv(
+            os.path.join(args['workspace_dir'], 'monthly_quickflow_baseflow.csv'))
+        expected_qf = [
+            2.116894e-5, 2.541212e-5, 2.514573e-5, 2.805605e-5, 2.916745e-5, 3.223471e-5,
+            3.323317e-5, 3.528216e-5, 3.858638e-5, 3.941151e-5, 4.287460e-5, 4.358156e-5]
+        expected_b = [
+            2.276286e-5, 2.498535e-5, 2.276286e-5, 2.352162e-5, 2.276286e-5, 2.352162e-5,
+            2.276286e-5, 2.276286e-5, 2.352162e-5, 2.276286e-5, 2.352162e-5, 2.276286e-5]
+        expected_p = [
+            4.106930e-5, 4.917726e-5, 4.853644e-5, 5.401235e-5, 5.600358e-5, 6.172840e-5,
+            6.347073e-5, 6.720430e-5, 7.330247e-5, 7.467145e-5, 8.101852e-5, 8.213859e-5]
+        for expected_val, col_name in [(expected_qf, 'quickflow'),
+                (expected_b, 'baseflow'), (expected_p, 'precipitation')]:
+            numpy.testing.assert_allclose(expected_val, actual_result_df[col_name],
+                                          rtol=1e-6)
+
     def test_base_regression_d8(self):
         """SWY base regression test on sample data in D8 mode.
 
@@ -668,6 +685,23 @@ def test_base_regression_d8(self):
         if mismatch_list:
             raise RuntimeError(f'results not expected: {mismatch_list}')
 
+        # check the values in the avg monthly quickflow baseflow precip csv
+        actual_result_df = pandas.read_csv(
+            os.path.join(args['workspace_dir'], 'monthly_quickflow_baseflow.csv'))
+        expected_qf = [
+            2.083911e-5, 2.501863e-5, 2.475883e-5, 2.762715e-5, 2.872444e-5, 3.174831e-5,
+            3.273498e-5, 3.475670e-5, 3.801546e-5, 3.883215e-5, 4.224839e-5, 4.294909e-5]
+        expected_b = [
+            2.313658e-5, 2.539555e-5, 2.313658e-5, 2.390779e-5, 2.313658e-5, 2.390779e-5,
+            2.313658e-5, 2.313658e-5, 2.390779e-5, 2.313658e-5, 2.390779e-5, 2.313658e-5]
+        expected_p = [
+            4.106930e-5, 4.917726e-5, 4.853644e-5, 5.401235e-5, 5.600358e-5, 6.172840e-5,
+            6.347073e-5, 6.720430e-5, 7.330247e-5, 7.467145e-5, 8.101852e-5, 8.213859e-5]
+        for expected_val, col_name in [(expected_qf, 'quickflow'),
+                (expected_b, 'baseflow'), (expected_p, 'precipitation')]:
+            numpy.testing.assert_allclose(expected_val, actual_result_df[col_name],
+                                          rtol=1e-6)
+
     def test_base_regression_nodata_inf(self):
         """SWY base regression test on sample data with really small nodata.
 
@@ -746,6 +780,23 @@ def test_base_regression_nodata_inf(self):
             os.path.join(args['workspace_dir'], 'aggregated_results_swy.shp'),
             agg_results_csv_path)
 
+        # check the values in the avg monthly quickflow baseflow precip csv
+        actual_result_df = pandas.read_csv(
+            os.path.join(args['workspace_dir'], 'monthly_quickflow_baseflow.csv'))
+        expected_qf = [
+            2.102520e-5, 2.524000e-5, 2.497585e-5, 2.786701e-5, 2.897143e-5, 3.201866e-5,
+            3.301102e-5, 3.504694e-5, 3.832983e-5, 3.915016e-5, 4.259103e-5, 4.329407e-5]
+        expected_b = [
+            2.242553e-5, 2.461508e-5, 2.242553e-5, 2.317305e-5, 2.242553e-5, 2.317305e-5,
+            2.242553e-5, 2.242553e-5, 2.317305e-5, 2.242553e-5, 2.317305e-5, 2.242553e-5]
+        expected_p = [
+            4.065860e-5, 4.868549e-5, 4.805108e-5, 5.347222e-5, 5.544355e-5, 6.111111e-5,
+            6.283602e-5, 6.653226e-5, 7.256944e-5, 7.392473e-5, 8.020833e-5, 8.131720e-5]
+        for expected_val, col_name in [(expected_qf, 'quickflow'),
+                (expected_b, 'baseflow'), (expected_p, 'precipitation')]:
+            numpy.testing.assert_allclose(expected_val, actual_result_df[col_name],
+                                          rtol=1e-6)
+
     def test_bad_biophysical_table(self):
         """SWY bad biophysical table with non-numerical values."""
         from natcap.invest.seasonal_water_yield import seasonal_water_yield
@@ -799,6 +850,23 @@ def test_monthly_alpha_regression(self):
             os.path.join(args['workspace_dir'], 'aggregated_results_swy.shp'),
             agg_results_csv_path)
 
+        # check the values in the avg monthly quickflow baseflow precip csv
+        actual_result_df = pandas.read_csv(
+            os.path.join(args['workspace_dir'], 'monthly_quickflow_baseflow.csv'))
+        expected_qf = [
+            2.116894e-5, 2.541212e-5, 2.514573e-5, 2.805605e-5, 2.916745e-5, 3.223471e-5,
+            3.323317e-5, 3.528216e-5, 3.858638e-5, 3.941151e-5, 4.287460e-5, 4.358156e-5]
+        expected_b = [
+            2.276286e-5, 2.498535e-5, 2.276286e-5, 2.352162e-5, 2.276286e-5, 2.352162e-5,
+            2.276286e-5, 2.276286e-5, 2.352162e-5, 2.276286e-5, 2.352162e-5, 2.276286e-5]
+        expected_p = [
+            4.106930e-5, 4.917726e-5, 4.853644e-5, 5.401235e-5, 5.600358e-5, 6.172840e-5,
+            6.347073e-5, 6.720430e-5, 7.330247e-5, 7.467145e-5, 8.101852e-5, 8.213859e-5]
+        for expected_val, col_name in [(expected_qf, 'quickflow'),
+                (expected_b, 'baseflow'), (expected_p, 'precipitation')]:
+            numpy.testing.assert_allclose(expected_val, actual_result_df[col_name],
+                                          rtol=1e-6)
+
     def test_climate_zones_missing_cz_id(self):
         """SWY climate zone regression test fails on bad cz table data.
 
@@ -932,7 +1000,14 @@ def test_user_recharge(self):
         make_recharge_raster(recharge_ras_path)
         args['l_path'] = recharge_ras_path
 
-        seasonal_water_yield.execute(args)
+        execute_kwargs = {
+            'generate_report': bool(seasonal_water_yield.MODEL_SPEC.reporter),
+            'save_file_registry': True,
+            'check_outputs': True
+        }
+        seasonal_water_yield.MODEL_SPEC.execute(args, **execute_kwargs)
+        assert_complete_execute(
+            args, seasonal_water_yield.MODEL_SPEC, **execute_kwargs)
 
         # generate aggregated results csv table for assertion
         agg_results_csv_path = os.path.join(args['workspace_dir'],

From 26753e342a7c66e0bf806607a2a51be0dc870e8d Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Wed, 18 Mar 2026 15:53:33 -0400
Subject: [PATCH 20/37] always pass model spec to `raster_inputs_summary`
 (#2321)

---
 src/natcap/invest/carbon/reporter.py                  | 2 +-
 src/natcap/invest/reports/raster_utils.py             | 5 ++---
 src/natcap/invest/reports/sdr_ndr_report_generator.py | 2 +-
 src/natcap/invest/seasonal_water_yield/reporter.py    | 2 +-
 4 files changed, 5 insertions(+), 6 deletions(-)

diff --git a/src/natcap/invest/carbon/reporter.py b/src/natcap/invest/carbon/reporter.py
index 52e01df3bb..115c773e09 100644
--- a/src/natcap/invest/carbon/reporter.py
+++ b/src/natcap/invest/carbon/reporter.py
@@ -201,7 +201,7 @@ def report(file_registry: dict, args_dict: dict, model_spec: ModelSpec,
     ]
 
     input_raster_stats_table = raster_utils.raster_inputs_summary(
-        args_dict).to_html(na_rep='')
+        args_dict, model_spec).to_html(na_rep='')
 
     output_raster_stats_table = raster_utils.raster_workspace_summary(
         file_registry).to_html(na_rep='')
diff --git a/src/natcap/invest/reports/raster_utils.py b/src/natcap/invest/reports/raster_utils.py
index 5986cbe1fa..4f6ab68f86 100644
--- a/src/natcap/invest/reports/raster_utils.py
+++ b/src/natcap/invest/reports/raster_utils.py
@@ -674,15 +674,14 @@ def _summarize_output(path):
     return pandas.DataFrame(raster_summary).T
 
 
-def raster_inputs_summary(args_dict, model_spec=None, check_csv_for_rasters=False):
+def raster_inputs_summary(args_dict, model_spec):
     """Create a table of stats for all rasters in an args_dict."""
     raster_summary = {}
 
     paths_to_check = [v for v in args_dict.values()
                       if isinstance(v, str) and os.path.isfile(v)]
 
-    if model_spec and check_csv_for_rasters:
-        paths_to_check.extend(_parse_csv_paths_from_spec(args_dict, model_spec))
+    paths_to_check.extend(_parse_csv_paths_from_spec(args_dict, model_spec))
 
     for v in paths_to_check:
         resource = geometamaker.describe(v, compute_stats=True)
diff --git a/src/natcap/invest/reports/sdr_ndr_report_generator.py b/src/natcap/invest/reports/sdr_ndr_report_generator.py
index 4057db96c3..5a709171a4 100644
--- a/src/natcap/invest/reports/sdr_ndr_report_generator.py
+++ b/src/natcap/invest/reports/sdr_ndr_report_generator.py
@@ -69,7 +69,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath,
         file_registry).to_html(na_rep='')
 
     input_raster_stats_table = raster_utils.raster_inputs_summary(
-        args_dict).to_html(na_rep='')
+        args_dict, model_spec).to_html(na_rep='')
 
     model_description = model_spec.about
     model_description += gettext(
diff --git a/src/natcap/invest/seasonal_water_yield/reporter.py b/src/natcap/invest/seasonal_water_yield/reporter.py
index c1c2506ca1..114ecee2e0 100644
--- a/src/natcap/invest/seasonal_water_yield/reporter.py
+++ b/src/natcap/invest/seasonal_water_yield/reporter.py
@@ -328,7 +328,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
             file_registry).to_html(na_rep='')
 
     input_raster_stats_table = raster_utils.raster_inputs_summary(
-            args_dict, model_spec, check_csv_for_rasters=True).to_html(na_rep='')
+            args_dict, model_spec).to_html(na_rep='')
 
     inputs_img_src = raster_utils.plot_and_base64_encode_rasters(
             input_raster_config_list)

From cecba2eb4bd4b10f4bc8e9ae3406f284e292ddb7 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Thu, 19 Mar 2026 10:21:51 -0400
Subject: [PATCH 21/37] swy report template testing (#2321)

---
 tests/reports/test_swy_template.py | 105 +++++++++++++++++++++++++++++
 1 file changed, 105 insertions(+)
 create mode 100644 tests/reports/test_swy_template.py

diff --git a/tests/reports/test_swy_template.py b/tests/reports/test_swy_template.py
new file mode 100644
index 0000000000..24f4d24b24
--- /dev/null
+++ b/tests/reports/test_swy_template.py
@@ -0,0 +1,105 @@
+import time
+import unittest
+
+import pandas
+from bs4 import BeautifulSoup
+
+from natcap.invest.seasonal_water_yield import MODEL_SPEC
+from natcap.invest.reports import jinja_env
+
+TEMPLATE = jinja_env.get_template('models/seasonal_water_yield.html')
+
+BSOUP_HTML_PARSER = 'html.parser'
+
+
+def _get_render_args(model_spec):
+    report_filepath = 'swy_report_test.html'
+    invest_version = '987.65.0'
+    timestamp = time.strftime('%Y-%m-%d %H:%M')
+    args_dict = {'suffix': 'test'}
+    img_src = 'bAse64eNcoDEdIMagE'
+    output_stats_table = '<table class="test__output-stats-table"></table>'
+    input_stats_table = '<table class="test__input-stats-table"></table>'
+    stats_table_note = 'This is a test!'
+    raster_group_caption = 'This is another test!'
+    stream_caption = ['stream.tif:Stream map.']
+    heading = 'Test heading'
+    outputs_caption = ['results.tif:Results map.']
+    inputs_caption = ['input.tif:Input map.']
+    vegalite_json = '{}'
+    caption = 'figure caption'
+    agg_map_source_list = ['/source/file.shp']
+
+    return {
+        'report_script': model_spec.reporter,
+        'invest_version': invest_version,
+        'report_filepath': report_filepath,
+        'model_id': model_spec.model_id,
+        'model_name': model_spec.model_title,
+        'model_description': model_spec.about,
+        'userguide_page': model_spec.userguide,
+        'timestamp': timestamp,
+        'args_dict': args_dict,
+        'raster_group_caption': raster_group_caption,
+        'stats_table_note': stats_table_note,
+        'stream_img_src': img_src,
+        'stream_caption': stream_caption,
+        'stream_outputs_heading': heading,
+        'outputs_heading': heading,
+        'outputs_img_src': img_src,
+        'outputs_caption': outputs_caption,
+        'qf_rasters': None,
+        'output_raster_stats_table': output_stats_table,
+        'input_raster_stats_table': input_stats_table,
+        'inputs_img_src': img_src,
+        'inputs_caption': inputs_caption,
+        'qf_b_charts': None,
+        'qb_map_json': vegalite_json,
+        'qb_map_caption': caption,
+        'vri_sum_map_json': vegalite_json,
+        'vri_sum_map_caption': caption,
+        'aggregate_map_source_list': agg_map_source_list,
+        'model_spec_outputs': model_spec.outputs
+    }
+
+class SeasonalWaterYieldTemplateTests(unittest.TestCase):
+    """Unit tests for SWY template."""
+
+    def test_render_without_user_defined_recharge(self):
+        """Test report rendering when user_defined_local_recharge is False."""
+
+        render_args = _get_render_args(MODEL_SPEC)
+
+        render_args['qf_rasters'] = {
+            'annual_qf_img_src': 'bAse64eNcoDEdIMagE',
+            'monthly_qf_img_src': 'bAse64eNcoDEdIMagE',
+            'qf_caption': 'test caption'
+        }
+        render_args['qf_b_charts'] = {
+            'json': '{}',
+            'caption': 'test caption',
+            'sources': ['/src/path.shp', '/src/path.csv']
+        }
+        html = TEMPLATE.render(render_args)
+        soup = BeautifulSoup(html, BSOUP_HTML_PARSER)
+
+        sections = soup.find_all(class_='accordion-section')
+        # Includes quickflow raster section and monthly qf + b charts section
+        self.assertEqual(len(sections), 10)
+
+        self.assertEqual(
+            soup.h1.string, f'InVEST Results: {MODEL_SPEC.model_title}')
+
+    def test_render_with_user_defined_recharge(self):
+        """Test report rendering when user_defined_local_recharge is True."""
+
+        render_args = _get_render_args(MODEL_SPEC)
+        html = TEMPLATE.render(render_args)
+        soup = BeautifulSoup(html, BSOUP_HTML_PARSER)
+
+        sections = soup.find_all(class_='accordion-section')
+        # No quickflow raster section or monthly qf + b charts section
+        self.assertEqual(len(sections), 8)
+
+        self.assertEqual(
+            soup.h1.string, f'InVEST Results: {MODEL_SPEC.model_title}')

From 375e96c691e4b9ccb7d7d0831100a89da63f1d1f Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Thu, 19 Mar 2026 16:38:54 -0400
Subject: [PATCH 22/37] accept either str or colormap object for
 RasterPlotConfig.colormap; unified colorbar for raster facet plots (#2321)

---
 src/natcap/invest/reports/raster_utils.py | 36 +++++++++--------------
 1 file changed, 14 insertions(+), 22 deletions(-)

diff --git a/src/natcap/invest/reports/raster_utils.py b/src/natcap/invest/reports/raster_utils.py
index 4f6ab68f86..9008d18d3e 100644
--- a/src/natcap/invest/reports/raster_utils.py
+++ b/src/natcap/invest/reports/raster_utils.py
@@ -6,6 +6,7 @@
 import os
 from io import BytesIO
 from enum import Enum
+from typing import Union, Optional
 
 import distinctipy
 import geometamaker
@@ -158,8 +159,8 @@ class RasterPlotConfig:
     """For highly skewed data, a transformation can help reveal variation."""
     title: str | None = None
     """An optional plot title. If ``None``, the filename is used."""
-    custom_colormap: str | None = None
-    """The name of a custom matplotlib colormap."""
+    colormap: Optional[Union[str, object]] = None
+    """The string name of a registered matplotlib colormap or a colormap object."""
 
     def __post_init__(self):
         if self.title is None:
@@ -396,7 +397,7 @@ def plot_raster_list(raster_list: list[RasterPlotConfig]):
         colorbar_kwargs = {}
         imshow_kwargs['norm'] = transform
         imshow_kwargs['interpolation'] = 'none'
-        cmap = config.custom_colormap if config.custom_colormap else COLORMAPS[dtype]
+        cmap = plt.get_cmap(config.colormap if config.colormap else COLORMAPS[dtype])
         if dtype == 'divergent':
             if transform == 'log':
                 transform = matplotlib.colors.SymLogNorm(linthresh=0.03)
@@ -508,7 +509,7 @@ def plot_and_base64_encode_rasters(raster_list: list[RasterPlotConfig]) -> str:
 
 
 def plot_raster_facets(tif_list, datatype, transform=None, title_list=None,
-                       small_plots=False, custom_colormap=None, supertitle=None):
+                       small_plots=False, colormap=None, supertitle=None):
     """Plot a list of rasters that will all share a fixed colorscale.
 
     When all the rasters have the same shape and represent the same variable,
@@ -527,8 +528,9 @@ def plot_raster_facets(tif_list, datatype, transform=None, title_list=None,
         small_plots (bool): Defaults to False. If True, the typical number of
             columns calculated for plotting facets will be increased by 1,
             making the plots smaller so more can be viewed side-by-side.
-        custom_colormap (str): Optional name of a matplotlib colormap to use
-            in place of the default derived from the raster datatype.
+        colormap (str): Optional string name of a registered matplotlib
+            colormap or a colormap object to use in place of the default
+            derived from the raster datatype.
         supertitle (str): Optional title to use for the entire group of
             raster facets.
 
@@ -539,7 +541,6 @@ def plot_raster_facets(tif_list, datatype, transform=None, title_list=None,
     xy_ratio = _get_aspect_ratio(bbox)
     fig, axes = _figure_subplots(xy_ratio, n_plots, small_plots=small_plots)
 
-    cmap_str = custom_colormap if custom_colormap else COLORMAPS[datatype]
     if transform is None:
         transform = 'linear'
     if title_list is None:
@@ -563,7 +564,7 @@ def plot_raster_facets(tif_list, datatype, transform=None, title_list=None,
     # instead of storing all arrays in memory
     vmin = numpy.nanmin(ndarray)
     vmax = numpy.nanmax(ndarray)
-    cmap = plt.get_cmap(cmap_str)
+    cmap = plt.get_cmap(colormap if colormap else COLORMAPS[datatype])
     if datatype == 'divergent':
         if transform == 'log':
             normalizer = matplotlib.colors.SymLogNorm(linthresh=0.03, vmin=vmin, vmax=vmax)
@@ -581,22 +582,13 @@ def plot_raster_facets(tif_list, datatype, transform=None, title_list=None,
         # all rasters are identical size; `resampled` will be the same for all
         title_line_width = _get_title_line_width(n_plots, xy_ratio)
         ax.set_title(**_get_title_kwargs(title, resampled, title_line_width))
-        if not small_plots:
-            units = _get_raster_units(raster_path)
-            if units:
-                (ylim_kwargs,
-                 text_kwargs) = _get_units_text_kwargs(units, len(arr))
-                ax.set_ylim(**ylim_kwargs)
-                ax.text(**text_kwargs)
-            fig.colorbar(mappable, ax=ax)
 
-    if small_plots:
-        units = _get_raster_units(tif_list[0])
-        legend_label = f"{UNITS_TEXT}: {units}" if units else None
-        fig.colorbar(mappable, ax=axes.ravel().tolist(), label=legend_label)
+    units = _get_raster_units(tif_list[0])
+    legend_label = f"{UNITS_TEXT}: {units}" if units else None
+    fig.colorbar(mappable, ax=axes.ravel().tolist(), label=legend_label)
 
     if supertitle:
-        fig.suptitle(f"{supertitle}", fontsize=TITLE_FONT_SIZE)
+        fig.suptitle(supertitle, fontsize=TITLE_FONT_SIZE)
 
     [ax.set_axis_off() for ax in axes.flatten()]
     return fig
@@ -702,7 +694,7 @@ def _parse_csv_paths_from_spec(args_dict, spec):
     for input_ in spec.inputs:
         if isinstance(input_, CSVInput):
             table_map_inputs.extend([
-                (input_.id, col.id) for col in spec.get_input(input_.id).columns
+                (input_.id, col.id) for col in input_.columns
                 if isinstance(col, SingleBandRasterInput)])
 
     paths_to_check = []

From 7e51bdd94c89b7c2491716de6140abed9ed5b75d Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Thu, 19 Mar 2026 16:41:27 -0400
Subject: [PATCH 23/37] improvements to SWY reporter, including better axis and
 legend config for linked chart+map (#2321)

---
 .../models/seasonal_water_yield.html          |  6 +-
 .../invest/seasonal_water_yield/reporter.py   | 79 +++++++++++--------
 2 files changed, 47 insertions(+), 38 deletions(-)

diff --git a/src/natcap/invest/reports/templates/models/seasonal_water_yield.html b/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
index 4367d097ba..ceb96ba281 100644
--- a/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
+++ b/src/natcap/invest/reports/templates/models/seasonal_water_yield.html
@@ -26,10 +26,8 @@ <h2 class="section-header">Results</h2>
     {{ accordion_section(
       'Average Quickflow and Baseflow by Month',
       content_grid([
-        (content_grid([
-          ('<div id="qf_b_charts"></div>', 100),
-          (caption(qf_b_charts['caption'], qf_b_charts['sources']), 100)
-        ]), 100)
+        ('<div id="qf_b_charts"></div>', 100),
+        (caption(qf_b_charts['caption'], qf_b_charts['sources']), 100)
       ])
     ) }}
   {% endif %}
diff --git a/src/natcap/invest/seasonal_water_yield/reporter.py b/src/natcap/invest/seasonal_water_yield/reporter.py
index 114ecee2e0..1a201313b8 100644
--- a/src/natcap/invest/seasonal_water_yield/reporter.py
+++ b/src/natcap/invest/seasonal_water_yield/reporter.py
@@ -32,7 +32,7 @@
 MAP_WIDTH = 450 # pixels
 
 qf_label_month_map = {
-    f"qf_{month_index+1}": str(month_index+1) for month_index in range(12)
+    f"qf_{month_index}": str(month_index) for month_index in range(1, 13)
 }
 
 
@@ -40,16 +40,14 @@ def _label_to_month(row):
     return qf_label_month_map[row['MonthLabel']]
 
 
-def _create_aggregate_map(geodataframe, extent_feature, xy_ratio, attribute,
+def _create_aggregate_map(geodataframe, xy_ratio, attribute,
                           title):
     attr_map = altair.Chart(geodataframe).mark_geoshape(
-        clip=True,
         stroke="white",
         strokeWidth=0.5
     ).project(
         type='identity',
-        reflectY=True,
-        fit=extent_feature
+        reflectY=True
     ).encode(
         color=attribute,
         tooltip=[altair.Tooltip(attribute, title=attribute)]
@@ -62,31 +60,29 @@ def _create_aggregate_map(geodataframe, extent_feature, xy_ratio, attribute,
     return attr_map.to_json()
 
 
-def _create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
+def _create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path, xy_ratio):
     map_df = geopandas.read_file(aoi_vector_path)
     values_df = pandas.read_csv(aggregate_csv_path)
     values_df.month = values_df.month.astype(str)
 
-    extent_feature, xy_ratio = vector_utils.get_geojson_bbox(map_df)
-
     feat_select = altair.selection_point(fields=["geom_id"], name="feat_select", value=0)
 
     attr_map = altair.Chart(map_df).mark_geoshape(
-        clip=True, stroke="white", strokeWidth=0.5
-        ).project(
-            type='identity',
-            reflectY=True,
-            fit=extent_feature
+        stroke="white",
+        strokeWidth=0.5
+    ).project(
+        type='identity',
+        reflectY=True
     ).encode(
         color=altair.condition(
             feat_select,
             altair.value("seagreen"),
             altair.value("lightgray")
         ),
-        tooltip=[altair.Tooltip("geom_id", title="FID")]
+        tooltip=[altair.Tooltip("geom_id", title="Feature")]
     ).properties(
-        width=MAP_WIDTH,
-        height=MAP_WIDTH / xy_ratio,
+        width=MAP_WIDTH*1.25,
+        height=MAP_WIDTH*1.25 / xy_ratio,
         title="AOI"
     ).add_params(
         feat_select
@@ -124,11 +120,23 @@ def _create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path):
         feat_select
     ).properties(
         title=altair.Title(altair.expr(
-            f'"Mean Quickflow + Baseflow for Feature, FID " + {feat_select.name}.geom_id')
+            f'"Mean Quickflow + Baseflow for Feature " + {feat_select.name}.geom_id')
         )
     )
 
-    chart = attr_map | combined_chart
+    legend_config = vector_utils.LEGEND_CONFIG
+    legend_config['orient'] = 'right'
+
+    chart = altair.hconcat(
+        attr_map, combined_chart,
+        spacing=30
+    ).configure_legend(
+        **legend_config
+    ).configure_axis(
+        **vector_utils.AXIS_CONFIG
+    ).configure_view(
+        discreteWidth=300
+    )
     return chart.to_json()
 
 
@@ -151,10 +159,10 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
 
     # qb and vri_sum plots from the output aggregate vector
     aggregated_results = geopandas.read_file(file_registry['aggregate_vector'])
-    extent_feature, xy_ratio = vector_utils.get_geojson_bbox(aggregated_results)
+    _, xy_ratio = vector_utils.get_geojson_bbox(aggregated_results)
 
     qb_map_json = _create_aggregate_map(
-        aggregated_results, extent_feature, xy_ratio, 'qb',
+        aggregated_results, xy_ratio, 'qb',
         gettext("Mean local recharge value within the watershed "
                 f"({model_spec.get_output('aggregate_vector').get_field('qb').units})"))
     qb_map_caption = [
@@ -163,7 +171,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
                 'relative values, not absolute values.')]
 
     vri_sum_map_json = _create_aggregate_map(
-        aggregated_results, extent_feature, xy_ratio, 'vri_sum',
+        aggregated_results, xy_ratio, 'vri_sum',
         gettext("Total recharge contribution of the watershed "
                 f"({model_spec.get_output('aggregate_vector').get_field('vri_sum').units})"))
     vri_sum_map_caption = [
@@ -180,13 +188,15 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
         # Monthly quickflow + baseflow plots and map
         qf_b_charts_json = _create_linked_monthly_plots(
                 file_registry['aggregate_vector'],
-                file_registry['monthly_qf_table'])
+                file_registry['monthly_qf_table'],
+                xy_ratio)
         qf_b_charts_caption = gettext(
             """
             This chart displays the monthly combined average baseflow + quickflow for
-            each feature within the AOI. Select a feature on the AOI map to see the
-            values for that feature. Shift+Click to select multiple features; the chart
-            will display the sum of their values.
+            each feature within the AOI, as well as the monthly average precipitation.
+            Select a feature on the AOI map to see the values for that feature.
+            Shift+Click to select multiple features; the chart will display the sum of
+            their values.
             """
         )
         qf_b_charts_source_list = [
@@ -231,7 +241,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
         input_raster_config_list.append(
             RasterPlotConfig(
                 raster_path=args_dict['l_path'],
-                datatype=RasterDatatype.continuous,
+                datatype=RasterDatatype.divergent,
                 spec=model_spec.get_input('l_path')))
         qf_rasters = None
         raster_outputs_heading = 'Annual Baseflow'
@@ -252,7 +262,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
                 spec=model_spec.get_output('cn')),
             RasterPlotConfig(
                 raster_path=file_registry['l'],
-                datatype=RasterDatatype.continuous,
+                datatype=RasterDatatype.divergent,
                 spec=model_spec.get_output('l'))])
         raster_outputs_heading = 'Additional Raster Outputs'
 
@@ -269,7 +279,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
                 spec=model_spec.get_output('qf'),
                 title=gettext("Annual Quickflow"),
                 transform=RasterTransform.log,
-                custom_colormap='GnBu')
+                colormap='GnBu')
 
         monthly_qf_raster_config_list = [
             RasterPlotConfig(
@@ -278,7 +288,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
                 spec=model_spec.get_output('qf_[MONTH]'),
                 title=gettext(f"{calendar.month_name[month]}"),
                 transform=RasterTransform.log,
-                custom_colormap='GnBu'
+                colormap='GnBu'
             ) for month in range(1, 13)]
 
         annual_qf_img_src = raster_utils.plot_and_base64_encode_rasters(
@@ -291,16 +301,17 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
                 title_list=[raster_config.title for raster_config
                             in monthly_qf_raster_config_list],
                 small_plots=True,
-                custom_colormap='GnBu',
+                colormap='GnBu',
                 supertitle=gettext("Monthly Quickflow"))
+        annual_qf_displayname = os.path.basename(
+                annual_qf_raster_config.raster_path)
         monthly_qf_img_src = raster_utils.base64_encode(monthly_qf_plots)
         monthly_qf_displayname = os.path.basename(
                 monthly_qf_raster_config_list[0].raster_path).replace('1', '[MONTH]')
         qf_raster_caption = [
-            (f'{annual_qf_raster_config.title}:'
-             f'{annual_qf_raster_config.spec.about}'),
-            (f'{monthly_qf_displayname}:'
-             f'{monthly_qf_raster_config_list[0].spec.about}')
+                gettext(f'Map of Annual Quickflow: {annual_qf_displayname}'),
+                gettext(f'Maps of Monthly Quickflow: {monthly_qf_displayname}'
+                        ' (1 = January… 12 = December)')
         ]
 
         qf_rasters = {

From 93f9c736568950c7f4e69e5d0bb264776f623908 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Thu, 19 Mar 2026 16:45:52 -0400
Subject: [PATCH 24/37] document geom_id in aggregate vector output spec;
 adjust gemon_id description in csv output spec (#2321)

---
 .../seasonal_water_yield/seasonal_water_yield.py  | 15 +++++++++++----
 1 file changed, 11 insertions(+), 4 deletions(-)

diff --git a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
index 406ac684c1..a400e90e53 100644
--- a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
+++ b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
@@ -452,11 +452,19 @@
                 spec.NumberOutput(
                     id="vri_sum",
                     about=gettext(
-                        "Total recharge contribution, (positive or negative) within the"
+                        "Total recharge contribution (positive or negative) within the"
                         " watershed."
                     ),
                     units=u.millimeter / u.year
+                ),
+                spec.IntegerOutput(
+                    id="geom_id",
+                    about=gettext(
+                        "A unique ID for the watershed."
+                    ),
+                    units=u.none
                 )
+
             ]
         ),
         spec.CSVOutput(
@@ -470,8 +478,8 @@
                 spec.IntegerOutput(
                     id="geom_id",
                     about=gettext(
-                        "The Feature ID for the watershed. This will correspond to"
-                        " the FIDs in the Aggregate Results shapefile."
+                        "A unique ID for the watershed. This will correspond to"
+                        " the 'geom_id' column in the Aggregate Results shapefile."
                     ),
                     units=u.none
                 ),
@@ -1514,7 +1522,6 @@ def _aggregate_recharge(
     for raster_path, aggregate_field_id, op_type in [
             (l_path, 'qb', 'mean'), (vri_path, 'vri_sum', 'sum')]:
 
-        # aggregate carbon stocks by the new ID field
         aggregate_stats = pygeoprocessing.zonal_statistics(
             (raster_path, 1), aggregate_vector_path)
 

From 8c45a1a2403eef1ff344951083191c2a3cf74dfc Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Thu, 19 Mar 2026 16:46:55 -0400
Subject: [PATCH 25/37] update SWY model to index all monthly outputs
 consistently by month, where January=1 (#2321)

---
 .../seasonal_water_yield.py                   | 70 +++++++++----------
 tests/test_seasonal_water_yield_regression.py |  7 +-
 2 files changed, 38 insertions(+), 39 deletions(-)

diff --git a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
index a400e90e53..203b1e449a 100644
--- a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
+++ b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
@@ -22,6 +22,7 @@
 
 TARGET_NODATA = -1
 N_MONTHS = 12
+MONTH_RANGE = range(1, N_MONTHS+1)
 MONTH_ID_TO_LABEL = [
     'jan', 'feb', 'mar', 'apr', 'may', 'jun', 'jul', 'aug', 'sep', 'oct',
     'nov', 'dec']
@@ -754,7 +755,7 @@ def execute(args):
         # make all 12 entries equal to args['alpha_m']
         alpha_m = float(fractions.Fraction(args['alpha_m']))
         alpha_month_map = dict(
-            (month_index+1, alpha_m) for month_index in range(N_MONTHS))
+            (month_index, alpha_m) for month_index in MONTH_RANGE)
 
     beta_i = float(fractions.Fraction(args['beta_i']))
     gamma = float(fractions.Fraction(args['gamma']))
@@ -763,23 +764,22 @@ def execute(args):
         args['dem_raster_path'])['pixel_size']
 
     LOGGER.info('Checking that the AOI is not the output aggregate vector')
-    LOGGER.debug("aoi_path: %s", args['aoi_path'])
-    LOGGER.debug("aggregate_vector_path: %s",
-                 os.path.normpath(file_registry['aggregate_vector']))
+    LOGGER.debug(f"aoi_path: {args['aoi_path']}")
+    LOGGER.debug("aggregate_vector_path: "
+                 f"{os.path.normpath(file_registry['aggregate_vector'])}")
     if (os.path.normpath(args['aoi_path']) ==
             os.path.normpath(file_registry['aggregate_vector'])):
         raise ValueError(
             "The input AOI is the same as the output aggregate vector, "
             "please choose a different workspace or move the AOI file "
-            "out of the current workspace %s" %
-            file_registry['aggregate_vector'])
+            f"out of the current workspace {file_registry['aggregate_vector']}")
 
     LOGGER.info('Aligning and clipping dataset list')
     input_align_list = [args['lulc_raster_path'], args['dem_raster_path']]
     output_align_list = [
         file_registry['lulc_aligned'], file_registry['dem_aligned']]
     if not args['user_defined_local_recharge']:
-        month_indexes = [m+1 for m in range(N_MONTHS)]
+        month_indexes = [m for m in MONTH_RANGE]
 
         precip_df = MODEL_SPEC.get_input(
             'precip_raster_table').get_validated_dataframe(
@@ -804,9 +804,9 @@ def execute(args):
             precip_path_list + [args['soil_group_path']] + et0_path_list +
             input_align_list)
         output_align_list = (
-            [file_registry['prcp_a[MONTH]', month] for month in range(12)] +
+            [file_registry['prcp_a[MONTH]', month] for month in MONTH_RANGE] +
             [file_registry['soil_group_aligned']] +
-            [file_registry['et0_a[MONTH]', month] for month in range(12)] +
+            [file_registry['et0_a[MONTH]', month] for month in MONTH_RANGE] +
             output_align_list)
 
     align_index = len(input_align_list) - 1  # this aligns with the DEM
@@ -918,7 +918,7 @@ def execute(args):
         reclass_error_details = {
             'raster_name': 'Climate Zone', 'column_name': 'cz_id',
             'table_name': 'Climate Zone'}
-        for month_id in range(N_MONTHS):
+        for month_id in MONTH_RANGE:
             if args['user_defined_climate_zones']:
                 cz_rain_events_df = MODEL_SPEC.get_input(
                     'climate_zone_table_path').get_validated_dataframe(
@@ -936,7 +936,7 @@ def execute(args):
                     target_path_list=[
                         file_registry['n_events[MONTH]', month_id]],
                     dependent_task_list=[align_task],
-                    task_name='n_events for month %d' % month_id)
+                    task_name=f'n_events for month {month_id}')
                 reclassify_n_events_task_list.append(n_events_task)
             else:
                 n_events_task = task_graph.add_task(
@@ -946,12 +946,12 @@ def execute(args):
                         file_registry['n_events[MONTH]', month_id],
                         gdal.GDT_Float32, [TARGET_NODATA]),
                     kwargs={'fill_value_list': (
-                        rain_events_df['events'][month_id+1],)},
+                        rain_events_df['events'][month_id],)},
                     target_path_list=[
                         file_registry['n_events[MONTH]', month_id]],
                     dependent_task_list=[align_task],
                     task_name=(
-                        'n_events as a constant raster month %d' % month_id))
+                        f'n_events as a constant raster month {month_id}'))
                 reclassify_n_events_task_list.append(n_events_task)
 
         curve_number_task = task_graph.add_task(
@@ -975,8 +975,8 @@ def execute(args):
             task_name='calculate Si raster')
 
         quick_flow_task_list = []
-        for month_index in range(N_MONTHS):
-            LOGGER.info('calculate quick flow for month %d', month_index+1)
+        for month_index in MONTH_RANGE:
+            LOGGER.info(f'calculate quick flow for month {month_index}')
             monthly_quick_flow_task = task_graph.add_task(
                 func=_calculate_monthly_quick_flow,
                 args=(
@@ -984,21 +984,20 @@ def execute(args):
                     file_registry['n_events[MONTH]', month_index],
                     file_registry['stream'],
                     file_registry['si'],
-                    file_registry['qf_[MONTH]', month_index + 1]),
+                    file_registry['qf_[MONTH]', month_index]),
                 target_path_list=[
-                    file_registry['qf_[MONTH]', month_index + 1]],
+                    file_registry['qf_[MONTH]', month_index]],
                 dependent_task_list=[
-                    align_task, reclassify_n_events_task_list[month_index],
+                    align_task, reclassify_n_events_task_list[month_index-1],
                     si_task, stream_threshold_task],
-                task_name='calculate quick flow for month %d' % (
-                    month_index+1))
+                task_name=f'calculate quick flow for month {month_index}')
             quick_flow_task_list.append(monthly_quick_flow_task)
 
         qf_task = task_graph.add_task(
             func=pygeoprocessing.raster_map,
             kwargs=dict(
                 op=qfi_sum_op,
-                rasters=[file_registry['qf_[MONTH]', month] for month in range(1, 13)],
+                rasters=[file_registry['qf_[MONTH]', month] for month in MONTH_RANGE],
                 target_path=file_registry['qf']),
             target_path_list=[file_registry['qf']],
             dependent_task_list=quick_flow_task_list,
@@ -1009,8 +1008,8 @@ def execute(args):
         reclass_error_details = {
             'raster_name': 'LULC', 'column_name': 'lucode',
             'table_name': 'Biophysical'}
-        for month_index in range(N_MONTHS):
-            kc_lookup = biophysical_df['kc_%d' % (month_index+1)].to_dict()
+        for month_index in MONTH_RANGE:
+            kc_lookup = biophysical_df[f'kc_{month_index}'].to_dict()
             kc_task = task_graph.add_task(
                 func=utils.reclassify_raster,
                 args=(
@@ -1019,7 +1018,7 @@ def execute(args):
                     gdal.GDT_Float32, TARGET_NODATA, reclass_error_details),
                 target_path_list=[file_registry['kc_[MONTH]', month_index]],
                 dependent_task_list=[align_task],
-                task_name='classify kc month %d' % month_index)
+                task_name=f'classify kc month {month_index}')
             kc_task_list.append(kc_task)
 
         # call through to a cython function that does the necessary routing
@@ -1027,11 +1026,11 @@ def execute(args):
         calculate_local_recharge_task = task_graph.add_task(
             func=seasonal_water_yield_core.calculate_local_recharge,
             args=(
-                [file_registry['prcp_a[MONTH]', month] for month in range(12)],
-                [file_registry['et0_a[MONTH]', month] for month in range(12)],
-                [file_registry['qf_[MONTH]', month] for month in range(1, 13)],
+                [file_registry['prcp_a[MONTH]', month] for month in MONTH_RANGE],
+                [file_registry['et0_a[MONTH]', month] for month in MONTH_RANGE],
+                [file_registry['qf_[MONTH]', month] for month in MONTH_RANGE],
                 file_registry['flow_dir'],
-                [file_registry['kc_[MONTH]', month] for month in range(12)],
+                [file_registry['kc_[MONTH]', month] for month in MONTH_RANGE],
                 alpha_month_map,
                 beta_i, gamma, file_registry['stream'],
                 file_registry['l'],
@@ -1126,8 +1125,8 @@ def execute(args):
             args=(
                 file_registry['aggregate_vector'],
                 file_registry['b'],
-                [file_registry['qf_[MONTH]', month] for month in range(1, 13)],
-                [file_registry['prcp_a[MONTH]', month] for month in range(12)],
+                [file_registry['qf_[MONTH]', month] for month in MONTH_RANGE],
+                [file_registry['prcp_a[MONTH]', month] for month in MONTH_RANGE],
                 file_registry['monthly_qf_table']
             ),
             target_path_list=[file_registry['monthly_qf_table']],
@@ -1506,8 +1505,7 @@ def _aggregate_recharge(
     """
     if os.path.exists(aggregate_vector_path):
         LOGGER.warning(
-            '%s exists, deleting and writing new output',
-            aggregate_vector_path)
+            f'{aggregate_vector_path} exists, deleting and writing new output')
         os.remove(aggregate_vector_path)
 
     original_aoi_vector = gdal.OpenEx(aoi_path, gdal.OF_VECTOR)
@@ -1606,10 +1604,10 @@ def _generate_monthly_qf_b_p_csv(
     raster_path_tuples = [(annual_baseflow_path, 0, "baseflow")]
     raster_path_tuples.extend([
             (qf_path, month, "quickflow") for qf_path, month
-            in zip(monthly_quickflow_path_list, range(1, 13))])
+            in zip(monthly_quickflow_path_list, MONTH_RANGE)])
     raster_path_tuples.extend([
             (precip_path, month, "precipitation") for precip_path, month
-            in zip(monthly_precip_path_list, range(1, 13))])
+            in zip(monthly_precip_path_list, MONTH_RANGE)])
 
     raster_stats = pygeoprocessing.zonal_statistics(
             [(path_tuple[0], 1) for path_tuple in raster_path_tuples],
@@ -1628,11 +1626,11 @@ def _generate_monthly_qf_b_p_csv(
                                 for k, v in stats_by_field['baseflow'][0].items()}
 
     values_dict = {fid: {month: {'baseflow': b_val / seconds_per_month[month]}
-                         for month in range(1, 13)}
+                         for month in MONTH_RANGE}
                    for fid, b_val in b_avg_per_feat_per_month.items()}
 
     for value_name in ['quickflow', 'precipitation']:
-        for month in range(1, 13):
+        for month in MONTH_RANGE:
             avg_per_feat_per_month = {
                     k: v['sum'] * 0.001 * pixel_area_m2
                     for k, v in stats_by_field[value_name][month].items()}
diff --git a/tests/test_seasonal_water_yield_regression.py b/tests/test_seasonal_water_yield_regression.py
index c75b7f2bc5..5d44e08437 100644
--- a/tests/test_seasonal_water_yield_regression.py
+++ b/tests/test_seasonal_water_yield_regression.py
@@ -1337,10 +1337,11 @@ def test_local_recharge_undefined_nodata(self):
         target_aet_path = os.path.join(self.workspace_dir,
                                        'target_aet_path.tif')
 
+        month_range = range(1, 13)
         seasonal_water_yield_core.calculate_local_recharge(
-            [precip_path for i in range(12)], [et0_path for i in range(12)],
-            [quickflow_path for i in range(12)], flow_dir_path,
-            [kc_path for i in range(12)], alpha_month_map, beta,
+            [precip_path for i in month_range], [et0_path for i in month_range],
+            [quickflow_path for i in month_range], flow_dir_path,
+            [kc_path for i in month_range], alpha_month_map, beta,
             gamma, stream_path, target_li_path, target_li_avail_path,
             target_l_sum_avail_path, target_aet_path,
             os.path.join(self.workspace_dir, 'target_precip_path.tif'),

From 4c4a78969618923551e988022ac88c53660745fe Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Thu, 19 Mar 2026 16:48:15 -0400
Subject: [PATCH 26/37] reduce tests checking quickflow/baseflow CSV values to
 the base regression mfd and d8 cases (#2321)

---
 tests/test_seasonal_water_yield_regression.py | 36 -------------------
 1 file changed, 36 deletions(-)

diff --git a/tests/test_seasonal_water_yield_regression.py b/tests/test_seasonal_water_yield_regression.py
index 5d44e08437..3084f018de 100644
--- a/tests/test_seasonal_water_yield_regression.py
+++ b/tests/test_seasonal_water_yield_regression.py
@@ -505,12 +505,10 @@ class SeasonalWaterYieldRegressionTests(unittest.TestCase):
     def setUp(self):
         """Create temporary workspace."""
         self.workspace_dir = tempfile.mkdtemp()
-#        self.workspace_dir = "/Users/megannissel/dev/debugging/swy_reports/"
 
     def tearDown(self):
         """Remove temporary workspace."""
         shutil.rmtree(self.workspace_dir)
-#        pass
 
     @staticmethod
     def generate_base_args(workspace_dir):
@@ -780,23 +778,6 @@ def test_base_regression_nodata_inf(self):
             os.path.join(args['workspace_dir'], 'aggregated_results_swy.shp'),
             agg_results_csv_path)
 
-        # check the values in the avg monthly quickflow baseflow precip csv
-        actual_result_df = pandas.read_csv(
-            os.path.join(args['workspace_dir'], 'monthly_quickflow_baseflow.csv'))
-        expected_qf = [
-            2.102520e-5, 2.524000e-5, 2.497585e-5, 2.786701e-5, 2.897143e-5, 3.201866e-5,
-            3.301102e-5, 3.504694e-5, 3.832983e-5, 3.915016e-5, 4.259103e-5, 4.329407e-5]
-        expected_b = [
-            2.242553e-5, 2.461508e-5, 2.242553e-5, 2.317305e-5, 2.242553e-5, 2.317305e-5,
-            2.242553e-5, 2.242553e-5, 2.317305e-5, 2.242553e-5, 2.317305e-5, 2.242553e-5]
-        expected_p = [
-            4.065860e-5, 4.868549e-5, 4.805108e-5, 5.347222e-5, 5.544355e-5, 6.111111e-5,
-            6.283602e-5, 6.653226e-5, 7.256944e-5, 7.392473e-5, 8.020833e-5, 8.131720e-5]
-        for expected_val, col_name in [(expected_qf, 'quickflow'),
-                (expected_b, 'baseflow'), (expected_p, 'precipitation')]:
-            numpy.testing.assert_allclose(expected_val, actual_result_df[col_name],
-                                          rtol=1e-6)
-
     def test_bad_biophysical_table(self):
         """SWY bad biophysical table with non-numerical values."""
         from natcap.invest.seasonal_water_yield import seasonal_water_yield
@@ -850,23 +831,6 @@ def test_monthly_alpha_regression(self):
             os.path.join(args['workspace_dir'], 'aggregated_results_swy.shp'),
             agg_results_csv_path)
 
-        # check the values in the avg monthly quickflow baseflow precip csv
-        actual_result_df = pandas.read_csv(
-            os.path.join(args['workspace_dir'], 'monthly_quickflow_baseflow.csv'))
-        expected_qf = [
-            2.116894e-5, 2.541212e-5, 2.514573e-5, 2.805605e-5, 2.916745e-5, 3.223471e-5,
-            3.323317e-5, 3.528216e-5, 3.858638e-5, 3.941151e-5, 4.287460e-5, 4.358156e-5]
-        expected_b = [
-            2.276286e-5, 2.498535e-5, 2.276286e-5, 2.352162e-5, 2.276286e-5, 2.352162e-5,
-            2.276286e-5, 2.276286e-5, 2.352162e-5, 2.276286e-5, 2.352162e-5, 2.276286e-5]
-        expected_p = [
-            4.106930e-5, 4.917726e-5, 4.853644e-5, 5.401235e-5, 5.600358e-5, 6.172840e-5,
-            6.347073e-5, 6.720430e-5, 7.330247e-5, 7.467145e-5, 8.101852e-5, 8.213859e-5]
-        for expected_val, col_name in [(expected_qf, 'quickflow'),
-                (expected_b, 'baseflow'), (expected_p, 'precipitation')]:
-            numpy.testing.assert_allclose(expected_val, actual_result_df[col_name],
-                                          rtol=1e-6)
-
     def test_climate_zones_missing_cz_id(self):
         """SWY climate zone regression test fails on bad cz table data.
 

From 5e78639f586ba195ad833295dd4fe22398367b66 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Fri, 20 Mar 2026 11:25:31 -0400
Subject: [PATCH 27/37] snapshot tests for plot_raster_facets small plots
 (#2321)

---
 src/natcap/invest/reports/raster_utils.py |  9 ++--
 tests/reports/test_raster_utils.py        | 52 +++++++++++++++++++++++
 2 files changed, 58 insertions(+), 3 deletions(-)

diff --git a/src/natcap/invest/reports/raster_utils.py b/src/natcap/invest/reports/raster_utils.py
index 9008d18d3e..3be4bede49 100644
--- a/src/natcap/invest/reports/raster_utils.py
+++ b/src/natcap/invest/reports/raster_utils.py
@@ -313,16 +313,19 @@ def _get_title_line_width(n_plots: int, xy_ratio: float) -> int:
         return 31  # 3-column layout
 
 
-def _get_title_kwargs(title: str, resampled: bool, line_width: int):
+def _get_title_kwargs(title: str, resampled: bool, line_width: int, facets=False):
     label = f"{title}{' (resampled)' if resampled else ''}"
     label = textwrap.fill(label, width=line_width)
+    padding = 1.5
+    if not facets:
+        padding *= SUBTITLE_FONT_SIZE
     return {
         'fontfamily': 'monospace',
         'fontsize': TITLE_FONT_SIZE,
         'fontweight': 700,
         'label': label,
         'loc': 'left',
-        'pad': 1.5 * SUBTITLE_FONT_SIZE,
+        'pad': padding,
         'verticalalignment': 'bottom',
     }
 
@@ -581,7 +584,7 @@ def plot_raster_facets(tif_list, datatype, transform=None, title_list=None,
         mappable = ax.imshow(arr, cmap=cmap, norm=normalizer)
         # all rasters are identical size; `resampled` will be the same for all
         title_line_width = _get_title_line_width(n_plots, xy_ratio)
-        ax.set_title(**_get_title_kwargs(title, resampled, title_line_width))
+        ax.set_title(**_get_title_kwargs(title, resampled, title_line_width, facets=True))
 
     units = _get_raster_units(tif_list[0])
     legend_label = f"{UNITS_TEXT}: {units}" if units else None
diff --git a/tests/reports/test_raster_utils.py b/tests/reports/test_raster_utils.py
index 66469782e6..20462e63e3 100644
--- a/tests/reports/test_raster_utils.py
+++ b/tests/reports/test_raster_utils.py
@@ -339,6 +339,24 @@ def tearDown(self):
         """Override tearDown function to remove temporary directory."""
         shutil.rmtree(self.workspace_dir)
 
+    @staticmethod
+    @patch('natcap.invest.reports.raster_utils._get_raster_units')
+    def create_small_plots_grid(workspace_dir, shape, mock_get_raster_units,
+                                supertitle=None):
+        raster_paths = [os.path.join(workspace_dir, f'{s}.tif')
+                        for s in ['a', 'b', 'c', 'd']]
+        arrays = [numpy.linspace(
+            i, i+1, num=numpy.multiply(*shape)).reshape(*shape) for i in range(4)]
+        for raster_path, array in zip(raster_paths, arrays):
+            pygeoprocessing.numpy_array_to_raster(
+                array, target_nodata=None, pixel_size=(1, 1), origin=(0, 0),
+                projection_wkt=PROJ_WKT, target_path=raster_path)
+
+        mock_get_raster_units.return_value = 'flux capacitrons'
+        return raster_utils.plot_raster_facets(
+            raster_paths, RasterDatatype.continuous, small_plots=True,
+            supertitle=supertitle)
+
     def test_plot_raster_facets(self):
         """Test rasters share a common colorscale."""
         figname = 'plot_raster_facets.png'
@@ -363,6 +381,40 @@ def test_plot_raster_facets(self):
         save_figure(fig, actual_png)
         compare_snapshots(reference, actual_png)
 
+    def test_plot_raster_facets_small_plots(self):
+        """Test small plots: standard AOI width should have 4 columns."""
+        figname = 'plot_raster_facets_small_plots.png'
+        reference = os.path.join(REFS_DIR, figname)
+        shape = (4, 4)
+        fig = self.create_small_plots_grid(self.workspace_dir, shape)
+
+        actual_png = os.path.join(self.workspace_dir, figname)
+        save_figure(fig, actual_png)
+        compare_snapshots(reference, actual_png)
+
+    def test_plot_raster_facets_small_plots_wide_aoi(self):
+        """Test small plots: wide AOI width should have 3 columns."""
+        figname = 'plot_raster_facets_small_plots_wide_aoi.png'
+        reference = os.path.join(REFS_DIR, figname)
+        shape = (6, 12)
+        fig = self.create_small_plots_grid(self.workspace_dir, shape)
+
+        actual_png = os.path.join(self.workspace_dir, figname)
+        save_figure(fig, actual_png)
+        compare_snapshots(reference, actual_png)
+
+    def test_plot_raster_facets_small_plots_supertitle(self):
+        """Test small plots with optional supertitle."""
+        figname = 'plot_raster_facets_small_plots_supertitle.png'
+        reference = os.path.join(REFS_DIR, figname)
+        shape = (2, 10)
+        fig = self.create_small_plots_grid(self.workspace_dir, shape,
+                                           supertitle="Custom Title")
+
+        actual_png = os.path.join(self.workspace_dir, figname)
+        save_figure(fig, actual_png)
+        compare_snapshots(reference, actual_png)
+
 
 class RasterPlotTitleTests(unittest.TestCase):
     """Snapshot tests for plotting rasters with various titles."""

From 68f4641a51bbd657a9967c8b4c13e8a5e8821840 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Fri, 20 Mar 2026 12:37:40 -0400
Subject: [PATCH 28/37] raster_workspace_summary: recursively handle nested
 file_registry (#2321)

---
 src/natcap/invest/reports/raster_utils.py | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/src/natcap/invest/reports/raster_utils.py b/src/natcap/invest/reports/raster_utils.py
index 3be4bede49..2cde5bceb0 100644
--- a/src/natcap/invest/reports/raster_utils.py
+++ b/src/natcap/invest/reports/raster_utils.py
@@ -651,20 +651,20 @@ def raster_workspace_summary(file_registry):
     """Create a table of stats for all rasters in a file_registry."""
     raster_summary = {}
 
-    def _summarize_output(path):
-        resource = _get_raster_metadata(path)
-        band = resource.get_band_description(1) if resource else None
-        if band:
-            filename = os.path.basename(resource.path)
-            raster_summary[filename] = _build_stats_table_row(
-                resource, band)
-
-    for item in file_registry.values():
+    def _summarize_output(item):
         if isinstance(item, collections.abc.Mapping):
             for path in item.values():
                 _summarize_output(path)
         else:
-            _summarize_output(item)
+            resource = _get_raster_metadata(item)
+            band = resource.get_band_description(1) if resource else None
+            if band:
+                filename = os.path.basename(resource.path)
+                raster_summary[filename] = _build_stats_table_row(
+                    resource, band)
+
+    for item in file_registry.values():
+        _summarize_output(item)
 
     return pandas.DataFrame(raster_summary).T
 

From 3903ad4cbc01ca0e6723ca2749b42b5e89a70be4 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Fri, 20 Mar 2026 12:45:55 -0400
Subject: [PATCH 29/37] provide quickflow, baseflow, and precip CSV values in
 m3/month (#2321)

---
 .../seasonal_water_yield.py                   | 39 ++++++-------------
 tests/test_seasonal_water_yield_regression.py | 30 +++++---------
 2 files changed, 22 insertions(+), 47 deletions(-)

diff --git a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
index 203b1e449a..59ba7cee76 100644
--- a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
+++ b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
@@ -496,29 +496,29 @@
                     id="quickflow",
                     about=gettext(
                         "The average quickflow value for the month in the watershed,"
-                        " expressed in cubic meters per second."
+                        " expressed in cubic meters."
                     ),
-                    units=u.meter ** 3 / u.second
+                    units=u.meter ** 3 / u.month
                 ),
                 spec.NumberOutput(
                     id="baseflow",
                     about=gettext(
                         "The average baseflow value for the month in the watershed,"
-                        " expressed in cubic meters per second. Since baseflow is"
-                        " calculated on an annual scale, the values for each watershed"
-                        " have been distributed evenly across the year (annual average"
-                        " divided by 12)."
+                        " expressed in cubic meters. Since baseflow is calculated on"
+                        " an annual scale, the values for each watershed have been"
+                        " distributed evenly across the year (annual average divided"
+                        " by 12)."
                     ),
-                    units=u.meter ** 3 / u.second
+                    units=u.meter ** 3 / u.month
                 ),
                 spec.NumberOutput(
                     id="precipitation",
                     about=gettext(
                         "The average precipitation value for the month in the watershed,"
-                        " expressed in cubic meters per second. Values are based on"
-                        " the aligned input monthly precipitation rasters."
+                        " expressed in cubic meters. Values are based on the aligned"
+                        " input monthly precipitation rasters."
                     ),
-                    units=u.meter ** 3 / u.second
+                    units=u.meter ** 3 / u.month
                 )
             ],
             index_col="geom_id",
@@ -1582,20 +1582,6 @@ def _generate_monthly_qf_b_p_csv(
         LOGGER.warning(f'{target_csv_path} exists, deleting and writing new output')
         os.remove(target_csv_path)
 
-    seconds_per_month = {
-        1: 2678400,
-        2: 2440152,
-        3: 2678400,
-        4: 2592000,
-        5: 2678400,
-        6: 2592000,
-        7: 2678400,
-        8: 2678400,
-        9: 2592000,
-        10: 2678400,
-        11: 2592000,
-        12: 2678400}
-
     # Use the baseflow raster to get the pixel_size; all rasters should be aligned
     raster_info = pygeoprocessing.get_raster_info(annual_baseflow_path)
     pixel_area_m2 = numpy.prod([abs(x) for x in raster_info['pixel_size']])
@@ -1625,7 +1611,7 @@ def _generate_monthly_qf_b_p_csv(
     b_avg_per_feat_per_month = {k: v['sum'] * 0.001 * pixel_area_m2 / 12
                                 for k, v in stats_by_field['baseflow'][0].items()}
 
-    values_dict = {fid: {month: {'baseflow': b_val / seconds_per_month[month]}
+    values_dict = {fid: {month: {'baseflow': b_val}
                          for month in MONTH_RANGE}
                    for fid, b_val in b_avg_per_feat_per_month.items()}
 
@@ -1635,8 +1621,7 @@ def _generate_monthly_qf_b_p_csv(
                     k: v['sum'] * 0.001 * pixel_area_m2
                     for k, v in stats_by_field[value_name][month].items()}
             for fid, value in avg_per_feat_per_month.items():
-                value_per_second = value / seconds_per_month[month]
-                values_dict[fid][month][value_name] = value_per_second
+                values_dict[fid][month][value_name] = value
 
     with open(target_csv_path, 'w', newline='') as csvfile:
         writer = csv.writer(csvfile, delimiter=',')
diff --git a/tests/test_seasonal_water_yield_regression.py b/tests/test_seasonal_water_yield_regression.py
index 3084f018de..37354dc031 100644
--- a/tests/test_seasonal_water_yield_regression.py
+++ b/tests/test_seasonal_water_yield_regression.py
@@ -617,19 +617,14 @@ def test_base_regression(self):
         # check the values in the avg monthly quickflow baseflow precip csv
         actual_result_df = pandas.read_csv(
             os.path.join(args['workspace_dir'], 'monthly_quickflow_baseflow.csv'))
-        expected_qf = [
-            2.116894e-5, 2.541212e-5, 2.514573e-5, 2.805605e-5, 2.916745e-5, 3.223471e-5,
-            3.323317e-5, 3.528216e-5, 3.858638e-5, 3.941151e-5, 4.287460e-5, 4.358156e-5]
-        expected_b = [
-            2.276286e-5, 2.498535e-5, 2.276286e-5, 2.352162e-5, 2.276286e-5, 2.352162e-5,
-            2.276286e-5, 2.276286e-5, 2.352162e-5, 2.276286e-5, 2.352162e-5, 2.276286e-5]
-        expected_p = [
-            4.106930e-5, 4.917726e-5, 4.853644e-5, 5.401235e-5, 5.600358e-5, 6.172840e-5,
-            6.347073e-5, 6.720430e-5, 7.330247e-5, 7.467145e-5, 8.101852e-5, 8.213859e-5]
+        expected_qf = [56.69889, 62.00944, 67.35032, 72.72129, 78.12209, 83.55236,
+                       89.01173, 94.49973, 100.01591, 105.55979, 111.13096, 116.72885]
+        expected_b = [60.96804 for i in range(12)]
+        expected_p = [110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 220]
         for expected_val, col_name in [(expected_qf, 'quickflow'),
                 (expected_b, 'baseflow'), (expected_p, 'precipitation')]:
             numpy.testing.assert_allclose(expected_val, actual_result_df[col_name],
-                                          rtol=1e-6)
+                                          rtol=1e-5)
 
     def test_base_regression_d8(self):
         """SWY base regression test on sample data in D8 mode.
@@ -686,19 +681,14 @@ def test_base_regression_d8(self):
         # check the values in the avg monthly quickflow baseflow precip csv
         actual_result_df = pandas.read_csv(
             os.path.join(args['workspace_dir'], 'monthly_quickflow_baseflow.csv'))
-        expected_qf = [
-            2.083911e-5, 2.501863e-5, 2.475883e-5, 2.762715e-5, 2.872444e-5, 3.174831e-5,
-            3.273498e-5, 3.475670e-5, 3.801546e-5, 3.883215e-5, 4.224839e-5, 4.294909e-5]
-        expected_b = [
-            2.313658e-5, 2.539555e-5, 2.313658e-5, 2.390779e-5, 2.313658e-5, 2.390779e-5,
-            2.313658e-5, 2.313658e-5, 2.390779e-5, 2.313658e-5, 2.390779e-5, 2.313658e-5]
-        expected_p = [
-            4.106930e-5, 4.917726e-5, 4.853644e-5, 5.401235e-5, 5.600358e-5, 6.172840e-5,
-            6.347073e-5, 6.720430e-5, 7.330247e-5, 7.467145e-5, 8.101852e-5, 8.213859e-5]
+        expected_qf = [55.81547, 61.04926, 66.31405, 71.60957, 76.93555, 82.29161,
+                       87.67738, 93.09235, 98.53606, 104.00803, 109.50784, 115.03485]
+        expected_b = [61.969 for i in range(12)]
+        expected_p = [110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 220]
         for expected_val, col_name in [(expected_qf, 'quickflow'),
                 (expected_b, 'baseflow'), (expected_p, 'precipitation')]:
             numpy.testing.assert_allclose(expected_val, actual_result_df[col_name],
-                                          rtol=1e-6)
+                                          rtol=1e-5)
 
     def test_base_regression_nodata_inf(self):
         """SWY base regression test on sample data with really small nodata.

From 47b69cb993c7507628a2294b7e2adb06d5a38584 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Fri, 20 Mar 2026 13:05:22 -0400
Subject: [PATCH 30/37] update qf+b+precip chart axis labels to reflect
 m3/month (#2321)

---
 src/natcap/invest/seasonal_water_yield/reporter.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/natcap/invest/seasonal_water_yield/reporter.py b/src/natcap/invest/seasonal_water_yield/reporter.py
index 1a201313b8..3bea7ffa25 100644
--- a/src/natcap/invest/seasonal_water_yield/reporter.py
+++ b/src/natcap/invest/seasonal_water_yield/reporter.py
@@ -94,7 +94,7 @@ def _create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path, xy_ratio):
         ['baseflow', 'quickflow']
     ).encode(
         altair.X("month(month):O").title("Month"),
-        altair.Y("sum(value):Q").title("Quickflow + Baseflow (m3/s)"),
+        altair.Y("sum(value):Q").title("Quickflow + Baseflow (m3/month)"),
         altair.Order(field='key', sort='ascending'),
         color=altair.Color('key:N').scale(
             domain=['quickflow', "baseflow", "precipitation"],
@@ -110,7 +110,7 @@ def _create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path, xy_ratio):
         altair.Y(
             "sum(precipitation)",
             axis=altair.Axis(orient="right")
-        ).title("Precipitation (m3/s)"),
+        ).title("Precipitation (m3/month)"),
         color=altair.value('#0500a3')
     )
 

From 6645ee4ce44127329ad468cda29edc133fee00ce Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Fri, 20 Mar 2026 14:51:17 -0400
Subject: [PATCH 31/37] update MONTH_ID_TO_LABEL so that it is indexed in line
 with month numbers (to support consistent month indexing) (#2321)

---
 .../invest/seasonal_water_yield/seasonal_water_yield.py     | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
index 59ba7cee76..f52d2a2af8 100644
--- a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
+++ b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
@@ -23,9 +23,9 @@
 TARGET_NODATA = -1
 N_MONTHS = 12
 MONTH_RANGE = range(1, N_MONTHS+1)
-MONTH_ID_TO_LABEL = [
-    'jan', 'feb', 'mar', 'apr', 'may', 'jun', 'jul', 'aug', 'sep', 'oct',
-    'nov', 'dec']
+MONTH_ID_TO_LABEL = {
+    1: 'jan', 2: 'feb', 3: 'mar', 4: 'apr', 5: 'may', 6: 'jun',
+    7: 'jul', 8: 'aug', 9: 'sep', 10: 'oct',  11: 'nov', 12: 'dec'}
 
 _model_description = gettext(
     """

From 7150fca8a90c2dbd9160240bf47bd1361fd0c7a6 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Fri, 20 Mar 2026 15:16:54 -0400
Subject: [PATCH 32/37] add regression test for user_defined_climate_zones base
 case (#2321)

---
 tests/test_seasonal_water_yield_regression.py | 43 +++++++++++++++++++
 1 file changed, 43 insertions(+)

diff --git a/tests/test_seasonal_water_yield_regression.py b/tests/test_seasonal_water_yield_regression.py
index 37354dc031..bd37b16871 100644
--- a/tests/test_seasonal_water_yield_regression.py
+++ b/tests/test_seasonal_water_yield_regression.py
@@ -972,6 +972,49 @@ def test_user_recharge(self):
             os.path.join(args['workspace_dir'], 'aggregated_results_swy.shp'),
             agg_results_csv_path)
 
+    def test_climate_zones(self):
+        """SWY user recharge regression test on sample data.
+
+        Executes SWY in user defined local recharge mode and checks that the
+        output files are generated and that the aggregate shapefile fields
+        are the same as the regression case.
+        """
+        from natcap.invest.seasonal_water_yield import seasonal_water_yield
+
+        # use predefined directory so test can clean up files during teardown
+        workspace_dir = os.path.join(self.workspace_dir, 'workspace')
+        os.mkdir(workspace_dir)
+        args = SeasonalWaterYieldRegressionTests.generate_base_args(
+            workspace_dir)
+        args['monthly_alpha'] = False
+        args['results_suffix'] = ''
+
+        cz_csv_path = os.path.join(self.workspace_dir, 'cz.csv')
+        make_climate_zone_csv(cz_csv_path)
+        cz_ras_path = os.path.join(args['workspace_dir'], 'dem.tif')
+        make_gradient_raster(cz_ras_path)
+        args['climate_zone_raster_path'] = cz_ras_path
+        args['climate_zone_table_path'] = cz_csv_path
+        args['user_defined_climate_zones'] = True
+
+        execute_kwargs = {
+            'generate_report': bool(seasonal_water_yield.MODEL_SPEC.reporter),
+            'save_file_registry': True,
+            'check_outputs': True
+        }
+        seasonal_water_yield.MODEL_SPEC.execute(args, **execute_kwargs)
+        assert_complete_execute(
+            args, seasonal_water_yield.MODEL_SPEC, **execute_kwargs)
+
+        # generate aggregated results csv table for assertion
+        agg_results_csv_path = os.path.join(args['workspace_dir'],
+                                            'agg_results_cz.csv')
+        make_agg_results_csv(agg_results_csv_path, climate_zones=True)
+
+        SeasonalWaterYieldRegressionTests._assert_regression_results_equal(
+            os.path.join(args['workspace_dir'], 'aggregated_results_swy.shp'),
+            agg_results_csv_path)
+
     @staticmethod
     def _assert_regression_results_equal(
             result_vector_path, agg_results_path):

From 24b9f968abd9c1658106536d2de34c52937313d3 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Fri, 20 Mar 2026 15:22:36 -0400
Subject: [PATCH 33/37] fix test docstring (#2321)

---
 tests/test_seasonal_water_yield_regression.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/tests/test_seasonal_water_yield_regression.py b/tests/test_seasonal_water_yield_regression.py
index bd37b16871..fb54023707 100644
--- a/tests/test_seasonal_water_yield_regression.py
+++ b/tests/test_seasonal_water_yield_regression.py
@@ -972,10 +972,10 @@ def test_user_recharge(self):
             os.path.join(args['workspace_dir'], 'aggregated_results_swy.shp'),
             agg_results_csv_path)
 
-    def test_climate_zones(self):
-        """SWY user recharge regression test on sample data.
+    def test_user_climate_zones(self):
+        """SWY user climate zones test on sample data.
 
-        Executes SWY in user defined local recharge mode and checks that the
+        Executes SWY in user defined climate zones mode and checks that the
         output files are generated and that the aggregate shapefile fields
         are the same as the regression case.
         """

From 779f84f8015fdf92f88b342cd2417e29a9edc5ce Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Fri, 20 Mar 2026 16:47:38 -0400
Subject: [PATCH 34/37] fixes and improvements to utils; different color scheme
 for qb, vri_sum plots; point Makefile to invest-test-data revision (#2321)

---
 Makefile                                          |  2 +-
 src/natcap/invest/reports/raster_utils.py         | 15 ++++++++++-----
 src/natcap/invest/reports/vector_utils.py         |  5 +++--
 .../invest/seasonal_water_yield/reporter.py       |  7 +++++--
 4 files changed, 19 insertions(+), 10 deletions(-)

diff --git a/Makefile b/Makefile
index bc895612fb..8b3e5c108e 100644
--- a/Makefile
+++ b/Makefile
@@ -6,7 +6,7 @@ GIT_SAMPLE_DATA_REPO_REV    := cfd1f07673e66823fd22989a2b87afb017aac447
 
 GIT_TEST_DATA_REPO          := https://bitbucket.org/natcap/invest-test-data.git
 GIT_TEST_DATA_REPO_PATH     := $(DATA_DIR)/invest-test-data
-GIT_TEST_DATA_REPO_REV      := c791f2b50e67680832054536899efacbc72e9e0b
+GIT_TEST_DATA_REPO_REV      := 94c4bc9f0f22082d2251b6b14063eb0ff4094451
 
 GIT_UG_REPO                 := https://github.com/natcap/invest.users-guide
 GIT_UG_REPO_PATH            := doc/users-guide
diff --git a/src/natcap/invest/reports/raster_utils.py b/src/natcap/invest/reports/raster_utils.py
index 2cde5bceb0..3bb8152d92 100644
--- a/src/natcap/invest/reports/raster_utils.py
+++ b/src/natcap/invest/reports/raster_utils.py
@@ -6,7 +6,6 @@
 import os
 from io import BytesIO
 from enum import Enum
-from typing import Union, Optional
 
 import distinctipy
 import geometamaker
@@ -14,12 +13,13 @@
 import pygeoprocessing
 import matplotlib
 import matplotlib.colors
-from matplotlib.colors import ListedColormap
+from matplotlib.colors import Colormap, ListedColormap
 import matplotlib.patches
 import matplotlib.pyplot as plt
 import pandas
 import yaml
 from osgeo import gdal
+from pydantic import ConfigDict
 from pydantic.dataclasses import dataclass
 
 from natcap.invest import gettext
@@ -145,7 +145,7 @@ class RasterTransform(str, Enum):
     log = 'log'
 
 
-@dataclass
+@dataclass(config=ConfigDict(arbitrary_types_allowed=True))
 class RasterPlotConfig:
     """A definition for how to plot a raster."""
 
@@ -159,7 +159,7 @@ class RasterPlotConfig:
     """For highly skewed data, a transformation can help reveal variation."""
     title: str | None = None
     """An optional plot title. If ``None``, the filename is used."""
-    colormap: Optional[Union[str, object]] = None
+    colormap: str | Colormap | None = None
     """The string name of a registered matplotlib colormap or a colormap object."""
 
     def __post_init__(self):
@@ -167,6 +167,9 @@ def __post_init__(self):
             self.title = os.path.basename(self.raster_path)
         self.caption = f'{self.title}:{self.spec.about}'
 
+        self.colormap = plt.get_cmap(self.colormap if self.colormap
+                                     else COLORMAPS[self.datatype])
+
 
 def build_raster_plot_configs(id_lookup_table, raster_plot_tuples):
     """Build RasterPlotConfigs for use in plotting input or output rasters.
@@ -318,6 +321,8 @@ def _get_title_kwargs(title: str, resampled: bool, line_width: int, facets=False
     label = textwrap.fill(label, width=line_width)
     padding = 1.5
     if not facets:
+        # Faceted plots don't need extra padding for title because their units
+        # label appears with the legend instead of under the title
         padding *= SUBTITLE_FONT_SIZE
     return {
         'fontfamily': 'monospace',
@@ -400,7 +405,7 @@ def plot_raster_list(raster_list: list[RasterPlotConfig]):
         colorbar_kwargs = {}
         imshow_kwargs['norm'] = transform
         imshow_kwargs['interpolation'] = 'none'
-        cmap = plt.get_cmap(config.colormap if config.colormap else COLORMAPS[dtype])
+        cmap = config.colormap
         if dtype == 'divergent':
             if transform == 'log':
                 transform = matplotlib.colors.SymLogNorm(linthresh=0.03)
diff --git a/src/natcap/invest/reports/vector_utils.py b/src/natcap/invest/reports/vector_utils.py
index 5dca731962..f53abf0776 100644
--- a/src/natcap/invest/reports/vector_utils.py
+++ b/src/natcap/invest/reports/vector_utils.py
@@ -21,12 +21,13 @@ def get_geojson_bbox(geodataframe):
 
     Args:
         geodataframe (geopandas.GeoDataFrame):
+
     Returns:
         tuple: A 2-tuple containing:
             - extent_feature (dict): A GeoJSON feature representing the bounding
-            box of the input GeoDataFrame.
+              box of the input GeoDataFrame.
             - xy_ratio (float): The aspect ratio of the bounding box
-            (width/height).
+              (width/height).
 
     """
     xmin, ymin, xmax, ymax = geodataframe.total_bounds
diff --git a/src/natcap/invest/seasonal_water_yield/reporter.py b/src/natcap/invest/seasonal_water_yield/reporter.py
index 3bea7ffa25..31c3433f4a 100644
--- a/src/natcap/invest/seasonal_water_yield/reporter.py
+++ b/src/natcap/invest/seasonal_water_yield/reporter.py
@@ -49,7 +49,10 @@ def _create_aggregate_map(geodataframe, xy_ratio, attribute,
         type='identity',
         reflectY=True
     ).encode(
-        color=attribute,
+        color=altair.Color(
+            attribute,
+            scale=altair.Scale(domainMid=0, scheme="brownbluegreen")
+        ),
         tooltip=[altair.Tooltip(attribute, title=attribute)]
     ).properties(
         width=MAP_WIDTH,
@@ -124,7 +127,7 @@ def _create_linked_monthly_plots(aoi_vector_path, aggregate_csv_path, xy_ratio):
         )
     )
 
-    legend_config = vector_utils.LEGEND_CONFIG
+    legend_config = vector_utils.LEGEND_CONFIG.copy()
     legend_config['orient'] = 'right'
 
     chart = altair.hconcat(

From c3a7f952a239ed44a4ea1f58350c0d5388e73082 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Fri, 20 Mar 2026 17:18:14 -0400
Subject: [PATCH 35/37] custom matplotlib colormap for local recharge raster
 that matches the qb + vri_sum plot color scheme (#2321)

---
 .../invest/seasonal_water_yield/reporter.py       | 15 +++++++++++++--
 1 file changed, 13 insertions(+), 2 deletions(-)

diff --git a/src/natcap/invest/seasonal_water_yield/reporter.py b/src/natcap/invest/seasonal_water_yield/reporter.py
index 31c3433f4a..43410e9fb0 100644
--- a/src/natcap/invest/seasonal_water_yield/reporter.py
+++ b/src/natcap/invest/seasonal_water_yield/reporter.py
@@ -9,6 +9,7 @@
 import numpy
 import pandas
 import pygeoprocessing
+from matplotlib.colors import LinearSegmentedColormap
 from osgeo import gdal
 from osgeo import ogr
 
@@ -31,6 +32,14 @@
 
 MAP_WIDTH = 450 # pixels
 
+# Custom matplotlib colormap that matches the Altair color scheme
+# used for the qb and vri_sum maps
+brownbluegreen = LinearSegmentedColormap.from_list(
+    'brownbluegreen', (
+        (0.000, (0.627, 0.396, 0.102)),
+        (0.500, (0.933, 0.945, 0.918)),
+        (1.000, (0.094, 0.443, 0.447))))
+
 qf_label_month_map = {
     f"qf_{month_index}": str(month_index) for month_index in range(1, 13)
 }
@@ -245,7 +254,8 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
             RasterPlotConfig(
                 raster_path=args_dict['l_path'],
                 datatype=RasterDatatype.divergent,
-                spec=model_spec.get_input('l_path')))
+                spec=model_spec.get_input('l_path'),
+                colormap=brownbluegreen))
         qf_rasters = None
         raster_outputs_heading = 'Annual Baseflow'
 
@@ -266,7 +276,8 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
             RasterPlotConfig(
                 raster_path=file_registry['l'],
                 datatype=RasterDatatype.divergent,
-                spec=model_spec.get_output('l'))])
+                spec=model_spec.get_output('l'),
+                colormap=brownbluegreen)])
         raster_outputs_heading = 'Additional Raster Outputs'
 
         input_raster_config_list.append(

From 30c8535b734e19e15d3e9ae006bce5c0d381f873 Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Mon, 23 Mar 2026 09:13:05 -0400
Subject: [PATCH 36/37] use built-in matplotlib BrBG instead of custom version
 (#2321)

---
 src/natcap/invest/seasonal_water_yield/reporter.py | 12 ++----------
 1 file changed, 2 insertions(+), 10 deletions(-)

diff --git a/src/natcap/invest/seasonal_water_yield/reporter.py b/src/natcap/invest/seasonal_water_yield/reporter.py
index 43410e9fb0..44ad7cb9af 100644
--- a/src/natcap/invest/seasonal_water_yield/reporter.py
+++ b/src/natcap/invest/seasonal_water_yield/reporter.py
@@ -32,14 +32,6 @@
 
 MAP_WIDTH = 450 # pixels
 
-# Custom matplotlib colormap that matches the Altair color scheme
-# used for the qb and vri_sum maps
-brownbluegreen = LinearSegmentedColormap.from_list(
-    'brownbluegreen', (
-        (0.000, (0.627, 0.396, 0.102)),
-        (0.500, (0.933, 0.945, 0.918)),
-        (1.000, (0.094, 0.443, 0.447))))
-
 qf_label_month_map = {
     f"qf_{month_index}": str(month_index) for month_index in range(1, 13)
 }
@@ -255,7 +247,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
                 raster_path=args_dict['l_path'],
                 datatype=RasterDatatype.divergent,
                 spec=model_spec.get_input('l_path'),
-                colormap=brownbluegreen))
+                colormap='BrBG'))
         qf_rasters = None
         raster_outputs_heading = 'Annual Baseflow'
 
@@ -277,7 +269,7 @@ def report(file_registry, args_dict, model_spec, target_html_filepath):
                 raster_path=file_registry['l'],
                 datatype=RasterDatatype.divergent,
                 spec=model_spec.get_output('l'),
-                colormap=brownbluegreen)])
+                colormap='BrBG')])
         raster_outputs_heading = 'Additional Raster Outputs'
 
         input_raster_config_list.append(

From b86deabe40ec350bc73d830a89cdebad327c78ca Mon Sep 17 00:00:00 2001
From: Megan Nissel <megnissel@gmail.com>
Date: Mon, 23 Mar 2026 09:49:21 -0400
Subject: [PATCH 37/37] update HISTORY.rst to reflect updates to the SWY model
 and report (#2321)

---
 HISTORY.rst | 24 +++++++++++++++++++++---
 1 file changed, 21 insertions(+), 3 deletions(-)

diff --git a/HISTORY.rst b/HISTORY.rst
index b578021b2d..6f9207b1b7 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -61,10 +61,28 @@
   7. InVEST model Z (model names should be sorted A-Z)
 
 
+Unreleased Changes
+------------------
+
+Seasonal Water Yield
+====================
+* The model now generates a report, a visual summary of results, available in
+  the output workspace and also viewable from the Workbench after the model run
+  completes. (`#2321 <https://github.com/natcap/invest/issues/2321>`_)
+* The model now generates an additional output, a CSV containing average monthly
+  quickflow, baseflow, and precipitation values, in cubic meters per month, for
+  each feature in the AOI. This output is used by the report to generate some
+  plots. Note that this CSV is only created when the model is run without
+  inputting a Local Recharge raster.
+  (`#2321 <https://github.com/natcap/invest/issues/2321>`_)
+* Various updates to model output data metadata, including correcting the
+  units of some outputs.
+  (`#2450 <https://github.com/natcap/invest/issues/2450>`_)
+* Updated the naming convention of several monthly intermediate outputs to be
+  1-indexed rather than 0-indexed. This makes filenames consistent throughout
+  the model, where 1=January and 12=December.
+  (`#2451 <https://github.com/natcap/invest/issues/2451>`_)
 
-..
-  Unreleased Changes
-  ------------------
 
 3.18.0 (2026-02-25)
 -------------------