diff --git a/pygem/__init__.py b/pygem/__init__.py index 65c2ce17..725509ee 100755 --- a/pygem/__init__.py +++ b/pygem/__init__.py @@ -3,7 +3,7 @@ copyright © 2018 David Rounce David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license compress OGGM glacier directories """ diff --git a/pygem/bin/op/duplicate_gdirs.py b/pygem/bin/op/duplicate_gdirs.py index af236931..1a3e2bb1 100644 --- a/pygem/bin/op/duplicate_gdirs.py +++ b/pygem/bin/op/duplicate_gdirs.py @@ -3,7 +3,7 @@ copyright © 2024 Brandon Tober David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license duplicate OGGM glacier directories """ diff --git a/pygem/bin/op/initialize.py b/pygem/bin/op/initialize.py index c6f052e1..614a7447 100644 --- a/pygem/bin/op/initialize.py +++ b/pygem/bin/op/initialize.py @@ -3,7 +3,7 @@ copyright © 2024 Brandon Tober David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license initialization script (ensure config.yaml and get sample datasets) """ diff --git a/pygem/bin/op/list_failed_simulations.py b/pygem/bin/op/list_failed_simulations.py index 29e98b22..757cb103 100644 --- a/pygem/bin/op/list_failed_simulations.py +++ b/pygem/bin/op/list_failed_simulations.py @@ -3,7 +3,7 @@ copyright © 2018 David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license script to check for failed glaciers for a given simulation and export a pickle file containing a list of said glacier numbers to be reprocessed """ diff --git a/pygem/bin/postproc/postproc_binned_monthly_mass.py b/pygem/bin/postproc/postproc_binned_monthly_mass.py index cf57db07..b123bced 100644 --- a/pygem/bin/postproc/postproc_binned_monthly_mass.py +++ b/pygem/bin/postproc/postproc_binned_monthly_mass.py @@ -3,7 +3,7 @@ copyright © 2024 Brandon Tober David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license derive binned monthly ice thickness and mass from PyGEM simulation """ diff --git a/pygem/bin/postproc/postproc_compile_simulations.py b/pygem/bin/postproc/postproc_compile_simulations.py index 0de32b39..a9665f87 100644 --- a/pygem/bin/postproc/postproc_compile_simulations.py +++ b/pygem/bin/postproc/postproc_compile_simulations.py @@ -3,7 +3,7 @@ copyright © 2024 Brandon Tober David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license compile individual glacier simulations to the regional level """ diff --git a/pygem/bin/postproc/postproc_distribute_ice.py b/pygem/bin/postproc/postproc_distribute_ice.py index cfbab6fc..ec909397 100644 --- a/pygem/bin/postproc/postproc_distribute_ice.py +++ b/pygem/bin/postproc/postproc_distribute_ice.py @@ -3,7 +3,7 @@ copyright © 2024 Brandon Tober David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license """ # Built-in libraries diff --git a/pygem/bin/postproc/postproc_monthly_mass.py b/pygem/bin/postproc/postproc_monthly_mass.py index 267da544..d3a3678f 100644 --- a/pygem/bin/postproc/postproc_monthly_mass.py +++ b/pygem/bin/postproc/postproc_monthly_mass.py @@ -3,7 +3,7 @@ copyright © 2024 Brandon Tober David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license derive monthly glacierwide mass for PyGEM simulation using annual glacier mass and monthly total mass balance """ diff --git a/pygem/bin/preproc/preproc_fetch_mbdata.py b/pygem/bin/preproc/preproc_fetch_mbdata.py index 70d1f2b5..a10063a7 100644 --- a/pygem/bin/preproc/preproc_fetch_mbdata.py +++ b/pygem/bin/preproc/preproc_fetch_mbdata.py @@ -3,7 +3,7 @@ copyright © 2018 David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license Fetch filled Hugonnet reference mass balance data """ diff --git a/pygem/bin/preproc/preproc_wgms_estimate_kp.py b/pygem/bin/preproc/preproc_wgms_estimate_kp.py index 841cdbf3..f9047452 100644 --- a/pygem/bin/preproc/preproc_wgms_estimate_kp.py +++ b/pygem/bin/preproc/preproc_wgms_estimate_kp.py @@ -3,7 +3,7 @@ copyright © 2018 David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license Process the WGMS data to connect with RGIIds and evaluate potential precipitation biases diff --git a/pygem/bin/run/__init__.py b/pygem/bin/run/__init__.py index e7b24f6e..fefe01a7 100755 --- a/pygem/bin/run/__init__.py +++ b/pygem/bin/run/__init__.py @@ -3,5 +3,5 @@ copyright © 2018 David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license Run model calibration """ @@ -37,15 +37,10 @@ # read the config pygem_prms = config_manager.read_config() +import pygem.oggm_compat as oggm_compat import pygem.pygem_modelsetup as modelsetup from pygem import class_climate, mcmc from pygem.massbalance import PyGEMMassBalance - -# from pygem.glacierdynamics import MassRedistributionCurveModel -from pygem.oggm_compat import ( - single_flowline_glacier_directory, - single_flowline_glacier_directory_with_calving, -) from pygem.utils.stats import mcmc_stats # from oggm.core import climate @@ -175,6 +170,11 @@ def getparser(): action='store_true', help='Flag to keep glacier lists ordered (default is false)', ) + parser.add_argument( + '-spinup', + action='store_true', + help='Flag to use spinup flowlines (default is false)', + ) parser.add_argument( '-p', '--progress_bar', action='store_true', help='Flag to show progress bar' ) @@ -605,11 +605,13 @@ def run(list_packed_vars): glacier_rgi_table['TermType'] not in [1, 5] or not pygem_prms['setup']['include_frontalablation'] ): - gdir = single_flowline_glacier_directory(glacier_str) + gdir = oggm_compat.single_flowline_glacier_directory(glacier_str) gdir.is_tidewater = False else: # set reset=True to overwrite non-calving directory that may already exist - gdir = single_flowline_glacier_directory_with_calving(glacier_str) + gdir = oggm_compat.single_flowline_glacier_directory_with_calving( + glacier_str + ) gdir.is_tidewater = True fls = gdir.read_pickle('inversion_flowlines') @@ -704,6 +706,10 @@ def run(list_packed_vars): 'Mass balance data missing. Check dataset and column names' ) + # if spinup, grab appropriate flowlines + if args.spinup: + fls = oggm_compat.get_spinup_flowlines(gdir, y0=args.ref_startyear) + # ----- CALIBRATION OPTIONS ------ if (fls is not None) and (gdir.mbdata is not None) and (glacier_area.sum() > 0): modelprms = { diff --git a/pygem/bin/run/run_calibration_frontalablation.py b/pygem/bin/run/run_calibration_frontalablation.py index 9523174b..efe2084b 100644 --- a/pygem/bin/run/run_calibration_frontalablation.py +++ b/pygem/bin/run/run_calibration_frontalablation.py @@ -3,7 +3,7 @@ copyright © 2018 David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license Calibrate frontal ablation parameters for tidewater glaciers """ @@ -260,10 +260,10 @@ def reg_calving_flux( cfg.PARAMS['calving_k'] = calving_k cfg.PARAMS['inversion_calving_k'] = cfg.PARAMS['calving_k'] - if pygem_prms['sim']['oggm_dynamics']['use_reg_glena']: + if pygem_prms['sim']['oggm_dynamics']['use_regional_glen_a']: glena_df = pd.read_csv( pygem_prms['root'] - + pygem_prms['sim']['oggm_dynamics']['glena_reg_relpath'] + + pygem_prms['sim']['oggm_dynamics']['glen_a_regional_relpath'] ) glena_idx = np.where(glena_df.O1Region == glacier_rgi_table.O1Region)[ 0 diff --git a/pygem/bin/run/run_calibration_reg_glena.py b/pygem/bin/run/run_calibration_reg_glena.py index e73963b9..5efd9dae 100644 --- a/pygem/bin/run/run_calibration_reg_glena.py +++ b/pygem/bin/run/run_calibration_reg_glena.py @@ -3,7 +3,7 @@ copyright © 2018 David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license Find the optimal values of glens_a_multiplier to match the consensus ice thickness estimates """ @@ -559,7 +559,7 @@ def to_minimize(a_multiplier): try: glena_df = pd.read_csv( - f'{pygem_prms["root"]}/{pygem_prms["out"]["glena_reg_relpath"]}' + f'{pygem_prms["root"]}/{pygem_prms["out"]["glen_a_regional_relpath"]}' ) # Add or overwrite existing file @@ -578,7 +578,7 @@ def to_minimize(a_multiplier): glena_df = glena_df.sort_values('O1Region', ascending=True) glena_df.reset_index(inplace=True, drop=True) glena_df.to_csv( - f'{pygem_prms["root"]}/{pygem_prms["out"]["glena_reg_relpath"]}', + f'{pygem_prms["root"]}/{pygem_prms["out"]["glen_a_regional_relpath"]}', index=False, ) diff --git a/pygem/bin/run/run_inversion.py b/pygem/bin/run/run_inversion.py index d118ad5a..323224fc 100644 --- a/pygem/bin/run/run_inversion.py +++ b/pygem/bin/run/run_inversion.py @@ -1,4 +1,5 @@ import argparse +import json import os from functools import partial @@ -21,6 +22,7 @@ # from pygem.glacierdynamics import MassRedistributionCurveModel from pygem.oggm_compat import update_cfg from pygem.shop import debris, mbdata +from pygem.utils._funcs import str2bool cfg.initialize() cfg.PATHS['working_dir'] = ( @@ -28,28 +30,35 @@ ) -def run(glac_no, ncores=1, debug=False): +def run( + glac_no, ncores=1, calibrate_regional_glen_a=False, reset_gdirs=False, debug=False +): """ Run OGGM's bed inversion for a list of RGI glacier IDs using PyGEM's mass balance model. """ - update_cfg({'continue_on_error': True}, 'PARAMS') + update_cfg({'continue_on_error': False}, 'PARAMS') if ncores > 1: update_cfg({'use_multiprocessing': True}, 'PARAMS') update_cfg({'mp_processes': ncores}, 'PARAMS') + if not isinstance(glac_no, list): + glac_no = [glac_no] main_glac_rgi = modelsetup.selectglaciersrgitable(glac_no=glac_no) # get list of RGIId's for each rgitable being run rgiids = main_glac_rgi['RGIId'].tolist() # initialize glacier directories - gdirs = workflow.init_glacier_directories( - rgiids, - from_prepro_level=2, - prepro_border=cfg.PARAMS['border'], - prepro_base_url=pygem_prms['oggm']['base_url'], - prepro_rgi_version='62', - ) + if reset_gdirs: + gdirs = workflow.init_glacier_directories( + rgiids, + from_prepro_level=2, + prepro_border=cfg.PARAMS['border'], + prepro_base_url=pygem_prms['oggm']['base_url'], + prepro_rgi_version='62', + ) + else: + gdirs = workflow.init_glacier_directories(rgiids) # PyGEM setup - model datestable, climate data import, prior model parameters # model dates @@ -134,12 +143,10 @@ def run(glac_no, ncores=1, debug=False): overwrite_gdir=True, ) - ############################## - ### INVERSION - no calving ### - ############################## - if debug: - print('Running initial inversion') - # note, PyGEMMassBalance_wrapper is passed to `tasks.apparent_mb_from_any_mb` as the `mb_model_class` so that PyGEMs mb model is used for inversion + ####################################### + ### CALCULATE APPARENT MASS BALANCE ### + ####################################### + # note, PyGEMMassBalance_wrapper is passed to `tasks.apparent_mb_from_any_mb` as the `mb_model_class` so that PyGEMs mb model is used for apparent mb workflow.execute_entity_task( tasks.apparent_mb_from_any_mb, gdirs, @@ -147,7 +154,6 @@ def run(glac_no, ncores=1, debug=False): PyGEMMassBalance_wrapper, fl_str='inversion_flowlines', option_areaconstant=True, - inversion_filter=True, ), ) # add debris data to flowlines @@ -159,69 +165,112 @@ def run(glac_no, ncores=1, debug=False): ### CALIBRATE GLEN'S A ### ########################## # fit ice thickness to consensus estimates to find "best" Glen's A - if debug: - print("Calibrating Glen's A") - workflow.calibrate_inversion_from_consensus( - gdirs, - apply_fs_on_mismatch=True, - error_on_mismatch=False, # if you running many glaciers some might not work - filter_inversion_output=True, # this partly filters the overdeepening due to - # the equilibrium assumption for retreating glaciers (see. Figure 5 of Maussion et al. 2019) - volume_m3_reference=None, # here you could provide your own total volume estimate in m3 - ) + # note, this runs task.mass_conservation_inversion internally + if calibrate_regional_glen_a: + if debug: + print("Calibrating Glen's A") + workflow.calibrate_inversion_from_consensus( + gdirs, + apply_fs_on_mismatch=True, + error_on_mismatch=False, # if you running many glaciers some might not work + filter_inversion_output=True, # this partly filters the overdeepening due to + # the equilibrium assumption for retreating glaciers (see. Figure 5 of Maussion et al. 2019) + volume_m3_reference=None, # here you could provide your own total volume estimate in m3 + ) - ################################ - ### INVERSION - with calving ### - ################################ - cfg.PARAMS['use_kcalving_for_inversion'] = True for gdir in gdirs: - # note, tidewater glacier inversion is not done in parallel since there's not currently a way to pass a different inversion_calving_k to each gdir + if calibrate_regional_glen_a: + glen_a = gdir.get_diagnostics()['inversion_glen_a'] + fs = gdir.get_diagnostics()['inversion_fs'] + else: + # get glen_a and fs values from prior calibration or manual entry + if pygem_prms['sim']['oggm_dynamics']['use_regional_glen_a']: + glen_a_df = pd.read_csv( + f'{pygem_prms["root"]}/{pygem_prms["sim"]["oggm_dynamics"]["glen_a_regional_relpath"]}' + ) + glen_a_O1regions = [int(x) for x in glen_a_df.O1Region.values] + assert gdir.glacier_rgi_table.O1Region in glen_a_O1regions, ( + '{0:0.5f}'.format(gd.glacier_rgi_table['RGIId_float']) + + ' O1 region not in glen_a_df' + ) + glen_a_idx = np.where( + glen_a_O1regions == gdir.glacier_rgi_table.O1Region + )[0][0] + glen_a_multiplier = glen_a_df.loc[glen_a_idx, 'glens_a_multiplier'] + fs = glen_a_df.loc[glen_a_idx, 'fs'] + else: + glen_a_multiplier = pygem_prms['sim']['oggm_dynamics'][ + 'glen_a_multiplier' + ] + fs = pygem_prms['sim']['oggm_dynamics']['fs'] + glen_a = cfg.PARAMS['glen_a'] * glen_a_multiplier + + # non-tidewater if ( gdir.glacier_rgi_table['TermType'] not in [1, 5] or not pygem_prms['setup']['include_frontalablation'] ): - continue - if debug: - print(f'Running inversion for {gdir.rgi_id} with calving') - - # Load quality controlled frontal ablation data - fp = f'{pygem_prms["root"]}/{pygem_prms["calib"]["data"]["frontalablation"]["frontalablation_relpath"]}/analysis/{pygem_prms["calib"]["data"]["frontalablation"]["frontalablation_cal_fn"]}' - assert os.path.exists(fp), 'Calibrated calving dataset does not exist' - calving_df = pd.read_csv(fp) - calving_rgiids = list(calving_df.RGIId) - - # Use calibrated value if individual data available - if gdir.rgi_id in calving_rgiids: - calving_idx = calving_rgiids.index(gdir.rgi_id) - calving_k = calving_df.loc[calving_idx, 'calving_k'] - # Otherwise, use region's median value - else: - calving_df['O1Region'] = [ - int(x.split('-')[1].split('.')[0]) for x in calving_df.RGIId.values - ] - calving_df_reg = calving_df.loc[ - calving_df['O1Region'] == int(gdir.rgi_id[6:8]), : - ] - calving_k = np.median(calving_df_reg.calving_k) - - # increase calving line for inversion so that later spinup will work - cfg.PARAMS['calving_line_extension'] = 120 - # set inversion_calving_k - cfg.PARAMS['inversion_calving_k'] = calving_k - if debug: - print(f'inversion_calving_k = {calving_k}') + if calibrate_regional_glen_a: + # nothing else to do here - already ran inversion when calibrating Glen's A + continue + + # run inversion using regionally calibrated Glen's A values + cfg.PARAMS['use_kcalving_for_inversion'] = False + + ############################### + ### INVERSION - no calving ### + ################################ + tasks.prepare_for_inversion(gdir) + tasks.mass_conservation_inversion( + gdir, + glen_a=glen_a, + fs=fs, + ) - tasks.find_inversion_calving_from_any_mb( - gdir, - mb_model=PyGEMMassBalance_wrapper( + # tidewater + else: + cfg.PARAMS['use_kcalving_for_inversion'] = True + + # Load quality controlled frontal ablation data + fp = f'{pygem_prms["root"]}/{pygem_prms["calib"]["data"]["frontalablation"]["frontalablation_relpath"]}/analysis/{pygem_prms["calib"]["data"]["frontalablation"]["frontalablation_cal_fn"]}' + assert os.path.exists(fp), 'Calibrated calving dataset does not exist' + calving_df = pd.read_csv(fp) + calving_rgiids = list(calving_df.RGIId) + + # Use calibrated value if individual data available + if gdir.rgi_id in calving_rgiids: + calving_idx = calving_rgiids.index(gdir.rgi_id) + calving_k = calving_df.loc[calving_idx, 'calving_k'] + # Otherwise, use region's median value + else: + calving_df['O1Region'] = [ + int(x.split('-')[1].split('.')[0]) for x in calving_df.RGIId.values + ] + calving_df_reg = calving_df.loc[ + calving_df['O1Region'] == int(gdir.rgi_id[6:8]), : + ] + calving_k = np.median(calving_df_reg.calving_k) + + # increase calving line for inversion so that later spinup will work + cfg.PARAMS['calving_line_extension'] = 120 + # set inversion_calving_k + cfg.PARAMS['inversion_calving_k'] = calving_k + if debug: + print(f'inversion_calving_k = {calving_k}') + + ################################ + ### INVERSION - with calving ### + ################################ + tasks.find_inversion_calving_from_any_mb( gdir, - fl_str='inversion_flowlines', - option_areaconstant=True, - inversion_filter=True, - ), - glen_a=gdir.get_diagnostics()['inversion_glen_a'], - fs=gdir.get_diagnostics()['inversion_fs'], - ) + mb_model=PyGEMMassBalance_wrapper( + gdir, + fl_str='inversion_flowlines', + option_areaconstant=True, + ), + glen_a=glen_a, + fs=fs, + ) ###################### ### POSTPROCESSING ### @@ -236,7 +285,7 @@ def run(glac_no, ncores=1, debug=False): def main(): # define ArgumentParser parser = argparse.ArgumentParser( - description="perform glacier bed inversion (defaults to find best Glen's A for each RGI order 01 region)" + description="Perform glacier bed inversion (defaults to find best Glen's A for each RGI order 01 region)" ) # add arguments parser.add_argument( @@ -246,31 +295,76 @@ def main(): help='Randoph Glacier Inventory region (can take multiple, e.g. `-run_region01 1 2 3`)', nargs='+', ) + parser.add_argument( + '-rgi_glac_number', + action='store', + type=float, + default=pygem_prms['setup']['glac_no'], + nargs='+', + help='Randoph Glacier Inventory glacier number (can take multiple)', + ) + ( + parser.add_argument( + '-rgi_glac_number_fn', + action='store', + type=str, + default=None, + help='Filepath containing list of rgi_glac_number, helpful for running batches on spc', + ), + ) + parser.add_argument( + '-calibrate_regional_glen_a', + type=str2bool, + default=True, + help="If True (False) run ice thickness inversion and regionally calibrate (use previously calibrated or user-input) Glen's A values. Default is True", + ) parser.add_argument( '-ncores', action='store', type=int, default=1, - help='number of simultaneous processes (cores) to use', + help='Number of simultaneous processes (cores) to use', + ) + parser.add_argument( + '-reset_gdirs', + action='store_true', + help='If True (False) reset OGGM glacier directories. Default is False', ) parser.add_argument('-v', '--debug', action='store_true', help='Flag for debugging') args = parser.parse_args() - # RGI glacier number - batches = [ - modelsetup.selectglaciersrgitable( - rgi_regionsO1=[r01], - rgi_regionsO2='all', - include_landterm=pygem_prms['setup']['include_landterm'], - include_laketerm=pygem_prms['setup']['include_laketerm'], - include_tidewater=pygem_prms['setup']['include_tidewater'], - min_glac_area_km2=pygem_prms['setup']['min_glac_area_km2'], - )['rgino_str'].values.tolist() - for r01 in args.rgi_region01 - ] + # RGI glacier batches + if args.rgi_region01: + batches = [ + modelsetup.selectglaciersrgitable( + rgi_regionsO1=[r01], + rgi_regionsO2='all', + include_landterm=pygem_prms['setup']['include_landterm'], + include_laketerm=pygem_prms['setup']['include_laketerm'], + include_tidewater=pygem_prms['setup']['include_tidewater'], + min_glac_area_km2=pygem_prms['setup']['min_glac_area_km2'], + )['rgino_str'].values.tolist() + for r01 in args.rgi_region01 + ] + else: + batches = None + if args.rgi_glac_number: + glac_no = args.rgi_glac_number + # format appropriately + glac_no = [float(g) for g in glac_no] + batches = [f'{g:.5f}' if g >= 10 else f'0{g:.5f}' for g in glac_no] + elif args.rgi_glac_number_fn is not None: + with open(args.rgi_glac_number_fn, 'r') as f: + batches = json.load(f) # set up partial function with common arguments - run_partial = partial(run, ncores=args.ncores, debug=args.debug) + run_partial = partial( + run, + ncores=args.ncores, + calibrate_regional_glen_a=args.calibrate_regional_glen_a, + reset_gdirs=args.reset_gdirs, + debug=args.debug, + ) for i, batch in enumerate(batches): run_partial(batch) diff --git a/pygem/bin/run/run_simulation.py b/pygem/bin/run/run_simulation.py index 7bf357d2..0f365d7c 100755 --- a/pygem/bin/run/run_simulation.py +++ b/pygem/bin/run/run_simulation.py @@ -3,7 +3,7 @@ copyright © 2018 David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license Run a model simulation """ @@ -269,10 +269,10 @@ def getparser(): help='glacier dynamics scheme (options: ``OGGM`, `MassRedistributionCurves`, `None`)', ) parser.add_argument( - '-use_reg_glena', + '-use_regional_glen_a', action='store', type=bool, - default=pygem_prms['sim']['oggm_dynamics']['use_reg_glena'], + default=pygem_prms['sim']['oggm_dynamics']['use_regional_glen_a'], help='Take the glacier flow parameterization from regionally calibrated priors (boolean: `0` or `1`, `True` or `False`)', ) parser.add_argument( @@ -335,6 +335,12 @@ def getparser(): action='store_true', help='Flag to keep glacier lists ordered (default is off)', ) + parser.add_argument( + '-spinup', + action='store_true', + default=False, + help='Flag to perform dynamical spinup before calibration', + ) parser.add_argument('-v', '--debug', action='store_true', help='Flag for debugging') return parser @@ -839,9 +845,9 @@ def run(list_packed_vars): if debug: print('cfl number:', cfg.PARAMS['cfl_number']) - if args.use_reg_glena: + if args.use_regional_glen_a: glena_df = pd.read_csv( - f'{pygem_prms["root"]}/{pygem_prms["sim"]["oggm_dynamics"]["glena_reg_relpath"]}' + f'{pygem_prms["root"]}/{pygem_prms["sim"]["oggm_dynamics"]["glen_a_regional_relpath"]}' ) glena_O1regions = [int(x) for x in glena_df.O1Region.values] assert glacier_rgi_table.O1Region in glena_O1regions, ( @@ -860,6 +866,20 @@ def run(list_packed_vars): glen_a_multiplier = pygem_prms['sim']['oggm_dynamics'][ 'glen_a_multiplier' ] + glen_a = cfg.PARAMS['glen_a'] * glen_a_multiplier + + # spinup + if args.spinup: + try: + # see if model_flowlines from spinup exist + nfls = gdir.read_pickle( + 'model_flowlines', + filesuffix=f'_dynamic_spinup_pygem_mb_{args.sim_startyear}', + ) + except: + raise + glen_a = gdir.get_diagnostics()['inversion_glen_a'] + fs = gdir.get_diagnostics()['inversion_fs'] # Time attributes and values if pygem_prms['climate']['sim_wateryear'] == 'hydro': @@ -978,65 +998,59 @@ def run(list_packed_vars): else: inversion_filter = False - # Perform inversion based on PyGEM MB using reference directory - mbmod_inv = PyGEMMassBalance( - gdir_ref, - modelprms, - glacier_rgi_table, - fls=fls, - option_areaconstant=True, - inversion_filter=inversion_filter, - ) - # if debug: - # h, w = gdir.get_inversion_flowline_hw() - # mb_t0 = (mbmod_inv.get_annual_mb(h, year=0, fl_id=0, fls=fls) * cfg.SEC_IN_YEAR * - # pygem_prms['constants']['density_ice'] / pygem_prms['constants']['density_water']) - # plt.plot(mb_t0, h, '.') - # plt.ylabel('Elevation') - # plt.xlabel('Mass balance (mwea)') - # plt.show() - - # Non-tidewater glaciers - if ( - not gdir.is_tidewater - or not pygem_prms['setup']['include_frontalablation'] - ): - # Arbitrariliy shift the MB profile up (or down) until mass balance is zero (equilibrium for inversion) - apparent_mb_from_any_mb(gdir, mb_model=mbmod_inv) - tasks.prepare_for_inversion(gdir) - tasks.mass_conservation_inversion( - gdir, - glen_a=cfg.PARAMS['glen_a'] * glen_a_multiplier, - fs=fs, + # run inversion + if not args.spinup: + # Perform inversion based on PyGEM MB using reference directory + mbmod_inv = PyGEMMassBalance( + gdir_ref, + modelprms, + glacier_rgi_table, + fls=fls, + option_areaconstant=True, + inversion_filter=inversion_filter, ) - # Tidewater glaciers - else: - cfg.PARAMS['use_kcalving_for_inversion'] = True - cfg.PARAMS['use_kcalving_for_run'] = True - tasks.find_inversion_calving_from_any_mb( - gdir, - mb_model=mbmod_inv, - glen_a=cfg.PARAMS['glen_a'] * glen_a_multiplier, - fs=fs, - ) + # Non-tidewater glaciers + if ( + not gdir.is_tidewater + or not pygem_prms['setup']['include_frontalablation'] + ): + # Arbitrariliy shift the MB profile up (or down) until mass balance is zero (equilibrium for inversion) + apparent_mb_from_any_mb(gdir, mb_model=mbmod_inv) + tasks.prepare_for_inversion(gdir) + tasks.mass_conservation_inversion( + gdir, + glen_a=cfg.PARAMS['glen_a'] * glen_a_multiplier, + fs=fs, + ) - # ----- INDENTED TO BE JUST WITH DYNAMICS ----- - tasks.init_present_time_glacier(gdir) # adds bins below + # Tidewater glaciers + else: + cfg.PARAMS['use_kcalving_for_inversion'] = True + cfg.PARAMS['use_kcalving_for_run'] = True + tasks.find_inversion_calving_from_any_mb( + gdir, + mb_model=mbmod_inv, + glen_a=cfg.PARAMS['glen_a'] * glen_a_multiplier, + fs=fs, + ) - if not os.path.isfile(gdir.get_filepath('model_flowlines')): - tasks.compute_downstream_line(gdir) - tasks.compute_downstream_bedshape(gdir) + # ----- INDENTED TO BE JUST WITH DYNAMICS ----- tasks.init_present_time_glacier(gdir) # adds bins below - try: - if pygem_prms['mb']['include_debris']: - debris.debris_binned( - gdir, fl_str='model_flowlines' - ) # add debris enhancement factors to flowlines - nfls = gdir.read_pickle('model_flowlines') - except: - raise + if not os.path.isfile(gdir.get_filepath('model_flowlines')): + tasks.compute_downstream_line(gdir) + tasks.compute_downstream_bedshape(gdir) + tasks.init_present_time_glacier(gdir) # adds bins below + + try: + if pygem_prms['mb']['include_debris']: + debris.debris_binned( + gdir, fl_str='model_flowlines' + ) # add debris enhancement factors to flowlines + nfls = gdir.read_pickle('model_flowlines') + except: + raise # Water Level # Check that water level is within given bounds diff --git a/pygem/bin/run/run_spinup.py b/pygem/bin/run/run_spinup.py index 2ae7eb24..102e93a8 100644 --- a/pygem/bin/run/run_spinup.py +++ b/pygem/bin/run/run_spinup.py @@ -27,12 +27,13 @@ ) -def run(glacno_list, spinup_start_yr, **kwargs): +def run(glacno_list, mb_model_params, debug=False, **kwargs): + # remove None kwargs + kwargs = {k: v for k, v in kwargs.items() if v is not None} + main_glac_rgi = modelsetup.selectglaciersrgitable(glac_no=glacno_list) # model dates - dt = modelsetup.datesmodelrun( - startyear=spinup_start_yr, endyear=2019 - ) # will have to cover the time period of inversion (2000-2019) and spinup (1979-~2010 by default) + dt = modelsetup.datesmodelrun(startyear=1979, endyear=2019) # load climate data ref_clim = class_climate.GCM(name='ERA5') @@ -91,22 +92,44 @@ def run(glacno_list, spinup_start_yr, **kwargs): } gd.dates_table = dt - # get modelprms from regional priors - priors_idx = np.where( - (priors_df.O1Region == gd.glacier_rgi_table['O1Region']) - & (priors_df.O2Region == gd.glacier_rgi_table['O2Region']) - )[0][0] - tbias_mu = float(priors_df.loc[priors_idx, 'tbias_mean']) - kp_mu = float(priors_df.loc[priors_idx, 'kp_mean']) - gd.modelprms = { - 'kp': kp_mu, - 'tbias': tbias_mu, - 'ddfsnow': pygem_prms['calib']['MCMC_params']['ddfsnow_mu'], - 'ddfice': pygem_prms['calib']['MCMC_params']['ddfsnow_mu'] - / pygem_prms['sim']['params']['ddfsnow_iceratio'], - 'precgrad': pygem_prms['sim']['params']['precgrad'], - 'tsnow_threshold': pygem_prms['sim']['params']['tsnow_threshold'], - } + if mb_model_params == 'regional_priors': + # get modelprms from regional priors + priors_idx = np.where( + (priors_df.O1Region == gd.glacier_rgi_table['O1Region']) + & (priors_df.O2Region == gd.glacier_rgi_table['O2Region']) + )[0][0] + tbias_mu = float(priors_df.loc[priors_idx, 'tbias_mean']) + kp_mu = float(priors_df.loc[priors_idx, 'kp_mean']) + gd.modelprms = { + 'kp': kp_mu, + 'tbias': tbias_mu, + 'ddfsnow': pygem_prms['calib']['MCMC_params']['ddfsnow_mu'], + 'ddfice': pygem_prms['calib']['MCMC_params']['ddfsnow_mu'] + / pygem_prms['sim']['params']['ddfsnow_iceratio'], + 'precgrad': pygem_prms['sim']['params']['precgrad'], + 'tsnow_threshold': pygem_prms['sim']['params']['tsnow_threshold'], + } + elif mb_model_params == 'emulator': + # get modelprms from emulator mass balance calibration + modelprms_fn = glacier_str + '-modelprms_dict.json' + modelprms_fp = ( + pygem_prms['root'] + + '/Output/calibration/' + + glacier_str.split('.')[0].zfill(2) + + '/' + ) + modelprms_fn + with open(modelprms_fp, 'r') as f: + modelprms_dict = json.load(f) + + modelprms_all = modelprms_dict['emulator'] + gd.modelprms = { + 'kp': modelprms_all['kp'][0], + 'tbias': modelprms_all['tbias'][0], + 'ddfsnow': modelprms_all['ddfsnow'][0], + 'ddfice': modelprms_all['ddfice'][0], + 'tsnow_threshold': modelprms_all['tsnow_threshold'][0], + 'precgrad': modelprms_all['precgrad'][0], + } # update cfg.PARAMS update_cfg({'continue_on_error': True}, 'PARAMS') @@ -116,10 +139,10 @@ def run(glacno_list, spinup_start_yr, **kwargs): workflow.execute_entity_task( tasks.run_dynamic_spinup, gd, - spinup_start_yr=spinup_start_yr, # When to start the spinup + # spinup_start_yr=spinup_start_yr, # When to start the spinup minimise_for='area', # what target to match at the RGI date # target_yr=target_yr, # The year at which we want to match area or volume. If None, gdir.rgi_date + 1 is used (the default) - # ye=, # When the simulation should stop + # ye=ye, # When the simulation should stop output_filesuffix='_dynamic_spinup_pygem_mb', store_fl_diagnostics=True, store_model_geometry=True, @@ -169,10 +192,15 @@ def main(): action='store', type=str, default=None, - help='filepath containing list of rgi_glac_number, helpful for running batches on spc', + help='Filepath containing list of rgi_glac_number, helpful for running batches on spc', ), ) - parser.add_argument('-spinup_start_yr', type=int, default=1979) + parser.add_argument( + '-spinup_period', + type=int, + default=None, + help='Fixed spinup period (years). If not provided, OGGM default is used.', + ) parser.add_argument('-target_yr', type=int, default=None) parser.add_argument('-ye', type=int, default=2020) parser.add_argument( @@ -180,7 +208,14 @@ def main(): action='store', type=int, default=1, - help='number of simultaneous processes (cores) to use', + help='Number of simultaneous processes (cores) to use', + ) + parser.add_argument( + '-mb_model_params', + type=str, + default='regional_priors', + choices=['regional_priors', 'emulator'], + help='Which mass balance model parameters to use ("regional_priors" or "emulator")', ) args = parser.parse_args() @@ -221,7 +256,10 @@ def main(): # set up partial function with debug argument run_partial = partial( - run, spinup_start_yr=args.spinup_start_yr, target_yr=args.target_yr, ye=args.ye + run, + mb_model_params=args.mb_model_params, + target_yr=args.target_yr, + spinup_period=args.spinup_period, ) # parallel processing print(f'Processing with {ncores} cores... \n{glac_no_lsts}') diff --git a/pygem/class_climate.py b/pygem/class_climate.py index e39e19dd..300d13d5 100755 --- a/pygem/class_climate.py +++ b/pygem/class_climate.py @@ -3,7 +3,7 @@ copyright © 2018 David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license Run bias adjustments a given climate dataset """ diff --git a/pygem/glacierdynamics.py b/pygem/glacierdynamics.py index 5ba4b8ad..02b607df 100755 --- a/pygem/glacierdynamics.py +++ b/pygem/glacierdynamics.py @@ -3,7 +3,7 @@ copyright © 2018 David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license """ diff --git a/pygem/setup/config.py b/pygem/setup/config.py index 33cdc1d0..86d28094 100644 --- a/pygem/setup/config.py +++ b/pygem/setup/config.py @@ -3,7 +3,7 @@ copyright © 2018 David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license """ import os @@ -129,8 +129,8 @@ def _validate_config(self, config): 'user.institution': (str, type(None)), 'user.email': (str, type(None)), 'setup': dict, - 'setup.rgi_region01': list, - 'setup.rgi_region02': str, + 'setup.rgi_region01': (list, type(None)), + 'setup.rgi_region02': (str, type(None)), 'setup.glac_no_skip': (list, type(None)), 'setup.glac_no': (list, type(None)), 'setup.min_glac_area_km2': int, @@ -276,8 +276,8 @@ def _validate_config(self, config): 'sim.oggm_dynamics': dict, 'sim.oggm_dynamics.cfl_number': float, 'sim.oggm_dynamics.cfl_number_calving': float, - 'sim.oggm_dynamics.glena_reg_relpath': str, - 'sim.oggm_dynamics.use_reg_glena': bool, + 'sim.oggm_dynamics.glen_a_regional_relpath': str, + 'sim.oggm_dynamics.use_regional_glen_a': bool, 'sim.oggm_dynamics.fs': int, 'sim.oggm_dynamics.glen_a_multiplier': int, 'sim.icethickness_advancethreshold': int, diff --git a/pygem/setup/config.yaml b/pygem/setup/config.yaml index 57b6ae27..6bb8e578 100644 --- a/pygem/setup/config.yaml +++ b/pygem/setup/config.yaml @@ -226,8 +226,8 @@ sim: oggm_dynamics: cfl_number: 0.02 # Time step threshold (seconds) cfl_number_calving: 0.01 # Time step threshold for marine-terimating glaciers (seconds) - use_reg_glena: true - glena_reg_relpath: /Output/calibration/glena_region.csv + use_regional_glen_a: true + glen_a_regional_relpath: /Output/calibration/glena_region.csv # glen_a multiplier if not using regionally calibrated glens_a fs: 0 glen_a_multiplier: 1 diff --git a/pygem/shop/debris.py b/pygem/shop/debris.py index 0c4c922b..66c548f4 100755 --- a/pygem/shop/debris.py +++ b/pygem/shop/debris.py @@ -3,7 +3,7 @@ copyright © 2018 David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license """ import logging diff --git a/pygem/shop/icethickness.py b/pygem/shop/icethickness.py index 256bf526..a17cfbb4 100755 --- a/pygem/shop/icethickness.py +++ b/pygem/shop/icethickness.py @@ -3,7 +3,7 @@ copyright © 2018 David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license """ import logging diff --git a/pygem/shop/mbdata.py b/pygem/shop/mbdata.py index b62dc856..9b23897b 100755 --- a/pygem/shop/mbdata.py +++ b/pygem/shop/mbdata.py @@ -3,7 +3,7 @@ copyright © 2018 David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license """ # Built-in libaries diff --git a/pygem/shop/oib.py b/pygem/shop/oib.py index b7be9611..a35ce02d 100644 --- a/pygem/shop/oib.py +++ b/pygem/shop/oib.py @@ -3,7 +3,7 @@ copyright © 2024 Brandon Tober , David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license NASA Operation IceBridge data and processing class """ diff --git a/pygem/tests/__init__.py b/pygem/tests/__init__.py index e7b24f6e..fefe01a7 100755 --- a/pygem/tests/__init__.py +++ b/pygem/tests/__init__.py @@ -3,5 +3,5 @@ copyright © 2018 David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license Functions that didn't fit into other modules """ +import argparse import json import numpy as np @@ -20,6 +21,25 @@ pygem_prms = config_manager.read_config() +def str2bool(v): + """ + Convert a string to a boolean. + + Accepts: "yes", "true", "t", "1" → True; + "no", "false", "f", "0" → False. + + Raises an error if input is unrecognized. + """ + if isinstance(v, bool): + return v + if v.lower() in ('yes', 'true', 't', '1', 'y'): + return True + elif v.lower() in ('no', 'false', 'f', '0', 'n'): + return False + else: + raise argparse.ArgumentTypeError('Boolean value expected.') + + def annualweightedmean_array(var, dates_table): """ Calculate annual mean of variable according to the timestep. diff --git a/pygem/utils/_funcs_selectglaciers.py b/pygem/utils/_funcs_selectglaciers.py index 742d9314..86fcff02 100644 --- a/pygem/utils/_funcs_selectglaciers.py +++ b/pygem/utils/_funcs_selectglaciers.py @@ -3,7 +3,7 @@ copyright © 2018 David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license Functions of different ways to select glaciers """ diff --git a/pygem/utils/stats.py b/pygem/utils/stats.py index a09b15a2..8f61f98d 100644 --- a/pygem/utils/stats.py +++ b/pygem/utils/stats.py @@ -3,7 +3,7 @@ copyright © 2018 David Rounce -Distrubted under the MIT lisence +Distributed under the MIT license Model statistics module """