Commit last-minute

hydroshoot/example/figs/paper_figs.py

@@ -0,0 +1,995 @@
+#!/usr/bin/env python2
+# -*- coding: utf-8 -*-
+
+
+def plot_figure_6():
+    """Generates figure 6 of the paper. This figure traces the hourly values of absorbed irradiance, leaf temeprature,
+    net plant carbon assimilation rate, and net plant transpiration rate.
+    """
+
+    training_color = {'gdc': 'blue', 'vsp': 'red', 'lyre': 'green'}
+
+    beg_date = datetime(2009, 7, 29, 00, 00, 0, )
+    end_date = datetime(2009, 7, 29, 23, 00, 0, )
+    datet = pd.date_range(beg_date, end_date, freq='H')
+
+    meteo_df = pd.read_csv(example_pth / 'virtual_canopies' / 'gdc' / 'meteo.input',
+                           sep=';', decimal='.', index_col='time')  # all simus have the same meteo data
+    meteo_df.index = pd.DatetimeIndex(meteo_df.index)
+    meteo_df = meteo_df.loc[datet]
+
+    fig, axs = pyplot.subplots(nrows=2, ncols=2, sharex=True, figsize=(6.69, 6))
+    [ax.grid() for ax in axs.flatten()]
+
+    axs[1, 0].plot(datet, meteo_df['Tac'], 'k--')
+
+    for training in ('gdc', 'vsp', 'lyre'):
+        pth = example_pth / 'virtual_canopies' / training
+
+        sims_df = pd.read_csv(pth / 'output' / 'time_series.output', sep=';', decimal='.', index_col='time')
+        sims_df.index = [datetime.strptime(s, "%Y-%m-%d %H:%M:%S") for s in sims_df.index]
+
+        axs[0, 0].plot(datet, sims_df['Rg'], label=training, color=training_color[training])
+        axs[0, 1].plot(datet, sims_df['An'], label=training, color=training_color[training])
+        axs[1, 0].plot(datet, sims_df['Tleaf'], label=training, color=training_color[training])
+        axs[1, 1].plot(datet, sims_df['E'], label=training, color=training_color[training])
+
+    # some layout
+    for iax, ax in enumerate(axs.flatten()):
+        ax.text(0.875, 0.9, ('(a)', '(b)', '(c)', '(d)')[iax],
+                transform=ax.transAxes, fontdict={'size': 11})
+
+    axs[0, 0].set(xlim=(beg_date, end_date), ylim=(0, 600),
+                  ylabel='$\mathregular{\Phi_{R_g, plant}\/[W\/m^{-2}_{ground}]}$')
+    axs[1, 0].set(xlabel='hour', ylim=(10, 40),
+                  ylabel='$\mathregular{Temperature\/[^\circ C]}$')
+    axs[0, 1].set(ylim=(-10, 35),
+                  ylabel='$\mathregular{A_{n, plant}\/[\mu mol\/s^{-1}]}$')
+    axs[1, 1].set(xlabel='hour', ylim=(0, 400),
+                  ylabel='$\mathregular{E_{plant}\/[g\/h^{-1}]}$')
+
+    for ax in axs[1, :]:
+        ax.set(xlim=(beg_date, end_date))
+        ax.set_xticklabels(datet, rotation=90)
+        ax.xaxis.set_major_locator(dates.HourLocator())
+        ax.xaxis.set_major_formatter(dates.DateFormatter('%H'))
+
+    axs[0, 0].legend(loc='upper left', prop={'size': 11})
+    handles, _ = axs[1, 0].get_legend_handles_labels()
+    axs[1, 0].legend(handles=[handles[0], ], labels=('$\mathregular{T_{air}}$',), loc='upper left',
+                     prop={'size': 11})
+
+    axs[1, 0].set_xticks(axs[1, 0].get_xticks()[:-1:2])
+
+    fig.tight_layout()
+
+    fig.savefig('fig_6.png', dpi=600.)
+    pyplot.close(fig)
+
+
+def plot_figure_7():
+    """Generates figure 7 of the paper. This figure shows a snapshot at solar midday of water potential distribution
+    across three virtual canopies.
+    """
+    training_color = {'gdc': 'blue', 'vsp': 'red', 'lyre': 'green'}
+
+    obs_date = datetime(2009, 7, 29, 14, 00, 0, )
+
+    fig, axs = pyplot.subplots(nrows=3, ncols=2, sharex=True, figsize=(6.69, 6),
+                               gridspec_kw={'width_ratios': [0.8, 0.2]})
+
+    for i, training in enumerate(('vsp', 'gdc', 'lyre')):
+        pth = example_pth / 'virtual_canopies' / training
+
+        g, _ = mtg_load(pth / 'output' / f"mtg{obs_date.strftime('%Y%m%d%H%M%S')}")
+
+        axs[i, 0] = display.property_map(g, 'psi_head', color=training_color[training],
+                                         ax=axs[i, 0])
+        axs[i, 0].set(ylim=(0, 250), xlim=(-1.7, -1.2))
+        axs[i, 0].legend([training], prop={'size': 11})
+        axs[i, 0].xaxis.labelpad = 5
+
+        axs[i, 1].boxplot(g.property('psi_head').values(), vert=False, sym='')
+
+    # some layout
+    [ax.set_xlabel('') for ax in axs[:2, 0]]
+
+    for ax in axs[:, 1]:
+        ax.tick_params(axis='y', which='both', labelleft='off')
+        ax.set_xticklabels(np.arange(-1.7, -1.1, 0.1), rotation=90)
+
+    axs[2, 1].set_xlabel('$\mathregular{\Psi\/[MPa]}$')
+
+    fig.tight_layout()
+
+    fig.savefig('fig_7.png', dpi=600.)
+    pyplot.close(fig)
+
+
+def plot_figure_8():
+    """Generates figure 8 of the paper. This figure shows a snapshot at solar midday of leaf temperature distribution
+    across three virtual canopies.
+    """
+    training_color = {'gdc': 'blue', 'vsp': 'red', 'lyre': 'green'}
+
+    obs_date = datetime(2009, 7, 29, 14, 00, 0, )
+
+    fig, axs = pyplot.subplots(nrows=3, ncols=2, sharex=True, figsize=(6.69, 6),
+                               gridspec_kw={'width_ratios': [0.8, 0.2]})
+
+    for i, training in enumerate(('vsp', 'gdc', 'lyre')):
+        pth = example_pth / 'virtual_canopies' / training
+
+        g, _ = mtg_load(pth / 'output' / f"mtg{obs_date.strftime('%Y%m%d%H%M%S')}")
+
+        axs[i, 0] = display.property_map(g, 'Tlc', color=training_color[training],
+                                         ax=axs[i, 0])
+        axs[i, 0].set(ylim=(0, 250))
+        axs[i, 0].legend([training], prop={'size': 11})
+        axs[i, 0].xaxis.labelpad = 5
+
+        axs[i, 1].boxplot(g.property('Tlc').values(), vert=False, sym='')
+
+    # some layout
+    [ax.set_xlabel('') for ax in axs[:2, 0]]
+
+    for ax in axs[:, 1]:
+        ax.tick_params(axis='y', which='both', labelleft='off')
+        ax.set_xticklabels(range(32, 42), rotation=90)
+
+    axs[2, 1].set_xlabel('$\mathregular{T_{leaf}\/[^{\circ}C]}$')
+
+    fig.tight_layout()
+
+    fig.savefig('fig_8.png', dpi=600.)
+    pyplot.close(fig)
+
+
+def plot_figure_9():
+    """Generates figure 9 of the paper. This figure compares, simulated to observed xylem water potential.
+    """
+
+    beg_date = datetime(2012, 8, 1, 00, 00, 0, )
+    end_date = datetime(2012, 8, 3, 00, 00, 0, )
+    datet = pd.date_range(beg_date, end_date, freq='D')
+
+    fig, axs = pyplot.subplots(nrows=3, ncols=3, sharex='all', sharey='all', figsize=(6.69, 6))
+
+    for it, training in enumerate(('gdc_can1_grapevine', 'gdc_can2_grapevine', 'gdc_can3_grapevine')):
+
+        pth = example_pth / training
+
+        for i_day, date in enumerate(datet):
+            ax = axs[it, i_day]
+            obs_date = date + pd.Timedelta(hours=13)
+            g, _ = mtg_load(pth / 'output' / f"mtg{obs_date.strftime('%Y%m%d%H%M%S')}")
+            ax = display.property_map(g, 'psi_head', ax=ax, prop2='Eabs', color='grey',
+                                      colormap='autumn', add_color_bar=False)
+
+            obs_df = pd.read_csv(pth / 'var.obs', sep=';', index_col='date')
+            obs_df.index = pd.DatetimeIndex(obs_df.index)
+            psi_leaf = obs_df.loc[date, 'psi_leaf']
+            ax.add_patch(patches.Rectangle((max(psi_leaf), 50),
+                                           min(psi_leaf) - max(psi_leaf), 250,
+                                           color='0.8', zorder=-1))
+            if i_day == 0:
+                ax.text(0.05, 0.075, 'Canopy%d' % (it + 1),
+                        transform=ax.transAxes, fontdict={'size': 11})
+
+            if it == 0:
+                ax.text(0.675, 0.825, obs_date.strftime('%-d %b'),
+                        transform=ax.transAxes, fontdict={'size': 11})
+
+    [axi.set_xticklabels(axi.get_xticks(), rotation=90) for axi in axs[2]]
+
+    for ax, axi in enumerate(axs):
+        if ax < 2:
+            [axi[i_day].set_xlabel('') for i_day in range(3)]
+        [axi[i_day].set_ylabel('') for i_day in (1, 2)]
+        [axi[i_day].legend_.remove() for i_day in range(3)]
+
+    fig.tight_layout()
+
+    fig.subplots_adjust(bottom=0.3)
+    cbar_ax = fig.add_axes([0.395, 0.1, 0.30, 0.04])
+    norm = colors.Normalize(0, vmax=2000)
+    cbar = colorbar.ColorbarBase(cbar_ax, cmap='autumn', orientation='horizontal', norm=norm)
+    cbar.ax.set_xticklabels(cbar.ax.get_xticklabels(), rotation=90)
+    cbar.set_label('$\mathregular{PPFD\/[\mu mol\/m^{-2}_{leaf}\/s^{-1}]}$', labelpad=-20, x=-0.4)
+
+    fig.savefig('fig_9.png', dpi=600.)
+    pyplot.close(fig)
+
+
+def plot_figure_10():
+    """Generates figure 10 of the paper. This figure compares, simulated to observed stomatal conductance rates.
+    """
+
+    beg_date = datetime(2012, 8, 1, 00, 00, 0, )
+    end_date = datetime(2012, 8, 3, 00, 00, 0, )
+    datet = pd.date_range(beg_date, end_date, freq='D')
+
+    fig, axs = pyplot.subplots(nrows=3, ncols=3, sharex='all', sharey='all', figsize=(6.69, 6))
+
+    for it, training in enumerate(('gdc_can1_grapevine', 'gdc_can2_grapevine', 'gdc_can3_grapevine')):
+
+        pth = example_pth / training
+
+        for i_day, date in enumerate(datet):
+            ax = axs[it, i_day]
+            obs_date = date + pd.Timedelta(hours=13)
+            g, _ = mtg_load(pth / 'output' / f"mtg{obs_date.strftime('%Y%m%d%H%M%S')}")
+            ax = display.property_map(g, 'gs', ax=ax, prop2='Eabs', color='grey',
+                                      colormap='autumn', add_color_bar=False)
+
+            obs_df = pd.read_csv(pth / 'var.obs', sep=';', index_col='date')
+            obs_df.index = pd.DatetimeIndex(obs_df.index)
+            gs_leaf = obs_df.loc[date, 'gs'] / 1000.
+            ax.add_patch(patches.Rectangle((max(gs_leaf), 50),
+                                           min(gs_leaf) - max(gs_leaf), 250,
+                                           color='0.8', zorder=-1))
+            if i_day == 0:
+                ax.text(0.05, 0.075, 'Canopy%d' % (it + 1),
+                        transform=ax.transAxes, fontdict={'size': 11})
+
+            if it == 0:
+                ax.text(0.675, 0.825, obs_date.strftime('%-d %b'),
+                        transform=ax.transAxes, fontdict={'size': 11})
+
+    [axi.set_xticklabels(axi.get_xticks(), rotation=90) for axi in axs[2]]
+
+    for ax, axi in enumerate(axs):
+        if ax < 2:
+            [axi[i_day].set_xlabel('') for i_day in range(3)]
+        [axi[i_day].set_ylabel('') for i_day in (1, 2)]
+        [axi[i_day].legend_.remove() for i_day in range(3)]
+
+    fig.tight_layout()
+
+    fig.subplots_adjust(bottom=0.3)
+    cbar_ax = fig.add_axes([0.395, 0.1, 0.30, 0.04])
+    norm = colors.Normalize(0, vmax=2000)
+    cbar = colorbar.ColorbarBase(cbar_ax, cmap='autumn', orientation='horizontal', norm=norm)
+    cbar.ax.set_xticklabels(cbar.ax.get_xticklabels(), rotation=90)
+    cbar.set_label('$\mathregular{PPFD\/[\mu mol\/m^{-2}_{leaf}\/s^{-1}]}$', labelpad=-20, x=-0.4)
+
+    fig.savefig('fig_10.png', dpi=600.)
+    pyplot.close(fig)
+
+
+def plot_figure_11():
+    """Generates figure 11 of the paper. This figure compares simulated to observed leaf temperatures.
+    """
+    fig, axs = pyplot.subplots(nrows=2, ncols=2, sharey='row', figsize=(6.69, 6))
+
+    for iax, training in enumerate(('vsp_ww_grapevine', 'vsp_ws_grapevine')):
+        pth = example_pth / training
+        axs[0, iax].grid()
+        axs[1, iax].grid()
+
+        # set dates ************************
+        beg_date = datetime(2009, 7, 29, 00, 00, 0, )
+        end_date = datetime(2009, 8, 1, 23, 00, 0, )
+        datet = pd.date_range(beg_date, end_date, freq='H')
+
+        # read observations
+        obs = pd.read_csv(pth / 'temp.obs', sep=';', decimal='.', index_col='date')
+        obs.index = [datetime.strptime(s, "%d/%m/%Y %H:%M") for s in obs.index]
+
+        # read simulations
+        sims = pd.read_csv(pth / 'output' / 'time_series.output',
+                           sep=';', decimal='.', index_col='time')
+        sims.index = [datetime.strptime(s, "%Y-%m-%d %H:%M:%S") for s in sims.index]
+        sims.index = pd.DatetimeIndex(sims.index)
+
+        # plot median simulated leaf temperature
+        axs[0, iax].plot(sims['Tleaf'], '-k', label='$\mathregular{T_{leaf}}$', linewidth=1)
+
+        # plot simulated temperature magnitude
+        med_sim = np.array([])
+        q1_sim = np.array([])
+        q3_sim = np.array([])
+        for date in datet:
+            g, _ = mtg_load(pth / 'output' / f"mtg{date.strftime('%Y%m%d%H%M%S')}")
+            leaf_temp_sim = g.property('Tlc').values()
+            q1_sim = np.append(q1_sim, min(leaf_temp_sim))
+            q3_sim = np.append(q3_sim, max(leaf_temp_sim))
+            med_sim = np.append(med_sim, np.median(leaf_temp_sim))
+            print(date)
+        axs[0, iax].fill_between(datet, q1_sim, q3_sim, color='red', alpha=0.5, zorder=0)
+
+        # plot observed temperature magnitude
+        med_obs = np.array([])
+        q1_obs = np.array([])
+        q3_obs = np.array([])
+        for row, date in enumerate(datet):
+            pos = date.toordinal() + date.hour / 24.
+            leaf_temp_obs = obs.loc[date, ['Tleaf%d' % d for d in (2, 3, 4, 5, 6, 8, 9, 10)]]
+            leaf_temp_obs = leaf_temp_obs[~np.isnan(leaf_temp_obs)]
+            axs[0, iax].boxplot(leaf_temp_obs, positions=[pos], widths=0.01)
+            q1_obs = np.append(q1_obs, min(leaf_temp_obs))
+            q3_obs = np.append(q3_obs, max(leaf_temp_obs))
+            med_obs = np.append(med_obs, np.median(leaf_temp_obs))
+            print(date)
+
+        # compare obs to sim temperature on 1:1 plot
+        axs[1, iax].plot((10, 50), (10, 50), 'k--')
+        axs[1, iax].errorbar(x=med_obs, y=med_sim,
+                             xerr=(np.nan_to_num(med_obs - q1_obs), np.nan_to_num(q3_obs - med_obs)),
+                             yerr=(np.nan_to_num(med_sim - q1_sim), np.nan_to_num(q3_sim - med_sim)),
+                             fmt='ro', ecolor='0.5', capthick=1)
+
+        # plot MBE and RMSE results on it
+        diff_obs_sim = (med_sim - med_obs)[~np.isnan(med_sim - med_obs)]
+        bias = diff_obs_sim.mean()
+        rmse = np.sqrt(diff_obs_sim ** 2).mean()
+
+        axs[1, iax].text(0.50, 0.20, '$\mathregular{MBE\/=\/%.3f}$' % bias,
+                         transform=axs[1, iax].transAxes, fontdict={'size': 11})
+        axs[1, iax].text(0.50, 0.10, '$\mathregular{RMSE\/=\/%.1f}$' % rmse,
+                         transform=axs[1, iax].transAxes, fontdict={'size': 11})
+
+    # some layout
+
+    axs[0, 0].set_ylabel('$\mathregular{Temperature\/[^{\circ}C]}$')
+    axs[1, 0].set_ylabel('$\mathregular{T_{leaf, sim}\/[^{\circ}C]}$')
+
+    for ax_upper, ax_lower in zip(axs[0, :], axs[1, :]):
+        ax_upper.set(xlim=(beg_date, end_date), ylim=(5, 50), xlabel='')
+        ax_upper.set_xticklabels(datet, rotation=45)
+        ax_upper.xaxis.set_major_locator(dates.DayLocator())
+        ax_upper.xaxis.set_major_formatter(dates.DateFormatter('%-d %b'))
+        ax_lower.set(xlim=(10, 50), ylim=(10, 50),
+                     xlabel='$\mathregular{T_{leaf, obs}\/[^{\circ}C]}$')
+
+    for iax, ax in enumerate(axs.flatten()):
+        ax.text(0.05, 0.875, ('(a)', '(c)', '(b)', '(d)')[iax],
+                transform=ax.transAxes, fontdict={'size': 11})
+
+    fig.tight_layout()
+    fig.savefig('fig_11.png', dpi=600.)
+    pyplot.close(fig)
+
+
+def plot_figure_12():
+    """Generates figure 12 of the paper. This figure compares, leaf-to-leaf, simulated to observed leaf and leaf-to-air
+    temperature.
+    """
+
+    fig, axs = pyplot.subplots(ncols=3, figsize=(6.69, 2.7))
+    [ax.grid() for ax in axs]
+
+    daily_temp_obs = np.array([])
+    daily_temp_sim = np.array([])
+    daily_temp_air = np.array([])
+
+    for iax, training in enumerate(('vsp_ww_grapevine', 'vsp_ws_grapevine')):
+
+        pth = example_pth / training
+
+        # set dates ************************
+        beg_date = datetime(2009, 7, 29, 00, 00, 0, )
+        end_date = datetime(2009, 8, 1, 23, 00, 0, )
+        datet = pd.date_range(beg_date, end_date, freq='H')
+
+        # read observations
+        obs = pd.read_csv(pth / 'temp.obs', sep=';', decimal='.', index_col='date')
+        obs.index = [datetime.strptime(s, "%d/%m/%Y %H:%M") for s in obs.index]
+
+        # read simulations
+        sims = pd.read_csv(pth / 'output' / 'time_series.output',
+                           sep=';', decimal='.', index_col='time')
+        sims.index = [datetime.strptime(s, "%Y-%m-%d %H:%M:%S") for s in sims.index]
+        sims.index = pd.DatetimeIndex(sims.index)
+
+        # plot simulated temperature magnitude
+        med_sim = np.array([])
+        q_1_sim = np.array([])
+        q_3_sim = np.array([])
+        for date in datet:
+            g, _ = mtg_load(pth / 'output' / f"mtg{date.strftime('%Y%m%d%H%M%S')}")
+            leaf_temp_sim = g.property('Tlc').values()
+            q_1_sim = np.append(q_1_sim, min(leaf_temp_sim))
+            q_3_sim = np.append(q_3_sim, max(leaf_temp_sim))
+            med_sim = np.append(med_sim, np.median(leaf_temp_sim))
+            print(date)
+
+        # get observed temperature magnitude
+        med_obs = np.array([])
+        q1_obs = np.array([])
+        q3_obs = np.array([])
+        for row, date in enumerate(datet):
+            pos = date.toordinal() + date.hour / 24.
+            leaf_temp_obs = obs.loc[date, ['Tleaf%d' % d for d in (2, 3, 4, 5, 6, 8, 9, 10)]]
+            leaf_temp_obs = leaf_temp_obs[~np.isnan(leaf_temp_obs)]
+            q1_obs = np.append(q1_obs, min(leaf_temp_obs))
+            q3_obs = np.append(q3_obs, max(leaf_temp_obs))
+            med_obs = np.append(med_obs, np.median(leaf_temp_obs))
+            print(date)
+
+        # read meteo data
+        meteo_df = pd.read_csv(pth / 'meteo.input', sep=';', decimal='.', index_col='time')
+        meteo_df.index = pd.DatetimeIndex(meteo_df.index)
+
+        air_temperature = meteo_df.loc[datet, 'Tac']
+
+        daily_temp_obs = np.append(daily_temp_obs, q1_obs)
+        daily_temp_obs = np.append(daily_temp_obs, q3_obs)
+
+        daily_temp_sim = np.append(daily_temp_sim, q_1_sim)
+        daily_temp_sim = np.append(daily_temp_sim, q_3_sim)
+
+        daily_temp_air = np.append(daily_temp_air, air_temperature)
+        daily_temp_air = np.append(daily_temp_air, air_temperature)
+
+        # minimum temperature
+        axs[0].plot(q1_obs - air_temperature, q_1_sim - air_temperature,
+                    ['b^', 'r^'][iax], label=['WW', 'WD'][iax])
+        # maximum temperature
+        axs[0].plot(q3_obs - air_temperature, q_3_sim - air_temperature,
+                    ['bo', 'ro'][iax], label=['WW', 'WD'][iax])
+
+        axs[1].plot(q1_obs, q_1_sim, ['b^', 'r^'][iax], label=['WW', 'WD'][iax])
+        axs[1].plot(q3_obs, q_3_sim, ['bo', 'ro'][iax], label=['WW', 'WD'][iax])
+
+        axs[2].plot(q3_obs - q1_obs, q_3_sim - q_1_sim, ['bs', 'rs'][iax], label=['WW', 'WD'][iax])
+
+    # some layout
+    axs[0].plot((-8, 12), (-8, 12), 'k--')
+    axs[0].set(xlabel='$\mathregular{T_{leaf, obs}-T_{air}\/[^\circ C]}$',
+               ylabel='$\mathregular{T_{leaf, sim}-T_{air}\/[^\circ C]}$',
+               xlim=(-8, 12), ylim=(-8, 12))
+
+    axs[1].plot((10, 45), (10, 45), 'k--')
+    axs[1].set(xlabel='$\mathregular{T_{leaf, obs}\/[^\circ C]}$',
+               ylabel='$\mathregular{T_{leaf, sim}\/[^\circ C]}$')
+    axs[1].xaxis.set_major_locator(ticker.MultipleLocator(20))
+
+    axs[2].plot((-2, 14), (-2, 14), 'k--')
+    axs[2].set(xlabel='$\mathregular{\Delta T_{leaf, obs}\/[^\circ C]}$',
+               ylabel='$\mathregular{\Delta T_{leaf, sim}\/[^\circ C]}$')
+    axs[2].xaxis.set_major_locator(ticker.IndexLocator(base=2, offset=2))
+    axs[2].xaxis.set_major_locator(ticker.MultipleLocator(5))
+
+    for i in range(3):
+        if i == 0:
+            x, y = daily_temp_obs - daily_temp_air, daily_temp_sim - daily_temp_air
+        elif i == 1:
+            x, y = daily_temp_obs, daily_temp_sim
+        else:
+            x = np.array([])
+            x = np.append(x, daily_temp_obs[96:2 * 96] - daily_temp_obs[:96])
+            x = np.append(x, daily_temp_obs[3 * 96:4 * 96] - daily_temp_obs[2 * 96:3 * 96])
+            y = np.array([])
+            y = np.append(y, daily_temp_sim[96:2 * 96] - daily_temp_sim[:96])
+            y = np.append(y, daily_temp_sim[3 * 96:4 * 96] - daily_temp_sim[2 * 96:3 * 96])
+
+        diff_y_x = (y - x)[~np.isnan(y - x)]
+        bias = (diff_y_x).mean()
+        rmse = np.sqrt((diff_y_x ** 2).mean())
+
+        axs[i].text(0.05, -0.6, '$\mathregular{MBE\/=\/%.3f}$' % bias,
+                    transform=axs[i].transAxes, fontdict={'size': 11})
+        axs[i].text(0.05, -0.7, '$\mathregular{RMSE\/=\/%.1f}$' % rmse,
+                    transform=axs[i].transAxes, fontdict={'size': 11})
+
+        axs[i].text(0.05, 0.875, ['(a)', '(b)', '(c)'][i],
+                    transform=axs[i].transAxes, fontdict={'size': 11})
+
+    fig.tight_layout()
+    fig.subplots_adjust(bottom=0.4)
+    fig.savefig('fig_12.png', dpi=600.)
+    pyplot.close(fig)
+
+
+def plot_figure_13():
+    """Generates figure 13 of the paper.
+    """
+    # set dates
+    beg_date = datetime(2009, 7, 29, 00, 00, 0, )
+    end_date = datetime(2009, 8, 1, 23, 00, 0, )
+    datet = pd.date_range(beg_date, end_date, freq='H')
+
+    fig = pyplot.figure()
+    gs1 = gridspec.GridSpec(2, 2)
+    gs1.update(left=0.135, right=0.65, top=0.975, bottom=0.125, wspace=0.1, hspace=0.35)
+    ax1 = pyplot.subplot(gs1[0, 0])
+    ax2 = pyplot.subplot(gs1[0, 1], sharex=ax1)
+    ax3 = pyplot.subplot(gs1[1, 0], sharex=ax1)
+    ax4 = pyplot.subplot(gs1[1, 1], sharex=ax1)
+
+    gs2 = gridspec.GridSpec(2, 1)
+    gs2.update(left=0.67, right=0.865, top=0.975, bottom=0.125, hspace=0.35)
+    ax5 = pyplot.subplot(gs2[0])
+    ax6 = pyplot.subplot(gs2[1])
+
+    fig.set_figheight(6)
+    fig.set_figwidth(6.69)
+    axs = np.array([[ax1, ax2, ax5],
+                    [ax3, ax4, ax6]])
+
+    [ax.grid() for ax in axs.flatten()]
+
+    for iax, training in enumerate(('vsp_ww_grapevine', 'vsp_ws_grapevine')):
+        pth = example_pth / training
+
+        # read observations
+        obs_tab = pd.read_csv(pth / 'gas.obs', sep=';', decimal='.', header=0)
+        obs_tab.time = pd.DatetimeIndex(obs_tab.time)
+        obs_tab = obs_tab.set_index(obs_tab.time)
+
+        axs[0, iax].plot(obs_tab['time'], obs_tab['An_plante'],
+                         color='0.8', marker='o', markeredgecolor='none', label='$\mathregular{A_{n,\/obs}}$')
+        axs[1, iax].plot(obs_tab['time'], obs_tab['E_plante'],
+                         color='0.8', marker='o', markeredgecolor='none', label='$\mathregular{E_{obs}}$')
+
+        # read simulations
+        sims = pd.read_csv(pth / 'output' / 'time_series.output',
+                           sep=';', decimal='.', index_col='time')
+        sims.index = [datetime.strptime(s, "%Y-%m-%d %H:%M:%S") for s in sims.index]
+        sims.index = pd.DatetimeIndex(sims.index)
+
+        axs[0, iax].plot(sims['An'],
+                         'b-', label='$\mathregular{A_{n,\/sim}}$', linewidth=1)
+        axs[1, iax].plot(sims['E'],
+                         'b-', label='$\mathregular{E_{sim}}$', linewidth=1)
+
+        sims['itime'] = sims.index
+        obs_tab['time'] = pd.DatetimeIndex(obs_tab.time)
+
+        time_group = []
+        for itime in obs_tab.time:
+            if itime.minute < 30:
+                time = itime - pd.Timedelta(minutes=itime.minute)
+            else:
+                time = itime + pd.Timedelta(minutes=60 - itime.minute)
+            time_group.append(time)
+        obs_tab['itime'] = time_group
+
+        obs_grouped = obs_tab.groupby('itime').aggregate(np.mean)
+        obs_grouped['itime'] = obs_grouped.index
+
+        m_df = pd.merge(obs_grouped, sims)
+
+        axs[0, 2].plot(m_df['An_plante'], m_df['An'], ('bo', 'ro')[iax])
+        axs[1, 2].plot(m_df['E_plante'], m_df['E'], ('bo', 'ro')[iax])
+
+        for egas, igas in enumerate(('An', 'E')):
+            x = m_df[igas + '_plante'].values
+            y = m_df[igas].values
+
+            xindex = np.isfinite(x)
+            yindex = np.isfinite(y)
+            xyindex = xindex * yindex
+            x, y = x[xyindex], y[xyindex]
+
+            bias = (y - x).mean()
+            rmse = np.sqrt(((x - y) ** 2).mean())
+
+            axs[egas, iax].text(0.25, 0.90,
+                                '$\mathregular{MBE\/=\/%.1f}$' % bias,
+                                transform=axs[egas, iax].transAxes,
+                                fontdict={'size': 11})
+            axs[egas, iax].text(0.25, 0.8,
+                                '$\mathregular{RMSE\/=\/%.1f}$' % rmse,
+                                transform=axs[egas, iax].transAxes,
+                                fontdict={'size': 11})
+
+    axs[0, 2].plot((-20, 50), (-20, 50), 'k--')
+    axs[1, 2].plot((-200, 1000), (-200, 1000), 'k--')
+
+    for i, ax in enumerate(axs.flatten()):
+        ax.text(0.05, 0.9, ('(a)', '(b)', '(e)', '(c)', '(d)', '(f)')[i],
+                transform=ax.transAxes, fontdict={'size': 11})
+
+    # some layout
+    axs[0, 0].set(
+        ylabel='$\mathregular{A_{n, plant}\/[\mu mol\/s^{-1}]}$',
+        xlim=(beg_date, end_date),
+        ylim=(-20, 50))
+    axs[1, 0].set(
+        ylabel='$\mathregular{E_{plant}\/[g\/h^{-1}]}$',
+        xlim=(beg_date, end_date),
+        ylim=(-200, 1000))
+    axs[0, 1].set(
+        xlim=(beg_date, end_date),
+        ylim=(-20, 50))
+    axs[1, 1].set(
+        xlim=(beg_date, end_date),
+        ylim=(-200, 1000))
+    axs[0, 2].set(
+        xlabel='$\mathregular{A_{n, plant, obs}\/[\mu mol\/s^{-1}]}$',
+        ylabel='$\mathregular{A_{n, plant, sim}\/[\mu mol\/s^{-1}]}$',
+        xlim=(-20, 50),
+        ylim=(-20, 50))
+    axs[1, 2].set(
+        xlabel='$\mathregular{E_{plant, obs}\/[g\/h^{-1}]}$',
+        ylabel='$\mathregular{E_{plant, sim}\/[g\/h^{-1}]}$',
+        xlim=(-200, 1000),
+        ylim=(-200, 1000))
+
+    axs[0, 1].get_yaxis().set_ticklabels('')
+    axs[1, 1].get_yaxis().set_ticklabels('')
+    axs[1, 1].get_yaxis().set_ticklabels('')
+    axs[0, 2].yaxis.tick_right()
+    axs[1, 2].yaxis.tick_right()
+    axs[0, 2].yaxis.set_label_position('right')
+    axs[1, 2].yaxis.set_label_position('right')
+    axs[0, 2].xaxis.set_major_locator(ticker.MultipleLocator(25))
+    axs[1, 2].xaxis.set_major_locator(ticker.MultipleLocator(500))
+
+    for ax in axs[:, :2].flatten():
+        ax.set_xticklabels(datet, rotation=90)
+        ax.xaxis.set_major_locator(dates.DayLocator())
+        ax.xaxis.set_major_formatter(dates.DateFormatter('%-d %b'))
+
+    fig.savefig('fig_13.png', dpi=600.)
+    pyplot.close(fig)
+
+
+def plot_figure_14():
+    """Generates figure 14 of the paper. This figure compares, simulated to observed plant transpiration rates.
+    """
+
+    fig, axs = pyplot.subplots(nrows=3, ncols=4, sharey='row', figsize=(6.69, 6))
+
+    for i_treat, training in enumerate(('gdc_can1_grapevine', 'gdc_can2_grapevine', 'gdc_can3_grapevine')):
+
+        pth = example_pth / training
+
+        # read observations
+        obs_df = pd.read_csv(pth / 'sapflow.obs', sep=';', decimal='.', index_col='date')
+        obs_df.index = [datetime.strptime(s, "%d/%m/%Y %H:%M") for s in obs_df.index]
+
+        time_group = []
+        for itime in obs_df.index:
+            if itime.minute < 30:
+                time = itime - pd.Timedelta(minutes=itime.minute)
+            else:
+                time = itime + pd.Timedelta(minutes=60 - itime.minute)
+            time_group.append(time)
+        obs_df['itime'] = time_group
+
+        obs_grouped = obs_df.groupby('itime').aggregate(np.mean)
+        obs_grouped['itime'] = obs_grouped.index
+
+        # read simulations
+        sims = pd.read_csv(pth / 'output' / 'time_series.output',
+                           sep=';', decimal='.', index_col='time')
+        sims.index = [datetime.strptime(s, "%Y-%m-%d %H:%M:%S") for s in sims.index]
+        sims.index = pd.DatetimeIndex(sims.index)
+
+        sims['itime'] = sims.index
+        m_df = pd.merge(obs_grouped, sims)
+
+        x = m_df['west'].values + m_df['east'].values
+        y = m_df['sapWest'].values + m_df['sapEast'].values
+
+        x_index = np.isfinite(x)
+        y_index = np.isfinite(y)
+        xy_index = x_index * y_index
+        x, y = x[xy_index], y[xy_index]
+
+        for iax, ax in enumerate(axs[i_treat, :]):
+            ax.xaxis.grid(which='minor', zorder=0)
+            ax.yaxis.grid(which='major', zorder=0)
+
+            ax.plot(sims.index, sims.sapEast, c='r', linewidth=1)
+            ax.plot(sims.index, sims.sapWest, c='b', linewidth=1)
+
+            ax.plot(obs_df.index, obs_df['east'], label='East', color='red',
+                    marker='o', markeredgecolor='none', alpha=0.1)
+            ax.plot(obs_df.index, obs_df['west'], label='West', color='blue',
+                    marker='o', markeredgecolor='none', alpha=0.1)
+
+            x_t = x[iax * 24:iax * 24 + 23]
+            y_t = y[iax * 24:iax * 24 + 23]
+
+            x_index = np.isfinite(x_t)
+            y_index = np.isfinite(y_t)
+            xy_index = x_index * y_index
+            x_t, y_t = x_t[xy_index], y_t[xy_index]
+
+            bias = (y_t - x_t).mean()
+            rmse = np.sqrt(((x_t - y_t) ** 2).mean())
+
+            ax.text(0.45, 0.675,
+                    '$\mathregular{MBE\/=\/%.1f}$' % bias,
+                    transform=ax.transAxes, fontdict={'size': 7})
+            ax.text(0.45, 0.575,
+                    '$\mathregular{RMSE\/=\/%.1f}$' % rmse,
+                    transform=ax.transAxes, fontdict={'size': 7})
+
+    # some layout
+    for day in range(4):
+        day_sdate = datetime(2012, 8, 1, 00, 00, 0, ) + timedelta(days=day)
+        day_edate = day_sdate + timedelta(hours=23)
+        for can in range(3):
+            axs[can, day].set_xlim(day_sdate, day_edate)
+            axs[can, day].yaxis.set_major_locator(ticker.MultipleLocator(100))
+            if can != 2:
+                axs[can, day].xaxis.set_ticklabels('')
+                axs[can, day].xaxis.set_minor_locator(dates.HourLocator(interval=5))
+            else:
+                axs[can, day].set_xticklabels(pd.date_range(day_sdate, day_edate, freq='H'),
+                                              rotation=90)
+                axs[can, day].xaxis.set_major_locator(dates.DayLocator())
+                axs[can, day].xaxis.set_major_formatter(dates.DateFormatter('%-d %b'))
+                axs[can, day].xaxis.set_minor_locator(dates.HourLocator(interval=5))
+                axs[can, day].xaxis.set_minor_formatter(dates.DateFormatter('%H'))
+                axs[can, day].tick_params(axis='x', which='major', pad=15)
+
+    axs[1, 0].set_ylabel('$\mathregular{E\/[g\/h^{-1}]}$')
+
+    axs[0, 0].set_ylim(0, 800)
+    axs[1, 0].set_ylim(0, 450)
+    axs[2, 0].set_ylim(0, 500)
+
+    for can in range(3):
+        axs[can, 0].text(0.05, 0.85, 'Canopy%s' % str(can + 1),
+                         transform=axs[can, 0].transAxes, fontdict={'size': 11})
+
+    fig.tight_layout()
+    fig.subplots_adjust(left=0.13, bottom=0.3)
+
+    h1, l1 = axs[2, 2].get_legend_handles_labels()
+
+    axs[2, 0].legend(h1, ('$\mathregular{SapEast_{sim}}$', '$\mathregular{SapWest_{sim}}$',
+                          '$\mathregular{SapEast_{obs}}$', '$\mathregular{SapWest_{obs}}$'),
+                     frameon=True, bbox_to_anchor=(-0.36, -1.2, 2, .102),
+                     loc=3, ncol=8, prop={'size': 11})
+    fig.subplots_adjust(wspace=0.075)
+
+    fig.savefig('fig_14.png', dpi=600.)
+    pyplot.close(fig)
+
+
+def plot_figure_15():
+    """Generates figure 15 of the paper. This figure compares simulated to observed plant photosynthesis and
+    transpiration rates using 4 simulation parameters' sets (modularity test)
+    """
+    # set dates
+    beg_date = datetime(2009, 7, 29, 00, 00, 0, )
+    end_date = datetime(2009, 8, 1, 23, 00, 0, )
+    datet = pd.date_range(beg_date, end_date, freq='H')
+
+    style = ('b-', 'k--', 'k-.', 'k-', 'k:')
+    vpd_air = np.vectorize(VPDa)
+
+    fig, axs = pyplot.subplots(nrows=2, ncols=2, sharey='row', sharex='all', figsize=(6.69, 6))
+    [ax.grid() for ax in axs.flatten()]
+
+    for i_treat, treat in enumerate(('ww', 'ws')):
+
+        pth = example_pth / ('vsp_%s_grapevine' % treat)
+
+        # read and plot observations
+        obs_df = pd.read_csv(pth / 'gas.obs', sep=';', decimal='.')
+        obs_df.time = pd.DatetimeIndex(obs_df.time)
+        obs_df = obs_df.set_index(obs_df.time)
+
+        axs[0, i_treat].plot(obs_df['time'], obs_df['An_plante'],
+                             color='0.8', marker='o', markeredgecolor='none', label='$\mathregular{A_{n,\/obs}}$')
+        axs[1, i_treat].plot(obs_df['time'], obs_df['E_plante'],
+                             color='0.8', marker='o', markeredgecolor='none', label='$\mathregular{E_{obs}}$')
+
+        # read and plot vpd
+        meteo_df = pd.read_csv(pth / 'meteo.input', sep=';', decimal='.', index_col='time')
+        meteo_df.index = pd.DatetimeIndex(meteo_df.index)
+        meteo_df['vpd'] = vpd_air(meteo_df.Tac, meteo_df.Tac, meteo_df.hs)
+        axt = axs[1, i_treat].twinx()
+        axt.plot(datet, meteo_df.loc[datet, 'vpd'], 'r-')
+
+        # read and plot simulations
+        for case in range(5):
+            if case == 0:
+                sim_pth = pth / 'output' / 'time_series.output'
+            else:
+                sim_pth = example_pth / 'modularity' / treat / ('sim_%d' % case) / 'output' / 'time_series.output'
+
+            sims = pd.read_csv(sim_pth, sep=';', decimal='.', index_col='time')
+            sims.index = [datetime.strptime(s, "%Y-%m-%d %H:%M:%S") for s in sims.index]
+            sims.index = pd.DatetimeIndex(sims.index)
+
+            axs[0, i_treat].plot(sims['An'], style[case], label='sim0', linewidth=1)
+            axs[1, i_treat].plot(sims['E'], style[case], label='sim0', linewidth=1)
+
+    for i, ax in enumerate(axs.flatten()):
+        ax.text(0.05, 0.9, ('(a)', '(b)', '(c)', '(d)')[i], transform=ax.transAxes,
+                fontdict={'size': 11})
+
+    axs[0, 0].set(
+        ylabel='$\mathregular{A_{n, plant}\/[\mu mol\/s^{-1}]}$',
+        ylim=(-20, 50))
+    axs[1, 0].set(
+        ylabel='$\mathregular{E_{plant}\/[g\/h^{-1}]}$',
+        xlim=(beg_date, end_date),
+        ylim=(-200, 1600))
+
+    [ax.set_ylim(0, 5) for ax in fig.get_axes()[-2:]]
+    fig.get_axes()[-1].set(ylabel='VPD [kPa]')
+
+    for ax in axs[1, :].flatten():
+        ax.set_xticklabels(datet, rotation=90)
+        ax.xaxis.set_major_locator(dates.DayLocator())
+        ax.xaxis.set_major_formatter(dates.DateFormatter('%-d %b'))
+
+    fig.tight_layout()
+    fig.subplots_adjust(bottom=0.25)
+
+    h1, l1 = axs[1, 1].get_legend_handles_labels()
+    h2, l2 = fig.get_axes()[-1].get_legend_handles_labels()
+
+    axs[1, 0].legend(h1 + h2, ['obs'] + ['sim%d' % d for d in range(5)] + ['VPD'], frameon=True,
+                     bbox_to_anchor=(0.15, -0.75, 2, .102), loc=3, ncol=4,
+                     prop={'size': 11})
+
+    fig.savefig('fig_15.png', dpi=600.)
+    pyplot.close(fig)
+
+
+def write_table_1():
+    """Generates table 1 of the paper. This table compares simulated to observed plant photosynthesis and
+    transpiration rates using 4 simulation parameters' sets (modularity test)
+    """
+    indices = [['ww'] * 5 + ['ws'] * 5,
+               ['sim%d' % d for d in range(5)] * 2]
+    multiple_index = pd.MultiIndex.from_tuples(list(zip(*indices)), names=['treat', 'sim'])
+
+    df = pd.DataFrame(columns=('an_mbe', 'an_rmse', 'e_mbe', 'e_rmse'),
+                      index=multiple_index)
+
+    for i_treat, treat in enumerate(('ww', 'ws')):
+
+        pth = example_pth / ('vsp_%s_grapevine' % treat)
+
+        # read observations
+        obs_df = pd.read_csv(pth / 'gas.obs', sep=';', decimal='.')
+        obs_df.time = pd.DatetimeIndex(obs_df.time)
+        obs_df = obs_df.set_index(obs_df.time)
+
+        # read simulations
+        for case in range(5):
+            if case == 0:
+                sim_pth = pth / 'output' / 'time_series.output'
+            else:
+                sim_pth = example_pth / 'modularity' / treat / ('sim_%d' % case) / 'output' / 'time_series.output'
+
+            sims_df = pd.read_csv(sim_pth, sep=';', decimal='.', index_col='time')
+            sims_df.index = [datetime.strptime(s, "%Y-%m-%d %H:%M:%S") for s in sims_df.index]
+            sims_df.index = pd.DatetimeIndex(sims_df.index)
+
+            # match sims to obs
+            sims_df['itime'] = sims_df.index
+            obs_df['time'] = pd.DatetimeIndex(obs_df.time)
+
+            time_group = []
+            for itime in obs_df.index:
+                if itime.minute < 30:
+                    time = itime - pd.Timedelta(minutes=itime.minute)
+                else:
+                    time = itime + pd.Timedelta(minutes=60 - itime.minute)
+                time_group.append(time)
+            obs_df['itime'] = time_group
+
+            obs_grouped = obs_df.groupby('itime').aggregate(np.mean)
+            obs_grouped['itime'] = obs_grouped.index
+
+            m_df = pd.merge(obs_grouped, sims_df)
+
+            for gas in ('An', 'E'):
+                rate_obs = m_df['%s_plante' % gas].values
+                rate_sim = m_df[gas].values
+
+                x_index = np.isfinite(rate_obs)
+                y_index = np.isfinite(rate_sim)
+                xy_index = x_index * y_index
+                rate_obs, rate_sim = rate_obs[xy_index], rate_sim[xy_index]
+
+                bias = (rate_sim - rate_obs).mean()
+                rmse = np.sqrt(((rate_obs - rate_sim) ** 2).mean())
+
+                df.loc[(treat, 'sim%d' % case), '%s_mbe' % gas.lower()] = bias
+                df.loc[(treat, 'sim%d' % case), '%s_rmse' % gas.lower()] = rmse
+
+    print(df)
+    df.to_csv('table_1.csv')
+
+
+def estimate_energy_balance_contribution():
+    """Generates table 1 of the paper. This table compares simulated to observed plant photosynthesis and
+    transpiration rates using 4 simulation parameters' sets (modularity test)
+    """
+    indices = [['ww'] * 4 + ['ws'] * 4,
+               ['day%d' % d for d in range(1, 5)] * 2]
+    multiple_index = pd.MultiIndex.from_tuples(list(zip(*indices)), names=['treat', 'day'])
+
+    df = pd.DataFrame(columns=('e_obs', 'e_sim0', 'e_sim3'),
+                      index=multiple_index)
+
+    for i_treat, treat in enumerate(('ww', 'ws')):
+
+        pth = example_pth / ('vsp_%s_grapevine' % treat)
+
+        # read observations
+        obs_df = pd.read_csv(pth / 'gas.obs', sep=';', decimal='.')
+        obs_df.time = pd.DatetimeIndex(obs_df.time)
+        obs_df = obs_df.set_index(obs_df.time)
+
+        # read simulations
+        for case in (0, 3):
+            if case == 0:
+                sim_pth = pth / 'output' / 'time_series.output'
+            else:
+                sim_pth = example_pth / 'modularity' / treat / ('sim_%d' % case) / 'output' / 'time_series.output'
+
+            sims_df = pd.read_csv(sim_pth, sep=';', decimal='.', index_col='time')
+            sims_df.index = [datetime.strptime(s, "%Y-%m-%d %H:%M:%S") for s in sims_df.index]
+            sims_df.index = pd.DatetimeIndex(sims_df.index)
+
+            # match sims to obs
+            sims_df['itime'] = sims_df.index
+            obs_df['time'] = pd.DatetimeIndex(obs_df.time)
+
+            time_group = []
+            for itime in obs_df.index:
+                if itime.minute < 30:
+                    time = itime - pd.Timedelta(minutes=itime.minute)
+                else:
+                    time = itime + pd.Timedelta(minutes=60 - itime.minute)
+                time_group.append(time)
+            obs_df['itime'] = time_group
+
+            obs_grouped = obs_df.groupby('itime').aggregate(np.mean)
+            obs_grouped['itime'] = obs_grouped.index
+
+            m_df = pd.merge(obs_grouped, sims_df)
+
+            rate_obs = m_df['E_plante'].values
+            rate_sim = m_df['E'].values
+
+            for day in range(4):
+                rate_obs_daily = rate_obs[day * 24: day * 24 + 24]
+                rate_sim_daily = rate_sim[day * 24: day * 24 + 24]
+
+                x_index = np.isfinite(rate_obs_daily)
+                y_index = np.isfinite(rate_sim_daily)
+                xy_index = x_index * y_index
+                rate_obs_daily, rate_sim_daily = rate_obs_daily[xy_index], rate_sim_daily[xy_index]
+
+                df.loc[(treat, 'day%d' % (day + 1)), 'e_obs'] = sum(rate_obs_daily)
+                df.loc[(treat, 'day%d' % (day + 1)), 'e_sim%d' % case] = sum(rate_sim_daily)
+
+    df['energy_effect'] = (df['e_sim3'] - df['e_sim0']) / df['e_sim0']
+
+    print(df)
+
+
+if __name__ == '__main__':
+    import pandas as pd
+    import numpy as np
+    from pathlib import Path
+    from datetime import datetime, timedelta
+    from matplotlib import dates, pyplot, patches, colors, colorbar, rcParams, ticker, gridspec
+
+    from hydroshoot.architecture import mtg_load
+    from hydroshoot.utilities import vapor_pressure_deficit as VPDa
+    from hydroshoot import display
+
+    rcParams.update({'font.size': 11})
+    pyplot.style.use('seaborn-ticks')
+
+    example_pth = Path(__file__).parents[2] / 'example'
+
+    plot_figure_6()
+    plot_figure_7()
+    plot_figure_8()
+    plot_figure_9()
+    plot_figure_10()
+    plot_figure_11()
+    plot_figure_12()
+    plot_figure_13()
+    plot_figure_14()
+    plot_figure_15()
+    write_table_1()
+    estimate_energy_balance_contribution()

hydroshoot/example/figs/table_1.csv

@@ -0,0 +1,11 @@
+treat,sim,an_mbe,an_rmse,e_mbe,e_rmse
+ww,sim0,-0.051909718477230367,5.4855457305525137,89.466239528284021,141.66446008942523
+ww,sim1,1.259010978690791,7.4301344033772123,225.00334182464707,361.18357575613237
+ww,sim2,1.2014013791105675,5.2713845531137755,170.15288191369015,247.4531022336021
+ww,sim3,-0.25619092657093223,5.5829424838741728,109.17880562355603,154.60276044804414
+ww,sim4,0.50396645797126982,5.3957488856738012,100.03530354699046,156.68214220659749
+ws,sim0,1.010557860506649,3.7581259962191074,27.650360054899654,67.570546085109285
+ws,sim1,3.0316942947008552,6.9534095582668689,193.89281738817965,315.4147399719518
+ws,sim2,1.363692781794585,3.8564779052298519,35.407294393274391,71.131528719311277
+ws,sim3,1.0450312828027237,3.8226861782122543,34.120378493320707,70.808298231510008
+ws,sim4,3.1606701990475412,6.0690600845156348,127.39656220185651,196.01912731881146