import logging from datetime import timedelta import os import numpy as np import xarray as xr # OGGM core import oggm.cfg as cfg from oggm import utils, tasks, entity_task from oggm.shop.gcm_climate import process_gcm_data log = logging.getLogger(__name__) import seaborn as sns ### copied from Fabien SCENARIO_PARENTS: dict[str, dict] = { # Leaf entries have no parent (or parent = None / piControl not loaded here) "piControl": {"parent": None, "branch_year": None}, "hist": {"parent": None, "branch_year": None}, # standalone 1850–2014 "up2p0": {"parent": None, "branch_year": None}, # standalone from 1850, parallel to hist # GWL1.5 branch and its children "up2p0-gwl1p5": {"parent": "up2p0", "branch_year": 1919}, "up2p0-gwl1p5-50y-dn1p0": {"parent": "up2p0-gwl1p5", "branch_year": 1969}, "up2p0-gwl1p5-200y-dn1p0": {"parent": "up2p0-gwl1p5", "branch_year": 2119}, # GWL2.0 branch and its children "up2p0-gwl2p0": {"parent": "up2p0", "branch_year": 1944}, "up2p0-gwl2p0-50y-dn0p5": {"parent": "up2p0-gwl2p0", "branch_year": 1994}, "up2p0-gwl2p0-50y-dn1p0": {"parent": "up2p0-gwl2p0", "branch_year": 1994}, "up2p0-gwl2p0-50y-dn2p0": {"parent": "up2p0-gwl2p0", "branch_year": 1994}, "up2p0-gwl2p0-200y-dn0p5": {"parent": "up2p0-gwl2p0", "branch_year": 2144}, "up2p0-gwl2p0-200y-dn1p0": {"parent": "up2p0-gwl2p0", "branch_year": 2144}, # GWL3.0 branch and its children "up2p0-gwl3p0": {"parent": "up2p0", "branch_year": 1992}, "up2p0-gwl3p0-50y-dn0p5": {"parent": "up2p0-gwl3p0", "branch_year": 2043}, "up2p0-gwl3p0-50y-dn1p0": {"parent": "up2p0-gwl3p0", "branch_year": 2043}, "up2p0-gwl3p0-50y-dn2p0": {"parent": "up2p0-gwl3p0", "branch_year": 2044}, "up2p0-gwl3p0-200y-dn0p5": {"parent": "up2p0-gwl3p0", "branch_year": 2193}, "up2p0-gwl3p0-200y-dn1p0": {"parent": "up2p0-gwl3p0", "branch_year": 2182}, # GWL4.0 branch and its children "up2p0-gwl4p0": {"parent": "up2p0", "branch_year": 2044}, "up2p0-gwl4p0-50y-dn0p5": {"parent": "up2p0-gwl4p0", "branch_year": 2095}, "up2p0-gwl4p0-50y-dn2p0": {"parent": "up2p0-gwl4p0", "branch_year": 2094}, "up2p0-gwl4p0-50y-dn2p0-gwl2p0": {"parent": "up2p0-gwl4p0-50y-dn2p0", "branch_year": 2233}, "up2p0-gwl4p0-200y-dn1p0": {"parent": "up2p0-gwl4p0", "branch_year": 2245}, # GWL5.0 branch and its children "up2p0-gwl5p0": {"parent": "up2p0", "branch_year": 2082}, "up2p0-gwl5p0-50y-dn2p0": {"parent": "up2p0-gwl5p0", "branch_year": 2132}, "up2p0-gwl5p0-200y-dn0p5": {"parent": "up2p0-gwl5p0", "branch_year": 2282}, "up2p0-gwl5p0-200y-dn1p0": {"parent": "up2p0-gwl5p0", "branch_year": 2282}, "up2p0-gwl5p0-200y-dn2p0": {"parent": "up2p0-gwl5p0", "branch_year": 2282}, # GWL6.0 branch and its children "up2p0-gwl6p0": {"parent": "up2p0", "branch_year": 2137}, "up2p0-gwl6p0-50y-dn1p0": {"parent": "up2p0-gwl6p0", "branch_year": 2187}, "up2p0-gwl6p0-50y-dn2p0": {"parent": "up2p0-gwl6p0", "branch_year": 2187}, "up2p0-gwl6p0-200y-dn0p5": {"parent": "up2p0-gwl6p0", "branch_year": 2337}, "up2p0-gwl6p0-200y-dn1p0": {"parent": "up2p0-gwl6p0", "branch_year": 2337}, "up2p0-gwl6p0-200y-dn2p0": {"parent": "up2p0-gwl6p0", "branch_year": 2337}, } @entity_task(log, writes=['gcm_data']) def process_terrafirma_ukesm_flattened_concat_data( gdir, scenario='up2p0-gwl2p0', y0 = None, y1=None, filesuffix=None, year_range=('1975', '2014'), shift_timeseries_by_year_range=True, **kwargs ): """ Process Terrafirma UKESM1-2-LL GCM data for a single glacier. Extracts monthly temperature and precipitation time series from preprocessed flattened and concatenated Terrafirma-provided UKESM climate projection files at the nearest gridpoint to the glacier, applies unit conversions and finally stores the result in OGGM's standard ``gcm_data`` format. Per default, the timeseries is overwritten to start in 1975. Parameters ---------- gdir : oggm.GlacierDirectory Glacier directory to process. scenario : str Climate scenario identifier. Supported scenarios in SCENARIO_PARENTS y0,y1 : int or str, optional Start, End year for subsetting the scenario data (start year is always the beginning of simulation data). Based on the non-shifted timeseries! filesuffix : str, optional Explicit output filesuffix for the generated ``gcm_data`` file. If not provided, a suffix is constructed from the scenario name and ``year_range``. year_range : tuple(str, str) Reference period used for bias correction (passed through to OGGM). shift_timeseries_by_year_range : bool If True, overwrite the time coordinate such that the climate time series starts at ``year_range[0]``, i.e. per default 1975 for Terrafirma-style scenarios. **kwargs Additional keyword arguments passed to :func:`oggm.shop.gcm_climate.process_gcm_data`. Writes ------ gcm_data : NetCDF Glacier-specific monthly temperature and precipitation time series in OGGM GCM format. """ # here y0 should always be start of the simulation period... assert y0 == None if filesuffix is not None: output_filesuffix = filesuffix else: output_filesuffix = ( f'_terrafirma_UKESM1-2-LL_esm_{scenario}_bc_' f'{year_range[0]}_{year_range[1]}' ) glon = gdir.cenlon glat = gdir.cenlat file_path = ( 'https://cluster.klima.uni-bremen.de/~fmaussion/' f'terrafirma/summary/flattened/ukesm_{scenario}_flat_glaciers.nc' ) time_coder = xr.coders.CFDatetimeCoder(use_cftime=True) # concatenation was done in preprocessing step by Fabien Maussion with utils.get_lock(): f_clim = utils.file_downloader(file_path) with xr.open_dataset(f_clim, decode_times=time_coder) as ds: if ds.lon.min() >= 0 and glon <= 0: glon += 360 # computing all the distances and choose the nearest gridpoint c = ((ds.lon - glon) ** 2 + (ds.lat - glat) ** 2) ds_sel = ds.isel(points=np.argmin(c.data)) temp = ds_sel.tas.sel(time=slice(y0,y1)) precip = ds_sel.pr.sel(time=slice(y0,y1)) assert 'K' in temp.units # Convert kg m-2 s-1 to mm mth-1 => 1 kg m-2 = 1 mm !!! assert 'kg m-2 s-1' in precip.units, 'Precip units not understood' ny, r = divmod(len(temp), 12) assert r == 0 dimo = [cfg.DAYS_IN_MONTH[m - 1] for m in temp['time.month']] precip = precip * dimo * (60 * 60 * 24) # Back to [-180, 180] for OGGM precip.lon.values = precip.lon if precip.lon <= 180 else precip.lon - 360 temp.lon.values = temp.lon if temp.lon <= 180 else temp.lon - 360 if shift_timeseries_by_year_range and scenario != 'hist': new_time = xr.date_range( start=f"{year_range[0]}-01-16", periods=temp.sizes["time"], freq="30D", calendar="360_day", use_cftime=True, ) temp = temp.assign_coords(time=new_time) precip = precip.assign_coords(time=new_time) process_gcm_data( gdir, output_filesuffix=output_filesuffix, prcp=precip, temp=temp, source=output_filesuffix, year_range=year_range, **kwargs, ) import json with open('/home/www/fmaussion/terrafirma/summary/oggm_scenarios_run_until.json', 'r') as f: yr_stop_scenarios = json.load(f) @entity_task(log) def run_with_terrafirma_data( gdir, scenarios=SCENARIOS_sel, stop_early=False, ye=None ): """ Run OGGM forward simulations using Terrafirma UKESM climate forcing. For each specified scenario, the glacier model is initialized from the historical spinup geometry (1979) and forced with Terrafirma-processed UKESM GCM data. The starting year is always the first year of the simulations which we interpret as 1975. Each scenario is run independently so that failures in one scenario do not affect the others. Parameters ---------- gdir : oggm.GlacierDirectory Glacier directory to run. scenarios : iterable of str Climate scenarios to run. Each scenario must correspond to an existing ``gcm_data`` file created by :func:`process_terrafirma_ukesm_data`. stop_early : boolean if True, simulation stopped early (according to https://cluster.klima.uni-bremen.de/~fmaussion/terrafirma/summary/ oggm_scenarios_run_until.json). Otherwise, simulation run until ye. ye : int Default until the end of the simulation period. but some have to end earlier because it gets too cold otherwise. Attention, here ye corresponds to the shifted year ... Writes ------ model_run : NetCDF Glacier evolution outputs for each scenario, stored with a scenario-specific filesuffix. """ #curently stopping critetion independent of region #if stop_early: # assert ye is None # gdir_region = gdir.rgi_region for scenario in scenarios: rid = f'_terrafirma_UKESM1-2-LL_esm_{scenario}_bc_1975_2014' add = '' if stop_early: if scenario == 'hist': ye_stop_early = None else: ye_stop_early = yr_stop_scenarios[scenario] # we still have to shift the estimated year! if ye_stop_early is None: ye = None else: ye = ye_stop_early + (1975-1850) add = '_stop_early' tasks.run_from_climate_data( gdir, ys=1975, ye=ye, init_model_yr=1979, climate_filename='gcm_data', climate_input_filesuffix=rid, init_model_filesuffix='_spinup_historical', output_filesuffix=f'_future_run{add}', )