Skip to content

Commit 1c57bd9

Browse files
JohnVolkJohnVolk
committed
* add check in optimizer constructor for original model's
 statvar.dat file, if not found run original model in location of other input files for later plot comparisons * improve optimizer before and after plots * add monthly aggregate plots, ability to show simulation swrad results versus original model and/or measured swrad
1 parent 21b5c43 commit 1c57bd9

1 file changed

Lines changed: 86 additions & 26 deletions

File tree

‎prms_python/optimizer.py‎

Lines changed: 86 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@
1515
from .data import Data
1616
from .parameters import Parameters
1717
from .simulation import Simulation, SimulationSeries
18-
18+
from .util import load_statvar
1919

2020
OPJ = os.path.join
2121

@@ -58,9 +58,18 @@ def __init__(self, parameters, data, control_file, working_dir,
5858
else:
5959
raise TypeError('data must be instance of Data')
6060

61+
input_dir = '{}'.format(os.sep).join(control_file.split(os.sep)[:-1])
62+
if not os.path.isfile(OPJ(input_dir, 'statvar.dat')):
63+
print('You have no statvar.dat file in your current model directory')
64+
print('Running PRMS on original data in {} for later comparison'\
65+
.format(input_dir))
66+
sim = Simulation(input_dir)
67+
sim.run()
68+
6169
if not os.path.isdir(working_dir):
6270
os.mkdir(working_dir)
63-
71+
72+
self.input_dir = input_dir
6473
self.control_file = control_file
6574
self.working_dir = working_dir
6675
self.title = title
@@ -142,7 +151,7 @@ def _error(x, y):
142151
'start_time' : str(srad_start_time),
143152
'end_time' : str(srad_end_time),
144153
'measured_swrad' : reference_srad_path,
145-
'sim_dirs' : [],
154+
'sim_dirs' : [],
146155
'original_params' : self.parameters.base_file,
147156
'n_sims' : n_sims
148157
}
@@ -159,7 +168,8 @@ def _error(x, y):
159168

160169
print('{0}\nOutput information sent to {1}\n'.format('-' * 80, json_outfile))
161170

162-
def plot_srad_optimization(self, freq='daily', method='time_series'):
171+
def plot_srad_optimization(self, freq='daily', method='time_series',\
172+
plot_vars='both', return_fig=False):
163173
"""
164174
Basic plotting of current srad optimization results with
165175
limited options for quick viewing, measured, original,
@@ -171,13 +181,17 @@ def plot_srad_optimization(self, freq='daily', method='time_series'):
171181
172182
Kwargs:
173183
freq (str): frequency of time series plots, value can be 'daily'
174-
or 'monthly' for solar radiation TODO:need to finish monthly!!!
184+
or 'monthly' for solar radiation !!!need to finish monthly!!!
175185
method (str): 'time_series' for time series sub plot of each
176186
simulation alongside measured radiation. Other choice is
177187
'correlation' which plots each measured daily solar radiation
178188
value versus the corresponding simulated variable as subplots
179189
one for each simulation in the optimization. With coefficients
180190
of determiniationi i.e. square of pearson correlation coef.
191+
plot_vars (str): what to plot alongside simulated srad:
192+
'meas': plot simulated along with measured swrad
193+
'orig': plot simulated along with the original simulated swrad
194+
'both': plot simulated, with original simulation and measured
181195
"""
182196
if not self.srad_outputs:
183197
raise ValueError('You have not run any srad optimizations')
@@ -187,32 +201,75 @@ def plot_srad_optimization(self, freq='daily', method='time_series'):
187201
idx = X.index.intersection(self.srad_outputs[0]['statvar']\
188202
['swrad_{}'.format(self.srad_hru)].index)
189203
X = X[idx]
204+
orig = load_statvar(OPJ(self.input_dir, 'statvar.dat'))['swrad_{}'\
205+
.format(self.srad_hru)][idx]
206+
meas = self.measured_srad[idx]
207+
# styles for each plot
208+
ms = 4 # markersize for all points
209+
orig_sty = dict(linestyle='none',markersize=ms,\
210+
markerfacecolor='none', marker='s',\
211+
markeredgecolor='royalblue', color='royalblue')
212+
meas_sty = dict(linestyle='none',markersize=ms+1,\
213+
markerfacecolor='none', marker='1',\
214+
markeredgecolor='k', color='k')
215+
sim_sty = dict(linestyle='none',markersize=ms,\
216+
markerfacecolor='none', marker='o',\
217+
markeredgecolor='r', color='r')
190218
n = len(self.srad_outputs) # number of simulations to plot
191-
219+
## number of subplots and rows (two plots per row)
220+
nrow = n//2 # round down if odd n
221+
ncol = 2
222+
odd_n = False
223+
if n/2. - nrow == 0.5:
224+
nrow+=1 # odd number need extra row
225+
odd_n = True
226+
########
227+
## Start plots depnding on key word arguments
228+
########
192229
if freq == 'daily' and method == 'time_series':
193230
fig, ax = plt.subplots(n, sharex=True, sharey=True,\
194231
figsize=(12,n*3.5))
195232
axs = ax.ravel()
196233
for i,out in enumerate(self.srad_outputs):
234+
if plot_vars in ('meas', 'both'):
235+
axs[i].plot(meas, label='Measured', **meas_sty)
236+
if plot_vars in ('orig', 'both'):
237+
axs[i].plot(orig, label='Original sim.', **orig_sty)
197238
axs[i].plot(out['statvar']['swrad_{}'.format(self.srad_hru)]\
198-
[idx], 'r.', markersize=3, label='Simulated')
199-
axs[i].plot(self.measured_srad[idx], 'k.', markersize=3,\
200-
label='Measured')
239+
[idx], **sim_sty)
201240
axs[i].set_ylabel('sim: {}'.format(out['simulation_dir'].\
202-
split(os.sep)[-1].replace('_', ' ')))
241+
split(os.sep)[-1].replace('_', ' ')), fontsize=10)
203242
if i == 0: axs[i].legend(markerscale=5, loc='best')
204243
fig.subplots_adjust(hspace=0)
205244
fig.autofmt_xdate()
206-
plt.show()
245+
246+
elif freq == 'monthly' and method == 'time_series':
247+
# compute monthly means
248+
meas = meas.groupby(meas.index.month).mean()
249+
orig = orig.groupby(orig.index.month).mean()
250+
# change line styles for monthly plots to lines not points
251+
for d in (orig_sty, meas_sty, sim_sty):
252+
d['linestyle'] = '-'
253+
d['marker'] = None
254+
255+
fig, ax = plt.subplots(nrows=nrow, ncols=ncol, figsize=(12,n*3.5))
256+
axs = ax.ravel()
257+
for i,out in enumerate(self.srad_outputs):
258+
if plot_vars in ('meas', 'both'):
259+
axs[i].plot(meas, label='Measured', **meas_sty)
260+
if plot_vars in ('orig', 'both'):
261+
axs[i].plot(orig, label='Original sim.', **orig_sty)
262+
sim = out['statvar']['swrad_{}'.format(self.srad_hru)][idx]
263+
sim = sim.groupby(sim.index.month).mean()
264+
axs[i].plot(sim, **sim_sty)
265+
axs[i].set_ylabel('sim: {}\nmean swrad'.format(out['simulation_dir'].\
266+
split(os.sep)[-1].replace('_', ' ')), fontsize=10)
267+
axs[i].set_xlim(0.5,12.5)
268+
if i == 0: axs[i].legend(markerscale=5, loc='best')
269+
if odd_n: # empty subplot if odd number of simulations
270+
fig.delaxes(axs[n])
207271

208272
elif method == 'correlation':
209-
## number of subplots and rows (two plots per row)
210-
nrow = n//2 # round down if odd n
211-
ncol = 2
212-
odd_n = False
213-
if n/2. - nrow == 0.5:
214-
nrow+=1 # odd number need extra row
215-
odd_n = True
216273
## figure
217274
fig, ax = plt.subplots(nrows=nrow, ncols=ncol, figsize=(12,n*3))
218275
axs = ax.ravel()
@@ -228,7 +285,7 @@ def plot_srad_optimization(self, freq='daily', method='time_series'):
228285
axs[i].plot([0, m], [0, m], 'k--', lw=2) ## one to one line
229286
axs[i].set_xlim(meas_min,meas_max)
230287
axs[i].set_ylim(sim_min, sim_max)
231-
axs[i].scatter(X, Y, facecolors='none', edgecolor='r', s=3)
288+
axs[i].plot(X, Y, **sim_sty)
232289
axs[i].set_ylabel('sim: {}'.format(out['simulation_dir']\
233290
.split(os.sep)[-1].replace('_', ' ')))
234291
axs[i].set_xlabel('Measured shortwave radiation')
@@ -237,6 +294,9 @@ def plot_srad_optimization(self, freq='daily', method='time_series'):
237294
ha='left', va='center', transform=axs[i].transAxes)
238295
if odd_n: # empty subplot if odd number of simulations
239296
fig.delaxes(axs[n])
297+
298+
if return_fig:
299+
return fig
240300

241301
def _create_metafile_name(out_dir, opt_title, stage):
242302
"""
@@ -252,14 +312,14 @@ def _create_metafile_name(out_dir, opt_title, stage):
252312
location where simulation series outputs and optimization
253313
json files are located, aka Optimizer.working_dir
254314
opt_title (str): optimization instance title for file search
255-
stage (str): stage of optimization, e.g. 'swrad', 'pet'
315+
stage (str): stage of optimization, e.g. 'swrad', 'pet'
256316
257317
Returns:
258-
name (str): file name for the current optimization simulation series
259-
metadata json file. E.g 'dry_creek_swrad_opt.json', or if
260-
this is the second time you have run an optimization titled
261-
'dry_creek' the next json file will be returned as
262-
'dry_creek_swrad_opt1.json' and so on with integer increments
318+
name (str): file name for the current optimization simulation series
319+
metadata json file. E.g 'dry_creek_swrad_opt.json', or if
320+
this is the second time you have run an optimization titled
321+
'dry_creek' the next json file will be returned as
322+
'dry_creek_swrad_opt1.json' and so on with integer increments
263323
"""
264324
swrad_meta_re = re.compile(r'^{}_{}_opt(\d*)\.json'.format(opt_title, stage))
265325
reps = []
@@ -315,7 +375,7 @@ def _resample_param(param, p_min, p_max, noise_factor=0.1 ):
315375
return tmp
316376

317377
def _mod_params(parameters, intcp, slope):
318-
# deepcopy on parameters was rasining:
378+
# deepcopy was crashing, raising:
319379
# TypeError: cannot serialize '_io.TextIOWrapper' object
320380
ret = copy(parameters)
321381
#print (intcp, slope)

0 commit comments

Comments
 (0)