Reading model output

Opening output files

Let us open demo data. This will download a model output data file from the web, and store it into a ~/.cache/hyoga directory so that it can be reused the next time.

In [1]: import xarray as xr

In [2]: import hyoga.demo

In [3]: ds = xr.open_dataset(hyoga.demo.get('pism.alps.out.2d.nc'))

Alternatively, hyoga.open provides thin wrappers around xarray functions to open a single-file dataset, a multi-file dataset, and a single time slice within a multi-file dataset.

hyoga.open.dataset

Open a single-file model output dataset.

hyoga.open.mfdataset

Open a multi-file model output dataset.

hyoga.open.subdataset

Open a single file in a multi-file dataset.

Warning

These functions convert model time to an age coordinate in ka before the present. This functionality is not always useful and will be reworked in future versions using Climate and Forecast (CF) time conventions. If you do not need an age coordinate, please use regular xarray.open_dataset() for the time being.

Selecting variables

Please refer to the excellent xarray docs to explore the many other ways of indexing and subsetting data. As long as hyoga (or any submodule) has been imported, new functionality will be available in a so-called “dataset accessor” attribute ds.hyoga:

In [4]: ds.hyoga
Out[4]: <hyoga.hyoga.HyogaDataset at 0x7f61f45ff7d0>

Note

I am thinking about renaming ds.hyoga to ds.ice starting from v0.2.0. The name may be used in other projects though, I am not sure. In any case ds.hyoga would remain backward-compatible for a while. If you have an opinion about that, feel free to share it on Github (#13).

In particular, hyoga never accesses model variables by their “short names”. While thk, for instance refers to ice thickness in PISM, it may refer to a different quantity, or to nothing at all, in another ice-sheet model. This where CF standard names come into play. To access a variable by standard name you may use:

In [5]: var = ds.hyoga.getvar('land_ice_thickness');

In [6]: var.max()
Out[6]: 
<xarray.DataArray 'thk' ()>
array(2571.5942, dtype=float32)
Coordinates:
    time     object ...

In [7]: var.units
Out[7]: 'm'

If a particular variable is missing, hyoga will additionally try to reconstruct it from others, such as the sum of bedrock altitude and ice thickness for surface altitude, or the norm of velocity components for its magnitude.

In [8]: ds.hyoga.getvar('surface_altitude');

In [9]: ds.hyoga.getvar('magnitude_of_land_ice_surface_velocity');

This mechanism can be disabled using infer=False. Because surface altitude is not actually present in the demo dataset, the following will raise an exception:

In [10]: ds.hyoga.getvar('surface_altitude', infer=False)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-10-ffd77ed650e4> in <module>
----> 1 ds.hyoga.getvar('surface_altitude', infer=False)

~/checkouts/readthedocs.org/user_builds/hyoga/conda/stable/lib/python3.7/site-packages/hyoga/hyoga.py in getvar(self, standard_name, infer, directions)
    313         # really nothing worked, give up
    314         raise ValueError(
--> 315             f"No variable found with standard name {standard_name}.")
    316 
    317     def interp(self, datasource, ax=None, sigma=None, threshold=None):

ValueError: No variable found with standard name surface_altitude.

Because CF standard names for land ice variables are relatively recent, older ice sheet models may not include them in output metadata. For PISM, a mechanism has been implemented to fill (some of) these missing standard names during initialization.

Note

While hyoga has only been tested with PISM so far, it is my hope that it will become compatible with some other glacier and ice sheet models in the future. If you want to make your glacier model compatible with hyoga, please consider implementing CF standard names.

Adding new variables

New variables can be added using using xarray’s dictionary interface or methods such as xarray.Dataset.assign(). Besides, hyoga provides a dataset method to assign new variables by their standard name.

In [11]: bedrock = ds.hyoga.getvar('bedrock_altitude')

In [12]: thickness = ds.hyoga.getvar('land_ice_thickness')

In [13]: surface = bedrock + thickness

In [14]: new = ds.hyoga.assign(surface_altitude=surface)

This returns a new dataset including the surface altitude variable. Some control on the variable (short) name can be achieved by preceding the assign call with xarray.DataArray.rename().

In [15]: ds = ds.hyoga.assign(surface_altitude=surface.rename('usurf'))

In [16]: 'usurf' in ds
Out[16]: True

However, this only works if the data does not already contain a variable with the standard name surface_altitude. In that case, that variable’s data is quietly replaced, and the variable is not renamed.

In [17]: ds = ds.hyoga.assign(surface_altitude=surface.rename('newsurf'))

In [18]: 'newsurf' in ds
Out[18]: False