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)