Compute area and thickness changes on a grid (experimental!)#
from oggm import cfg
from oggm import tasks, utils, workflow, graphics
import xarray as xr
import matplotlib.pyplot as plt
This notebook proposes a method for redistributing glacier ice that has been simulated along the flowline after a glacier retreat simulation. Extrapolating the glacier ice onto a map involves certain assumptions and trade-offs. Depending on the purpose, different choices may be preferred. For example, higher resolution may be desired for visualization compared to using the output for a hydrological model. It is possible to add different options to the final function to allow users to select the option that best suits their needs.
This notebook demonstrates the redistribution process using a single glacier. Its purpose is to initiate further discussion before incorporating it into the main OGGM code base (currently in the sandbox).
Pick a glacier#
# Initialize OGGM and set up the default run parameters
cfg.initialize(logging_level='WARNING')
# Local working directory (where OGGM will write its output)
# WORKING_DIR = utils.gettempdir('OGGM_distr4')
cfg.PATHS['working_dir'] = utils.get_temp_dir('OGGM_distributed', reset=True)
2023-05-22 07:59:44: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2023-05-22 07:59:44: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2023-05-22 07:59:44: oggm.cfg: Multiprocessing: using all available processors (N=2)
rgi_ids = ['RGI60-11.01450'] # This is Aletsch
base_url = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/2023.1/elev_bands/W5E5'
gdir = workflow.init_glacier_directories(rgi_ids, prepro_base_url=base_url, from_prepro_level=3, prepro_border=80)[0]
2023-05-22 07:59:46: oggm.workflow: init_glacier_directories from prepro level 3 on 1 glaciers.
2023-05-22 07:59:46: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers
Experiment: currently only a random warming simulation#
There is a little bit more work needed on OGGM for more complex workflows including past and future simulations (not much! coming soon).
# Do a random run with a bit of warming
tasks.run_random_climate(gdir, nyears=100,
y0=2000, # Climate of 1985-2015
seed=1, # Change for another randomness
temperature_bias=1, # casual warming - change for other scenarios
store_fl_diagnostics=True, # important! This will be needed for the redistribution
output_filesuffix='_rdn_1', # optional - here I just want to make things explicit as to which run we are using afterwards
);
Redistribute: preprocessing#
The required tasks can be found in the distribute_2d
module of the sandbox:
from oggm.sandbox import distribute_2d
# This is to add a new topography to the file (smoothed differently)
distribute_2d.add_smoothed_glacier_topo(gdir)
# This is to get the bed map at the start of the simulation
tasks.distribute_thickness_per_altitude(gdir)
# This is to prepare the glacier directory for the interpolation (needs to be done only once)
distribute_2d.assign_points_to_band(gdir)
Let’s have a look at what we just did:
with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
ds = ds.load()
# Inititial glacier thickness
f, ax = plt.subplots()
ds.distributed_thickness.plot(ax=ax);
ax.axis('equal');

# Which points belongs to which band, and then within one band which are the first to melt
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
ds.band_index.plot(ax=ax1);
ds.rank_per_band.plot(ax=ax2);
ax1.axis('equal'); ax2.axis('equal'); plt.tight_layout();

Redistribute simulation#
The tasks above need to be run only once. The next one however should be done for each simulation:
distribute_2d.distribute_thickness_from_simulation(gdir, input_filesuffix='_rdn_1')
2023-05-22 07:59:50: oggm.sandbox.distribute_2d: TypeError occurred during task distribute_thickness_from_simulation on RGI60-11.01450: distribute_thickness_from_simulation() got an unexpected keyword argument 'input_filesuffix'
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[10], line 1
----> 1 distribute_2d.distribute_thickness_from_simulation(gdir, input_filesuffix='_rdn_1')
File /usr/local/pyenv/versions/3.10.11/lib/python3.10/site-packages/oggm/utils/_workflow.py:491, in entity_task.__call__.<locals>._entity_task(gdir, reset, print_log, return_value, continue_on_error, add_to_log_file, **kwargs)
489 signal.alarm(cfg.PARAMS['task_timeout'])
490 ex_t = time.time()
--> 491 out = task_func(gdir, **kwargs)
492 ex_t = time.time() - ex_t
493 if cfg.PARAMS['task_timeout'] > 0:
TypeError: distribute_thickness_from_simulation() got an unexpected keyword argument 'input_filesuffix'
Plot#
Let’s have a look!
with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
ds = ds.load()
ds.time
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[12], line 1
----> 1 ds.time
File /usr/local/pyenv/versions/3.10.11/lib/python3.10/site-packages/xarray/core/common.py:276, in AttrAccessMixin.__getattr__(self, name)
274 with suppress(KeyError):
275 return source[name]
--> 276 raise AttributeError(
277 f"{type(self).__name__!r} object has no attribute {name!r}"
278 )
AttributeError: 'Dataset' object has no attribute 'time'
f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(14, 4))
ds.simulation_distributed_thickness_rdn_1.sel(time=0).plot(ax=ax1, vmax=400);
ds.simulation_distributed_thickness_rdn_1.sel(time=40).plot(ax=ax2, vmax=400);
ds.simulation_distributed_thickness_rdn_1.sel(time=80).plot(ax=ax3, vmax=400);
ax1.axis('equal'); ax2.axis('equal'); plt.tight_layout();
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[13], line 2
1 f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(14, 4))
----> 2 ds.simulation_distributed_thickness_rdn_1.sel(time=0).plot(ax=ax1, vmax=400);
3 ds.simulation_distributed_thickness_rdn_1.sel(time=40).plot(ax=ax2, vmax=400);
4 ds.simulation_distributed_thickness_rdn_1.sel(time=80).plot(ax=ax3, vmax=400);
File /usr/local/pyenv/versions/3.10.11/lib/python3.10/site-packages/xarray/core/common.py:276, in AttrAccessMixin.__getattr__(self, name)
274 with suppress(KeyError):
275 return source[name]
--> 276 raise AttributeError(
277 f"{type(self).__name__!r} object has no attribute {name!r}"
278 )
AttributeError: 'Dataset' object has no attribute 'simulation_distributed_thickness_rdn_1'

Animation!#
from matplotlib import animation
from IPython.display import HTML, display
# Get a handle on the figure and the axes
fig, ax = plt.subplots()
thk = ds['simulation_distributed_thickness_rdn_1']
# Plot the initial frame.
cax = thk.isel(time=0).plot(ax=ax,
add_colorbar=True,
cmap='viridis',
vmin=0, vmax=350,
cbar_kwargs={
'extend':'neither'
}
)
ax.axis('equal')
def animate(frame):
ax.set_title(f'Year {int(frame)}')
cax.set_array(thk.values[frame, :].flatten())
ani_glacier = animation.FuncAnimation(fig, animate, frames=len(thk.time), interval=200);
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
File /usr/local/pyenv/versions/3.10.11/lib/python3.10/site-packages/xarray/core/dataset.py:1353, in Dataset._construct_dataarray(self, name)
1352 try:
-> 1353 variable = self._variables[name]
1354 except KeyError:
KeyError: 'simulation_distributed_thickness_rdn_1'
During handling of the above exception, another exception occurred:
KeyError Traceback (most recent call last)
Cell In[14], line 7
4 # Get a handle on the figure and the axes
5 fig, ax = plt.subplots()
----> 7 thk = ds['simulation_distributed_thickness_rdn_1']
9 # Plot the initial frame.
10 cax = thk.isel(time=0).plot(ax=ax,
11 add_colorbar=True,
12 cmap='viridis',
(...)
16 }
17 )
File /usr/local/pyenv/versions/3.10.11/lib/python3.10/site-packages/xarray/core/dataset.py:1444, in Dataset.__getitem__(self, key)
1442 return self.isel(**key)
1443 if utils.hashable(key):
-> 1444 return self._construct_dataarray(key)
1445 if utils.iterable_of_hashable(key):
1446 return self._copy_listed(key)
File /usr/local/pyenv/versions/3.10.11/lib/python3.10/site-packages/xarray/core/dataset.py:1355, in Dataset._construct_dataarray(self, name)
1353 variable = self._variables[name]
1354 except KeyError:
-> 1355 _, name, variable = _get_virtual_variable(self._variables, name, self.dims)
1357 needed_dims = set(variable.dims)
1359 coords: dict[Hashable, Variable] = {}
File /usr/local/pyenv/versions/3.10.11/lib/python3.10/site-packages/xarray/core/dataset.py:185, in _get_virtual_variable(variables, key, dim_sizes)
183 split_key = key.split(".", 1)
184 if len(split_key) != 2:
--> 185 raise KeyError(key)
187 ref_name, var_name = split_key
188 ref_var = variables[ref_name]
KeyError: 'simulation_distributed_thickness_rdn_1'

HTML(ani_glacier.to_jshtml())
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[15], line 1
----> 1 HTML(ani_glacier.to_jshtml())
NameError: name 'ani_glacier' is not defined
# Write to mp4?
# FFwriter = animation.FFMpegWriter(fps=10)
# ani2.save('animation.mp4', writer=FFwriter)