Gridded Datasets II#

import numpy as np
import xarray as xr
import holoviews as hv
import geoviews as gv
import geoviews.feature as gf

from geoviews import opts
from cartopy import crs as ccrs

gv.extension('matplotlib', 'bokeh')

gv.output(size=200)

The main strength of HoloViews and its extensions (like GeoViews) is the ability to quickly explore complex datasets by declaring lower-dimensional views into a higher-dimensional space. In HoloViews we refer to the interface that allows you to do this as the conversion API. To begin with we will load a multi-dimensional dataset of surface temperatures for different “realizations” (modelling parameters) using XArray:

xr_ensembles = xr.open_dataset('../data/ensembles.nc')
xr_ensembles
<xarray.Dataset> Size: 9MB
Dimensions:                  (realization: 13, time: 6, latitude: 145,
                              longitude: 192, bnds: 2)
Coordinates:
  * realization              (realization) int32 52B 0 1 2 3 4 ... 9 10 11 12 13
  * time                     (time) datetime64[ns] 48B 2011-08-16T12:00:00 .....
    forecast_reference_time  (realization, time) datetime64[ns] 624B ...
  * latitude                 (latitude) float32 580B -90.0 -88.75 ... 88.75 90.0
  * longitude                (longitude) float32 768B 0.0 1.875 ... 356.2 358.1
Dimensions without coordinates: bnds
Data variables:
    surface_temperature      (realization, time, latitude, longitude) float32 9MB ...
    latitude_longitude       int32 4B ...
    time_bnds                (time, bnds) datetime64[ns] 96B ...
Attributes:
    source:       Data from Met Office Unified Model
    um_version:   7.6
    Conventions:  CF-1.5

As we saw in the Gridded Datasets I Tutorial we can easily wrap this xarray data structure as a HoloViews Dataset:

dataset = gv.Dataset(xr_ensembles, vdims='surface_temperature')
dataset
:Dataset   [realization,time,latitude,longitude]   (surface_temperature)

From the repr we can immediately see the list of key dimensions (time, realization, longitude and latitude) and the value dimension of the cube (surface_temperature). However, unlike most other HoloViews Elements, the Dataset Element does not display itself visually. This is because it can be n-dimensional and therefore does not have any specific straightforward visual representation on a 2D display. To view the cube, we first have to transform it into individually visualizable chunks. Before doing so, we will want to supply a custom value formatter for the time dimension so that it is readable by humans:

hv.Dimension.type_formatters[np.datetime64] = '%Y-%m-%d'

Conversions#

A HoloViews Dataset is a wrapper around a complex multi-dimensional datastructure which allows the user to convert their data into individually visualizable views, each usually of lower dimensionality. This is done by grouping the data by some dimension and then casting it to a specific Element type, which visualizes itself.

The dataset.to interface makes this especially easy. To use it, you supply the Element type that you want to view the data as and the key dimensions of that view and it will figure out the rest. Depending on the type of Element, you can specify one or more dimensions to be displayed. GeoViews provides a set of GeoElements that allow you to display geographic data on a cartographic projection, but you can use any Elements from HoloViews for non-geographic plots.

Recall that the cube we are working with has 4 coordinate dimensions (or key dimensions as they are known in HoloViews) – time, realization, longitude, and latitude. For our purposes, a geographic plot is defined as a plot that has longitude along the x axis and latitude along the y axis. To declare a two-dimensional geographic plot, we therefore simply request a gv.Image plot with longitude and latitude as key dimensions. There is one value dimension (vdim) available, surface_temperature, and any remaining key dimensions (time and realization in this case) are assigned to a HoloMap data structure by default. The resulting HoloMap gives you widgets automatically to allow you to explore the data across the two “remaining” key dimensions (those not mapped onto axes of the image):

geo_dims = ['longitude', 'latitude']
(dataset.to(gv.Image, geo_dims) * gf.coastline)[::5, ::5]

In this way we can visualize the geographic data in a number of ways, currently either as an Image (as above) or as LineContours, FilledContours, or Points:

layout = hv.Layout([dataset.to(el, geo_dims)[::10, ::5] * gf.coastline
                    for el in [gv.FilledContours, gv.LineContours, gv.Points]]).cols(1)

layout.opts(opts.Points(color='surface_temperature', cmap='jet'))

Note that by default the conversion interface will automatically expand all the individual Elements, which can take some time if the data is very large. Instead we can also request the objects to be expanded dynamically using the dynamic keyword:

dataset.to(gv.Image, geo_dims, dynamic=True) * gf.coastline

Using dynamic mode means that the data for each frame is only extracted when you’re actually viewing that part of the data, which can have huge benefits in terms of speed and memory consumption. However, it relies on having a running Python process to render and serve each image, and so it cannot be used when generating static HTML output such as for the GeoViews web site.

Irregularly sampled data#

Often gridded datasets are not regularly sampled, instead providing irregularly sampled multi-dimensional coordinates. Such datasets can be easily visualized using the QuadMesh element.

xrds = xr.tutorial.open_dataset('rasm').load()
qmesh = gv.Dataset(xrds.Tair).to(gv.QuadMesh, ['xc', 'yc'], dynamic=True).redim.range(Tair=(-30, 30))
qmesh.opts(colorbar=True, cmap='RdBu_r', projection=ccrs.Robinson()) * gf.coastline

Non-geographic views#

So far we have focused entirely on geographic views of the data, plotting the data on a projection. However the conversion interface is completely general, allowing us to slice and dice the data in any way we like. For these views we will switch to the bokeh plotting extension:

hv.output(backend='bokeh')

The simplest example of this capability is simply a view showing the temperature over time for each realization, longitude, and latitude coordinate:

curves = dataset.to(hv.Curve, 'time', dynamic=True).overlay('realization')

curves.opts(
    opts.Curve(xrotation=25, width=600, height=400, framewise=True),
    opts.NdOverlay(legend_position='right', toolbar='right'))

Note that the longitude slider will have no effect, if latitude is -90 or +90, since there is only one data point for the North or South poles (regardless of the declared longitude). Here the .overlay gives a different curve for each realization; without it all realization values would be pooled together.

We can also make non-geographic 2D plots, for instance as a HeatMap over time and realization, at a specified longitude and latitude:

dataset.to(hv.HeatMap, ['realization', 'time'], dynamic=True).opts(width=600, colorbar=True, tools=['hover'])
---------------------------------------------------------------------------
UFuncTypeError                            Traceback (most recent call last)
File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/IPython/core/formatters.py:1036, in MimeBundleFormatter.__call__(self, obj, include, exclude)
   1033     method = get_real_method(obj, self.print_method)
   1035     if method is not None:
-> 1036         return method(include=include, exclude=exclude)
   1037     return None
   1038 else:

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/holoviews/core/dimension.py:1439, in Dimensioned._repr_mimebundle_(self, include, exclude)
   1432 def _repr_mimebundle_(self, include=None, exclude=None):
   1433     """Resolves the class hierarchy for the class rendering the
   1434     object using any display hooks registered on Store.display
   1435     hooks.  The output of all registered display_hooks is then
   1436     combined and returned.
   1437 
   1438     """
-> 1439     return Store.render(self)

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/holoviews/core/options.py:1500, in Store.render(cls, obj)
   1498 data, metadata = {}, {}
   1499 for hook in hooks:
-> 1500     ret = hook(obj)
   1501     if ret is None:
   1502         continue

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/holoviews/ipython/display_hooks.py:328, in pprint_display(obj)
    326 if not ip.display_formatter.formatters["text/plain"].pprint:
    327     return None
--> 328 return display(obj, raw_output=True)

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/holoviews/ipython/display_hooks.py:302, in display(obj, raw_output, **kwargs)
    300 elif isinstance(obj, (HoloMap, DynamicMap)):
    301     with option_state(obj):
--> 302         output = map_display(obj)
    303 elif isinstance(obj, Plot):
    304     output = render(obj)

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/holoviews/ipython/display_hooks.py:191, in display_hook.<locals>.wrapped(element)
    189 try:
    190     max_frames = OutputSettings.options["max_frames"]
--> 191     mimebundle = fn(element, max_frames=max_frames)
    192     if mimebundle is None:
    193         return {}, {}

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/holoviews/ipython/display_hooks.py:244, in map_display(vmap, max_frames)
    241     max_frame_warning(max_frames)
    242     return None
--> 244 return render(vmap)

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/holoviews/ipython/display_hooks.py:60, in render(obj, **kwargs)
     57 if renderer.fig == "pdf":
     58     renderer = renderer.instance(fig="png")
---> 60 return renderer.components(obj, **kwargs)

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/holoviews/plotting/renderer.py:449, in Renderer.components(self, obj, fmt, comm, **kwargs)
    447 embed = not (dynamic or streams or self.widget_mode == "live") or config.embed
    448 if embed or config.comms == "default":
--> 449     return self._render_panel(plot, embed, comm)
    450 return self._render_ipywidget(plot)

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/holoviews/plotting/renderer.py:456, in Renderer._render_panel(self, plot, embed, comm)
    454 doc = Document()
    455 with config.set(embed=embed):
--> 456     model = plot.layout._render_model(doc, comm)
    457 if embed:
    458     return render_model(model, comm)

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/panel/viewable.py:769, in Viewable._render_model(self, doc, comm)
    767 if comm is None:
    768     comm = state._comm_manager.get_server_comm()
--> 769 model = self.get_root(doc, comm)
    771 if self._design and self._design.theme.bokeh_theme:
    772     doc.theme = self._design.theme.bokeh_theme

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/panel/layout/base.py:329, in Panel.get_root(self, doc, comm, preprocess)
    325 def get_root(
    326     self, doc: Document | None = None, comm: Comm | None = None,
    327     preprocess: bool = True
    328 ) -> Model:
--> 329     root = super().get_root(doc, comm, preprocess)
    330     # ALERT: Find a better way to handle this
    331     if hasattr(root, 'styles') and 'overflow-x' in root.styles:

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/panel/viewable.py:699, in Renderable.get_root(self, doc, comm, preprocess)
    697 wrapper = self._design._wrapper(self)
    698 if wrapper is self:
--> 699     root = self._get_model(doc, comm=comm)
    700     if preprocess:
    701         self._preprocess(root)

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/panel/layout/base.py:313, in Panel._get_model(self, doc, root, parent, comm)
    311 root = root or model
    312 self._models[root.ref['id']] = (model, parent)
--> 313 objects, _ = self._get_objects(model, [], doc, root, comm)
    314 props = self._get_properties(doc)
    315 props[self._property_mapping['objects']] = objects

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/panel/layout/base.py:295, in Panel._get_objects(self, model, old_objects, doc, root, comm)
    293 else:
    294     try:
--> 295         child = pane._get_model(doc, root, model, comm)
    296     except RerenderError as e:
    297         if e.layout is not None and e.layout is not self:

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/panel/pane/holoviews.py:488, in HoloViews._get_model(self, doc, root, parent, comm)
    486     plot = self.object
    487 else:
--> 488     plot = self._render(doc, comm, root)
    490 plot.pane = self
    491 backend = plot.renderer.backend

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/panel/pane/holoviews.py:582, in HoloViews._render(self, doc, comm, root)
    579     if comm:
    580         kwargs['comm'] = comm
--> 582 return renderer.get_plot(self.object, **kwargs)

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/holoviews/plotting/bokeh/renderer.py:75, in BokehRenderer.get_plot(self_or_cls, obj, doc, renderer, **kwargs)
     68 @bothmethod
     69 def get_plot(self_or_cls, obj, doc=None, renderer=None, **kwargs):
     70     """Given a HoloViews Viewable return a corresponding plot instance.
     71     Allows supplying a document attach the plot to, useful when
     72     combining the bokeh model with another plot.
     73 
     74     """
---> 75     plot = super().get_plot(obj, doc, renderer, **kwargs)
     76     if plot.document is None:
     77         plot.document = Document() if self_or_cls.notebook_context else curdoc()

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/holoviews/plotting/renderer.py:290, in Renderer.get_plot(self_or_cls, obj, doc, renderer, comm, **kwargs)
    286     defaults = [kd.default for kd in plot.dimensions]
    287     init_key = tuple(
    288         v if d is None else d for v, d in zip(plot.keys[0], defaults, strict=None)
    289     )
--> 290     plot.update(init_key)
    291 else:
    292     plot = obj

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/holoviews/plotting/plot.py:1064, in DimensionedPlot.update(self, key)
   1062 def update(self, key):
   1063     if len(self) == 1 and key in (0, self.keys[0]) and not self.drawn:
-> 1064         return self.initialize_plot()
   1065     item = self.__getitem__(key)
   1066     self.traverse(lambda x: setattr(x, "_updated", True))

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/holoviews/plotting/bokeh/element.py:2490, in ElementPlot.initialize_plot(self, ranges, plot, plots, source)
   2488 if self.autorange:
   2489     self._setup_autorange()
-> 2490 self._init_glyphs(plot, element, ranges, source)
   2491 if not self.overlaid:
   2492     self._update_plot(key, plot, style_element)

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/holoviews/plotting/bokeh/heatmap.py:341, in HeatMapPlot._init_glyphs(self, plot, element, ranges, source)
    340 def _init_glyphs(self, plot, element, ranges, source):
--> 341     super()._init_glyphs(plot, element, ranges, source)
    342     self._draw_markers(plot, element, self.xmarks, axis="x")
    343     self._draw_markers(plot, element, self.ymarks, axis="y")

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/holoviews/plotting/bokeh/element.py:2394, in ElementPlot._init_glyphs(self, plot, element, ranges, source)
   2392 else:
   2393     style = self.style[self.cyclic_index]
-> 2394     data, mapping, style = self.get_data(element, ranges, style)
   2395     current_id = element._plot_id
   2397 with abbreviated_exception():

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/holoviews/plotting/bokeh/heatmap.py:193, in HeatMapPlot.get_data(self, element, ranges, style)
    189         ys = element.dimension_values(y, expanded=False)
    190         ydiff = np.diff(ys)
    191         y_index = (
    192             ranges[y]["data"][0]
--> 193             if (not ydiff.size or np.allclose(ydiff, ydiff[0]))
    194             else None
    195         )
    197 self._is_contiguous_gridded = is_gridded and x_index is not None and y_index is not None
    198 if self.is_contiguous_gridded:

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/numpy/_core/numeric.py:2376, in allclose(a, b, rtol, atol, equal_nan)
   2290 @array_function_dispatch(_allclose_dispatcher)
   2291 def allclose(a, b, rtol=1.e-5, atol=1.e-8, equal_nan=False):
   2292     """
   2293     Returns True if two arrays are element-wise equal within a tolerance.
   2294 
   (...)   2374 
   2375     """
-> 2376     res = all(isclose(a, b, rtol=rtol, atol=atol, equal_nan=equal_nan))
   2377     return builtins.bool(res)

File ~/work/geoviews/geoviews/.pixi/envs/docs/lib/python3.11/site-packages/numpy/_core/numeric.py:2507, in isclose(a, b, rtol, atol, equal_nan)
   2503         print(err_msg)
   2505 with errstate(invalid='ignore'):
-> 2507     result = (less_equal(abs(x - y), atol + rtol * abs(y))
   2508               & isfinite(y)
   2509               | (x == y))
   2510     if equal_nan:
   2511         result |= isnan(x) & isnan(y)

UFuncTypeError: ufunc 'add' cannot use operands with types dtype('float64') and dtype('<m8[ns]')
:DynamicMap   [latitude,longitude]
   :HeatMap   [realization,time]   (surface_temperature)

Lower-dimensional views#

So far all the conversions shown have incorporated each of the available coordinate dimensions explicitly. However, often times we want to see the spread of values along one or more dimensions, pooling all the other dimensions together. A simple example of this is a box plot where we might want to see the spread of surface_temperature on each day, pooled across all latitude and longitude coordinates. To pool across particular dimensions, we can explicitly declare the “map” dimensions, which are the key dimensions of the HoloMap container rather than those of the Elements contained in the HoloMap. By explicitly declaring no dimensions to groupby, we can tell the conversion interface to pool across all dimensions except the particular key dimension(s) supplied, in this case the 'time' (plot A) and 'realization' (plot B):

hv.output(backend='matplotlib')

hv.Layout([dataset.to(hv.Violin, d, groupby=[], datatype=['dataframe']).opts(xrotation=25)
           for d in ['time', 'realization']])

Reducing the data#

So far all the examples we have seen have displayed all the data in some way or another. Another way to explore a dataset is to explicitly reduce the dimensionality or select subregions of a dataset. There are two main ways to do this—either we explicitly select a subset of the data, or we collapse a dimension using an aggregation function, e.g. by computing a mean along a particular dimension.

Selecting slices#

Using the select method we can easily select ranges of coordinates in the dataset. Unfortunately, the select method does not currently know that latitude and longitude are cyclic, so instead we have to select regions at both ends of the prime meridian (0$^\circ$ longitude) and overlay them. In this way we can stitch together multiple cubes or xarrays or simply view specific subregions:

northern = dataset.select(latitude=slice(25, 75))
(northern.select(longitude=slice(260, 305)).to(gv.Image, geo_dims) *
 northern.select(longitude=slice(330, 362)).to(gv.Image, geo_dims) *
 gf.coastline)[::5, ::5]

Selecting a particular coordinate#

To examine one particular coordinate, we can select it, cast the data to Curves, reindex the data to drop the now-constant latitude and longitude dimensions, and overlay the remaining ‘realization’ dimension:

hv.output(backend='bokeh')

curves = dataset.select(latitude=0, longitude=0).to(hv.Curve, ['time']).reindex().overlay()

curves.opts(width=600, height=400, legend_position='right', toolbar='above')

As you can see, with GeoViews and HoloViews it is very simple to select precisely which aspects of complex, multidimensional datasets that you want to focus on.

This web page was generated from a Jupyter notebook and not all interactivity will work on this website. Right click to download and run locally for full Python-backed interactivity.

Right click to download this notebook from GitHub.