diff --git a/postprocessing/compare_inputs.py b/postprocessing/compare_inputs.py
new file mode 100644
index 000000000..2e20d9485
--- /dev/null
+++ b/postprocessing/compare_inputs.py
@@ -0,0 +1,721 @@
+#%% Imports
+import argparse
+import os
+import re
+import sys
+import textwrap
+import traceback
+import cmocean
+
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+from matplotlib import patheffects as pe
+import pandas as pd
+
+sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '..')))
+import reeds
+from postprocessing import input_plots
+
+reeds_path = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+reeds.plots.plotparams()
+
+#%% Global Variables
+# Define regions for which hourly demand profiles will be plotted
+hourly_regions = {
+    'IL':['IL_MISO','IL_PJM'],
+    'TX':['TX_ERCOT_E', 'TX_ERCOT_W', 'TX_MISO', 'TX_SPP'],
+}
+cmap = cmocean.cm.tempo
+cmap_diff = plt.cm.RdBu_r
+
+#%% Plotting
+
+def _parse_weatheryears(weatheryear):
+    """
+    Parse a weather-year input into a sorted list of unique integer years.
+    
+    Parameters
+    ----------
+    weatheryear : int, str, or iterable of int
+        The weather year(s) to parse. Can be a single integer, a string with
+        years separated by underscores or spaces, or an iterable of integers.
+
+    Returns
+    -------
+    list of int
+        Sorted list of unique weather years.
+    """
+    if isinstance(weatheryear, int):
+        return [weatheryear]
+
+    if isinstance(weatheryear, str):
+        tokens = [tok.strip() for tok in re.split(r'[_\s]+', weatheryear) if tok.strip()]
+    else:
+        try:
+            tokens = list(weatheryear)
+        except TypeError as exc:
+            raise ValueError(
+                f'Invalid weather year input: {weatheryear}. Provide an integer year '
+                'or comma/space-delimited list of years.'
+            ) from exc
+
+    if not tokens:
+        raise ValueError('No weather years provided.')
+
+    try:
+        years = sorted(set(int(token) for token in tokens))
+    except (TypeError, ValueError) as exc:
+        raise ValueError(
+            f'Invalid weather year input: {weatheryear}. All weather years must be integers.'
+        ) from exc
+
+    return years
+
+def plot_hourly_demand_profiles(cases, colors, year='last', weatheryear=2012, region=None,):
+    '''
+    Plot hourly demand profiles for the specified cases and region
+    using reeds.input_plots.plot_profile().
+
+    Parameters
+    ----------
+    cases : dict
+        Dictionary of case names and their corresponding paths.
+    colors : dict
+        Dictionary of case names and their corresponding colors.
+    year : int or str, optional
+        Model year to plot. Default is 'last'.
+    weatheryear : int, str, or iterable of int, optional
+        Weather year(s) to plot. Default is 2012.
+    region : str, optional
+        Region for which to plot the demand profiles. Default is None.
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+        The figure object containing the plot.
+    matplotlib.axes.Axes
+        The axes object containing the plot.
+    '''
+    selected_weatheryears = _parse_weatheryears(weatheryear)
+
+    plt.close()
+    f, ax = plt.subplots(figsize=(13.33, 4))
+    
+    for idx, (casename, casepath) in enumerate(cases.items()):
+        color = colors.get(casename, f'C{idx}')
+
+        input_plots.plot_profile(
+            casepath, datum='demand', year=year, weatheryears=selected_weatheryears,
+            reg_sub=region, color=color, label=casename, f=f, ax=ax, hourly=True,
+        )
+    ax.margins(x=0)
+    ax.grid(True, which='major', axis='y', alpha=0.3)
+    
+    ax.legend(frameon=False, loc='lower center', ncol=len(cases))
+    plt.autoscale(enable=True, axis='y')
+    ax.set_ylim(0)
+
+    return f, ax
+
+def plot_demand_yearbymonth(cases, colors, year='last', weatheryear=2012):
+    '''
+    Plot demand profiles in rows by month for the specified cases.
+
+    Parameters
+    ----------
+    cases : dict
+        Dictionary of case names and their corresponding paths.
+    colors : dict
+        Dictionary of case names and their corresponding colors.
+    year : int or str, optional
+        Model year to plot. Default is 'last'.
+    weatheryear : int, str, or iterable of int, optional
+        Weather year(s) to plot. Default is 2012.
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+        The figure object containing the plot.
+    matplotlib.axes.Axes
+        The axes object containing the plot.
+    '''
+    selected_weatheryears = _parse_weatheryears(weatheryear)
+
+    plt.close()
+    for idx, (casename, casepath) in enumerate(cases.items()):
+        color = colors.get(casename, f'C{idx}')
+        
+        dfprofile = (
+            reeds.results.summarize_load_data(
+                casepath, 
+                use_run=True,
+                agg_reg_lvl='all',
+                model_year_sub=year if year in [0, 'last'] else [year], 
+                weather_year_sub=selected_weatheryears,
+            )
+            / 1e3
+        )
+        dfprofile.index = pd.to_datetime(dfprofile.index.get_level_values('datetime'))
+        dfprofile.rename(columns={'load_MWh': casename}, inplace=True)
+        if idx == 0:
+            f,ax = reeds.plots.plotyearbymonth(
+                dfprofile, style='line', colors=[color],
+            )
+        else:
+            reeds.plots.plotyearbymonth(
+                dfprofile, style='line', colors=[color], f=f, ax=ax,
+            )
+                   
+    ax[0].legend(loc='upper left', bbox_to_anchor=(1,1), frameon=False, fontsize='x-large')
+    
+    return f, ax
+
+def plot_peak_and_total_load(cases, colors, weatheryear=None):
+    """
+    Plots annual total load and annual peak load by model year
+    for the specified cases, including range across all weather years.
+    
+    Parameters
+    ----------
+    cases : dict
+        Dictionary of case names and their corresponding paths.
+    colors : dict
+        Dictionary of case names and their corresponding colors.
+    weatheryear : int, str, or iterable of int, optional
+        Weather year(s) to plot with a line. Default is None.
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+        The figure object containing the plot.
+    matplotlib.axes.Axes
+        The axes object containing the plot.
+    pandas.DataFrame
+        The statistics by case, indexed by model year.
+    """
+    selected_weatheryears = _parse_weatheryears(weatheryear)
+
+    f, (ax_total, ax_peak) = plt.subplots(
+        1, 2, figsize=(12, 4.5), constrained_layout=True
+    )
+    stats_by_case = {}
+    for idx, (casename, casepath) in enumerate(cases.items()):
+        print(f'  {casename}: loading annual and peak stats by model year...')
+        color = colors.get(casename, f'C{idx}')
+
+        sw = reeds.io.get_switches(casepath)
+        modelyears = reeds.io.get_years(casepath)
+        case_weatheryears = sw.resource_adequacy_years_list
+
+        dfprofile = (
+            reeds.results.summarize_load_data(
+                casepath, 
+                use_run=True,
+            )
+            / 1e3
+        )
+
+        rows = []
+        for t in modelyears:
+            profile = dfprofile.loc[t].sum(axis=1)
+            profile.index = pd.to_datetime(profile.index.get_level_values('datetime'))
+
+            # profile is in GW and timestamps are hourly, so annual sum is GWh.
+            # Convert to TWh for plotting.
+            annual_total_by_wy = profile.groupby(profile.index.year).sum() / 1e3
+            annual_peak_by_wy = profile.groupby(profile.index.year).max()
+
+            rows.append(
+                {
+                    't': t,
+                    'total_min': annual_total_by_wy.min(),
+                    'total_max': annual_total_by_wy.max(),
+                    'peak_min': annual_peak_by_wy.min(),
+                    'peak_max': annual_peak_by_wy.max(),
+                    **{
+                        f'total_selected_{wy}': annual_total_by_wy.loc[wy]
+                        for wy in selected_weatheryears
+                    },
+                    **{
+                        f'peak_selected_{wy}': annual_peak_by_wy.loc[wy]
+                        for wy in selected_weatheryears
+                    },
+                }
+            )
+
+        dfstats = pd.DataFrame(rows).set_index('t').sort_index()
+        stats_by_case[casename] = dfstats
+
+        x = dfstats.index.values
+
+        avg_total = dfstats[[f'total_selected_{wy}' for wy in selected_weatheryears]].mean(axis=1)
+        ax_total.plot(
+            x, avg_total.values, lw=1.2, color=color, label=casename,
+        )
+        ax_total.fill_between(x, dfstats['total_max'].values, dfstats['total_min'].values, lw=0, alpha=0.15, color=color)
+        ax_total.plot(x, dfstats['total_max'].values, lw=0.8, linestyle=':', color=color, alpha=0.7)
+        ax_total.plot(x, dfstats['total_min'].values, lw=0.8, linestyle=':', color=color, alpha=0.7)
+
+        avg_peak = dfstats[[f'peak_selected_{wy}' for wy in selected_weatheryears]].mean(axis=1)
+        ax_peak.plot(
+            x, avg_peak.values, lw=1.2, color=color, label=casename,
+        )
+        ax_peak.fill_between(x, dfstats['peak_max'].values, dfstats['peak_min'].values, lw=0, alpha=0.15, color=color)
+        ax_peak.plot(x, dfstats['peak_max'].values, lw=0.8, linestyle=':', color=color, alpha=0.7)
+        ax_peak.plot(x, dfstats['peak_min'].values, lw=0.8, linestyle=':', color=color, alpha=0.7)
+
+    ax_total.yaxis.set_minor_locator(mpl.ticker.AutoMinorLocator(2))
+    ax_total.set_ylabel('Annual total load [TWh]')
+    ax_total.legend(frameon=False, loc='best')
+    ax_total.margins(x=0)
+    ax_total.set_xlim(left=2020)
+    ax_total.grid(True, which='major', axis='y', alpha=0.3)
+    reeds.plots.despine(ax_total)
+
+    ax_peak.yaxis.set_minor_locator(mpl.ticker.AutoMinorLocator(2))
+    ax_peak.set_xlabel('Model year')
+    ax_peak.set_ylabel('Annual peak load [GW]')
+    ax_peak.legend(frameon=False, loc='best')
+    ax_peak.margins(x=0)
+    ax_peak.set_xlim(left=2020)
+    ax_peak.grid(True, which='major', axis='y', alpha=0.3)
+    reeds.plots.despine(ax_peak)
+
+    plt.draw()
+
+
+    return f, (ax_total, ax_peak), stats_by_case
+
+def plot_regional_peak_demand_maps(cases, year='last', weatheryear=None):
+    """
+    Maps peak demand by state across specified weather years (default all available),
+    comparing across cases.
+
+    Layout: 1 row x N-cases columns.
+      Base case (first case): absolute peak demand [GW] with a sequential colormap.
+      All other cases: percent difference vs. the base case with a diverging colormap.
+
+    Parameters
+    ----------
+    cases : dict  {casename: casepath}
+    year : int or 'last'
+        Model year to plot. Defaults to the last solved year of the base case.
+    weatheryear : int, list, or None
+        Specific weather year(s) to use. Defaults to None, which uses all available weather years.
+    Returns
+    -------
+    matplotlib.figure.Figure
+        The figure object containing the plot.
+    matplotlib.axes.Axes
+        The axes object containing the plot.
+    pandas.DataFrame
+        The peak demand by state for each case, indexed by state, values in GW.
+    """
+    basecasename = list(cases.keys())[0]
+    basecasepath = list(cases.values())[0]
+
+    selected_weatheryears = _parse_weatheryears(weatheryear) if weatheryear is not None else None
+
+    dfmap = reeds.io.get_dfmap(basecasepath)
+    dfstates = dfmap['st']
+
+    peak_st = {}
+    for casename, casepath in cases.items():
+        print(f'  {casename}: loading state-level peak demand...')
+        sw = reeds.io.get_switches(casepath)
+        case_weatheryears = sw.resource_adequacy_years_list
+        region_to_state = reeds.io.get_hierarchy(casepath)['st']
+
+        dfprofile = (
+            reeds.results.summarize_load_data(
+                casepath, 
+                use_run=True, 
+                model_year_sub=year if year in [0, None, 'last'] else [year], 
+                weather_year_sub=selected_weatheryears,
+            )
+            / 1e3
+        )
+
+        dfprofile.index = pd.to_datetime(dfprofile.index.get_level_values('datetime'))
+        # Aggregate BAs to states, then take the max simultaneous hour
+        df_t_st = dfprofile.T.groupby(dfprofile.columns.map(region_to_state)).sum().T
+        peak_st[casename] = df_t_st.max()
+
+    ncols = len(cases)
+    cmap_abs = cmap
+
+    plt.close()
+    f, ax = plt.subplots(
+        1, ncols,
+        figsize=(max(5 * ncols, 10), 5),
+        gridspec_kw={'wspace': 0.0},
+    )
+    if ncols == 1:
+        ax = [ax]
+
+    base_vals = peak_st[basecasename]
+    abs_vmax = max(
+        (peak_st[c].reindex(dfstates.index).fillna(0).max() for c in cases),
+        default=0.0,
+    )
+    pct_diffs = {
+        c: ((peak_st[c] - base_vals) / base_vals * 100)
+           .replace([float('inf'), float('-inf')], float('nan'))
+        for c in list(cases.keys())[1:]
+    }
+    diff_absmax = max(
+        (pct_diffs[c].abs().max()
+         for c in pct_diffs if not pct_diffs[c].dropna().empty),
+        default=1.0,
+    )
+
+    for col, casename in enumerate(cases):
+        _ax = ax[col]
+        dfplot = dfstates.copy()
+
+        if casename == basecasename:
+            dfplot['val'] = base_vals.reindex(dfstates.index).fillna(0)
+            _cmap, vmin, vmax = cmap_abs, 0, abs_vmax
+            cbar_label = 'Peak demand [GW]'
+        else:
+            dfplot['val'] = pct_diffs[casename].reindex(dfstates.index)
+            _cmap, vmin, vmax = cmap_diff, -diff_absmax, diff_absmax
+            cbar_label = f'Change in peak demand [%]'
+
+        dfstates.plot(ax=_ax, facecolor='none', edgecolor='k', lw=0.2, zorder=10000)
+        dfplot.plot(ax=_ax, column='val', cmap=_cmap, vmin=vmin, vmax=vmax, legend=False)
+
+        # Value labels for states
+        reeds.reedsplots.label_region_value(
+            dfplot,
+            ax=_ax,
+            column='val',
+            fmt='{:.0f}' if casename == basecasename else '{:+.0f}',
+            text_kwargs={'fontsize': 7, 'ha': 'center', 'va': 'center', 'color': 'k'},
+            pe_kwargs={'linewidth': 1.5, 'foreground': 'white'}
+        )
+
+        if col == 0 or col == 1:
+            reeds.plots.addcolorbarhist(
+                f=f, ax0=_ax, data=dfplot['val'].dropna().values,
+                title=cbar_label, cmap=_cmap, vmin=vmin, vmax=vmax,
+                orientation='horizontal', labelpad=2.25, histratio=0.,
+                cbarwidth=0.05, cbarheight=0.85,
+                cbarbottom=-0.05, cbarhoffset=0.,
+                title_fontsize=18, ticklabel_fontsize=16
+            )
+
+        _ax.annotate(casename, (0.1, 1), xycoords='axes fraction', fontsize=18)
+        _ax.axis('off')
+
+    df_peak = pd.DataFrame(peak_st).reindex(dfstates.index)
+    return f, ax, df_peak
+
+
+def plot_regional_total_demand_maps(cases, year='last', weatheryear=2012):
+    """
+    Maps total annual demand by state, comparing across cases.
+
+    For each state, total demand is summed across all BAs in that state over
+    the selected weather year(s) and averaged if more than one year is given.
+
+    Layout: 1 row x N-cases columns.
+      Base case (first case): absolute total demand [TWh] with a sequential colormap.
+      All other cases: percent difference vs. the base case with a diverging colormap.
+
+    Parameters
+    ----------
+    cases : dict  {casename: casepath}
+    year : int or 'last'
+        Model year to plot. Defaults to the last solved year of the base case.
+    weatheryear : int, str, or iterable
+        Weather year(s) to use. Values are averaged when multiple are given.
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+        The figure object containing the plot.
+    matplotlib.axes.Axes
+        The axes object containing the plot.
+    pandas.DataFrame
+        The total demand by state for each case, indexed by state, values in TWh.
+    """
+    selected_weatheryears = _parse_weatheryears(weatheryear)
+    basecasename = list(cases.keys())[0]
+    basecasepath = list(cases.values())[0]
+
+    if year in [0, None, 'last']:
+        t = reeds.io.get_years(basecasepath)[-1]
+    else:
+        t = int(year)
+
+    dfmap = reeds.io.get_dfmap(basecasepath)
+    dfstates = dfmap['st']
+
+    total_st = {}
+    for casename, casepath in cases.items():
+        print(f'  {casename}: loading state-level total demand...')
+        sw = reeds.io.get_switches(casepath)
+        case_weatheryears = sw.resource_adequacy_years_list
+        region_to_state = reeds.io.get_hierarchy(casepath)['st']
+
+        dfprofile = (
+            reeds.results.summarize_load_data(
+                casepath, 
+                use_run=True, 
+                model_year_sub=[t], 
+                weather_year_sub=selected_weatheryears,
+            )
+            / 1e3
+        )
+        dfprofile.index = pd.to_datetime(dfprofile.index.get_level_values('datetime'))
+        dfprofile = dfprofile.loc[
+            str(min(case_weatheryears)) : str(max(case_weatheryears))
+        ]
+
+        # Aggregate BAs to states, then average annual total across ALL available weather years
+        dfprofile_st = dfprofile.T.groupby(dfprofile.columns.map(region_to_state)).sum().T
+        wy_totals = pd.concat(
+            [dfprofile_st.loc[dfprofile_st.index.year == wy].sum() / 1e3  # GWh -> TWh
+             for wy in selected_weatheryears],
+            axis=1,
+        ).mean(axis=1)
+        total_st[casename] = wy_totals
+
+    ncols = len(cases)
+    cmap_abs = cmap
+
+    plt.close()
+    f, ax = plt.subplots(
+        1, ncols,
+        figsize=(max(5 * ncols, 10), 5),
+        gridspec_kw={'wspace': 0.0},
+    )
+    if ncols == 1:
+        ax = [ax]
+
+    base_vals = total_st[basecasename]
+    abs_vmax = max(
+        (total_st[c].reindex(dfstates.index).fillna(0).max() for c in cases),
+        default=0.0,
+    )
+    pct_diffs = {
+        c: ((total_st[c] - base_vals) / base_vals * 100)
+           .replace([float('inf'), float('-inf')], float('nan'))
+        for c in list(cases.keys())[1:]
+    }
+    diff_absmax = max(
+        (pct_diffs[c].abs().max()
+         for c in pct_diffs if not pct_diffs[c].dropna().empty),
+        default=1.0,
+    )
+
+    for col, casename in enumerate(cases):
+        _ax = ax[col]
+        dfplot = dfstates.copy()
+
+        if casename == basecasename:
+            dfplot['val'] = base_vals.reindex(dfstates.index).fillna(0)
+            _cmap, vmin, vmax = cmap_abs, 0, abs_vmax
+            cbar_label = 'Annual demand [TWh]'
+        else:
+            dfplot['val'] = pct_diffs[casename].reindex(dfstates.index)
+            _cmap, vmin, vmax = cmap_diff, -diff_absmax, diff_absmax
+            cbar_label = f'Change in annual demand [%]'
+
+        dfstates.plot(ax=_ax, facecolor='none', edgecolor='k', lw=0.2, zorder=10000)
+        dfplot.plot(ax=_ax, column='val', cmap=_cmap, vmin=vmin, vmax=vmax, legend=False)
+
+        # Value labels for states
+        reeds.reedsplots.label_region_value(
+            dfplot,
+            ax=_ax,
+            column='val',
+            fmt='{:.0f}' if casename == basecasename else '{:+.0f}',
+            text_kwargs={'fontsize': 7, 'ha': 'center', 'va': 'center', 'color': 'k'},
+            pe_kwargs={'linewidth': 1.5, 'foreground': 'white'}
+        )
+        
+        if col == 0 or col == 1:
+            reeds.plots.addcolorbarhist(
+                f=f, ax0=_ax, data=dfplot['val'].dropna().values,
+                title=cbar_label, cmap=_cmap, vmin=vmin, vmax=vmax,
+                orientation='horizontal', labelpad=2.25, histratio=0.,
+                cbarwidth=0.05, cbarheight=0.85,
+                cbarbottom=-0.05, cbarhoffset=0.,
+                title_fontsize=18, ticklabel_fontsize=16
+            )
+
+        _ax.annotate(casename, (0.1, 1), xycoords='axes fraction', fontsize=18)
+        _ax.axis('off')
+
+
+    df_total = pd.DataFrame(total_st).reindex(dfstates.index)
+    return f, ax, df_total
+
+#%% Main
+
+
+
+if __name__ == '__main__':
+    #%% Argument inputs
+    parser = argparse.ArgumentParser(
+        description='Compare ReEDS input data for multiple cases',
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+    )
+    parser.add_argument(
+        'caselist',
+        type=str,
+        nargs='+',
+        help=(
+            'Space-delimited list of case paths, OR shared case prefix, OR CSV '
+            'with casepath/casename/color columns (e.g., postprocessing/example.csv).'
+        ),
+    )
+    parser.add_argument(
+        '--casenames', '-n', type=str, default='',
+        help='Comma-delimited case names when not provided in CSV.',
+    )
+    parser.add_argument(
+        '--titleshorten', '-s', type=str, default='',
+        help='Characters to cut from start of case name (if not using CSV casenames).',
+    )
+    parser.add_argument(
+        '--basecase', '-b', type=str, default='',
+        help='Case-path suffix to force as base case.',
+    )
+    parser.add_argument(
+        '--output', '-o', choices=['png', 'ppt', 'pptx'], default='png',
+        help='Output format (png or pptx)')
+    parser.add_argument(
+        '--year', '-y', type=str, default='2050',
+        help='Model year for hourly demand plot (integer year or "last").',
+    )
+    parser.add_argument(
+        '--weatheryear', '-w', type=str, default=None,
+        help=(
+            'Weather year(s) to plot (single year or comma/space-delimited list, e.g., '
+            '"2012" or "2007,2012"). Defaults to GSw_HourlyWeatherYears from the first case.'
+        ),
+    )
+    args = parser.parse_args()
+
+    try:
+        titleshorten = int(args.titleshorten)
+    except ValueError:
+        titleshorten = len(args.titleshorten)
+
+    cases, colors, basecase, basemap = reeds.report_utils.parse_caselist(
+        args.caselist,
+        args.casenames,
+        args.basecase,
+        titleshorten,
+    )
+    if len(cases) < 2:
+        raise ValueError('Need at least 2 cases to compare inputs.')
+
+    firstcasepath = list(cases.values())[0]
+    outpath = os.path.join(firstcasepath, 'outputs', 'comparisons','inputs')
+    os.makedirs(outpath, exist_ok=True)
+
+    #%% Create output container
+    suffix = args.output.strip('.')
+    if suffix in ['pdf', 'png']:
+        savepath = os.path.join(firstcasepath, 'outputs', 'comparisons','inputs')
+        os.makedirs(savepath, exist_ok=True)
+
+        def saveit(savename):
+            outpath = os.path.join(savepath, savename.lower().replace(' ', '-') + f'.{suffix}')
+            plt.savefig(outpath)
+            print(os.path.basename(outpath))
+            if mpl.interactive:
+                plt.show()
+
+    elif suffix in ['ppt', 'pptx']:
+        savepath = os.path.join(firstcasepath, 'outputs', 'comparisons', 'inputs', 'inputs.pptx')
+        prs = reeds.report_utils.init_pptx()
+        def saveit(savename, **kwargs):
+            reeds.report_utils.add_to_pptx(savename, prs=prs, **kwargs)
+            if mpl.interactive:
+                plt.show()
+
+    try:
+        year = int(args.year)
+    except ValueError:
+        year = args.year
+
+    if args.weatheryear is None:
+        sw = reeds.io.get_switches(firstcasepath)
+        selected_weatheryears = sw.GSw_HourlyWeatherYears
+    else:
+        selected_weatheryears = _parse_weatheryears(args.weatheryear)
+
+    try:
+        f, ax = plot_hourly_demand_profiles(
+            cases,
+            colors,
+            year=year,
+            weatheryear=selected_weatheryears,
+        )
+        saveit(f'Demand hourly profile {selected_weatheryears}')
+    except Exception as e:
+        print(traceback.format_exc())
+    
+    for region_label, region_list in hourly_regions.items():
+        try:
+            f, ax = plot_hourly_demand_profiles(
+                cases,
+                colors,
+                year=year,
+                weatheryear=selected_weatheryears,
+                region=region_list,
+            )
+            saveit(f'Demand hourly profile {selected_weatheryears} {region_label}')
+        except Exception as e:
+            print(traceback.format_exc())
+    
+    try:
+        f, ax = plot_demand_yearbymonth(
+            cases,
+            colors,
+            year=year,
+            weatheryear=selected_weatheryears,
+        )
+        saveit(f'Demand by year and month {year}')
+    except Exception as e:
+        print(traceback.format_exc())
+    
+    try:
+        f_modelyear, _, _ = plot_peak_and_total_load(
+            cases,
+            colors,
+            weatheryear=selected_weatheryears,
+        )
+        saveit(f'Demand by year')
+    except Exception as e:
+        print(traceback.format_exc())
+    
+    try:
+        f_peak, _, _ = plot_regional_peak_demand_maps(
+            cases,
+            year=year,
+            weatheryear=selected_weatheryears,
+        )
+        saveit(f'Demand peak by region {year}')
+    except Exception as e:
+        print(traceback.format_exc())
+
+    try:
+        f_total, _, _ = plot_regional_total_demand_maps(
+            cases,
+            year=year,
+            weatheryear=selected_weatheryears,
+        )
+        saveit(f'Demand total by region {year}')
+    except Exception as e:
+        print(traceback.format_exc())
+
+    #%% Save the powerpoint file if necessary
+    if suffix in ['ppt', 'pptx']:
+        print(f'\ninput_plots.py results saved to:\n{savepath}')
+        prs.save(savepath)
\ No newline at end of file
diff --git a/postprocessing/input_plots.py b/postprocessing/input_plots.py
index 4f16b249d..804b14b51 100644
--- a/postprocessing/input_plots.py
+++ b/postprocessing/input_plots.py
@@ -100,6 +100,7 @@ def plot_profile(
     datum='demand',
     year=0,
     weatheryears=None,
+    reg_sub=None,
     color='k',
     hourly=False,
     f=None,
@@ -118,6 +119,8 @@ def plot_profile(
     ## Parse inputs
     sw = reeds.io.get_switches(case)
     t = reeds.io.get_years(case)[-1] if year in [0, None, 'last'] else year
+    if t not in reeds.io.get_years(case):
+        raise KeyError("Year {} not in the data.".format(t))
     rs = reeds.inputs.parse_regions(case)
     if weatheryears is None:
         weatheryears = sw.resource_adequacy_years_list
@@ -131,11 +134,21 @@ def plot_profile(
             os.path.join(case, 'inputs_case', 'load.h5'),
         ## Convert to GW
         ) / 1e3
-        dfprofile = (
-            dfprofile
-            .loc[t, [r for r in dfprofile if r in rs]]
-            .sum(axis=1)
-        )
+        if reg_sub is not None:
+            if not all(reg in dfprofile.columns for reg in reg_sub):
+                raise KeyError("Region(s) {} not in the data.".format(
+                    " ".join([reg for reg in reg_sub if reg not in dfprofile.columns])))
+            dfprofile = (
+                dfprofile
+                .loc[t, reg_sub]
+                .sum(axis=1)
+            )
+        else:
+            dfprofile = (
+                dfprofile
+                .loc[t, [r for r in dfprofile if r in rs]]
+                .sum(axis=1)
+            )
     elif datum in ['temperature']:
         ylabel = 'Temperature [°C]'
         ...
@@ -148,6 +161,9 @@ def plot_profile(
             ylabel = 'PV CF [%]'
             dfprofile = recf['upv'].squeeze(1)
 
+    if not all(wy in dfprofile.index.year.unique() for wy in weatheryears):
+                    raise KeyError("Weather year(s) {} not in the data.".format(
+                        " ".join([str(wy) for wy in weatheryears if wy not in dfprofile.index.year.unique()])))
     dfprofile = dfprofile.loc[str(min(weatheryears)):str(max(weatheryears))].copy()
     ## Use a continuous set of datetimes to avoid interpolating over missing years
     full_timeseries = pd.date_range(dfprofile.index[0], dfprofile.index[-1], freq='h')
diff --git a/reeds/results.py b/reeds/results.py
index 22743d651..34a6a9836 100644
--- a/reeds/results.py
+++ b/reeds/results.py
@@ -12,6 +12,7 @@
 import h5py
 import numpy as np
 import pandas as pd
+import time
 from pathlib import Path
 from itertools import product
 
@@ -808,3 +809,156 @@ def diff_outputs(
             reeds.io.write_to_h5(df.reset_index(), key, outpath)
         print(f'Difference written to {outpath}')
     return dictout
+
+def summarize_load_data(
+    scenario, 
+    use_run=False, 
+    run_path=None,
+    reg_sub=None,
+    agg_reg_lvl=None, 
+    model_year_sub=None,
+    weather_year_sub=None,
+    agg_by_year=False
+):
+    """helper function to read and summarize load data
+    Parameters
+    ----------
+    scenario
+        Either the name of a ReEDS run (if use_run=True) or the name of a load file in the repo
+    use_run, optional
+        If True looks for load data from a ReEDS run, by default False
+    run_path, optional
+        Folder to look for the runs (if none then uses the run folder in the ReEDS repo)
+        If use use_run=False and run_path is specified, use this as the ReEDS repo.
+    reg_sub, optional
+        List of 'r' regions to subset on
+    agg_reg_lvl, optional
+        Regionality to aggregate to using a column from hierarchy.csv; if 'all' then aggregates all
+        regions, by default None
+    model_year_sub, optional
+        Model years to subset on, specified either as a list of years or 'last' for the final year
+    weather_year_sub
+        Weather years to subset on, specified as a list of years, by default None 
+    Returns
+    -------
+        pd.DataFrame: load data summarized as specified
+    """
+
+    print(f"\nprocessing {scenario}")
+    start_time = time.perf_counter()
+
+    ## load data
+    if use_run:
+        # set path to load file in inputs case folder
+        if run_path is None:
+            run_path = os.path.join(reeds_path, 'runs', scenario)
+        else:
+            run_path = os.path.join(run_path, scenario)
+        load_file_path = os.path.join(run_path, 'inputs_case', 'load.h5')
+    else:
+        # use input load in repo
+        if run_path is None:
+            load_file_path = os.path.join(reeds_path,'inputs','load', scenario + "_load_hourly.h5")
+        else:
+            load_file_path = os.path.join(run_path,'inputs','load', scenario + "_load_hourly.h5")
+
+    # check that file exists; if not, return empty data frame
+    if not os.path.exists(load_file_path):
+        print(f"'{load_file_path}' does not exist.")
+        return pd.DataFrame()
+    print(f"...loading {load_file_path}")
+    dfin = reeds.io.read_file(os.path.join(load_file_path), parse_timestamps=True)
+
+    ## subset to desired regions (defined in terms of r regions)
+    if reg_sub is None:
+        dfout = dfin.copy()
+    else:
+        print(f"...subsetting to {' '.join([reg for reg in reg_sub])}")
+        if not all(reg in dfin.columns for reg in reg_sub):
+            raise KeyError("Region(s) {} not in the data.".format(
+                " ".join([reg for reg in reg_sub if reg not in dfin.columns])))
+        dfout = dfin.loc[:, [reg for reg in reg_sub if reg in dfin.columns]].copy()
+
+    ## aggregate to desired regionality
+    if agg_reg_lvl is None:
+    # no aggregation 
+        pass
+    elif agg_reg_lvl == 'all':
+    # sum all regions
+        print(f"...aggregating all regions")
+        dfout = dfout.sum(axis=1)
+        dfout.name = 'load_MWh'
+    else:
+    # sum to new mapping
+        print(f"...aggregating regionality to {agg_reg_lvl}")
+        # load hierarchy
+        if use_run:
+            hierarchy_file_path = os.path.join(run_path, 'inputs_case', 'hierarchy.csv')
+            reg_col = "*r"
+        else:
+            hierarchy_file_path = os.path.join(reeds_path, 'inputs', 'hierarchy.csv')
+            reg_col = "ba"
+        hierarchy = pd.read_csv(hierarchy_file_path)
+        # set up region mapping
+        if agg_reg_lvl not in hierarchy.columns:
+            raise UserWarning(f"{agg_reg_lvl} is not a column in hierarchy file. Skipping aggregation.")
+        else:
+            region_map = hierarchy[[reg_col, agg_reg_lvl]].set_index(reg_col).squeeze()
+
+        #TODO: add option to aggregate without region mapping
+        # map load to regionality in relevant hierarchy column
+        dfout.columns = dfout.columns.map(region_map)
+        # aggregate
+        dfout = dfout.groupby(level=0, axis=1).sum()
+
+    # subset to desired model years
+    if model_year_sub is not None:
+        print(f"...subsetting to {' '.join([str(my) for my in model_year_sub])} model years")
+        model_year_sub = reeds.io.get_years(casepath)[-1] if model_year_sub in [0,'last'] else model_year_sub
+        if isinstance(model_year_sub, int):
+            model_year_sub = [model_year_sub]
+        my_vals = dfout.index.get_level_values('year')
+        if not all(my in my_vals for my in model_year_sub):
+            raise KeyError("Model year(s) {} not in the data.".format(
+                " ".join([str(my) for my in model_year_sub if my not in my_vals])))
+        dfout = dfout[my_vals.isin(model_year_sub)]
+
+    ## subset to desired weather years
+    if weather_year_sub is not None:
+        print(f"...subsetting to {' '.join([str(wy) for wy in weather_year_sub])} weather years")
+        # temporarily reset index to slice
+        idx_names = dfout.index.names
+        dfout = dfout.reset_index()
+        # isolate weather year from datetime
+        dfout['weather_year'] = dfout['datetime'].dt.year
+        if isinstance(weather_year_sub, int):
+            weather_year_sub = [weather_year_sub]
+        elif not isinstance(weather_year_sub, list):
+            print("'weather_year_sub' only accepts int or list.")
+            raise Exception
+        if not all(wy in dfout['weather_year'].unique() for wy in weather_year_sub):
+            raise KeyError("Weather year(s) {} not in the data.".format(
+                " ".join([str(wy) for wy in weather_year_sub if wy not in dfout['weather_year'].unique()])))
+        dfout = dfout.loc[dfout.weather_year.isin(weather_year_sub)]
+        # reset index
+        dfout = dfout.drop('weather_year', axis=1).set_index(idx_names)
+
+    ## aggregate data
+    if agg_by_year:
+        print(f"...aggregating hourly data to annual")
+        # temporarily reset index to slice
+        idx_names = dfout.index.names
+        dfout = dfout.reset_index()
+        # isolate weather year from datetime
+        dfout['weather_year'] = dfout['datetime'].dt.year
+        dfout = dfout.drop('datetime', axis=1).groupby(['year','weather_year']).sum()
+
+    # ensure output is pd.DataFrame and not a Series
+    if isinstance(dfout, pd.Series):
+        dfout = dfout.to_frame()
+
+    end_time = time.perf_counter()
+    elapsed_time = end_time - start_time
+    print(f"\ndone processing {scenario} ({elapsed_time:.0f} seconds)")
+
+    return dfout
\ No newline at end of file