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 # should be the same everywhere ... pal_colorblind = sns.color_palette("colorblind") scenarios = ['up2p0-gwl4p0-50y-dn2p0', 'up2p0-gwl4p0', 'up2p0-gwl2p0-50y-dn2p0', 'up2p0-gwl2p0', 'up2p0', 'hist'] color_scenario = {'hist': 'black', 'up2p0': pal_colorblind[4], 'up2p0-gwl4p0': pal_colorblind[3], 'up2p0-gwl4p0-50y-dn2p0': pal_colorblind[1], 'up2p0-gwl2p0': pal_colorblind[2], 'up2p0-gwl2p0-50y-dn2p0': pal_colorblind[0]} ls_scenario_tuple = { 'hist': (None, None), # solid 'up2p0': (None, None), # solid 'up2p0-gwl4p0': (1, 2), # dotted 'up2p0-gwl4p0-50y-dn2p0': (6, 2), # dashed 'up2p0-gwl2p0': (1, 2), # dotted 'up2p0-gwl2p0-50y-dn2p0': (6, 2) # dashed } ls_scenario = {'hist': '-', 'up2p0': '-', 'up2p0-gwl4p0': ':', 'up2p0-gwl4p0-50y-dn2p0': '--', 'up2p0-gwl2p0': ':', 'up2p0-gwl2p0-50y-dn2p0': '--'} @entity_task(log, writes=['gcm_data']) def process_terrafirma_ukesm_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 Terrafirma-provided UKESM NetCDF files at the nearest gridpoint to the glacier, applies unit conversions and, optionally concatenates a historical / ramp-up period with a future scenario, and finally stores the result in OGGM's standard ``gcm_data`` format. Parameters ---------- gdir : oggm.GlacierDirectory Glacier directory to process. scenario : str Climate scenario identifier. Supported values include ``'hist'``, ``'up2p0'``, and derived scenarios such as ``'up2p0-gwl2p0'`` or ``'up2p0-gwl4p0'``. y0, y1 : int or str, optional Start and end years for subsetting the scenario data. The ramp-up or historical period is always fully included. 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 This is mainly used for Terrafirma-style stitched 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. """ 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 # We always want the full ramp-up / historical period y0 = None dir_path = ( 'https://cluster.klima.uni-bremen.de/~lschuster/' 'terrafirma_oggm_proj/terrafirma_data' ) time_coder = xr.coders.CFDatetimeCoder(use_cftime=True) if scenario in ('up2p0', 'hist'): fp = f'{dir_path}/pr_Amon_UKESM1-2-LL_esm-{scenario}_r1i1p1f1_gn.nc' ft = f'{dir_path}/tas_Amon_UKESM1-2-LL_esm-{scenario}_r1i1p1f1_gn.nc' with utils.get_lock(): ft = utils.file_downloader(ft) fp = utils.file_downloader(fp) y1_ramp_up_or_hist = y1 else: fp = f'{dir_path}/pr_Amon_UKESM1-2-LL_esm-up2p0_r1i1p1f1_gn.nc' ft = f'{dir_path}/tas_Amon_UKESM1-2-LL_esm-up2p0_r1i1p1f1_gn.nc' fp_scen = f'{dir_path}/pr_Amon_UKESM1-2-LL_esm-{scenario}_r1i1p1f1_gn.nc' ft_scen = f'{dir_path}/tas_Amon_UKESM1-2-LL_esm-{scenario}_r1i1p1f1_gn.nc' with utils.get_lock(): ft = utils.file_downloader(ft) fp = utils.file_downloader(fp) ft_scen = utils.file_downloader(ft_scen) fp_scen = utils.file_downloader(fp_scen) with xr.open_dataset(ft_scen, decode_times=time_coder) as tempds_scen: y0_scen = tempds_scen.time[0].values y1_ramp_up_or_hist = y0_scen - timedelta(days=30) if y0 is not None: y0 = str(y0) if y1 is not None: y1 = str(y1) if y1_ramp_up_or_hist is not None: y1_ramp_up_or_hist = str(y1_ramp_up_or_hist) with xr.open_dataset(ft, decode_times=time_coder) as tempds, \ xr.open_dataset(fp, decode_times=time_coder) as precipds: if (y0 is not None) or (y1_ramp_up_or_hist is not None): tempds = tempds.sel(time=slice(y0, y1_ramp_up_or_hist)) precipds = precipds.sel(time=slice(y0, y1_ramp_up_or_hist)) if tempds.lon.min() >= 0 and glon <= 0: glon += 360 temp = tempds.tas.sel(lat=glat, lon=glon, method='nearest') precip = precipds.pr.sel(lat=glat, lon=glon, method='nearest') if scenario not in ('up2p0', 'hist'): with xr.open_dataset(ft_scen, decode_times=time_coder) as tempds_scen, \ xr.open_dataset(fp_scen, decode_times=time_coder) as precipds_scen: if (y0 is not None) or (y1 is not None): tempds_scen = tempds_scen.sel(time=slice(y0, y1)) precipds_scen = precipds_scen.sel(time=slice(y0, y1)) temp_scen = tempds_scen.tas.sel(lat=glat, lon=glon, method='nearest') precip_scen = precipds_scen.pr.sel(lat=glat, lon=glon, method='nearest') temp = xr.merge([temp, temp_scen], join='outer', compat="no_conflicts").tas precip = xr.merge([precip, precip_scen], join='outer', compat="no_conflicts").pr temp.lon.values = temp.lon if temp.lon <= 180 else temp.lon - 360 precip.lon.values = precip.lon if precip.lon <= 180 else precip.lon - 360 assert 'kg m-2 s-1' in precip.units dimo = [cfg.DAYS_IN_MONTH[m - 1] for m in temp['time.month']] precip = precip * dimo * 24 * 3600 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, ) @entity_task(log) def run_with_terrafirma_data( gdir, scenarios=( 'hist', 'up2p0', 'up2p0-gwl2p0', 'up2p0-gwl2p0-50y-dn2p0', 'up2p0-gwl4p0', 'up2p0-gwl4p0-50y-dn2p0', ), ): """ 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 starting in 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`. Writes ------ model_run : NetCDF Glacier evolution outputs for each scenario, stored with a scenario-specific filesuffix. """ for scenario in scenarios: rid = f'_terrafirma_UKESM1-2-LL_esm_{scenario}_bc_1975_2014' tasks.run_from_climate_data( gdir, ys=1975, init_model_yr=1979, climate_filename='gcm_data', climate_input_filesuffix=rid, init_model_filesuffix='_spinup_historical', output_filesuffix=f'_future_run{rid}', )