Skip to content

Commit 5f3bacf

Browse files
albinwwellsalbinwwells
authored
Issue 128 -- importing melt extent and snowline data into OGGM gdirs (#139)
Closes #128
1 parent 5949dd4 commit 5f3bacf

File tree

14 files changed

+1085
-42
lines changed

14 files changed

+1085
-42
lines changed

‎docs/_static/example_meltextent_and_snowline_1d/sample_melt_extent_elev.csv‎

Lines changed: 336 additions & 0 deletions
Large diffs are not rendered by default.

‎docs/_static/example_meltextent_and_snowline_1d/sample_snowline_elev.csv‎

Lines changed: 336 additions & 0 deletions
Large diffs are not rendered by default.

‎pygem/bin/run/run_calibration.py‎

Lines changed: 9 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -565,25 +565,30 @@ def run(list_packed_vars):
565565
gcm = class_climate.GCM(name=ref_climate_name)
566566
# Air temperature [degC]
567567
gcm_temp, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
568-
gcm.temp_fn, gcm.temp_vn, main_glac_rgi, dates_table
568+
gcm.temp_fn, gcm.temp_vn, main_glac_rgi, dates_table, verbose=debug
569569
)
570570
if pygem_prms['mb']['option_ablation'] == 2 and ref_climate_name in ['ERA5']:
571571
gcm_tempstd, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
572-
gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi, dates_table
572+
gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi, dates_table, verbose=debug
573573
)
574574
else:
575575
gcm_tempstd = np.zeros(gcm_temp.shape)
576576
# Precipitation [m]
577577
gcm_prec, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
578-
gcm.prec_fn, gcm.prec_vn, main_glac_rgi, dates_table
578+
gcm.prec_fn, gcm.prec_vn, main_glac_rgi, dates_table, verbose=debug
579579
)
580580
# Elevation [m asl]
581581
gcm_elev = gcm.importGCMfxnearestneighbor_xarray(
582582
gcm.elev_fn, gcm.elev_vn, main_glac_rgi
583583
)
584584
# Lapse rate [degC m-1] (always monthly)
585585
gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
586-
gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table, upscale_var_timestep=True
586+
gcm.lr_fn,
587+
gcm.lr_vn,
588+
main_glac_rgi,
589+
dates_table,
590+
upscale_var_timestep=True,
591+
verbose=debug,
587592
)
588593

589594
# ===== LOOP THROUGH GLACIERS TO RUN CALIBRATION =====

‎pygem/bin/run/run_calibration_frontalablation.py‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -136,25 +136,25 @@ def reg_calving_flux(
136136
gcm = class_climate.GCM(name=args.ref_climate_name)
137137
# Air temperature [degC]
138138
gcm_temp, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
139-
gcm.temp_fn, gcm.temp_vn, main_glac_rgi, dates_table
139+
gcm.temp_fn, gcm.temp_vn, main_glac_rgi, dates_table, verbose=debug
140140
)
141141
if pygem_prms['mb']['option_ablation'] == 2 and args.ref_climate_name in ['ERA5']:
142142
gcm_tempstd, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
143-
gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi, dates_table
143+
gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi, dates_table, verbose=debug
144144
)
145145
else:
146146
gcm_tempstd = np.zeros(gcm_temp.shape)
147147
# Precipitation [m]
148148
gcm_prec, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
149-
gcm.prec_fn, gcm.prec_vn, main_glac_rgi, dates_table
149+
gcm.prec_fn, gcm.prec_vn, main_glac_rgi, dates_table, verbose=debug
150150
)
151151
# Elevation [m asl]
152152
gcm_elev = gcm.importGCMfxnearestneighbor_xarray(
153153
gcm.elev_fn, gcm.elev_vn, main_glac_rgi
154154
)
155155
# Lapse rate [degC m-1]
156156
gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
157-
gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table
157+
gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table, verbose=debug
158158
)
159159

160160
# ===== CALIBRATE ALL THE GLACIERS AT ONCE =====

‎pygem/bin/run/run_calibration_reg_glena.py‎

Lines changed: 8 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -341,25 +341,29 @@ def main():
341341
gcm = class_climate.GCM(name=sim_climate_name)
342342
# Air temperature [degC]
343343
gcm_temp, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
344-
gcm.temp_fn, gcm.temp_vn, main_glac_rgi_subset, dates_table
344+
gcm.temp_fn, gcm.temp_vn, main_glac_rgi_subset, dates_table, verbose=debug
345345
)
346346
if pygem_prms['mbmod']['option_ablation'] == 2 and sim_climate_name in ['ERA5']:
347347
gcm_tempstd, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
348-
gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi_subset, dates_table
348+
gcm.tempstd_fn,
349+
gcm.tempstd_vn,
350+
main_glac_rgi_subset,
351+
dates_table,
352+
verbose=debug,
349353
)
350354
else:
351355
gcm_tempstd = np.zeros(gcm_temp.shape)
352356
# Precipitation [m]
353357
gcm_prec, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
354-
gcm.prec_fn, gcm.prec_vn, main_glac_rgi_subset, dates_table
358+
gcm.prec_fn, gcm.prec_vn, main_glac_rgi_subset, dates_table, verbose=debug
355359
)
356360
# Elevation [m asl]
357361
gcm_elev = gcm.importGCMfxnearestneighbor_xarray(
358362
gcm.elev_fn, gcm.elev_vn, main_glac_rgi_subset
359363
)
360364
# Lapse rate [degC m-1]
361365
gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
362-
gcm.lr_fn, gcm.lr_vn, main_glac_rgi_subset, dates_table
366+
gcm.lr_fn, gcm.lr_vn, main_glac_rgi_subset, dates_table, verbose=debug
363367
)
364368

365369
# ===== RUN MASS BALANCE =====

‎pygem/bin/run/run_inversion.py‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -70,19 +70,19 @@ def run(
7070

7171
# Air temperature [degC]
7272
temp, _ = ref_clim.importGCMvarnearestneighbor_xarray(
73-
ref_clim.temp_fn, ref_clim.temp_vn, main_glac_rgi, dt
73+
ref_clim.temp_fn, ref_clim.temp_vn, main_glac_rgi, dt, verbose=debug
7474
)
7575
# Precipitation [m]
7676
prec, _ = ref_clim.importGCMvarnearestneighbor_xarray(
77-
ref_clim.prec_fn, ref_clim.prec_vn, main_glac_rgi, dt
77+
ref_clim.prec_fn, ref_clim.prec_vn, main_glac_rgi, dt, verbose=debug
7878
)
7979
# Elevation [m asl]
8080
elev = ref_clim.importGCMfxnearestneighbor_xarray(
8181
ref_clim.elev_fn, ref_clim.elev_vn, main_glac_rgi
8282
)
8383
# Lapse rate [degC m-1]
8484
lr, _ = ref_clim.importGCMvarnearestneighbor_xarray(
85-
ref_clim.lr_fn, ref_clim.lr_vn, main_glac_rgi, dt
85+
ref_clim.lr_fn, ref_clim.lr_vn, main_glac_rgi, dt, verbose=debug
8686
)
8787

8888
# load prior regionally averaged modelprms (from Rounce et al. 2023)

‎pygem/bin/run/run_simulation.py‎

Lines changed: 20 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -430,17 +430,17 @@ def run(list_packed_vars):
430430
# ----- Select Temperature and Precipitation Data -----
431431
# Air temperature [degC]
432432
gcm_temp, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
433-
gcm.temp_fn, gcm.temp_vn, main_glac_rgi, dates_table_full
433+
gcm.temp_fn, gcm.temp_vn, main_glac_rgi, dates_table_full, verbose=debug
434434
)
435435
ref_temp, ref_dates = ref_gcm.importGCMvarnearestneighbor_xarray(
436-
ref_gcm.temp_fn, ref_gcm.temp_vn, main_glac_rgi, dates_table_ref
436+
ref_gcm.temp_fn, ref_gcm.temp_vn, main_glac_rgi, dates_table_ref, verbose=debug
437437
)
438438
# Precipitation [m]
439439
gcm_prec, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
440-
gcm.prec_fn, gcm.prec_vn, main_glac_rgi, dates_table_full
440+
gcm.prec_fn, gcm.prec_vn, main_glac_rgi, dates_table_full, verbose=debug
441441
)
442442
ref_prec, ref_dates = ref_gcm.importGCMvarnearestneighbor_xarray(
443-
ref_gcm.prec_fn, ref_gcm.prec_vn, main_glac_rgi, dates_table_ref
443+
ref_gcm.prec_fn, ref_gcm.prec_vn, main_glac_rgi, dates_table_ref, verbose=debug
444444
)
445445
# Elevation [m asl]
446446
try:
@@ -548,13 +548,17 @@ def run(list_packed_vars):
548548
ref_tempstd = np.zeros((main_glac_rgi.shape[0], dates_table_ref.shape[0]))
549549
elif pygem_prms['mb']['option_ablation'] == 2 and sim_climate_name in ['ERA5']:
550550
gcm_tempstd, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
551-
gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi, dates_table
551+
gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi, dates_table, verbose=debug
552552
)
553553
ref_tempstd = gcm_tempstd
554554
elif pygem_prms['mb']['option_ablation'] == 2 and args.ref_climate_name in ['ERA5']:
555555
# Compute temp std based on reference climate data
556556
ref_tempstd, ref_dates = ref_gcm.importGCMvarnearestneighbor_xarray(
557-
ref_gcm.tempstd_fn, ref_gcm.tempstd_vn, main_glac_rgi, dates_table_ref
557+
ref_gcm.tempstd_fn,
558+
ref_gcm.tempstd_vn,
559+
main_glac_rgi,
560+
dates_table_ref,
561+
verbose=debug,
558562
)
559563
# Monthly average from reference climate data
560564
gcm_tempstd = gcmbiasadj.monthly_avg_array_rolled(
@@ -567,13 +571,18 @@ def run(list_packed_vars):
567571
# Lapse rate
568572
if sim_climate_name in ['ERA-Interim', 'ERA5']:
569573
gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
570-
gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table
574+
gcm.lr_fn,
575+
gcm.lr_vn,
576+
main_glac_rgi,
577+
dates_table,
578+
upscale_var_timestep=True,
579+
verbose=debug,
571580
)
572581
ref_lr = gcm_lr
573582
else:
574583
# Compute lapse rates based on reference climate data
575584
ref_lr, ref_dates = ref_gcm.importGCMvarnearestneighbor_xarray(
576-
ref_gcm.lr_fn, ref_gcm.lr_vn, main_glac_rgi, dates_table_ref
585+
ref_gcm.lr_fn, ref_gcm.lr_vn, main_glac_rgi, dates_table_ref, verbose=debug
577586
)
578587
# Monthly average from reference climate data
579588
gcm_lr = gcmbiasadj.monthly_avg_array_rolled(
@@ -591,12 +600,9 @@ def run(list_packed_vars):
591600
else:
592601
nsims = 1
593602

594-
# Number of years (for OGGM's run_until_and_store)
595-
if pygem_prms['time']['timestep'] == 'monthly':
596-
nyears = int(dates_table.shape[0] / 12)
597-
nyears_ref = int(dates_table_ref.shape[0] / 12)
598-
else:
599-
assert True == False, 'Adjust nyears for non-monthly timestep'
603+
# Number of years
604+
nyears = dates_table.year.unique()[-1] - dates_table.year.unique()[0] + 1
605+
nyears_ref = dates_table_ref.year.unique()[-1] - dates_table.year.unique()[0] + 1
600606

601607
for glac in range(main_glac_rgi.shape[0]):
602608
if glac == 0:

‎pygem/bin/run/run_spinup.py‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -39,19 +39,19 @@ def run(glacno_list, mb_model_params, debug=False, **kwargs):
3939

4040
# Air temperature [degC]
4141
temp, _ = ref_clim.importGCMvarnearestneighbor_xarray(
42-
ref_clim.temp_fn, ref_clim.temp_vn, main_glac_rgi, dt
42+
ref_clim.temp_fn, ref_clim.temp_vn, main_glac_rgi, dt, verbose=debug
4343
)
4444
# Precipitation [m]
4545
prec, _ = ref_clim.importGCMvarnearestneighbor_xarray(
46-
ref_clim.prec_fn, ref_clim.prec_vn, main_glac_rgi, dt
46+
ref_clim.prec_fn, ref_clim.prec_vn, main_glac_rgi, dt, verbose=debug
4747
)
4848
# Elevation [m asl]
4949
elev = ref_clim.importGCMfxnearestneighbor_xarray(
5050
ref_clim.elev_fn, ref_clim.elev_vn, main_glac_rgi
5151
)
5252
# Lapse rate [degC m-1]
5353
lr, _ = ref_clim.importGCMvarnearestneighbor_xarray(
54-
ref_clim.lr_fn, ref_clim.lr_vn, main_glac_rgi, dt
54+
ref_clim.lr_fn, ref_clim.lr_vn, main_glac_rgi, dt, verbose=debug
5555
)
5656

5757
# load prior regionally averaged modelprms (from Rounce et al. 2023)

‎pygem/class_climate.py‎

Lines changed: 39 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@ class of climate data and functions associated with manipulating the dataset to
99
"""
1010

1111
import os
12+
import warnings
1213

1314
import numpy as np
1415

@@ -192,10 +193,7 @@ def __init__(self, name=str(), sim_climate_scenario=str(), realization=None):
192193
self.elev_fn = pygem_prms['climate']['paths']['era5_elev_fn']
193194
self.lr_fn = pygem_prms['climate']['paths']['era5_lr_fn']
194195
# Variable filepaths
195-
if pygem_prms['climate']['paths']['era5_fullpath']:
196-
self.var_fp = ''
197-
self.fx_fp = ''
198-
else:
196+
if pygem_prms['climate']['paths']['era5_relpath']:
199197
self.var_fp = (
200198
pygem_prms['root']
201199
+ pygem_prms['climate']['paths']['era5_relpath']
@@ -204,6 +202,10 @@ def __init__(self, name=str(), sim_climate_scenario=str(), realization=None):
204202
pygem_prms['root']
205203
+ pygem_prms['climate']['paths']['era5_relpath']
206204
)
205+
else:
206+
self.var_fp = ''
207+
self.fx_fp = ''
208+
207209
# Extra information
208210
self.timestep = pygem_prms['time']['timestep']
209211
self.rgi_lat_colname = pygem_prms['rgi']['rgi_lat_colname']
@@ -232,6 +234,7 @@ def __init__(self, name=str(), sim_climate_scenario=str(), realization=None):
232234
pygem_prms['root']
233235
+ pygem_prms['climate']['paths']['eraint_relpath']
234236
)
237+
235238
# Extra information
236239
self.timestep = pygem_prms['time']['timestep']
237240
self.rgi_lat_colname = pygem_prms['rgi']['rgi_lat_colname']
@@ -454,6 +457,7 @@ def importGCMvarnearestneighbor_xarray(
454457
dates_table,
455458
realizations=['r1i1p1f1', 'r4i1p1f1'],
456459
upscale_var_timestep=False,
460+
verbose=False,
457461
):
458462
"""
459463
Import time series of variables and extract nearest neighbor.
@@ -691,6 +695,37 @@ def importGCMvarnearestneighbor_xarray(
691695
start_idx : end_idx + 1, latlon[0], latlon[1]
692696
].values
693697

698+
# Check all glacier use appropriate climate data
699+
for i, latlon in enumerate(latlon_nearidx):
700+
rgi_id = main_glac_rgi[
701+
pygem_prms['rgi']['rgi_glacno_float_colname']
702+
].values[i]
703+
if (len(data[vn][self.lat_vn].values) == 1) or (
704+
len(data[vn][self.lon_vn].values) == 1
705+
):
706+
if verbose:
707+
warnings.warn(
708+
f'{vn} data has only one latitude or longitude value; check that the correct data is being used',
709+
Warning,
710+
stacklevel=2,
711+
)
712+
else:
713+
lat_res = abs(np.diff(data[vn][self.lat_vn].values)[0])
714+
lon_res = abs(np.diff(data[vn][self.lon_vn].values)[0])
715+
lat_dd = abs(
716+
main_glac_rgi[self.rgi_lat_colname].values[i]
717+
- data[vn][self.lat_vn].values[latlon[0]]
718+
)
719+
lon_dd = abs(
720+
main_glac_rgi[self.rgi_lon_colname].values[i]
721+
- data[vn][self.lon_vn].values[latlon[1]]
722+
)
723+
724+
assert lat_dd <= lat_res and lon_dd <= lon_res, (
725+
f'Climate data pixel for {vn} too from glacier {rgi_id}: Δlat={lat_dd:.3f}, '
726+
+ f'Δlon={lon_dd:.3f}, res=({lat_res:.3f}, {lon_res:.3f})'
727+
)
728+
694729
# Convert to series
695730
glac_variable_series = np.array(
696731
[glac_variable_dict[x] for x in latlon_nearidx]

‎pygem/oggm_compat.py‎

Lines changed: 22 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,7 @@
2525
from pygem.setup.config import ConfigManager
2626

2727
# from oggm.shop import rgitopo
28-
from pygem.shop import debris, icethickness, mbdata
28+
from pygem.shop import debris, icethickness, mbdata, meltextent_and_snowline_1d
2929

3030
# instantiate ConfigManager
3131
config_manager = ConfigManager()
@@ -127,6 +127,16 @@ def single_flowline_glacier_directory(
127127
):
128128
workflow.execute_entity_task(debris.debris_to_gdir, gdir)
129129
workflow.execute_entity_task(debris.debris_binned, gdir)
130+
# 1d melt extent calibration data
131+
if not os.path.isfile(gdir.get_filepath('meltextent_1d')):
132+
workflow.execute_entity_task(
133+
meltextent_and_snowline_1d.meltextent_1d_to_gdir, gdir
134+
)
135+
# 1d snowline calibration data
136+
if not os.path.isfile(gdir.get_filepath('snowline_1d')):
137+
workflow.execute_entity_task(
138+
meltextent_and_snowline_1d.snowline_1d_to_gdir, gdir
139+
)
130140

131141
return gdir
132142

@@ -138,7 +148,7 @@ def single_flowline_glacier_directory_with_calving(
138148
k_calving=1,
139149
logging_level=pygem_prms['oggm']['logging_level'],
140150
has_internet=pygem_prms['oggm']['has_internet'],
141-
working_dir=pygem_prms['root'] + pygem_prms['oggm']['oggm_gdir_relpath'],
151+
working_dir=f'{pygem_prms["root"]}/{pygem_prms["oggm"]["oggm_gdir_relpath"]}',
142152
facorrected=pygem_prms['setup']['include_frontalablation'],
143153
):
144154
"""Prepare a GlacierDirectory for PyGEM (single flowline to start with)
@@ -223,6 +233,16 @@ def single_flowline_glacier_directory_with_calving(
223233
workflow.execute_entity_task(
224234
mbdata.mb_df_to_gdir, gdir, **{'facorrected': facorrected}
225235
)
236+
# 1d melt extent calibration data
237+
if not os.path.isfile(gdir.get_filepath('meltextent_1d')):
238+
workflow.execute_entity_task(
239+
meltextent_and_snowline_1d.meltextent_1d_to_gdir, gdir
240+
)
241+
# 1d snowline calibration data
242+
if not os.path.isfile(gdir.get_filepath('snowline_1d')):
243+
workflow.execute_entity_task(
244+
meltextent_and_snowline_1d.snowline_1d_to_gdir, gdir
245+
)
226246

227247
return gdir
228248

0 commit comments

Comments
 (0)