Masking and interpolation

Ice thickness threshold

Let’s get started with the usual imports and the demo data.

In [1]: import matplotlib.pyplot as plt

In [2]: import as ccrs

In [3]: import xarray as xr

In [4]: import hyoga.demo

In [5]: ds = xr.open_dataset(hyoga.demo.get(''))

Hyoga’s plot methods use an ice mask to determine which grid cells are glacierized and which are not. According to CF conventions, this is defined by the standard variable land_ice_area_fraction. There are several ways to affect the ice mask. The easiest way is to use the (currently single) parametre in hyoga.config.

In [6]: hyoga.config.glacier_masking_point
Out[6]: 1.0

If the land_ice_area_fraction variable is missing from the dataset, hyoga falls back to compute and ice mask from land_ice_thickness, using this parametre as an ice thickness threshold. The default value is 1 (metre). For PISM output files, a non-zero threshold may be advisable in case winter output files contain a thin cover of “seasonal ice” outside the glacier margin, as is the case in the demo files.

In [7]: plt.subplot(projection=ccrs.UTM(32))
Out[7]: <GeoAxesSubplot:>

In [8]: ds.hyoga.plot.bedrock_altitude(vmin=0, vmax=4500)
Out[8]: <matplotlib.image.AxesImage at 0x7f61f3abd050>

In [9]: for i, value in enumerate([0.1, 1, 500]):
   ...:     hyoga.config.glacier_masking_point = value
   ...:     ds.hyoga.plot.ice_margin(edgecolor=f'C{i}', linewidths=1)

In [10]: hyoga.config.glacier_masking_point = 1  # restore default

Custom glacier mask

For more control, on can set the land_ice_area_fraction variable using assign_icemask(). Suppose that we define glaciers as grid cells filled with ice at least a metre thick, and moving at least ten metres per year:

In [11]: plt.subplot(projection=ccrs.UTM(32))
Out[11]: <GeoAxesSubplot:>

In [12]: ds = ds.hyoga.assign_icemask(
   ....:     (ds.hyoga.getvar('land_ice_thickness') > 1) &
   ....:     (ds.hyoga.getvar('magnitude_of_land_ice_surface_velocity') > 10))

In [13]: ds.hyoga.plot.bedrock_altitude(vmin=0, vmax=4500)
Out[13]: <matplotlib.image.AxesImage at 0x7f61f714a090>

In [14]: ds.hyoga.plot.ice_margin(facecolor='tab:blue')
[<cartopy.mpl.contour.GeoContourSet at 0x7f61f4636590>,
 <cartopy.mpl.contour.GeoContourSet at 0x7f61f4636650>]

Note that the assign_icemask() method edits (or add) a land_ice_area_fraction variable without affecting the rest of the dataset. Such lossless masking is should be enough for internal use within Hyoga. This mask looks a bit strange, so let us get rid of it before we move on:

In [15]: ds = ds.drop_vars(ds.hyoga.getvar('land_ice_area_fraction').name)

However in some situations, a lossy (destructive) ice mask may be more useful. This includes exporting data to a compressed netCDF file for the web, where having homogeneous values outside the glacier mask can greatly reduce file size. This can be achieved with hyoga’s where(), where_icemask(), and where_thicker(). These methods behave like xarray’s where(): they replace data values with np.nan outside the where condition. However, they are meant to only affect “glacier variables” (currently any variable whose standard name does not start with bedrock_altitude).

Interpolation to another topography

For enhanced visuals, hyoga supports interpolating and remasking model output onto a different (usually higher-resolution) topography. This produces an image which deviates somewhat from the model results, but with detail potentially improving readability.

A necessary first step in many cases is to compute the bedrock deformation due to isostatic adjustment. If the dataset was produced by a model including bedrock deformation, there will be a systematic offset between the model bedrock topography and the new topography onto which data are interpolated.


Assigning isostasy is not strictly necessary for interpolation to succeed, but omitting it would produce errors in the new ice mask. Of course, it can be safely omitted if the model run does not include isostasy.

In this case, we compute the bedrock deformation by comparing bedrock altitude in the dataset with bedrock altitude in the initial state:

In [16]: ds = ds.hyoga.assign_isostasy(hyoga.demo.get(''))

The method assign_isostasy() assigns a new variable (standard name bedrock_altitude_change_due_to_isostatic_adjustment). Next we run interp() which interpolates all variables, and recalculates an ice mask based on the new topographies, corrected for bedrock depression in this case. This uses yet another demo file, which contains high-resolution topographic data over a small part of the model domain.

In [17]: ds = ds.hyoga.interp(hyoga.demo.get(''))

The new dataset can be plotted in the same way as any other hyoga dataset, only with a much higher resolution.

In [18]: plt.subplot(projection=ccrs.UTM(32))
Out[18]: <GeoAxesSubplot:>

In [19]: ds.hyoga.plot.bedrock_altitude(vmin=0, vmax=4500)
Out[19]: <matplotlib.image.AxesImage at 0x7f61f789f9d0>

In [20]: ds.hyoga.plot.surface_velocity(vmin=1e1, vmax=1e3)
Out[20]: <matplotlib.image.AxesImage at 0x7f61f77f9a10>

In [21]: ds.hyoga.plot.surface_altitude_contours()
[<cartopy.mpl.contour.GeoContourSet at 0x7f61f78a0990>,
 <cartopy.mpl.contour.GeoContourSet at 0x7f61f78a0e10>]

In [22]: ds.hyoga.plot.ice_margin(edgecolor='0.25')
Out[22]: <cartopy.mpl.contour.GeoContourSet at 0x7f61f4333550>