"""One chunk of a round 5 run, from the gdirs of one v7 experiment. Copy of round4a/run_dyn.py (R5_ variables), which keeps the calibration run itself and can shift the Hugonnet target. Do not run this directly - see run_dyn_slurm.sh and submit.sh. All the settings come from environment variables: R5_GDIR_URL base url of the v7 experiment (.../exps//) R5_PARAMS_FILE the experiment's params.cfg, the same the gdirs were made with R5_CHUNK_IDX, R5_CHUNK_SIZE R5_DYN_YS start year of the calibration runs R5_DYN_ERR_FACTOR the Hugonnet error is scaled by this (0.2 as oggm_prepro) R5_DYN_MAXITER max calibration iterations per glacier R5_DYN_STEP_MIN melt_f_max_step_length_minimum R5_DYN_ID name of the calibration run output (hist_dyn_calib) R5_HUG_SHIFT optional: added to every glacier's Hugonnet target (kg m-2 yr-1; > 0 = less mass loss). Empty or 0: none. R5_RUN_YS optional: a second run with the calibrated melt_f, from this year (R5_RUN_ID). Empty: none. R5_TEST_IDS optional: only these glaciers (smoke test) R5_CALIB 'scan' (scan_calib.py: scan + bisection, default) or 'oggm' (OGGM's run_dynamic_melt_f_calibration) R5_SELECT scan only: choose between roots by 'volume' (default) or 'static' R5_BORDER map border of the gdirs (the experiment's BORDER, default 160) OGGM_WORKDIR, OGGM_OUTDIR The dynamic calibration is OGGM's run_dynamic_melt_f_calibration with the plain dynamic_melt_f_run: no dynamic spinup, i.e. no temp_bias loop matching the RGI area. Each iteration runs the glacier from the RGI-date flowlines in R5_DYN_YS to the end of the climate, and compares its 2000-2020 mass balance with the Hugonnet value in the gdir (observations.yml). Only melt_f changes: prcp_fac and temp_bias keep their static calibration values. Every iteration writes its run under the same suffix, and the last one written is the run with the final melt_f (also for part successes and the fallback), so that is the run we keep. """ import os import logging import geopandas as gpd import pandas as pd import oggm.cfg as cfg from oggm import utils, workflow, tasks GDIR_URL = os.environ['R5_GDIR_URL'] PARAMS_FILE = os.environ['R5_PARAMS_FILE'] CHUNK_IDX = int(os.environ['R5_CHUNK_IDX']) CHUNK_SIZE = int(os.environ['R5_CHUNK_SIZE']) DYN_YS = int(os.environ.get('R5_DYN_YS', 1000)) DYN_ERR_FACTOR = float(os.environ.get('R5_DYN_ERR_FACTOR', 0.2)) DYN_MAXITER = int(os.environ.get('R5_DYN_MAXITER', 40)) DYN_STEP_MIN = float(os.environ.get('R5_DYN_STEP_MIN', 0.1)) DYN_ID = os.environ.get('R5_DYN_ID', 'hist_dyn_calib') HUG_SHIFT = float(os.environ.get('R5_HUG_SHIFT') or 0) RUN_YS = os.environ.get('R5_RUN_YS', '') RUN_ID = os.environ.get('R5_RUN_ID', 'lmr_mira_hist') RGI_REG = '11' BORDER = int(os.environ.get('R5_BORDER') or 160) CALIB = os.environ.get('R5_CALIB') or 'scan' SELECT = os.environ.get('R5_SELECT') or 'volume' # Same params as when the gdirs were made. The calibrated values (melt_f, # prcp_fac, temp_bias, glen_a...) live in each gdir's settings and take # precedence over these anyway. cfg.initialize(file=PARAMS_FILE, logging_level='WORKFLOW') cfg.PARAMS['border'] = BORDER WORKING_DIR = os.environ.get('OGGM_WORKDIR', '') if not WORKING_DIR: raise RuntimeError('Need a working dir') utils.mkdir(WORKING_DIR) cfg.PATHS['working_dir'] = WORKING_DIR OUTPUT_DIR = os.environ.get('OGGM_OUTDIR', '') if not OUTPUT_DIR: raise RuntimeError('Need an output dir') utils.mkdir(OUTPUT_DIR) cfg.PARAMS['use_multiprocessing'] = True cfg.PARAMS['continue_on_error'] = True # round3's variables, plus mass: the dynamic calibration needs it cfg.PARAMS['store_diagnostic_variables'] = ['volume', 'volume_bsl', 'area', 'area_min_h', 'length', 'mass'] cfg.PARAMS['min_ice_thick_for_area'] = 2. cfg.PARAMS['min_ice_thick_for_length'] = 2. workflow.init_mp_pool(True) log = logging.getLogger(__name__) # Our chunk of glaciers: the same blocks of the RGI id space as the chunked # prepro (see workflow.get_rgi_chunk) rgidf = gpd.read_file(utils.get_rgi_region_file(RGI_REG, version='62')) rgidf = workflow.get_rgi_chunk(rgidf, CHUNK_IDX, chunk_size=CHUNK_SIZE) # The Alps only: 11-02 (Pyrenees) is outside the LMR grid, and its gdirs fail # the mass balance calibration anyway rgidf = rgidf.loc[rgidf.O2Region == '1'] TEST_IDS = os.environ.get('R5_TEST_IDS', '').split() if TEST_IDS: rgidf = rgidf.loc[rgidf.RGIId.isin(TEST_IDS)] if len(rgidf) == 0: # The ids have gaps, so an empty chunk is a normal thing log.workflow(f'Chunk {CHUNK_IDX} is empty, nothing to do.') raise SystemExit(0) log.workflow(f'Starting round 5 for RGI reg {RGI_REG}, chunk {CHUNK_IDX} ' f'(size {CHUNK_SIZE}): {len(rgidf)} glaciers') # Our L3 is a delta: OGGM fetches the L1 and L2 it names as parents itself gdirs = workflow.init_glacier_directories(rgidf, from_prepro_level=3, prepro_base_url=GDIR_URL, prepro_rgi_version='62', prepro_border=BORDER) from oggm.core.dynamic_spinup import (run_dynamic_melt_f_calibration, dynamic_melt_f_run, dynamic_melt_f_run_fallback) # Per-glacier targets: the Hugonnet value in the gdir, shifted if asked. # Glaciers without an observation (failed static calibration) are left as # they are: the task fails for them and they get no output, as before. items = [] for gdir in gdirs: kw = {} if HUG_SHIFT: try: ref = dict(gdir.observations['ref_mb']) kw = dict(ref_mb=float(ref['value']) + HUG_SHIFT, ref_mb_err=float(ref['err']), ref_mb_period=ref['period'], overwrite_observations=True) except Exception as e: log.workflow(f'{gdir.rgi_id}: no Hugonnet target to shift ({e!r})') items.append((gdir, kw)) if CALIB == 'scan': from scan_calib import scan_calib, GRID log.workflow(f'Scan calibration, chunk {CHUNK_IDX}: runs from {DYN_YS}, ' f'grid {GRID.tolist()}, Hugonnet error x {DYN_ERR_FACTOR}, ' f'select {SELECT}, target shift {HUG_SHIFT:+g} kg m-2 yr-1') res = workflow.execute_entity_task(scan_calib, gdirs, ys=DYN_YS, err_scaling=DYN_ERR_FACTOR, hug_shift=HUG_SHIFT, select=SELECT, output_filesuffix='_' + DYN_ID) rows = {} for gdir, r in zip(gdirs, res): r = dict(r) if r else {'used_spinup_option': None} r['error'] = gdir.get_error_log() rows[gdir.rgi_id] = r else: log.workflow(f'Dynamic melt_f calibration, chunk {CHUNK_IDX}: runs from ' f'{DYN_YS}, Hugonnet error x {DYN_ERR_FACTOR}, ' f'maxiter {DYN_MAXITER}, min step {DYN_STEP_MIN}, ' f'target shift {HUG_SHIFT:+g} kg m-2 yr-1') workflow.execute_entity_task(run_dynamic_melt_f_calibration, items, ref_mb_err_scaling_factor=DYN_ERR_FACTOR, maxiter=DYN_MAXITER, melt_f_max_step_length_minimum=DYN_STEP_MIN, ys=DYN_YS, run_function=dynamic_melt_f_run, fallback_function=dynamic_melt_f_run_fallback, ignore_errors=True, output_filesuffix='_' + DYN_ID) # What the calibration did, per glacier keys = ['used_spinup_option', 'melt_f_before_dynamic_calibration', 'melt_f_dynamic_calibration', 'dmdtda_mismatch_dynamic_calibration_reference', 'dmdtda_dynamic_calibration_given_error', 'dmdtda_mismatch_with_initial_melt_f', 'dmdtda_mismatch_dynamic_calibration', 'run_dynamic_melt_f_calibration_iterations', 'run_dynamic_melt_f_calibration_stop_reason'] # Read from diagnostics.json: the calibration ran in the worker processes, # and the settings object of the gdirs in this process is a stale copy rows = {} for gdir in gdirs: diag = gdir.get_diagnostics() r = {k: diag.get(k) for k in keys} r['melt_f'] = diag.get('melt_f_dynamic_calibration') r['hug_shift'] = HUG_SHIFT r['error'] = gdir.get_error_log() rows[gdir.rgi_id] = r pd.DataFrame(rows).T.to_csv(os.path.join( OUTPUT_DIR, f'dyn_calib_reg{RGI_REG}_chunk{CHUNK_IDX:03d}.csv')) # The calibration run (DYN_YS -> end, final melt_f) log.workflow(f'Aggregating {DYN_ID}, chunk {CHUNK_IDX}') utils.compile_run_output(gdirs, input_filesuffix='_' + DYN_ID, path=os.path.join(OUTPUT_DIR, f'{DYN_ID}_reg{RGI_REG}_chunk{CHUNK_IDX:03d}.nc')) # Optional: a second run with the calibrated melt_f (it is in the gdir # settings), e.g. from year 1 as the attribution ensemble if RUN_YS: workflow.execute_entity_task(tasks.run_from_climate_data, gdirs, climate_filename='climate_historical', output_filesuffix=RUN_ID, return_value=False, ys=int(RUN_YS)) log.workflow(f'Aggregating {RUN_ID}, chunk {CHUNK_IDX}') utils.compile_run_output(gdirs, input_filesuffix=RUN_ID, path=os.path.join(OUTPUT_DIR, f'{RUN_ID}_reg{RGI_REG}_chunk{CHUNK_IDX:03d}.nc'))