{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gridded Datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import numpy as np\n",
    "import holoviews as hv\n",
    "from holoviews import opts\n",
    "hv.extension('matplotlib')\n",
    "\n",
    "opts.defaults(opts.Scatter3D(color='Value', cmap='fire', edgecolor='black', s=50))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the [Tabular Data](./08-Tabular_Datasets.ipynb) guide we covered how to work with columnar data in HoloViews. Apart from tabular or column based data there is another data format that is particularly common in the science and engineering contexts, namely multi-dimensional arrays. The gridded data interfaces allow working with grid-based datasets directly.\n",
    "\n",
    "Grid-based datasets have two types of dimensions:\n",
    "\n",
    "* they have coordinate or key dimensions, which describe the sampling of each dimension in the value arrays\n",
    "* they have value dimensions which describe the quantity of the multi-dimensional value arrays\n",
    "\n",
    "There are many different types of gridded datasets, which each approximate or measure a surface or space at discretely specified coordinates.  In HoloViews, gridded datasets are typically one of three possible types: Regular rectilinear grids, irregular rectilinear grids, and curvilinear grids.  Regular rectilinear grids can be defined by 1D coordinate arrays specifying the spacing along each dimension, while the other types require grid coordinates with the same dimensionality as the underlying value arrays, specifying the full n-dimensional coordinates of the corresponding array value. HoloViews provides many different elements supporting regularly spaced rectilinear grids, but currently only QuadMesh supports irregularly spaced rectilinear and curvilinear grids.\n",
    "\n",
    "The difference between uniform, rectilinear and curvilinear grids is best illustrated by the figure below:\n",
    "\n",
    "<figure>\n",
    "  <img src=\"http://earthsystemmodeling.org/docs/release/ESMF_8_1_1/ESMC_crefdoc/img11.png\" alt=\"grid-types\">\n",
    "  <figcaption>Types of logically rectangular grid tiles. Red circles show the values needed to specify grid coordinates for each type. Reproduced from <a href=\"http://earthsystemmodeling.org/docs/release/ESMF_8_1_1/ESMC_crefdoc/node5.html\">ESMF documentation</a></figcaption>\n",
    "</figure>\n",
    "\n",
    "\n",
    "In this section we will first discuss how to work with the simpler rectilinear grids and then describe how to define a curvilinear grid with 2D coordinate arrays.\n",
    "\n",
    "## Declaring gridded data\n",
    "\n",
    "All Elements that support a ColumnInterface also support the GridInterface. The simplest example of a multi-dimensional (or more precisely 2D) gridded dataset is an image, which has implicit or explicit x-coordinates, y-coordinates and an array representing the values for each combination of these coordinates. Let us start by declaring an Image with explicit x- and y-coordinates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = hv.Image((range(10), range(5), np.random.rand(5, 10)), datatype=['grid'])\n",
    "img"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above example we defined that there would be 10 samples along the x-axis, 5 samples along the y-axis and then defined a random ``5x10`` array, matching those dimensions. This follows the NumPy (row, column) indexing convention. When passing a tuple HoloViews will use the first gridded data interface, which stores the coordinates and value arrays as a dictionary mapping the dimension name to a NumPy array representing the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img.data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However HoloViews also ships with an interface for ``xarray`` and the [GeoViews](https://geoviews.org) library ships with an interface for ``iris`` objects, which are two common libraries for working with multi-dimensional datasets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr_img = img.clone(datatype=['image'])\n",
    "print(type(arr_img.data))\n",
    "\n",
    "try: \n",
    "    xr_img = img.clone(datatype=['xarray'])\n",
    "\n",
    "    print(type(xr_img.data))    \n",
    "except:\n",
    "    print('xarray interface could not be imported.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the case of an Image HoloViews also has a simple image representation which stores the data as a single array and converts the x- and y-coordinates to a set of bounds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Array type: %s with bounds %s\" % (type(arr_img.data), arr_img.bounds))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To summarize the constructor accepts a number of formats where the value arrays should always match the shape of the coordinate arrays:\n",
    "\n",
    "    1. A simple np.ndarray along with (l, b, r, t) bounds\n",
    "    2. A tuple of the coordinate and value arrays\n",
    "    3. A dictionary of the coordinate and value arrays indexed by their dimension names\n",
    "    3. XArray DataArray or XArray Dataset\n",
    "    4. An Iris cube"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Working with a multi-dimensional dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A gridded Dataset may have as many dimensions as desired, however individual Element types only support data of a certain dimensionality. Therefore we usually declare a ``Dataset`` to hold our multi-dimensional data and take it from there."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset3d = hv.Dataset((range(3), range(5), range(7), np.random.randn(7, 5, 3)),\n",
    "                       ['x', 'y', 'z'], 'Value')\n",
    "dataset3d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is because even a 3D multi-dimensional array represents volumetric data which we can display easily only if it contains few samples. In this simple case we can get an overview of what this data looks like by casting it to a ``Scatter3D`` Element (which will help us visualize the operations we are applying to the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.Scatter3D(dataset3d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Indexing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to explore the dataset we therefore often want to define a lower dimensional slice into the array and then convert the dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset3d.select(x=1).to(hv.Image, ['y', 'z']) + hv.Scatter3D(dataset3d.select(x=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Groupby\n",
    "\n",
    "Another common method to apply to our data is to facet or animate the data using ``groupby`` operations. HoloViews provides a convenient interface to apply ``groupby`` operations and select which dimensions to visualize. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(dataset3d.to(hv.Image, ['y', 'z'], 'Value', ['x']) +\n",
    "hv.HoloMap({x: hv.Scatter3D(dataset3d.select(x=x)) for x in range(3)}, kdims='x'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aggregating"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another common operation is to aggregate the data with a function thereby reducing a dimension. You can either ``aggregate`` the data by passing the dimensions to aggregate or ``reduce`` a specific dimension. Both have the same function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.Image(dataset3d.aggregate(['x', 'y'], np.mean)) + hv.Image(dataset3d.reduce(z=np.mean))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By aggregating the data we can reduce it to any number of dimensions we want. We can for example compute the spread of values for each z-coordinate and plot it using a ``Spread`` and ``Curve`` Element. We simply aggregate by that dimension and pass the aggregation functions we want to apply:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.Spread(dataset3d.aggregate('z', np.mean, np.std)) * hv.Curve(dataset3d.aggregate('z', np.mean))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible to generate lower-dimensional views into the dataset which can be useful to summarize the statistics of the data along a particular dimension. A simple example is a box-whisker of the ``Value`` for each x-coordinate. Using the ``.to`` conversion interface we declare that we want a ``BoxWhisker`` Element indexed by the ``x`` dimension showing the ``Value`` dimension. Additionally we have to ensure to set ``groupby`` to an empty list because by default the interface will group over any remaining dimension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset3d.to(hv.BoxWhisker, 'x', 'Value', groupby=[])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly we can generate a ``Distribution`` Element showing the ``Value`` dimension, group by the 'x' dimension and then overlay the distributions, giving us another statistical summary of the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset3d.to(hv.Distribution, 'Value', [], groupby='x').overlay()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Categorical dimensions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The key dimensions of the multi-dimensional arrays do not have to represent continuous values, we can display datasets with categorical variables as a ``HeatMap`` Element:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "heatmap = hv.HeatMap((['A', 'B', 'C'], ['a', 'b', 'c', 'd', 'e'], np.random.rand(5, 3)))\n",
    "heatmap + hv.Table(heatmap)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Non-uniform rectilinear grids\n",
    "\n",
    "As discussed above, there are two main types of grids handled by HoloViews.  So far, we have mainly dealt with uniform, rectilinear grids, but we can use the ``QuadMesh`` element to work with non-uniform rectilinear grids and curvilinear grids.\n",
    "\n",
    "In order to define a non-uniform, rectilinear grid we can declare explicit irregularly spaced x- and y-coordinates. In the example below we specify the x/y-coordinate bin edges of the grid as arrays of shape ``M+1`` and ``N+1`` and a value array (``zs``) of shape ``NxM``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 8  # Number of bins in each direction\n",
    "xs = np.logspace(1, 3, n)\n",
    "ys = np.linspace(1, 10, n)\n",
    "zs = np.arange((n-1)**2).reshape(n-1, n-1)\n",
    "print('Shape of x-coordinates:', xs.shape)\n",
    "print('Shape of y-coordinates:', ys.shape)\n",
    "print('Shape of value array:', zs.shape)\n",
    "hv.QuadMesh((xs, ys, zs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Curvilinear grids\n",
    "\n",
    "To define a curvilinear grid the x/y-coordinates of the grid should be defined as 2D arrays of shape ``NxM`` or ``N+1xM+1``, i.e. either as the bin centers or the bin edges of each 2D bin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n=20\n",
    "coords = np.linspace(-1.5,1.5,n)\n",
    "X,Y = np.meshgrid(coords, coords);\n",
    "Qx = np.cos(Y) - np.cos(X)\n",
    "Qy = np.sin(Y) + np.sin(X)\n",
    "Z = np.sqrt(X**2 + Y**2)\n",
    "\n",
    "print('Shape of x-coordinates:', Qx.shape)\n",
    "print('Shape of y-coordinates:', Qy.shape)\n",
    "print('Shape of value array:', Z.shape)\n",
    "\n",
    "qmesh = hv.QuadMesh((Qx, Qy, Z))\n",
    "qmesh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with xarray data types\n",
    "As demonstrated previously, `Dataset` comes with support for the `xarray` library, which offers a powerful way to work with multi-dimensional, regularly spaced data. In this example, we'll load an example dataset, turn it into a HoloViews `Dataset` and visualize it. First, let's have a look at the xarray dataset's contents:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xr_ds = xr.tutorial.open_dataset(\"air_temperature\").load()\n",
    "xr_ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is trivial to turn this xarray Dataset into a Holoviews `Dataset` (the same also works for DataArray):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv_ds = hv.Dataset(xr_ds)[:, :, \"2013-01-01\"]\n",
    "print(hv_ds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have used the usual slice notation in order to select one single day in the rather large dataset. Finally, let's visualize the dataset by converting it to a `HoloMap` of `Images` using the `to()` method. We need to specify which of the dataset's key dimensions will be consumed by the images (in this case \"lat\" and \"lon\"), where the remaining key dimensions will be associated with the HoloMap (here: \"time\"). We'll use the slice notation again to clip the longitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "airtemp = hv_ds.to(hv.Image, kdims=[\"lon\", \"lat\"], dynamic=False)\n",
    "airtemp[:, 220:320, :].opts(colorbar=True, fig_size=200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we have explicitly specified the default behaviour `dynamic=False`, which returns a HoloMap. Note, that this approach immediately converts all available data to images, which will take up a lot of RAM for large datasets. For these situations, use `dynamic=True` to generate a [DynamicMap](./07-Live_Data.ipynb) instead. Additionally, [xarray features dask support](http://xarray.pydata.org/en/stable/dask.html), which is helpful when dealing with large amounts of data.\n",
    "\n",
    "It is also possible to render curvilinear grids with xarray, and here we will load one such example. The dataset below defines a curvilinear grid of air temperatures varying over time. The curvilinear grid can be identified by the fact that the ``xc`` and ``yc`` coordinates are defined as two-dimensional arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rasm = xr.tutorial.open_dataset(\"rasm\").load()\n",
    "rasm.coords"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To simplify the example we will select a single timepoint and add explicit coordinates for the x and y dimensions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rasm = rasm.isel(time=0, x=slice(0, 200)).assign_coords(x=np.arange(200), y=np.arange(205))\n",
    "rasm.coords"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have defined both rectilinear and curvilinear coordinates we can visualize the difference between the two by explicitly defining which set of coordinates to use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.QuadMesh(rasm, ['x', 'y']) + hv.QuadMesh(rasm, ['xc', 'yc'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "Additional examples of visualizing xarrays in the context of geographical data can be found in the GeoViews documentation: [Gridded Datasets I](http://geoviews.org/user_guide/Gridded_Datasets_I.html) and\n",
    "[Gridded Datasets II](http://geoviews.org/user_guide/Gridded_Datasets_II.html). These guides also contain useful information on the interaction between xarray data structures and HoloViews Datasets in general."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# API"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Accessing the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to be able to work with data in different formats Holoviews defines a general interface to access the data. The dimension_values method allows returning underlying arrays.\n",
    "\n",
    "#### Key dimensions (coordinates)\n",
    "\n",
    "By default ``dimension_values`` will return the expanded columnar format of the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "heatmap.dimension_values('x')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To access just the unique coordinates along a dimension simply supply the ``expanded=False`` keyword:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "heatmap.dimension_values('x', expanded=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we can also get a non-flattened, expanded coordinate array returning a coordinate array of the same shape as the value arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "heatmap.dimension_values('x', flat=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Value dimensions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When accessing a value dimension the method will similarly return a flat view of the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "heatmap.dimension_values('z')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can pass the ``flat=False`` argument to access the multi-dimensional array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "heatmap.dimension_values('z', flat=False)"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}