Plotting glacier data

Plotting with xarray

Let us open the demo data again:

In [1]: import xarray as xr

In [2]: import hyoga.demo

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

Plotting is already quite convenient using xarray:

In [4]: ds.thk.plot.imshow()
Out[4]: <matplotlib.image.AxesImage at 0x7f61f4584ed0>

In [5]: plt.gca().set_aspect('equal')  # needed to avoid distortion

Plotting with hyoga

To make things even easier, hyoga provides wrappers around xarray and matplotlib methods to produce oft-used ice sheet model plots with a more practical default style.

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

In [7]: plt.gca().set_aspect('equal')

Now let us make our first composite plot. The previously defined ice mask allows us to plot an ice margin contour on top of the bedrock topography:

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

In [9]: ds.hyoga.plot.ice_margin(facecolor='tab:blue')
[<matplotlib.contour.QuadContourSet at 0x7f61f4528210>,
 <matplotlib.contour.QuadContourSet at 0x7f61f44b2450>]

In [10]: plt.gca().set_aspect('equal')


Hyoga makes choices regarding style (such as the grey colormap for topography and the half-transparent background above) that do not always correspond to the matplotlib defaults. However, these choices can always be overriden using the matplotlib keyword arguments.

As may be noticed in the above code, hyoga plot methods also do the job to look up required variables (by their standard names) and to infer missing variables when that is possible. For instance in the above example, the 'surface_altitude' was missing from the data and thus computed from 'bedrock_altitude' and 'land_ice_thickness'. Such computations are done on-the-fly though, and the results are not stored:

In [11]: for name, var in ds.items():
   ....:    print(name, var.attrs.get('standard_name', None))
mapping None
pism_config None
run_stats None
time_bounds None
topg bedrock_altitude
thk land_ice_thickness
uvelbase land_ice_basal_x_velocity
vvelbase land_ice_basal_y_velocity
uvelsurf land_ice_surface_x_velocity
vvelsurf land_ice_surface_y_velocity

Velocity maps are automatically log-scaled, but the limits can be customized:

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

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

In [14]: ds.hyoga.plot.ice_margin(edgecolor='0.25')
Out[14]: <matplotlib.contour.QuadContourSet at 0x7f61f3bde750>

In [15]: plt.gca().set_aspect('equal')

Plotting with cartopy

For enhanced visuals, hyoga plots can be georeferenced and combined with Natural Earth vector data shipped with cartopy.

In [16]: import matplotlib.pyplot as plt

In [17]: import as ccrs

In [18]: import cartopy.feature as cfeature

# initialize subplot with UTM projection
In [19]: ax = plt.subplot(projection=ccrs.UTM(32))

# add coastline and rivers
In [20]: ax.coastlines(edgecolor='0.25', linewidth=0.5)
Out[20]: <cartopy.mpl.feature_artist.FeatureArtist at 0x7f6207399790>

In [21]: ax.add_feature(
   ....:    cfeature.NaturalEarthFeature(
   ....:       category='physical', name='rivers_lake_centerlines', scale='10m'),
   ....:    edgecolor='0.25', facecolor='none', linewidth=0.5, zorder=0)
Out[21]: <cartopy.mpl.feature_artist.FeatureArtist at 0x7f61f4237a50>

# plot model output
In [22]: ds.hyoga.plot.bedrock_altitude(vmin=0, vmax=4500)
Out[22]: <matplotlib.image.AxesImage at 0x7f61f78a3410>

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

In [24]: ds.hyoga.plot.ice_margin()
Out[24]: <cartopy.mpl.contour.GeoContourSet at 0x7f61f429dad0>

More plotting methods are available. Please take a look at the Examples gallery gallery.