Quickstart of SuPy¶
This quickstart demonstrates the essential and simplest workflow of
supy
in SUEWS simulation:
More advanced use of supy
are available in the
tutorials
Before start, we need to load the following necessary packages.
In [1]:
import matplotlib.pyplot as plt
import supy as sp
import pandas as pd
import numpy as np
from pathlib import Path
get_ipython().run_line_magic('matplotlib', 'inline')
# produce high-quality figures, which can also be set as one of ['svg', 'pdf', 'retina', 'png']
# 'svg' produces high quality vector figures
%config InlineBackend.figure_format = 'svg'
Load input files¶
For existing SUEWS users:¶
First, a path to SUEWS RunControl.nml
should be specified, which
will direct supy
to locate input files.
In [2]:
path_runcontrol = Path('../sample_run') / 'RunControl.nml'
In [3]:
df_state_init = sp.init_supy(path_runcontrol)
A sample df_state_init
looks below (note that .T
is used here to
a nicer tableform view):
In [4]:
df_state_init.filter(like='method').T
Out[4]:
grid | 98 | |
---|---|---|
var | ind_dim | |
aerodynamicresistancemethod | 0 | 2 |
evapmethod | 0 | 2 |
emissionsmethod | 0 | 2 |
netradiationmethod | 0 | 3 |
roughlenheatmethod | 0 | 2 |
roughlenmommethod | 0 | 2 |
smdmethod | 0 | 0 |
stabilitymethod | 0 | 3 |
storageheatmethod | 0 | 1 |
waterusemethod | 0 | 0 |
Following the convention of SUEWS, supy
loads meteorological forcing
(met-forcing) files at the grid level.
In [5]:
grid = df_state_init.index[0]
df_forcing = sp.load_forcing_grid(path_runcontrol, grid)
For new users to SUEWS/SuPy:¶
To ease the input file preparation, a helper function
load_SampleData
is provided to get the sample input for SuPy
simulations
In [6]:
df_state_init, df_forcing = sp.load_SampleData()
Overview of SuPy input¶
df_state_init
¶
df_state_init
includes model Initial state consisting of:
- surface characteristics (e.g., albedo, emissivity, land cover fractions, etc.; full details refer to SUEWS documentation)
- model configurations (e.g., stability; full details refer to SUEWS documentation)
Detailed description of variables in df_state_init
refers to SuPy
input
Surface land cover fraction information in the sample input dataset:
In [30]:
df_state_init.loc[:,['bldgh','evetreeh','dectreeh']]
Out[30]:
var | bldgh | dectreeh | evetreeh |
---|---|---|---|
ind_dim | 0 | 0 | 0 |
grid | |||
98 | 22.0 | 13.1 | 13.1 |
In [7]:
df_state_init.filter(like='sfr')
Out[7]:
var | sfr | ||||||
---|---|---|---|---|---|---|---|
ind_dim | (0,) | (1,) | (2,) | (3,) | (4,) | (5,) | (6,) |
grid | |||||||
98 | 0.43 | 0.38 | 0.001 | 0.019 | 0.029 | 0.001 | 0.14 |
df_forcing
¶
df_forcing
includes meteorological and other external forcing
information.
Detailed description of variables in df_forcing
refers to SuPy
input.
Below is an overview of forcing variables of the sample data set used in the following simulations.
In [28]:
list_var_forcing = [
'kdown',
'Tair',
'RH',
'pres',
'U',
'rain',
]
dict_var_label = {
'kdown': 'Incoming Solar\n Radiation ($ \mathrm{W \ m^{-2}}$)',
'Tair': 'Air Temperature ($^{\circ}}$C)',
'RH': r'Relative Humidity (%)',
'pres': 'Air Pressure (hPa)',
'rain': 'Rainfall (mm)',
'U': 'Wind Speed (m $\mathrm{s^{-1}}$)'
}
df_plot_forcing_x = df_forcing.loc[:, list_var_forcing].copy().shift(
-1).dropna(how='any')
df_plot_forcing = df_plot_forcing_x.resample('1h').mean()
df_plot_forcing['rain'] = df_plot_forcing_x['rain'].resample('1h').sum()
axes = df_plot_forcing.plot(
subplots=True,
figsize=(8, 12),
legend=False,
)
fig = axes[0].figure
fig.tight_layout()
fig.autofmt_xdate(bottom=0.2, rotation=0, ha='center')
for ax, var in zip(axes, list_var_forcing):
ax.set_ylabel(dict_var_label[var])
Run simulations¶
Once met-forcing (via df_forcing
) and initial conditions (via
df_state_init
) are loaded in, we call sp.run_supy
to conduct a
SUEWS simulation, which will return two pandas
DataFrame
s:
df_output
and df_state
.
In [9]:
df_output, df_state_final = sp.run_supy(df_forcing, df_state_init)
df_output
¶
df_output
is an ensemble output collection of major SUEWS output
groups, including:
- SUEWS: the essential SUEWS output variables
- DailyState: variables of daily state information
- snow: snow output variables (effective when
snowuse = 1
set indf_state_init
)
Detailed description of variables in df_output
refers to SuPy
output
In [10]:
df_output.columns.levels[0]
Out[10]:
Index(['SUEWS', 'snow', 'DailyState'], dtype='object', name='group')
df_state_final
¶
df_state_final
is a DataFrame
for holding:
- all model states if
save_state
is set toTrue
when callingsp.run_supy
andsupy
may run significantly slower for a large simulation; - or, only the final state if
save_state
is set toFalse
(the default setting) in which modesupy
has a similar performance as the standalone compiled SUEWS executable.
Entries in df_state_final
have the same data structure as
df_state_init
and can thus be used for other SUEWS simulations
staring at the timestamp as in df_state_final
.
Detailed description of variables in df_state_final
refers to SuPy
output
In [11]:
df_state_final.T.head()
Out[11]:
grid | 98 | ||
---|---|---|---|
datetime | 2012-01-01 00:05:00 | 2013-01-01 00:05:00 | |
var | ind_dim | ||
aerodynamicresistancemethod | 0 | 2.0 | 2.0 |
ah_min | (0,) | 15.0 | 15.0 |
(1,) | 15.0 | 15.0 | |
ah_slope_cooling | (0,) | 2.7 | 2.7 |
(1,) | 2.7 | 2.7 |
Examine results¶
Thanks to the functionality inherited from pandas
and other packages
under the PyData stack, compared with the
standard SUEWS simulation workflow, supy
enables more convenient
examination of SUEWS results by statistics calculation, resampling,
plotting (and many more).
Ouptut structure¶
df_output
is organised with MultiIndex
(grid,timestamp)
and
(group,varaible)
as index
and columns
, respectively.
In [12]:
df_output.head()
Out[12]:
group | SUEWS | ... | DailyState | |||||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
var | Kdown | Kup | Ldown | Lup | Tsurf | QN | QF | QS | QH | QE | ... | DensSnow_Paved | DensSnow_Bldgs | DensSnow_EveTr | DensSnow_DecTr | DensSnow_Grass | DensSnow_BSoil | DensSnow_Water | a1 | a2 | a3 | |
grid | datetime | |||||||||||||||||||||
98 | 2012-01-01 00:05:00 | 0.153333 | 0.018279 | 344.310184 | 371.986259 | 11.775615 | -27.541021 | 40.574001 | -46.53243 | 62.420064 | 3.576493 | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
2012-01-01 00:10:00 | 0.153333 | 0.018279 | 344.310184 | 371.986259 | 11.775615 | -27.541021 | 39.724283 | -46.53243 | 61.654096 | 3.492744 | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | |
2012-01-01 00:15:00 | 0.153333 | 0.018279 | 344.310184 | 371.986259 | 11.775615 | -27.541021 | 38.874566 | -46.53243 | 60.885968 | 3.411154 | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | |
2012-01-01 00:20:00 | 0.153333 | 0.018279 | 344.310184 | 371.986259 | 11.775615 | -27.541021 | 38.024849 | -46.53243 | 60.115745 | 3.331660 | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | |
2012-01-01 00:25:00 | 0.153333 | 0.018279 | 344.310184 | 371.986259 | 11.775615 | -27.541021 | 37.175131 | -46.53243 | 59.343488 | 3.254200 | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
5 rows × 218 columns
Here we demonstrate several typical scenarios for SUEWS results examination.
The essential SUEWS
output collection is extracted as a separate
variable for easier processing in the following sections. More advanced
slicing
techniques
are available in pandas
documentation.
In [13]:
df_output_suews = df_output['SUEWS']
Statistics Calculation¶
We can use .describe()
method for a quick overview of the key
surface energy balance budgets.
In [14]:
df_output_suews.loc[:, ['QN', 'QS', 'QH', 'QE', 'QF']].describe()
Out[14]:
var | QN | QS | QH | QE | QF |
---|---|---|---|---|---|
count | 105408.000000 | 105408.000000 | 105408.000000 | 105408.000000 | 105408.000000 |
mean | 42.574626 | -2.128683 | 101.666596 | 22.880054 | 79.047033 |
std | 134.685026 | 83.616791 | 64.426005 | 28.535854 | 31.237533 |
min | -84.389073 | -81.747551 | -44.370665 | -0.649114 | 26.333882 |
25% | -41.096052 | -54.493597 | 52.976865 | 1.918672 | 50.066249 |
50% | -24.869018 | -43.957916 | 82.984095 | 13.057008 | 82.896007 |
75% | 77.405862 | 20.563903 | 140.933003 | 31.672138 | 104.858241 |
max | 689.067820 | 387.412149 | 338.167114 | 272.755143 | 160.076122 |
Plotting¶
Basic example¶
Plotting is very straightforward via the .plot
method bounded with
pandas.DataFrame
. Note the usage of loc
for to slices of the
output DataFrame
.
In [15]:
# a dict for better display variable names
dict_var_disp = {
'QN': '$Q^*$',
'QS': r'$\Delta Q_S$',
'QE': '$Q_E$',
'QH': '$Q_H$',
'QF': '$Q_F$',
'Kdown': r'$K_{\downarrow}$',
'Kup': r'$K_{\uparrow}$',
'Ldown': r'$L_{\downarrow}$',
'Lup': r'$L_{\uparrow}$',
'Rain': '$P$',
'Irr': '$I$',
'Evap': '$E$',
'RO': '$R$',
'TotCh': '$\Delta S$',
}
Quick look at the simulation results:
In [16]:
ax_output = df_output_suews\
.loc[grid]\
.loc['2012 6 1':'2012 6 7',
['QN', 'QS', 'QE', 'QH', 'QF']]\
.rename(columns=dict_var_disp)\
.plot()
ax_output.set_xlabel('Date')
ax_output.set_ylabel('Flux ($ \mathrm{W \ m^{-2}}$)')
ax_output.legend()
Out[16]:
<matplotlib.legend.Legend at 0x1263f4f98>
More examples¶
Below is a more complete example for examination of urban energy balance over the whole summer (June to August).
In [17]:
# energy balance
ax_output = df_output_suews.loc[grid]\
.loc['2012 6':'2012 8', ['QN', 'QS', 'QE', 'QH', 'QF']]\
.rename(columns=dict_var_disp)\
.plot(
figsize=(10, 3),
title='Surface Energy Balance',
)
ax_output.set_xlabel('Date')
ax_output.set_ylabel('Flux ($ \mathrm{W \ m^{-2}}$)')
ax_output.legend()
Out[17]:
<matplotlib.legend.Legend at 0x12642e8d0>
Resampling¶
The suggested runtime/simulation frequency of SUEWS is 300 s
, which
usually results a large output and may be over-weighted for storage and
analysis. Also, you may feel apparent slowdown in producing the above
figure as a large amount of data were used for the plotting. To slim
down the result size for analysis and output, we can resample
the
default output very easily.
In [18]:
rsmp_1d = df_output_suews.loc[grid].resample('1d')
# daily mean values
df_1d_mean = rsmp_1d.mean()
# daily sum values
df_1d_sum = rsmp_1d.sum()
We can then re-examine the above energy balance at hourly scale and plotting will be significantly faster.
In [19]:
# energy balance
ax_output = df_1d_mean\
.loc[:, ['QN', 'QS', 'QE', 'QH', 'QF']]\
.rename(columns=dict_var_disp)\
.plot(
figsize=(10, 3),
title='Surface Energy Balance',
)
ax_output.set_xlabel('Date')
ax_output.set_ylabel('Flux ($ \mathrm{W \ m^{-2}}$)')
ax_output.legend()
Out[19]:
<matplotlib.legend.Legend at 0x126dccb70>
Then we use the hourly results for other analyses.
In [20]:
# radiation balance
ax_output = df_1d_mean\
.loc[:, ['QN', 'Kdown', 'Kup', 'Ldown', 'Lup']]\
.rename(columns=dict_var_disp)\
.plot(
figsize=(10, 3),
title='Radiation Balance',
)
ax_output.set_xlabel('Date')
ax_output.set_ylabel('Flux ($ \mathrm{W \ m^{-2}}$)')
ax_output.legend()
Out[20]:
<matplotlib.legend.Legend at 0x1272d6358>
In [21]:
# water balance
ax_output = df_1d_sum\
.loc[:, ['Rain', 'Irr', 'Evap', 'RO', 'TotCh']]\
.rename(columns=dict_var_disp)\
.plot(
figsize=(10, 3),
title='Surface Water Balance',
)
ax_output.set_xlabel('Date')
ax_output.set_ylabel('Water amount (mm)')
ax_output.legend()
Out[21]:
<matplotlib.legend.Legend at 0x127610668>
Get an overview of partitioning in energy and water balance at monthly scales:
In [22]:
# get a monthly Resampler
df_plot=df_output_suews.loc[grid].copy()
df_plot.index=df_plot.index.set_names('Month')
rsmp_1M = df_plot\
.shift(-1)\
.dropna(how='all')\
.resample('1M', kind='period')
# mean values
df_1M_mean = rsmp_1M.mean()
# sum values
df_1M_sum = rsmp_1M.sum()
In [23]:
# month names
name_mon = [x.strftime('%b') for x in rsmp_1M.groups]
# create subplots showing two panels together
fig, axes = plt.subplots(2, 1, sharex=True)
# surface energy balance
df_1M_mean\
.loc[:, ['QN', 'QS', 'QE', 'QH', 'QF']]\
.rename(columns=dict_var_disp)\
.plot(
ax=axes[0], # specify the axis for plotting
figsize=(10, 6), # specify figure size
title='Surface Energy Balance',
kind='bar',
)
# surface water balance
df_1M_sum\
.loc[:, ['Rain', 'Irr', 'Evap', 'RO', 'TotCh']]\
.rename(columns=dict_var_disp)\
.plot(
ax=axes[1], # specify the axis for plotting
title='Surface Water Balance',
kind='bar'
)
# annotations
axes[0].set_ylabel('Mean Flux ($ \mathrm{W \ m^{-2}}$)')
axes[0].legend()
axes[1].set_xlabel('Month')
axes[1].set_ylabel('Total Water Amount (mm)')
axes[1].xaxis.set_ticklabels(name_mon, rotation=0)
axes[1].legend()
Out[23]:
<matplotlib.legend.Legend at 0x127add320>
Output¶
The resampled output can be outputed for a smaller file.
In [24]:
df_1d_mean.to_csv(
'suews_1d_mean.txt',
sep='\t',
float_format='%8.2f',
na_rep=-999,
)
For a justified format, we use the to_string
for better format
controlling and write the formatted string out to a file.
In [25]:
str_out = df_1d_mean.to_string(
float_format='%8.2f',
na_rep='-999',
justify='right',
)
with open('suews_sample.txt', 'w') as file_out:
print(str_out, file=file_out)