Skip to content

Commit ecc1bb9

Browse files
btobersbtobers
committed
handle failed spinup
1 parent 927d843 commit ecc1bb9

2 files changed

Lines changed: 94 additions & 83 deletions

File tree

‎pygem/bin/run/run_inversion.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -406,7 +406,7 @@ def main():
406406
type=str,
407407
default=None,
408408
help='Filepath containing list of rgi_glac_number, helpful for running batches on spc',
409-
),
409+
)
410410
parser.add_argument(
411411
'-calibrate_regional_glen_a',
412412
type=str2bool,

‎pygem/bin/run/run_spinup.py‎

Lines changed: 93 additions & 82 deletions
Original file line numberDiff line numberDiff line change
@@ -195,7 +195,9 @@ def run(glacno_list, mb_model_params, optimize=False, periods2try=[20], outdir=N
195195
gd = single_flowline_glacier_directory_with_calving(glacier_str, reset=False)
196196
gd.is_tidewater = True
197197
if kwargs['allow_calving']:
198-
kwargs['evolution_model'] = FluxBasedModel # use FluxBasedModel to allow for calving
198+
kwargs['evolution_model'] = partial(
199+
FluxBasedModel, water_level=gd.get_diagnostics().get('calving_water_level', None)
200+
) # use FluxBasedModel to allow for calving
199201
cfg.PARAMS['use_kcalving_for_inversion'] = True
200202
cfg.PARAMS['use_kcalving_for_run'] = True
201203
# Load quality controlled frontal ablation data
@@ -369,96 +371,105 @@ def _objective(**kwargs):
369371
print('All results:', {k: v[0] for k, v in results.items()})
370372
print(f'Best spinup_period = {best_period}, mismatch = {best_value}')
371373

374+
if all(v[1] is None for v in results.values()):
375+
raise ValueError('Spinup failed for all tested periods')
376+
372377
# find worst - ignore failed runs
373-
worst_period = max((k for k in results if results[k][0] != float('inf')), key=lambda k: results[k][0])
378+
worst_period = max(
379+
(k for k in results if results[k][0] != float('inf')),
380+
key=lambda k: results[k][0],
381+
default=best_period,
382+
)
374383
worst_value, worst_model = results[worst_period]
375384

376385
############################
377386
### diagnostics plotting ###
378387
############################
379-
# binned area
380-
ax.legend()
381-
ax.set_title(gd.rgi_id)
382-
ax.set_xlabel('distance along flowline (m)')
383-
ax.set_ylabel('surface area (m$^2$)')
384-
ax.set_xlim([0, ax.get_xlim()[1]])
385-
ax.set_ylim([0, ax.get_ylim()[1]])
386-
fig.tight_layout()
387-
if debug and ncores == 1:
388-
plt.show()
389-
if outdir:
390-
fig.savefig(f'{outdir}/{glac_no}-spinup_binned_area.png', dpi=300)
388+
if best_model is not None and worst_model is not None:
389+
# binned area
390+
ax.legend()
391+
ax.set_title(gd.rgi_id)
392+
ax.set_xlabel('distance along flowline (m)')
393+
ax.set_ylabel('surface area (m$^2$)')
394+
ax.set_xlim([0, ax.get_xlim()[1]])
395+
ax.set_ylim([0, ax.get_ylim()[1]])
396+
fig.tight_layout()
397+
if debug and ncores == 1:
398+
plt.show()
399+
if outdir:
400+
fig.savefig(f'{outdir}/{glac_no}-spinup_binned_area.png', dpi=300)
391401
plt.close()
392402

393-
# 1d elevation change
394-
labels = [
395-
(f'{start[:-2].replace("-", "")}:{end[:-3].replace("-", "")}')
396-
for start, end in gd.elev_change_1d['dates']
397-
]
398-
fig, ax = plt.subplots(figsize=(8, 5))
399-
400-
for t in range(gd.elev_change_1d['dhdt'].shape[1]):
401-
# plot Obs first, grab the color
402-
(line,) = ax.plot(
403-
gd.elev_change_1d['bin_centers'],
404-
gd.elev_change_1d['dhdt'][:, t],
405-
linestyle='-',
406-
marker='.',
407-
label=labels[t],
408-
)
409-
color = line.get_color()
410-
411-
# plot Best model with same color
412-
ax.plot(
413-
gd.elev_change_1d['bin_centers'],
414-
best_model[:, t],
415-
linestyle='--',
416-
marker='.',
417-
color=color,
418-
)
403+
if best_model is not None and worst_model is not None:
404+
# 1d elevation change
405+
labels = [
406+
(f'{start[:-2].replace("-", "")}:{end[:-3].replace("-", "")}')
407+
for start, end in gd.elev_change_1d['dates']
408+
]
409+
fig, ax = plt.subplots(figsize=(8, 5))
410+
411+
for t in range(gd.elev_change_1d['dhdt'].shape[1]):
412+
# plot Obs first, grab the color
413+
(line,) = ax.plot(
414+
gd.elev_change_1d['bin_centers'],
415+
gd.elev_change_1d['dhdt'][:, t],
416+
linestyle='-',
417+
marker='.',
418+
label=labels[t],
419+
)
420+
color = line.get_color()
421+
422+
# plot Best model with same color
423+
ax.plot(
424+
gd.elev_change_1d['bin_centers'],
425+
best_model[:, t],
426+
linestyle='--',
427+
marker='.',
428+
color=color,
429+
)
419430

420-
# plot Worst model with same color
421-
ax.plot(
422-
gd.elev_change_1d['bin_centers'],
423-
worst_model[:, t],
424-
linestyle=':',
425-
marker='.',
426-
color=color,
427-
)
428-
ax.axvline(gd.ela, c='grey', ls=':')
429-
ax.axhline(0, c='grey', ls='-')
430-
ax.plot([], [], 'k--', label=r'$\hat{best}$')
431-
ax.plot([], [], 'k:', label=r'$\hat{worst}$')
432-
ax.set_xlabel('elevation (m)')
433-
ax.set_ylabel(r'elevation change (m yr$^{-1}$)')
434-
ax.set_title(
435-
f'{glac_no}\nBest={best_period} (mismatch={best_value:.3f}), '
436-
f'Worst={worst_period} (mismatch={worst_value:.3f})'
437-
)
438-
ax.legend(handlelength=1, borderaxespad=0, fancybox=False)
439-
# plot area
440-
if 'bin_area' in gd.elev_change_1d:
441-
area = np.array(gd.elev_change_1d['bin_area'])
442-
area_mask = area > 0
443-
ax2 = ax.twinx() # shares x-axis
444-
ax2.fill_between(
445-
np.array(gd.elev_change_1d['bin_centers'])[area_mask],
446-
0,
447-
area[area_mask],
448-
color='gray',
449-
alpha=0.1,
431+
# plot Worst model with same color
432+
ax.plot(
433+
gd.elev_change_1d['bin_centers'],
434+
worst_model[:, t],
435+
linestyle=':',
436+
marker='.',
437+
color=color,
438+
)
439+
ax.axvline(gd.ela, c='grey', ls=':')
440+
ax.axhline(0, c='grey', ls='-')
441+
ax.plot([], [], 'k--', label=r'$\hat{best}$')
442+
ax.plot([], [], 'k:', label=r'$\hat{worst}$')
443+
ax.set_xlabel('elevation (m)')
444+
ax.set_ylabel(r'elevation change (m yr$^{-1}$)')
445+
ax.set_title(
446+
f'{glac_no}\nBest={best_period} (mismatch={best_value:.3f}), '
447+
f'Worst={worst_period} (mismatch={worst_value:.3f})'
450448
)
451-
ax2.set_ylim([0, ax2.get_ylim()[1]])
452-
ax2.set_ylabel(r'area (m $^{2}$)', color='gray')
453-
ax2.tick_params(axis='y', colors='gray')
454-
ax2.spines['right'].set_color('gray')
455-
ax2.yaxis.label.set_color('gray')
456-
fig.tight_layout()
457-
if debug and ncores == 1:
458-
plt.show()
459-
if outdir:
460-
fig.savefig(f'{outdir}/{glac_no}-spinup_optimization.png', dpi=300)
461-
plt.close()
449+
ax.legend(handlelength=1, borderaxespad=0, fancybox=False)
450+
# plot area
451+
if 'bin_area' in gd.elev_change_1d:
452+
area = np.array(gd.elev_change_1d['bin_area'])
453+
area_mask = area > 0
454+
ax2 = ax.twinx() # shares x-axis
455+
ax2.fill_between(
456+
np.array(gd.elev_change_1d['bin_centers'])[area_mask],
457+
0,
458+
area[area_mask],
459+
color='gray',
460+
alpha=0.1,
461+
)
462+
ax2.set_ylim([0, ax2.get_ylim()[1]])
463+
ax2.set_ylabel(r'area (m $^{2}$)', color='gray')
464+
ax2.tick_params(axis='y', colors='gray')
465+
ax2.spines['right'].set_color('gray')
466+
ax2.yaxis.label.set_color('gray')
467+
fig.tight_layout()
468+
if debug and ncores == 1:
469+
plt.show()
470+
if outdir:
471+
fig.savefig(f'{outdir}/{glac_no}-spinup_optimization.png', dpi=300)
472+
plt.close()
462473
############################
463474

464475
# update spinup_period if optimized or specified as CLI argument, else remove kwarg and use OGGM default
@@ -502,7 +513,7 @@ def main():
502513
type=str,
503514
default=None,
504515
help='Filepath containing list of rgi_glac_number, helpful for running batches on spc',
505-
),
516+
)
506517
parser.add_argument('-target_yr', type=int, default=None)
507518
parser.add_argument('-ye', type=int, default=None)
508519
parser.add_argument(

0 commit comments

Comments
 (0)