Issue 128 -- importing melt extent and snowline data into OGGM gdirs

pygem/bin/run/run_calibration.py

@@ -565,25 +565,30 @@ def run(list_packed_vars):
     gcm = class_climate.GCM(name=ref_climate_name)
     # Air temperature [degC]
     gcm_temp, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
-        gcm.temp_fn, gcm.temp_vn, main_glac_rgi, dates_table
+        gcm.temp_fn, gcm.temp_vn, main_glac_rgi, dates_table, verbose=debug
     )
     if pygem_prms['mb']['option_ablation'] == 2 and ref_climate_name in ['ERA5']:
         gcm_tempstd, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
-            gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi, dates_table
+            gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi, dates_table, verbose=debug
         )
     else:
         gcm_tempstd = np.zeros(gcm_temp.shape)
     # Precipitation [m]
     gcm_prec, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
-        gcm.prec_fn, gcm.prec_vn, main_glac_rgi, dates_table
+        gcm.prec_fn, gcm.prec_vn, main_glac_rgi, dates_table, verbose=debug
     )
     # Elevation [m asl]
     gcm_elev = gcm.importGCMfxnearestneighbor_xarray(
         gcm.elev_fn, gcm.elev_vn, main_glac_rgi
     )
     # Lapse rate [degC m-1] (always monthly)
     gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
-        gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table, upscale_var_timestep=True
+        gcm.lr_fn,
+        gcm.lr_vn,
+        main_glac_rgi,
+        dates_table,
+        upscale_var_timestep=True,
+        verbose=debug,
     )
 
     # ===== LOOP THROUGH GLACIERS TO RUN CALIBRATION =====

pygem/bin/run/run_calibration_frontalablation.py

@@ -136,25 +136,25 @@ def reg_calving_flux(
     gcm = class_climate.GCM(name=args.ref_climate_name)
     # Air temperature [degC]
     gcm_temp, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
-        gcm.temp_fn, gcm.temp_vn, main_glac_rgi, dates_table
+        gcm.temp_fn, gcm.temp_vn, main_glac_rgi, dates_table, verbose=debug
     )
     if pygem_prms['mb']['option_ablation'] == 2 and args.ref_climate_name in ['ERA5']:
         gcm_tempstd, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
-            gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi, dates_table
+            gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi, dates_table, verbose=debug
         )
     else:
         gcm_tempstd = np.zeros(gcm_temp.shape)
     # Precipitation [m]
     gcm_prec, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
-        gcm.prec_fn, gcm.prec_vn, main_glac_rgi, dates_table
+        gcm.prec_fn, gcm.prec_vn, main_glac_rgi, dates_table, verbose=debug
     )
     # Elevation [m asl]
     gcm_elev = gcm.importGCMfxnearestneighbor_xarray(
         gcm.elev_fn, gcm.elev_vn, main_glac_rgi
     )
     # Lapse rate [degC m-1]
     gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
-        gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table
+        gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table, verbose=debug
     )
 
     # ===== CALIBRATE ALL THE GLACIERS AT ONCE =====

pygem/bin/run/run_calibration_reg_glena.py

@@ -341,25 +341,29 @@ def main():
         gcm = class_climate.GCM(name=sim_climate_name)
         # Air temperature [degC]
         gcm_temp, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
-            gcm.temp_fn, gcm.temp_vn, main_glac_rgi_subset, dates_table
+            gcm.temp_fn, gcm.temp_vn, main_glac_rgi_subset, dates_table, verbose=debug
         )
         if pygem_prms['mbmod']['option_ablation'] == 2 and sim_climate_name in ['ERA5']:
             gcm_tempstd, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
-                gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi_subset, dates_table
+                gcm.tempstd_fn,
+                gcm.tempstd_vn,
+                main_glac_rgi_subset,
+                dates_table,
+                verbose=debug,
             )
         else:
             gcm_tempstd = np.zeros(gcm_temp.shape)
         # Precipitation [m]
         gcm_prec, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
-            gcm.prec_fn, gcm.prec_vn, main_glac_rgi_subset, dates_table
+            gcm.prec_fn, gcm.prec_vn, main_glac_rgi_subset, dates_table, verbose=debug
         )
         # Elevation [m asl]
         gcm_elev = gcm.importGCMfxnearestneighbor_xarray(
             gcm.elev_fn, gcm.elev_vn, main_glac_rgi_subset
         )
         # Lapse rate [degC m-1]
         gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
-            gcm.lr_fn, gcm.lr_vn, main_glac_rgi_subset, dates_table
+            gcm.lr_fn, gcm.lr_vn, main_glac_rgi_subset, dates_table, verbose=debug
         )
 
         # ===== RUN MASS BALANCE =====

pygem/bin/run/run_inversion.py

@@ -70,19 +70,19 @@ def run(
 
     # Air temperature [degC]
     temp, _ = ref_clim.importGCMvarnearestneighbor_xarray(
-        ref_clim.temp_fn, ref_clim.temp_vn, main_glac_rgi, dt
+        ref_clim.temp_fn, ref_clim.temp_vn, main_glac_rgi, dt, verbose=debug
     )
     # Precipitation [m]
     prec, _ = ref_clim.importGCMvarnearestneighbor_xarray(
-        ref_clim.prec_fn, ref_clim.prec_vn, main_glac_rgi, dt
+        ref_clim.prec_fn, ref_clim.prec_vn, main_glac_rgi, dt, verbose=debug
     )
     # Elevation [m asl]
     elev = ref_clim.importGCMfxnearestneighbor_xarray(
         ref_clim.elev_fn, ref_clim.elev_vn, main_glac_rgi
     )
     # Lapse rate [degC m-1]
     lr, _ = ref_clim.importGCMvarnearestneighbor_xarray(
-        ref_clim.lr_fn, ref_clim.lr_vn, main_glac_rgi, dt
+        ref_clim.lr_fn, ref_clim.lr_vn, main_glac_rgi, dt, verbose=debug
     )
 
     # load prior regionally averaged modelprms (from Rounce et al. 2023)

pygem/bin/run/run_simulation.py

@@ -430,17 +430,17 @@ def run(list_packed_vars):
     # ----- Select Temperature and Precipitation Data -----
     # Air temperature [degC]
     gcm_temp, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
-        gcm.temp_fn, gcm.temp_vn, main_glac_rgi, dates_table_full
+        gcm.temp_fn, gcm.temp_vn, main_glac_rgi, dates_table_full, verbose=debug
     )
     ref_temp, ref_dates = ref_gcm.importGCMvarnearestneighbor_xarray(
-        ref_gcm.temp_fn, ref_gcm.temp_vn, main_glac_rgi, dates_table_ref
+        ref_gcm.temp_fn, ref_gcm.temp_vn, main_glac_rgi, dates_table_ref, verbose=debug
     )
     # Precipitation [m]
     gcm_prec, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
-        gcm.prec_fn, gcm.prec_vn, main_glac_rgi, dates_table_full
+        gcm.prec_fn, gcm.prec_vn, main_glac_rgi, dates_table_full, verbose=debug
     )
     ref_prec, ref_dates = ref_gcm.importGCMvarnearestneighbor_xarray(
-        ref_gcm.prec_fn, ref_gcm.prec_vn, main_glac_rgi, dates_table_ref
+        ref_gcm.prec_fn, ref_gcm.prec_vn, main_glac_rgi, dates_table_ref, verbose=debug
     )
     # Elevation [m asl]
     try:
@@ -548,13 +548,17 @@ def run(list_packed_vars):
         ref_tempstd = np.zeros((main_glac_rgi.shape[0], dates_table_ref.shape[0]))
     elif pygem_prms['mb']['option_ablation'] == 2 and sim_climate_name in ['ERA5']:
         gcm_tempstd, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
-            gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi, dates_table
+            gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi, dates_table, verbose=debug
         )
         ref_tempstd = gcm_tempstd
     elif pygem_prms['mb']['option_ablation'] == 2 and args.ref_climate_name in ['ERA5']:
         # Compute temp std based on reference climate data
         ref_tempstd, ref_dates = ref_gcm.importGCMvarnearestneighbor_xarray(
-            ref_gcm.tempstd_fn, ref_gcm.tempstd_vn, main_glac_rgi, dates_table_ref
+            ref_gcm.tempstd_fn,
+            ref_gcm.tempstd_vn,
+            main_glac_rgi,
+            dates_table_ref,
+            verbose=debug,
         )
         # Monthly average from reference climate data
         gcm_tempstd = gcmbiasadj.monthly_avg_array_rolled(
@@ -567,13 +571,18 @@ def run(list_packed_vars):
     # Lapse rate
     if sim_climate_name in ['ERA-Interim', 'ERA5']:
         gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(
-            gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table
+            gcm.lr_fn,
+            gcm.lr_vn,
+            main_glac_rgi,
+            dates_table,
+            upscale_var_timestep=True,
+            verbose=debug,
         )
         ref_lr = gcm_lr
     else:
         # Compute lapse rates based on reference climate data
         ref_lr, ref_dates = ref_gcm.importGCMvarnearestneighbor_xarray(
-            ref_gcm.lr_fn, ref_gcm.lr_vn, main_glac_rgi, dates_table_ref
+            ref_gcm.lr_fn, ref_gcm.lr_vn, main_glac_rgi, dates_table_ref, verbose=debug
         )
         # Monthly average from reference climate data
         gcm_lr = gcmbiasadj.monthly_avg_array_rolled(
@@ -591,12 +600,9 @@ def run(list_packed_vars):
     else:
         nsims = 1
 
-    # Number of years (for OGGM's run_until_and_store)
-    if pygem_prms['time']['timestep'] == 'monthly':
-        nyears = int(dates_table.shape[0] / 12)
-        nyears_ref = int(dates_table_ref.shape[0] / 12)
-    else:
-        assert True == False, 'Adjust nyears for non-monthly timestep'
+    # Number of years
+    nyears = dates_table.year.unique()[-1] - dates_table.year.unique()[0] + 1
+    nyears_ref = dates_table_ref.year.unique()[-1] - dates_table.year.unique()[0] + 1
 
     for glac in range(main_glac_rgi.shape[0]):
         if glac == 0:

pygem/bin/run/run_spinup.py

@@ -39,19 +39,19 @@ def run(glacno_list, mb_model_params, debug=False, **kwargs):
 
     # Air temperature [degC]
     temp, _ = ref_clim.importGCMvarnearestneighbor_xarray(
-        ref_clim.temp_fn, ref_clim.temp_vn, main_glac_rgi, dt
+        ref_clim.temp_fn, ref_clim.temp_vn, main_glac_rgi, dt, verbose=debug
     )
     # Precipitation [m]
     prec, _ = ref_clim.importGCMvarnearestneighbor_xarray(
-        ref_clim.prec_fn, ref_clim.prec_vn, main_glac_rgi, dt
+        ref_clim.prec_fn, ref_clim.prec_vn, main_glac_rgi, dt, verbose=debug
     )
     # Elevation [m asl]
     elev = ref_clim.importGCMfxnearestneighbor_xarray(
         ref_clim.elev_fn, ref_clim.elev_vn, main_glac_rgi
     )
     # Lapse rate [degC m-1]
     lr, _ = ref_clim.importGCMvarnearestneighbor_xarray(
-        ref_clim.lr_fn, ref_clim.lr_vn, main_glac_rgi, dt
+        ref_clim.lr_fn, ref_clim.lr_vn, main_glac_rgi, dt, verbose=debug
     )
 
     # load prior regionally averaged modelprms (from Rounce et al. 2023)

pygem/class_climate.py

@@ -9,6 +9,7 @@ class of climate data and functions associated with manipulating the dataset to
 """
 
 import os
+import warnings
 
 import numpy as np
 
@@ -192,10 +193,7 @@ def __init__(self, name=str(), sim_climate_scenario=str(), realization=None):
                 self.elev_fn = pygem_prms['climate']['paths']['era5_elev_fn']
                 self.lr_fn = pygem_prms['climate']['paths']['era5_lr_fn']
                 # Variable filepaths
-                if pygem_prms['climate']['paths']['era5_fullpath']:
-                    self.var_fp = ''
-                    self.fx_fp = ''
-                else:
+                if pygem_prms['climate']['paths']['era5_relpath']:
                     self.var_fp = (
                         pygem_prms['root']
                         + pygem_prms['climate']['paths']['era5_relpath']
@@ -204,6 +202,10 @@ def __init__(self, name=str(), sim_climate_scenario=str(), realization=None):
                         pygem_prms['root']
                         + pygem_prms['climate']['paths']['era5_relpath']
                     )
+                else:
+                    self.var_fp = ''
+                    self.fx_fp = ''
+
                 # Extra information
                 self.timestep = pygem_prms['time']['timestep']
                 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):
                     pygem_prms['root']
                     + pygem_prms['climate']['paths']['eraint_relpath']
                 )
+
                 # Extra information
                 self.timestep = pygem_prms['time']['timestep']
                 self.rgi_lat_colname = pygem_prms['rgi']['rgi_lat_colname']
@@ -454,6 +457,7 @@ def importGCMvarnearestneighbor_xarray(
         dates_table,
         realizations=['r1i1p1f1', 'r4i1p1f1'],
         upscale_var_timestep=False,
+        verbose=False,
     ):
         """
         Import time series of variables and extract nearest neighbor.
@@ -691,6 +695,37 @@ def importGCMvarnearestneighbor_xarray(
                         start_idx : end_idx + 1, latlon[0], latlon[1]
                     ].values
 
+            # Check all glacier use appropriate climate data
+            for i, latlon in enumerate(latlon_nearidx):
+                rgi_id = main_glac_rgi[
+                    pygem_prms['rgi']['rgi_glacno_float_colname']
+                ].values[i]
+                if (len(data[vn][self.lat_vn].values) == 1) or (
+                    len(data[vn][self.lon_vn].values) == 1
+                ):
+                    if verbose:
+                        warnings.warn(
+                            f'{vn} data has only one latitude or longitude value; check that the correct data is being used',
+                            Warning,
+                            stacklevel=2,
+                        )
+                else:
+                    lat_res = abs(np.diff(data[vn][self.lat_vn].values)[0])
+                    lon_res = abs(np.diff(data[vn][self.lon_vn].values)[0])
+                    lat_dd = abs(
+                        main_glac_rgi[self.rgi_lat_colname].values[i]
+                        - data[vn][self.lat_vn].values[latlon[0]]
+                    )
+                    lon_dd = abs(
+                        main_glac_rgi[self.rgi_lon_colname].values[i]
+                        - data[vn][self.lon_vn].values[latlon[1]]
+                    )
+
+                    assert lat_dd <= lat_res and lon_dd <= lon_res, (
+                        f'Climate data pixel for {vn} too from glacier {rgi_id}: Δlat={lat_dd:.3f}, '
+                        + f'Δlon={lon_dd:.3f}, res=({lat_res:.3f}, {lon_res:.3f})'
+                    )
+
             # Convert to series
             glac_variable_series = np.array(
                 [glac_variable_dict[x] for x in latlon_nearidx]

pygem/oggm_compat.py

@@ -25,7 +25,7 @@
 from pygem.setup.config import ConfigManager
 
 # from oggm.shop import rgitopo
-from pygem.shop import debris, icethickness, mbdata
+from pygem.shop import debris, icethickness, mbdata, meltextent_and_snowline_1d
 
 # instantiate ConfigManager
 config_manager = ConfigManager()
@@ -127,6 +127,16 @@ def single_flowline_glacier_directory(
     ):
         workflow.execute_entity_task(debris.debris_to_gdir, gdir)
         workflow.execute_entity_task(debris.debris_binned, gdir)
+    # 1d melt extent calibration data
+    if not os.path.isfile(gdir.get_filepath('meltextent_1d')):
+        workflow.execute_entity_task(
+            meltextent_and_snowline_1d.meltextent_1d_to_gdir, gdir
+        )
+    # 1d snowline calibration data
+    if not os.path.isfile(gdir.get_filepath('snowline_1d')):
+        workflow.execute_entity_task(
+            meltextent_and_snowline_1d.snowline_1d_to_gdir, gdir
+        )
 
     return gdir
 
@@ -138,7 +148,7 @@ def single_flowline_glacier_directory_with_calving(
     k_calving=1,
     logging_level=pygem_prms['oggm']['logging_level'],
     has_internet=pygem_prms['oggm']['has_internet'],
-    working_dir=pygem_prms['root'] + pygem_prms['oggm']['oggm_gdir_relpath'],
+    working_dir=f'{pygem_prms["root"]}/{pygem_prms["oggm"]["oggm_gdir_relpath"]}',
     facorrected=pygem_prms['setup']['include_frontalablation'],
 ):
     """Prepare a GlacierDirectory for PyGEM (single flowline to start with)
@@ -223,6 +233,16 @@ def single_flowline_glacier_directory_with_calving(
         workflow.execute_entity_task(
             mbdata.mb_df_to_gdir, gdir, **{'facorrected': facorrected}
         )
+    # 1d melt extent calibration data
+    if not os.path.isfile(gdir.get_filepath('meltextent_1d')):
+        workflow.execute_entity_task(
+            meltextent_and_snowline_1d.meltextent_1d_to_gdir, gdir
+        )
+    # 1d snowline calibration data
+    if not os.path.isfile(gdir.get_filepath('snowline_1d')):
+        workflow.execute_entity_task(
+            meltextent_and_snowline_1d.snowline_1d_to_gdir, gdir
+        )
 
     return gdir
 

pygem/setup/config.py

@@ -157,8 +157,7 @@ def _validate_config(self, config):
         'climate.sim_wateryear': str,
         'climate.constantarea_years': int,
         'climate.paths': dict,
-        'climate.paths.era5_fullpath': bool,
-        'climate.paths.era5_relpath': str,
+        'climate.paths.era5_relpath': (str, type(None)),
         'climate.paths.era5_temp_fn': str,
         'climate.paths.era5_tempstd_fn': str,
         'climate.paths.era5_prec_fn': str,
@@ -223,6 +222,8 @@ def _validate_config(self, config):
         'calib.emulator_params.ftol_opt': float,
         'calib.emulator_params.eps_opt': float,
         'calib.MCMC_params': dict,
+        'calib.MCMC_params.option_calib_meltextent_1d': bool,
+        'calib.MCMC_params.option_calib_snowline_1d': bool,
         'calib.MCMC_params.option_use_emulator': bool,
         'calib.MCMC_params.emulator_sims': int,
         'calib.MCMC_params.tbias_step': float,
@@ -263,6 +264,8 @@ def _validate_config(self, config):
         'calib.data.icethickness': dict,
         'calib.data.icethickness.h_ref_relpath': str,
         'calib.icethickness_cal_frac_byarea': float,
+        'calib.data.meltextent_1d.meltextent_1d_relpath': (str, type(None)),
+        'calib.data.snowline_1d.snowline_1d_relpath': (str, type(None)),
         'sim': dict,
         'sim.option_dynamics': (str, type(None)),
         'sim.option_bias_adjustment': int,