diff --git a/doc/user_guide/Geometry_Data.rst b/doc/user_guide/Geometry_Data.rst
new file mode 100644
index 0000000000..219fe8eb90
--- /dev/null
+++ b/doc/user_guide/Geometry_Data.rst
@@ -0,0 +1,5 @@
+Geometry Data
+_____________
+
+.. notebook:: holoviews ../../examples/user_guide/Geometry_Data.ipynb
+    :offset: 1
diff --git a/doc/user_guide/index.rst b/doc/user_guide/index.rst
index 571b0b240a..ae7ba50b49 100644
--- a/doc/user_guide/index.rst
+++ b/doc/user_guide/index.rst
@@ -39,6 +39,9 @@ concepts in HoloViews:
 * `Gridded Datasets <Gridded_Datasets.html>`_
    Explore gridded data (n-dimensional arrays) with `NumPy <http://www.numpy.org/>`_ and `XArray <http://xarray.pydata.org/>`_.
 
+* `Geometry Data <Geometry_Data.html>`_
+   Working with and representing geometry data such as lines, multi-lines, polygons, multi-polygons and contours.
+
 * `Indexing and Selecting Data <Indexing_and_Selecting_Data.html>`_
    Select and index subsets of your data with HoloViews.
 
@@ -119,6 +122,7 @@ These guides provide detail about specific additional features in HoloViews:
     Live Data <Live_Data>
     Tabular Datasets <Tabular_Datasets>
     Gridded Datasets <Gridded_Datasets>
+    Geometry Data <Geometry_Data>
     Indexing and Selecting Data <Indexing_and_Selecting_Data>
     Transforming Elements <Transforming_Elements>
     Responding to Events <Responding_to_Events>
diff --git a/examples/reference/elements/bokeh/Contours.ipynb b/examples/reference/elements/bokeh/Contours.ipynb
index 309510c5d3..12024e3da5 100644
--- a/examples/reference/elements/bokeh/Contours.ipynb
+++ b/examples/reference/elements/bokeh/Contours.ipynb
@@ -28,7 +28,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "A ``Contours`` object is similar to a ``Path`` element but allows each individual path to be associated with one or more scalar values declared as value dimensions (``vdims``), which can be used to apply colormapping the ``Contours``. Just like the ``Path`` element ``Contours`` will accept a list of arrays, dataframes, a dictionaries of columns (or any of the other literal formats including tuples of columns and lists of tuples). In order to efficiently represent the scalar values associated with each path the dictionary format is preferable since it can store the scalar values without expanding them into a whole column.\n",
+    "A ``Contours`` object is similar to a ``Path`` element but allows each individual path to be associated with one or more scalar values declared as value dimensions (``vdims``), which can be used to apply colormapping the ``Contours``. Just like the ``Path`` element ``Contours`` will accept a list of arrays, dataframes, a dictionaries of columns (or any of the other literal formats including tuples of columns and lists of tuples). In order to efficiently represent the scalar values associated with each path the dictionary format is preferable since it can store the scalar values without expanding them into a whole column. For a full description of the path geometry data model see the [Geometry Data User Guide](../user_guide/Geometry_Data.ipynb).\n",
     "\n",
     "To see the effect we will create a number of concentric rings with increasing radii and define a colormap to apply color the circles: "
    ]
diff --git a/examples/reference/elements/bokeh/Path.ipynb b/examples/reference/elements/bokeh/Path.ipynb
index 221b45924c..021f2f0c65 100644
--- a/examples/reference/elements/bokeh/Path.ipynb
+++ b/examples/reference/elements/bokeh/Path.ipynb
@@ -28,7 +28,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "A ``Path`` object is actually a collection of lines, unlike ``Curve`` where the y-axis is the dependent variable, a ``Path`` consists of lines connecting arbitrary points in two-dimensional space. The individual subpaths should be supplied as a list and will be stored as NumPy arrays, DataFrames or dictionaries for each column, i.e. any of the formats accepted by columnar data formats.\n",
+    "A ``Path`` object is actually a collection of lines, unlike ``Curve`` where the y-axis is the dependent variable, a ``Path`` consists of lines connecting arbitrary points in two-dimensional space. The individual subpaths should be supplied as a list and will be stored as NumPy arrays, DataFrames or dictionaries for each column, i.e. any of the formats accepted by columnar data formats. For a full description of the path geometry data model see the [Geometry Data User Guide](../user_guide/Geometry_Data.ipynb). \n",
     "\n",
     "In this example we will create a Lissajous curve, which describe complex harmonic motion:"
    ]
diff --git a/examples/reference/elements/bokeh/Polygons.ipynb b/examples/reference/elements/bokeh/Polygons.ipynb
index de2fde5f58..a4a15e7a52 100644
--- a/examples/reference/elements/bokeh/Polygons.ipynb
+++ b/examples/reference/elements/bokeh/Polygons.ipynb
@@ -28,7 +28,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "A ``Polygons`` represents a contiguous filled area in a 2D space as a list of paths. Just like the ``Contours`` element additional scalar value dimensions maybe may be supplied, which can be used to color the ``Polygons`` with the defined ``cmap``. Like other ``Path`` types it accepts a list of arrays, dataframes, a dictionary of columns (or any of the other literal formats including tuples of columns and lists of tuples).\n",
+    "A ``Polygons`` represents a contiguous filled area in a 2D space as a list of polygon geometries. Just like the ``Contours`` element additional scalar value dimensions maybe may be supplied, which can be used to color the ``Polygons`` with the defined ``cmap``. Like other ``Path`` types it accepts a list of arrays, dataframes, a dictionary of columns (or any of the other literal formats including tuples of columns and lists of tuples), but also supports a special 'holes' key to represent empty interior regions. For a full description of the polygon geometry data model see the [Geometry Data User Guide](../user_guide/Geometry_Data.ipynb). \n",
     "\n",
     "In order to efficiently represent the scalar values associated with each path the dictionary format is preferable since it can store the scalar values without expanding them into a whole column. Additionally it allows passing multiple columns as a single array by specifying the dimension names as a tuple.\n",
     "\n",
diff --git a/examples/reference/elements/matplotlib/Contours.ipynb b/examples/reference/elements/matplotlib/Contours.ipynb
index bfe783d01f..3cedf8afa8 100644
--- a/examples/reference/elements/matplotlib/Contours.ipynb
+++ b/examples/reference/elements/matplotlib/Contours.ipynb
@@ -28,7 +28,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "A ``Contours`` object is similar to a ``Path`` element but allows each individual path to be associated with one or more scalar values declared as value dimensions (``vdims``), which can be used to apply colormapping the ``Contours``. Just like the ``Path`` element ``Contours`` will accept a list of arrays, dataframes, a dictionaries of columns (or any of the other literal formats including tuples of columns and lists of tuples). In order to efficiently represent the scalar values associated with each path the dictionary format is preferable since it can store the scalar values without expanding them into a whole column.\n",
+    "A ``Contours`` object is similar to a ``Path`` element but allows each individual path to be associated with one or more scalar values declared as value dimensions (``vdims``), which can be used to apply colormapping the ``Contours``. Just like the ``Path`` element ``Contours`` will accept a list of arrays, dataframes, a dictionaries of columns (or any of the other literal formats including tuples of columns and lists of tuples). In order to efficiently represent the scalar values associated with each path the dictionary format is preferable since it can store the scalar values without expanding them into a whole column. For a full description of the path geometry data model see the [Geometry Data User Guide](../user_guide/Geometry_Data.ipynb). \n",
     "\n",
     "To see the effect we will create a number of concentric rings with increasing radii and define a colormap to apply color the circles: "
    ]
diff --git a/examples/reference/elements/matplotlib/Path.ipynb b/examples/reference/elements/matplotlib/Path.ipynb
index 059d4bc9c2..df71f400be 100644
--- a/examples/reference/elements/matplotlib/Path.ipynb
+++ b/examples/reference/elements/matplotlib/Path.ipynb
@@ -28,7 +28,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "A ``Path`` object is actually a collection of lines, unlike ``Curve`` where the y-axis is the dependent variable, a ``Path`` consists of lines connecting arbitrary points in two-dimensional space. The individual subpaths should be supplied as a list and will be stored as NumPy arrays, DataFrames or dictionaries for each column, i.e. any of the formats accepted by columnar data formats.\n",
+    "A ``Path`` object is actually a collection of lines, unlike ``Curve`` where the y-axis is the dependent variable, a ``Path`` consists of lines connecting arbitrary points in two-dimensional space. The individual subpaths should be supplied as a list and will be stored as NumPy arrays, DataFrames or dictionaries for each column, i.e. any of the formats accepted by columnar data formats. For a full description of the path geometry data model see the [Geometry Data User Guide](../user_guide/Geometry_Data.ipynb). \n",
     "\n",
     "In this example we will create a Lissajous curve, which describe complex harmonic motion:"
    ]
diff --git a/examples/reference/elements/matplotlib/Polygons.ipynb b/examples/reference/elements/matplotlib/Polygons.ipynb
index 892398e763..da9c13bbaf 100644
--- a/examples/reference/elements/matplotlib/Polygons.ipynb
+++ b/examples/reference/elements/matplotlib/Polygons.ipynb
@@ -28,7 +28,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "A ``Polygons`` represents a contiguous filled area in a 2D space as a list of paths. Just like the ``Contours`` element additional scalar value dimensions maybe may be supplied, which can be used to color the ``Polygons`` with the defined ``cmap``. Like other ``Path`` types it accepts a list of arrays, dataframes, a dictionary of columns (or any of the other literal formats including tuples of columns and lists of tuples).\n",
+    "A ``Polygons`` represents a contiguous filled area in a 2D space as a list of polygon geometries. Just like the ``Contours`` element additional scalar value dimensions maybe may be supplied, which can be used to color the ``Polygons`` with the defined ``cmap``. Like other ``Path`` types it accepts a list of arrays, dataframes, a dictionary of columns (or any of the other literal formats including tuples of columns and lists of tuples), but also supports a special 'holes' key to represent empty interior regions. For a full description of the polygon geometry data model see the [Geometry Data User Guide](../user_guide/Geometry_Data.ipynb). \n",
     "\n",
     "In order to efficiently represent the scalar values associated with each path the dictionary format is preferable since it can store the scalar values without expanding them into a whole column. Additionally it allows passing multiple columns as a single array by specifying the dimension names as a tuple.\n",
     "\n",
diff --git a/examples/user_guide/Geometry_Data.ipynb b/examples/user_guide/Geometry_Data.ipynb
new file mode 100644
index 0000000000..a70ee7abea
--- /dev/null
+++ b/examples/user_guide/Geometry_Data.ipynb
@@ -0,0 +1,259 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In addition to the two main types of data, namely tabular/columnar and gridded data HoloViews also provide extensible interfaces to represent path geometry data. Specifically it has three main element types used to representing different types of geometries. In this section we will cover the HoloViews data model for representing different kinds of geometries.\n",
+    "\n",
+    "There are many different ways of representing path geometries but HoloViews' data model is oriented on GEOS geometry definitions and allows faithfully round-tripping data between its element types and GEOS geometry definitions such as ``LinearString``, ``Polygon``, ``MultiLineString`` and ``MultiPolygon`` geometries (even if this is not implemented in HoloViews itself). Since HoloViews interfaces are extensible many different formats for representing geometries could be supported (see [GeoViews](http://geoviews.org/user_guide/Geometries.html) for other representations) but here we will cover the native formats used by HoloViews to represent this data."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import holoviews as hv\n",
+    "\n",
+    "hv.extension('bokeh')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Representing paths\n",
+    "\n",
+    "The ``Path`` element represents a collection of path geometries with optional associated values. Each path geometry may be split into sub-geometries on NaN-values and may be associated with scalar values or array values varying along its length. In analogy to GEOS geometry types a Path is a collection of LineString and MultiLineString geometries with associated values.\n",
+    "\n",
+    "While many different formats are accepted in theory, natively HoloViews provides support for representing paths as lists of regular columnar data objects including arrays, dataframes and dictionaries of column arrays and scalars. A simple path geometry may therefore be drawn using:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hv.Path([{'x': [1, 2, 3, 4, 5], 'y': [0, 0, 1, 1, 2]}]).options(padding=0.1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Here the dictionary of x- and y-coordinates could also be an NumPy array with two columns or a dataframe with 'x' and 'y' columns. Since the format supports lists any number of geometries may be drawn in this way. Additionally, it is also possible to associate a value with each path by declaring it as a value dimension:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hv.Path([{'x': [1, 2, 3, 4, 5], 'y': [0, 0, 1, 1, 2], 'value': 0},\n",
+    "         {'x': [5, 4, 3, 2, 1], 'y': [2, 2, 1, 1, 0], 'value': 1}], vdims='value').options(padding=0.1, color_index='value')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Multi-geometry"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Splitting the geometries in this way allows assigning separate values to each geometry, however often multiple geometries share the same value in which case it may be desirable to represent them as a multi-geometry by combining the coordinates and separating them by a NaN value:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hv.Path([{'x': [1, 2, 3, 4, 5, np.nan, 5, 4, 3, 2, 1],\n",
+    "          'y': [0, 0, 1, 1, 2, np.nan, 2, 2, 1, 1, 0], 'value': 0}],\n",
+    "        vdims='value').options(padding=0.1, color_index='value')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This represents a more efficient format particularly when there are very many small geometries with the same value.\n",
+    "\n",
+    "#### Scalar vs. continuously varying value dimensions\n",
+    "\n",
+    "Unlike ``Contours`` which are limited to representing iso-contours or isoclines, i.e. a function of two variables which describes a curve along which the function has a constant value, a ``Path`` element may also have continuously varying values along its path. Below we will declare a path with a value that varies along its path:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "a, b, delta = 3, 5, np.pi/2.\n",
+    "\n",
+    "vs = np.linspace(0, np.pi*2, 200)\n",
+    "xs = np.sin(a * vs + delta)\n",
+    "ys = np.sin(b * vs)\n",
+    "\n",
+    "hv.Path([{'x': xs, 'y': ys, 'value': vs}], vdims='value').options(\n",
+    "    color_index='value', padding=0.1, cmap='hsv'\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Note that since not all data formats allow storing scalar values as actual scalars, 1D-arrays matching the length of the coordinates but with only one unique value are also considered scalar. For example the following is a valid ``Contours`` element despite the fact that the value dimension is not a scalar variable:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hv.Contours([{'x': xs, 'y': ys, 'value': np.ones(200)}], vdims='value').options(\n",
+    "    color_index='value', padding=0.1\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Representing Polygons\n",
+    "\n",
+    "The ``Polygons`` element represents a collection of polygon geometries with associated scalar values. Each polygon geometry may be split into sub-geometries on NaN-values and may be associated with scalar values. In analogy to GEOS geometry types a ``Polygons`` element is a collection of Polygon and MultiPolygon geometries. Polygon geometries are defined as a set of coordinates describing the exterior bounding ring and any number of interior holes.\n",
+    "\n",
+    "In summary ``Polygons`` can be represented in much the same way as ``Paths`` above but have a special reserved key to store the polygon interiors or 'holes'. The holes are stored as a list-of-lists of arrays. This nested format is necessary to unambiguously associate holes with the sub-geometries in a multi-geometry. In the simplest case of a single Polygon geometry the format looks like this:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "xs = [1, 2, 3]\n",
+    "ys = [2, 0, 7]\n",
+    "holes = [[[(1.5, 2), (2, 3), (1.6, 1.6)], [(2.1, 4.5), (2.5, 5), (2.3, 3.5)]]]\n",
+    "\n",
+    "hv.Polygons([{'x': xs, 'y': ys, 'holes': holes}]).options(padding=0.1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The 'x' and 'y' coordinates represent the exterior of the Polygon and the list-of-list of holes defines two interior regions inside the polygon.\n",
+    "\n",
+    "In a multi-Polygon arrangement where two Polygon geometries are separated by NaNs, the purpose of the nested format becomes a bit clearer. Here the polygon from above still has the two holes but the second polygon does not have any holes, which we declare with an empty list:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "xs = [1, 2, 3, np.nan, 6, 7, 3]\n",
+    "ys = [2, 0, 7, np.nan, 7, 5, 2]\n",
+    "\n",
+    "holes = [\n",
+    "    [[(1.5, 2), (2, 3), (1.6, 1.6)], [(2.1, 4.5), (2.5, 5), (2.3, 3.5)]],\n",
+    "    []\n",
+    "]\n",
+    "\n",
+    "hv.Polygons([{'x': xs, 'y': ys, 'holes': holes}]).options(padding=0.1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If a polygon has no holes at all the 'holes' key may be ommitted entirely:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hv.Polygons([{'x': xs, 'y': ys, 'holes': holes, 'value': 0},\n",
+    "             {'x': [4, 6, 6], 'y': [0, 2, 1], 'value': 1}, {'x': [-3, -1, -6], 'y': [3, 2, 1], 'value': 3}], vdims='value').options(padding=0.1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Accessing the data"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To access the underlying data the geometry elements (``Path``/``Contours``/``Polygons``) implement a ``split`` method. By default it simply returns a list of elements, where each contains only one geometry:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "poly = hv.Polygons([\n",
+    "    {'x': xs, 'y': ys, 'holes': holes, 'value': 0},\n",
+    "    {'x': [4, 6, 6], 'y': [0, 2, 1], 'value': 1}\n",
+    "], vdims='value')\n",
+    "\n",
+    "hv.Layout(poly.split())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Using the ``datatype`` argument the data may instead be returned in the desired format, e.g. a list of arrays:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "poly.split(datatype='array')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Note that this conversion may be lossy if the converted format has no way of representing 'holes' or other data."
+   ]
+  }
+ ],
+ "metadata": {
+  "language_info": {
+   "name": "python",
+   "pygments_lexer": "ipython3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/holoviews/core/data/dictionary.py b/holoviews/core/data/dictionary.py
index 933bcf1592..048392e14d 100644
--- a/holoviews/core/data/dictionary.py
+++ b/holoviews/core/data/dictionary.py
@@ -84,32 +84,36 @@ def init(cls, eltype, data, kdims, vdims):
 
         if not isinstance(data, cls.types):
             raise ValueError("DictInterface interface couldn't convert data.""")
-        elif isinstance(data, dict):
-            unpacked = []
-            for d, vals in data.items():
-                if isinstance(d, tuple):
-                    vals = np.asarray(vals)
-                    if vals.shape == (0,):
-                        for sd in d:
-                            unpacked.append((sd, np.array([], dtype=vals.dtype)))
-                    elif not vals.ndim == 2 and vals.shape[1] == len(d):
-                        raise ValueError("Values for %s dimensions did not have "
-                                         "the expected shape.")
-                    else:
-                        for i, sd in enumerate(d):
-                            unpacked.append((sd, vals[:, i]))
+
+        unpacked = []
+        for d, vals in data.items():
+            if isinstance(d, tuple):
+                vals = np.asarray(vals)
+                if vals.shape == (0,):
+                    for sd in d:
+                        unpacked.append((sd, np.array([], dtype=vals.dtype)))
+                elif not vals.ndim == 2 and vals.shape[1] == len(d):
+                    raise ValueError("Values for %s dimensions did not have "
+                                     "the expected shape.")
                 else:
-                    if not isscalar(vals):
-                        vals = np.asarray(vals)
-                        if not vals.ndim == 1 and d in dimensions:
-                            raise ValueError('DictInterface expects data for each column to be flat.')
-                    unpacked.append((d, vals))
-            if not cls.expanded([vs for d, vs in unpacked if d in dimensions and not isscalar(vs)]):
-                raise ValueError('DictInterface expects data to be of uniform shape.')
-            if isinstance(data, odict_types):
-                data.update(unpacked)
+                    for i, sd in enumerate(d):
+                        unpacked.append((sd, vals[:, i]))
+            elif d not in dimensions:
+                unpacked.append((d, vals))
             else:
-                data = OrderedDict(unpacked)
+                if not isscalar(vals):
+                    vals = np.asarray(vals)
+                    if not vals.ndim == 1 and d in dimensions:
+                        raise ValueError('DictInterface expects data for each column to be flat.')
+                unpacked.append((d, vals))
+
+        if not cls.expanded([vs for d, vs in unpacked if d in dimensions and not isscalar(vs)]):
+            raise ValueError('DictInterface expects data to be of uniform shape.')
+        if isinstance(data, odict_types):
+            data.update(unpacked)
+        else:
+            data = OrderedDict(unpacked)
+
         return data, {'kdims':kdims, 'vdims':vdims}, {}
 
 
@@ -151,7 +155,12 @@ def isscalar(cls, dataset, dim):
         values = dataset.data[name]
         if isscalar(values):
             return True
-        unique = set(values) if values.dtype.kind == 'O' else np.unique(values)
+        if values.dtype.kind == 'O':
+            unique = set(values)
+        else:
+            unique = np.unique(values)
+            if (~util.isfinite(unique)).all():
+                return True
         return len(unique) == 1
 
 
@@ -221,6 +230,15 @@ def sort(cls, dataset, by=[], reverse=False):
                             for d, v in dataset.data.items()])
 
 
+    @classmethod
+    def range(cls, dataset, dimension):
+        dim = dataset.get_dimension(dimension)
+        column = dataset.data[dim.name]
+        if isscalar(column):
+            return column, column
+        return Interface.range(dataset, dimension)
+
+
     @classmethod
     def values(cls, dataset, dim, expanded=True, flat=True):
         dim = dataset.get_dimension(dim).name
@@ -358,5 +376,20 @@ def iloc(cls, dataset, index):
             return arr if isscalar(arr) else arr[0]
         return new_data
 
+    @classmethod
+    def has_holes(cls, dataset):
+        from holoviews.element import Polygons
+        key = Polygons._hole_key
+        return key in dataset.data and isinstance(dataset.data[key], list)
+
+    @classmethod
+    def holes(cls, dataset):
+        from holoviews.element import Polygons
+        key = Polygons._hole_key
+        if key in dataset.data:
+            return [[[np.asarray(h) for h in hs] for hs in dataset.data[key]]]
+        else:
+            return super(DictInterface, cls).holes(dataset)
+
 
 Interface.register(DictInterface)
diff --git a/holoviews/core/data/interface.py b/holoviews/core/data/interface.py
index 2a9b3e5949..591f87c7ac 100644
--- a/holoviews/core/data/interface.py
+++ b/holoviews/core/data/interface.py
@@ -337,7 +337,9 @@ def range(cls, dataset, dimension):
         else:
             try:
                 assert column.dtype.kind not in 'SUO'
-                return (np.nanmin(column), np.nanmax(column))
+                with warnings.catch_warnings():
+                    warnings.filterwarnings('ignore', r'All-NaN (slice|axis) encountered')
+                    return (np.nanmin(column), np.nanmax(column))
             except (AssertionError, TypeError):
                 column = [v for v in util.python2sort(column) if v is not None]
                 if not len(column):
@@ -414,3 +416,13 @@ def nonzero(cls, dataset):
     @classmethod
     def redim(cls, dataset, dimensions):
         return dataset.data
+
+    @classmethod
+    def has_holes(cls, dataset):
+        return False
+
+    @classmethod
+    def holes(cls, dataset):
+        coords = cls.values(dataset, dataset.kdims[0])
+        splits = np.where(np.isnan(coords.astype('float')))[0]
+        return [[[]]*(len(splits)+1)]
diff --git a/holoviews/core/data/multipath.py b/holoviews/core/data/multipath.py
index 030d8663c1..b861174ba4 100644
--- a/holoviews/core/data/multipath.py
+++ b/holoviews/core/data/multipath.py
@@ -1,6 +1,7 @@
 import numpy as np
 
 from ..util import max_range
+from .dictionary import DictInterface
 from .interface import Interface, DataError
 
 
@@ -54,10 +55,22 @@ def init(cls, eltype, data, kdims, vdims):
     def validate(cls, dataset, vdims=True):
         if not dataset.data:
             return
+
+        from holoviews.element import Polygons
         ds = cls._inner_dataset_template(dataset)
         for d in dataset.data:
             ds.data = d
             ds.interface.validate(ds, vdims)
+            if isinstance(dataset, Polygons) and ds.interface is DictInterface:
+                holes = ds.interface.holes(ds)
+                if not isinstance(holes, list):
+                    raise DataError('Polygons holes must be declared as a list-of-lists.', cls)
+                subholes = holes[0]
+                coords = ds.data[ds.kdims[0].name]
+                splits = np.isnan(coords.astype('float')).sum()
+                if len(subholes) != (splits+1):
+                    raise DataError('Polygons with holes containing multi-geometries '
+                                    'must declare a list of holes for each geometry.', cls)
 
 
     @classmethod
@@ -100,6 +113,27 @@ def range(cls, dataset, dim):
         return max_range(ranges)
 
 
+    @classmethod
+    def has_holes(cls, dataset):
+        if not dataset.data:
+            return False
+        ds = cls._inner_dataset_template(dataset)
+        for d in dataset.data:
+            ds.data = d
+            if ds.interface.has_holes(ds):
+                return True
+        return False
+
+    @classmethod
+    def holes(cls, dataset):
+        holes = []
+        ds = cls._inner_dataset_template(dataset)
+        for d in dataset.data:
+            ds.data = d
+            holes += ds.interface.holes(ds)
+        return holes
+
+
     @classmethod
     def isscalar(cls, dataset, dim):
         """
diff --git a/holoviews/element/path.py b/holoviews/element/path.py
index 4d6413270c..2524b6d9c5 100644
--- a/holoviews/element/path.py
+++ b/holoviews/element/path.py
@@ -22,19 +22,34 @@
 
 class Path(Dataset, Element2D):
     """
-    The Path Element contains a list of Paths stored as tabular data
-    types including arrays, dataframes and dictionary of column
-    arrays. In addition a number of convenient constructors are
-    supported:
-
-    1) A list of lists containing x/y coordinate tuples.
-    2) A tuple containing an array of length N with the x-values and a
-       second array of shape NxP, where P is the number of paths.
-    3) A list of tuples each containing arrays x and y values.
-
-    A Path can be split into subpaths using the split method or combined
-    into a flat view using the dimension_values, table, and dframe methods,
-    where each path is separated by a NaN value.
+    The Path element represents a collection of path geometries with
+    associated values. Each path geometry may be split into
+    sub-geometries on NaN-values and may be associated with scalar
+    values or array values varying along its length. In analogy to
+    GEOS geometry types a Path is a collection of LineString and
+    MultiLineString geometries with associated values.
+
+    Like all other elements a Path may be defined through an
+    extensible list of interfaces. Natively, HoloViews provides the
+    MultiInterface which allows representing paths as lists of regular
+    columnar data objects including arrays, dataframes and
+    dictionaries of column arrays and scalars.
+
+    The canonical representation is a list of dictionaries storing the
+    x- and y-coordinates along with any other values:
+
+        [{'x': 1d-array, 'y': 1d-array, 'value': scalar, 'continuous': 1d-array}, ...]
+
+    Both scalar values and values continuously varying along the
+    geometries coordinates a Path may be used to color the geometry
+    by. Since not all formats allow storing scalar values as actual
+    scalars arrays which are the same length as the coordinates but
+    have only one unique value are also considered scalar.
+
+    The easiest way of accessing the individual geometries is using
+    the `Path.split` method, which returns each path geometry as a
+    separate entity, while the other methods assume a flattened
+    representation where all paths are separated by NaN values.
     """
 
     kdims = param.List(default=[Dimension('x'), Dimension('y')],
@@ -123,8 +138,29 @@ def split(self, start=None, end=None, datatype=None, **kwargs):
 
 class Contours(Path):
     """
-    Contours is a type of Path that is also associated with a value
-    (the contour level).
+    The Contours element is a subtype of a Path which is characterized
+    by the fact that each path geometry may only be associated with
+    scalar values. It supports all the same data formats as a `Path`
+    but does not allow continuously varying values along the path
+    geometry's coordinates. Conceptually Contours therefore represent
+    iso-contours or isoclines, i.e. a function of two variables which
+    describes a curve along which the function has a constant value.
+
+    The canonical representation is a list of dictionaries storing the
+    x- and y-coordinates along with any other (scalar) values:
+
+        [{'x': 1d-array, 'y': 1d-array, 'value': scalar}, ...]
+
+    Since not all formats allow storing scalar values as actual
+    scalars arrays which are the same length as the coordinates but
+    have only one unique value are also considered scalar. This is
+    strictly enforced, ensuring that each path geometry represents
+    a valid iso-contour.
+
+    The easiest way of accessing the individual geometries is using
+    the `Contours.split` method, which returns each path geometry as a
+    separate entity, while the other methods assume a flattened
+    representation where all paths are separated by NaN values.
     """
 
     level = param.Number(default=None, doc="""
@@ -169,8 +205,43 @@ def dimension_values(self, dim, expanded=True, flat=True):
 
 class Polygons(Contours):
     """
-    Polygons is a Path Element type that may contain any number of
-    closed paths with an associated value.
+    The Polygons element represents a collection of polygon geometries
+    with associated scalar values. Each polygon geometry may be split
+    into sub-geometries on NaN-values and may be associated with
+    scalar values. In analogy to GEOS geometry types a Polygons
+    element is a collection of Polygon and MultiPolygon
+    geometries. Polygon geometries are defined as a set of coordinates
+    describing the exterior bounding ring and any number of interior
+    holes.
+
+    Like all other elements a Polygons element may be defined through
+    an extensible list of interfaces. Natively HoloViews provides the
+    MultiInterface which allows representing paths as lists of regular
+    columnar data objects including arrays, dataframes and
+    dictionaries of column arrays and scalars.
+
+    The canonical representation is a list of dictionaries storing the
+    x- and y-coordinates, a list-of-lists of arrays representing the
+    holes, along with any other values:
+
+        [{'x': 1d-array, 'y': 1d-array, 'holes': list-of-lists-of-arrays, 'value': scalar}, ...]
+
+    The list-of-lists format of the holes corresponds to the potential
+    for each coordinate array to be split into a multi-geometry
+    through NaN-separators. Each sub-geometry separated by the NaNs
+    therefore has an unambiguous mapping to a list of holes. If a
+    (multi-)polygon has no holes, the 'holes' key may be ommitted.
+
+    Any value dimensions stored on a Polygons geometry must be scalar,
+    just like the Contours element. Since not all formats allow
+    storing scalar values as actual scalars arrays which are the same
+    length as the coordinates but have only one unique value are also
+    considered scalar.
+
+    The easiest way of accessing the individual geometries is using
+    the `Polygons.split` method, which returns each path geometry as a
+    separate entity, while the other methods assume a flattened
+    representation where all paths are separated by NaN values.
     """
 
     group = param.String(default="Polygons", constant=True)
@@ -181,6 +252,28 @@ class Polygons(Contours):
 
     _level_vdim = Dimension('Value')
 
+    # Defines which key the DictInterface uses to look for holes
+    _hole_key = 'holes'
+
+    @property
+    def has_holes(self):
+        """
+        Detects whether any polygon in the Polygons element defines
+        holes. Useful to avoid expanding Polygons unless necessary.
+        """
+        return self.interface.has_holes(self)
+
+    def holes(self):
+        """
+        Returns a list-of-lists-of-lists of hole arrays. The three levels
+        of nesting reflects the structure of the polygons:
+
+          1. The first level of nesting corresponds to the list of geometries
+          2. The second level corresponds to each Polygon in a MultiPolygon
+          3. The third level of nesting allows for multiple holes per Polygon
+        """
+        return self.interface.holes(self)
+
 
 class BaseShape(Path):
     """
diff --git a/holoviews/operation/element.py b/holoviews/operation/element.py
index c25926dfb5..117c376b4f 100644
--- a/holoviews/operation/element.py
+++ b/holoviews/operation/element.py
@@ -15,7 +15,7 @@
 from ..core.util import (group_sanitizer, label_sanitizer, pd,
                          basestring, datetime_types, isfinite, dt_to_int)
 from ..element.chart import Histogram, Scatter
-from ..element.raster import Raster, Image, RGB, QuadMesh
+from ..element.raster import Image, RGB
 from ..element.path import Contours, Polygons
 from ..element.util import categorical_aggregate2d # noqa (API import)
 from ..streams import RangeXY
@@ -429,15 +429,9 @@ def _process(self, element, key=None):
             raise ImportError("contours operation requires matplotlib.")
         extent = element.range(0) + element.range(1)[::-1]
 
-        if type(element) is Raster:
-            data = [np.flipud(element.data)]
-        elif isinstance(element, Image):
-            data = [np.flipud(element.dimension_values(2, flat=False))]
-        elif isinstance(element, QuadMesh):
-            data = (element.dimension_values(0, False, flat=False),
-                    element.dimension_values(1, False, flat=False),
-                    element.dimension_values(2, flat=False))
-
+        data = (element.dimension_values(0, False, flat=False),
+                element.dimension_values(1, False, flat=False),
+                element.dimension_values(2, flat=False))
         xdim, ydim = element.dimensions('key', label=True)
         if self.p.filled:
             contour_type = Polygons
@@ -445,35 +439,48 @@ def _process(self, element, key=None):
             contour_type = Contours
         vdims = element.vdims[:1]
 
+        levels = self.p.levels
+        zmin, zmax = element.range(2)
         if isinstance(self.p.levels, int):
-            levels = self.p.levels+2 if self.p.filled else self.p.levels+3
-            zmin, zmax = element.range(2)
-            levels = np.linspace(zmin, zmax, levels)
             if zmin == zmax:
                 contours = contour_type([], [xdim, ydim], vdims)
                 return (element * contours) if self.p.overlaid else contours
+            kwargs = {'N': self.p.levels}
         else:
-            levels = self.p.levels
+            kwargs = {'levels': levels}
 
         fig = Figure()
         ax = Axes(fig, [0, 0, 1, 1])
         contour_set = QuadContourSet(ax, *data, filled=self.p.filled,
-                                     extent=extent, levels=levels)
+                                     extent=extent, **kwargs)
+        levels = np.array(contour_set.get_array())
+        crange = levels.min(), levels.max()
+        if self.p.filled:
+            levels = levels[:-1] + np.diff(levels)/2.
+            vdims = [vdims[0].clone(range=crange)]
 
         paths = []
-        empty = np.full((1, 2), np.NaN)
-        for level, cset in zip(contour_set.get_array(), contour_set.collections):
-            subpaths = []
-            for path in cset.get_paths():
-                if path.codes is None:
-                    subpaths.append(path.vertices)
-                else:
-                    subpaths += np.split(path.vertices, np.where(path.codes==1)[0][1:])
-            if len(subpaths):
-                subpath = np.concatenate([p for sp in subpaths for p in (sp, empty)][:-1])
-            else:
-                subpath = np.array([])
-            paths.append({(xdim, ydim): subpath, element.vdims[0].name: level})
+        empty = np.array([[np.nan, np.nan]])
+        for level, cset in zip(levels, contour_set.collections):
+            exteriors = []
+            interiors = []
+            for geom in cset.get_paths():
+                interior = []
+                polys = geom.to_polygons(closed_only=False)
+                for ncp, cp in enumerate(polys):
+                    if ncp == 0:
+                        exteriors.append(cp)
+                        exteriors.append(empty)
+                    else:
+                        interior.append(cp)
+                if len(polys):
+                    interiors.append(interior)
+            if not exteriors:
+                continue
+            geom = {element.vdims[0].name: level, (xdim, ydim): np.concatenate(exteriors[:-1])}
+            if self.p.filled and interiors:
+                geom['holes'] = interiors
+            paths.append(geom)
         contours = contour_type(paths, label=element.label, kdims=element.kdims, vdims=vdims)
         if self.p.overlaid:
             contours = element * contours
diff --git a/holoviews/plotting/bokeh/path.py b/holoviews/plotting/bokeh/path.py
index 590a291176..335b28171f 100644
--- a/holoviews/plotting/bokeh/path.py
+++ b/holoviews/plotting/bokeh/path.py
@@ -4,8 +4,9 @@
 import numpy as np
 
 from ...core import util
+from ...element import Polygons
 from .element import ColorbarPlot, LegendPlot, line_properties, fill_properties
-from .util import expand_batched_style
+from .util import expand_batched_style, mpl_to_bokeh, bokeh_version, multi_polygons_data
 
 
 class PathPlot(ColorbarPlot):
@@ -65,12 +66,14 @@ def get_data(self, element, ranges, style):
             vals = {util.dimension_sanitizer(vd.name): [] for vd in element.vdims}
             for path in element.split():
                 cvals = path.dimension_values(cdim)
+                array = path.array(path.kdims)
                 splits = [0]+list(np.where(np.diff(cvals)!=0)[0]+1)
+                cols = {vd.name: path.dimension_values(vd) for vd in element.vdims}
                 if len(splits) == 1:
                     splits.append(len(path))
                 for (s1, s2) in zip(splits[:-1], splits[1:]):
                     for i, vd in enumerate(element.vdims):
-                        path_val = path.iloc[s1, i+2]
+                        path_val = cols[vd.name][s1]
                         vd_column = util.dimension_sanitizer(vd.name)
                         dt_column = vd_column+'_dt_strings'
                         vals[vd_column].append(path_val)
@@ -78,7 +81,7 @@ def get_data(self, element, ranges, style):
                             if dt_column not in vals:
                                 vals[dt_column] = []
                             vals[dt_column].append(vd.pprint_value(path_val))
-                    paths.append(path.iloc[s1:s2+1, :2].array())
+                    paths.append(array[s1:s2+1])
             xs, ys = ([path[:, idx] for path in paths] for idx in inds)
             data = dict(xs=xs, ys=ys, **{d: np.array(vs) for d, vs in vals.items()})
         cmapper = self._get_colormapper(cdim, element, ranges, style)
@@ -157,15 +160,20 @@ def _get_hover_data(self, data, element):
                 data[dim] = [v for _ in range(len(list(data.values())[0]))]
 
     def get_data(self, element, ranges, style):
-        paths = element.split(datatype='array', dimensions=element.kdims)
+        has_holes = isinstance(element, Polygons) and element.has_holes
         if self.static_source:
             data = dict()
         else:
-            inds = (1, 0) if self.invert_axes else (0, 1)
-            xs, ys = ([path[:, idx] for path in paths] for idx in inds)
+            if has_holes and bokeh_version >= '1.0':
+                xs, ys = multi_polygons_data(element)
+                style['has_holes'] = has_holes
+            else:
+                paths = element.split(datatype='array', dimensions=element.kdims)
+                xs, ys = ([path[:, idx] for path in paths] for idx in (0, 1))
+            if self.invert_axes:
+                xs, ys = ys, xs
             data = dict(xs=xs, ys=ys)
         mapping = dict(self._mapping)
-
         if None not in [element.level, self.color_index] and element.vdims:
             cdim = element.vdims[0]
         else:
@@ -174,7 +182,7 @@ def get_data(self, element, ranges, style):
         if cdim is None:
             return data, mapping, style
 
-        ncontours = len(paths)
+        ncontours = len(xs)
         dim_name = util.dimension_sanitizer(cdim.name)
         if element.level is not None:
             values = np.full(ncontours, float(element.level))
@@ -189,6 +197,24 @@ def get_data(self, element, ranges, style):
             mapping['legend'] = dim_name
         return data, mapping, style
 
+    def _init_glyph(self, plot, mapping, properties):
+        """
+        Returns a Bokeh glyph object.
+        """
+        has_holes = properties.pop('has_holes', False)
+        plot_method = properties.pop('plot_method', None)
+        properties = mpl_to_bokeh(properties)
+        data = dict(properties, **mapping)
+        if has_holes:
+            plot_method = 'multi_polygons'
+        elif plot_method is None:
+            plot_method = self._plot_methods.get('single')
+        renderer = getattr(plot, plot_method)(**data)
+        if self.colorbar and 'color_mapper' in self.handles:
+            self._draw_colorbar(plot, self.handles['color_mapper'])
+        return renderer, renderer.glyph
+
+
 
 class PolygonPlot(ContourPlot):
 
diff --git a/holoviews/plotting/bokeh/util.py b/holoviews/plotting/bokeh/util.py
index 437302bd6d..e07f42aaa4 100644
--- a/holoviews/plotting/bokeh/util.py
+++ b/holoviews/plotting/bokeh/util.py
@@ -675,3 +675,30 @@ def theme_attr_json(theme, attr):
         return theme._json['attrs'].get(attr, {})
     else:
         return {}
+
+
+def multi_polygons_data(element):
+    """
+    Expands polygon data which contains holes to a bokeh multi_polygons
+    representation. Multi-polygons split by nans are expanded and the
+    correct list of holes is assigned to each sub-polygon.
+    """
+    paths = element.split(datatype='array', dimensions=element.kdims)
+    xs, ys = ([path[:, idx] for path in paths] for idx in (0, 1))
+    holes = element.holes()
+    xsh, ysh = [], []
+    for x, y, multi_hole in zip(xs, ys, holes):
+        xhs = [[h[:, 0] for h in hole] for hole in multi_hole]
+        yhs = [[h[:, 1] for h in hole] for hole in multi_hole]
+        array = np.column_stack([x, y])
+        splits = np.where(np.isnan(array[:, :2].astype('float')).sum(axis=1))[0]
+        arrays = np.split(array, splits+1) if len(splits) else [array]
+        multi_xs, multi_ys = [], []
+        for i, (path, hx, hy) in enumerate(zip(arrays, xhs, yhs)):
+            if i != (len(arrays)-1):
+                path = path[:-1]
+            multi_xs.append([path[:, 0]]+hx)
+            multi_ys.append([path[:, 1]]+hy)
+        xsh.append(multi_xs)
+        ysh.append(multi_ys)
+    return xsh, ysh
diff --git a/holoviews/plotting/mpl/path.py b/holoviews/plotting/mpl/path.py
index 89fc171569..beabefba3d 100644
--- a/holoviews/plotting/mpl/path.py
+++ b/holoviews/plotting/mpl/path.py
@@ -1,9 +1,11 @@
-from matplotlib.collections import PolyCollection, LineCollection
-import numpy as np
 import param
+import numpy as np
+from matplotlib.collections import PatchCollection, LineCollection
 
 from ...core import util
+from ...element import Polygons
 from .element import ColorbarPlot
+from .util import polygons_to_path_patches
 
 
 class PathPlot(ColorbarPlot):
@@ -82,9 +84,17 @@ def get_data(self, element, ranges, style):
         else:
             cidx = self.color_index+2 if isinstance(self.color_index, int) else self.color_index
             cdim = element.get_dimension(cidx)
-        paths = element.split(datatype='array', dimensions=element.kdims)
-        if self.invert_axes:
-            paths = [p[:, ::-1] for p in paths]
+
+        if isinstance(element, Polygons):
+            subpaths = polygons_to_path_patches(element)
+            paths = [path for subpath in subpaths for path in subpath]
+            if self.invert_axes:
+                for p in paths:
+                    p._path.vertices = p._path.vertices[:, ::-1]
+        else:
+            paths = element.split(datatype='array', dimensions=element.kdims)
+            if self.invert_axes:
+                paths = [p[:, ::-1] for p in paths]
 
         if cdim is None:
             return (paths,), style, {}
@@ -93,9 +103,15 @@ def get_data(self, element, ranges, style):
             array = np.full(len(paths), element.level)
         else:
             array = element.dimension_values(cdim, expanded=False)
+            if len(paths) != len(array):
+                # If there are multi-geometries the list of scalar values
+                # will not match the list of paths and has to be expanded
+                array = np.array([v for v, sps in zip(array, subpaths)
+                                  for _ in range(len(sps))])
+
         if array.dtype.kind not in 'uif':
             array = np.searchsorted(np.unique(array), array)
-        style['array']= array
+        style['array'] = array
         self._norm_kwargs(element, ranges, style, cdim)
         style['clim'] = style.pop('vmin'), style.pop('vmax')
         return (paths,), style, {}
@@ -117,6 +133,6 @@ class PolygonPlot(ContourPlot):
                   'hatch', 'linestyle', 'joinstyle', 'fill', 'capstyle']
 
     def init_artists(self, ax, plot_args, plot_kwargs):
-        polys = PolyCollection(*plot_args, **plot_kwargs)
+        polys = PatchCollection(*plot_args, **plot_kwargs)
         ax.add_collection(polys)
         return {'artist': polys}
diff --git a/holoviews/plotting/mpl/util.py b/holoviews/plotting/mpl/util.py
index bd84eb284c..e87b397391 100644
--- a/holoviews/plotting/mpl/util.py
+++ b/holoviews/plotting/mpl/util.py
@@ -5,12 +5,13 @@
 import numpy as np
 import matplotlib
 from matplotlib import ticker
+from matplotlib.patches import Path, PathPatch
 from matplotlib.transforms import Bbox, TransformedBbox, Affine2D
 
 mpl_version = LooseVersion(matplotlib.__version__)  # noqa
 
 from ...core.util import basestring, _getargspec
-from ...element import Raster, RGB
+from ...element import Raster, RGB, Polygons
 
 
 def wrap_formatter(formatter):
@@ -185,3 +186,42 @@ def get_raster_array(image):
         else:
             data = np.flipud(data)
     return data
+
+
+def ring_coding(array):
+    """
+    Produces matplotlib Path codes for exterior and interior rings
+    of a polygon geometry.
+    """
+    # The codes will be all "LINETO" commands, except for "MOVETO"s at the
+    # beginning of each subpath
+    n = len(array)
+    codes = np.ones(n, dtype=Path.code_type) * Path.LINETO
+    codes[0] = Path.MOVETO
+    return codes
+
+
+def polygons_to_path_patches(element):
+    """
+    Converts Polygons into list of lists of matplotlib.patches.PathPatch
+    objects including any specified holes. Each list represents one
+    (multi-)polygon.
+    """
+    paths = element.split(datatype='array', dimensions=element.kdims)
+    has_holes = isinstance(element, Polygons) and element.interface.has_holes(element)
+    holes = element.interface.holes(element) if has_holes else None
+    mpl_paths = []
+    for i, path in enumerate(paths):
+        splits = np.where(np.isnan(path[:, :2].astype('float')).sum(axis=1))[0]
+        arrays = np.split(path, splits+1) if len(splits) else [path]
+        subpath = []
+        for j, array in enumerate(arrays):
+            if j != (len(arrays)-1):
+                array = array[:-1]
+            interiors = holes[i][j] if has_holes else []
+            vertices = np.concatenate([array]+interiors)
+            codes = np.concatenate([ring_coding(array)]+
+                                   [ring_coding(h) for h in interiors])
+            subpath.append(PathPatch(Path(vertices, codes)))
+        mpl_paths.append(subpath)
+    return mpl_paths
diff --git a/holoviews/tests/element/testpaths.py b/holoviews/tests/element/testpaths.py
index 2142816d71..86e25a14ce 100644
--- a/holoviews/tests/element/testpaths.py
+++ b/holoviews/tests/element/testpaths.py
@@ -2,10 +2,86 @@
 Unit tests of Path types.
 """
 import numpy as np
-from holoviews import Ellipse, Box
+from holoviews import Ellipse, Box, Polygons
+from holoviews.core.data.interface import DataError
 from holoviews.element.comparison import ComparisonTestCase
 
 
+class PolygonsTests(ComparisonTestCase):
+
+    def setUp(self):
+        xs = [1, 2, 3]
+        ys = [2, 0, 7]
+        holes = [[[(1.5, 2), (2, 3), (1.6, 1.6)], [(2.1, 4.5), (2.5, 5), (2.3, 3.5)]]]
+        self.single_poly = Polygons([{'x': xs, 'y': ys, 'holes': holes}])
+
+        xs = [1, 2, 3, np.nan, 6, 7, 3]
+        ys = [2, 0, 7, np.nan, 7, 5, 2]
+        holes = [
+            [[(1.5, 2), (2, 3), (1.6, 1.6)], [(2.1, 4.5), (2.5, 5), (2.3, 3.5)]],
+            []
+        ]
+        self.multi_poly = Polygons([{'x': xs, 'y': ys, 'holes': holes}])
+        self.multi_poly_no_hole = Polygons([{'x': xs, 'y': ys}])
+
+        self.distinct_polys = Polygons([
+            {'x': xs, 'y': ys, 'holes': holes, 'value': 0},
+            {'x': [4, 6, 6], 'y': [0, 2, 1], 'value': 1}], vdims='value')
+
+    def test_single_poly_holes_match(self):
+        self.assertTrue(self.single_poly.interface.has_holes(self.single_poly))
+        paths = self.single_poly.split(datatype='array')
+        holes = self.single_poly.interface.holes(self.single_poly)
+        self.assertEqual(len(paths), len(holes))
+        self.assertEqual(len(holes), 1)
+        self.assertEqual(len(holes[0]), 1)
+        self.assertEqual(len(holes[0][0]), 2)
+
+    def test_multi_poly_holes_match(self):
+        self.assertTrue(self.multi_poly.interface.has_holes(self.multi_poly))
+        paths = self.multi_poly.split(datatype='array')
+        holes = self.multi_poly.interface.holes(self.multi_poly)
+        self.assertEqual(len(paths), len(holes))
+        self.assertEqual(len(holes), 1)
+        self.assertEqual(len(holes[0]), 2)
+        self.assertEqual(len(holes[0][0]), 2)
+        self.assertEqual(len(holes[0][1]), 0)
+
+    def test_multi_poly_no_holes_match(self):
+        self.assertFalse(self.multi_poly_no_hole.interface.has_holes(self.multi_poly_no_hole))
+        paths = self.multi_poly_no_hole.split(datatype='array')
+        holes = self.multi_poly_no_hole.interface.holes(self.multi_poly_no_hole)
+        self.assertEqual(len(paths), len(holes))
+        self.assertEqual(len(holes), 1)
+        self.assertEqual(len(holes[0]), 2)
+        self.assertEqual(len(holes[0][0]), 0)
+        self.assertEqual(len(holes[0][1]), 0)
+
+    def test_distinct_multi_poly_holes_match(self):
+        self.assertTrue(self.distinct_polys.interface.has_holes(self.distinct_polys))
+        paths = self.distinct_polys.split(datatype='array')
+        holes = self.distinct_polys.interface.holes(self.distinct_polys)
+        self.assertEqual(len(paths), len(holes))
+        self.assertEqual(len(holes), 2)
+        self.assertEqual(len(holes[0]), 2)
+        self.assertEqual(len(holes[0][0]), 2)
+        self.assertEqual(len(holes[0][1]), 0)
+        self.assertEqual(len(holes[1]), 1)
+        self.assertEqual(len(holes[1][0]), 0)
+
+    def test_single_poly_hole_validation(self):
+        xs = [1, 2, 3]
+        ys = [2, 0, 7]
+        with self.assertRaises(DataError):
+            Polygons([{'x': xs, 'y': ys, 'holes': [[], []]}])
+
+    def test_multi_poly_hole_validation(self):
+        xs = [1, 2, 3, np.nan, 6, 7, 3]
+        ys = [2, 0, 7, np.nan, 7, 5, 2]
+        with self.assertRaises(DataError):
+            Polygons([{'x': xs, 'y': ys, 'holes': [[]]}])
+
+
 class EllipseTests(ComparisonTestCase):
 
     def setUp(self):
diff --git a/holoviews/tests/element/teststatselements.py b/holoviews/tests/element/teststatselements.py
index 4d988ad857..a2670a49fa 100644
--- a/holoviews/tests/element/teststatselements.py
+++ b/holoviews/tests/element/teststatselements.py
@@ -143,7 +143,7 @@ def test_distribution_composite_custom_vdim(self):
         self.assertEqual(area.vdims, [Dimension('Test')])
 
     def test_distribution_composite_not_filled(self):
-        dist = Distribution(np.array([0, 1, 2])).opts(plot=dict(filled=False))
+        dist = Distribution(np.array([0, 1, 2]), ).opts(plot=dict(filled=False))
         curve = Compositor.collapse_element(dist, backend='matplotlib')
         self.assertIsInstance(curve, Curve)
         self.assertEqual(curve.vdims, [Dimension(('Value_density', 'Value Density'))])
@@ -182,7 +182,7 @@ def test_bivariate_composite_filled(self):
         dist = Bivariate(np.random.rand(10, 2)).opts(plot=dict(filled=True))
         contours = Compositor.collapse_element(dist)
         self.assertIsInstance(contours, Polygons)
-        self.assertEqual(contours.vdims, [Dimension('Density')])
+        self.assertEqual(contours.vdims[0].name, 'Density')
 
     def test_bivariate_composite_empty_filled(self):
         dist = Bivariate([]).opts(plot=dict(filled=True))
diff --git a/holoviews/tests/operation/testoperation.py b/holoviews/tests/operation/testoperation.py
index ad8e9e114f..6438fe9b1e 100644
--- a/holoviews/tests/operation/testoperation.py
+++ b/holoviews/tests/operation/testoperation.py
@@ -61,9 +61,9 @@ def test_image_gradient(self):
     def test_image_contours(self):
         img = Image(np.array([[0, 1, 0], [3, 4, 5.], [6, 7, 8]]))
         op_contours = contours(img, levels=[0.5])
-        contour = Contours([[(-0.5,  0.416667, 0.5), (-0.25, 0.5, 0.5),
-                             (np.NaN, np.NaN, 0.5), (0.25, 0.5, 0.5),
-                             (0.5, 0.45, 0.5)]],
+        contour = Contours([[(-0.166667,  0.333333, 0.5), (-0.333333, 0.277778, 0.5),
+                             (np.NaN, np.NaN, 0.5), (0.333333, 0.3, 0.5),
+                             (0.166667, 0.333333, 0.5)]],
                             vdims=img.vdims)
         self.assertEqual(op_contours, contour)
 
@@ -101,9 +101,9 @@ def test_qmesh_curvilinear_contours(self):
     def test_image_contours_filled(self):
         img = Image(np.array([[0, 1, 0], [3, 4, 5.], [6, 7, 8]]))
         op_contours = contours(img, filled=True, levels=[2, 2.5])
-        data = [[(0., 0.333333, 2), (0.5, 0.3, 2), (0.5, 0.25, 2), (0., 0.25, 2),
-                 (-0.5, 0.08333333, 2), (-0.5, 0.16666667, 2), (0., 0.33333333, 2)]]
-        polys = Polygons(data, vdims=img.vdims)
+        data = [[(0., 0.166667, 2.25), (0.333333, 0.166667, 2.25), (0.333333, 0.2, 2.25), (0., 0.222222, 2.25),
+                 (-0.333333, 0.111111, 2.25), (-0.333333, 0.055556, 2.25), (0., 0.166667, 2.25)]]
+        polys = Polygons(data, vdims=img.vdims[0].clone(range=(2, 2.5)))
         self.assertEqual(op_contours, polys)
 
     def test_points_histogram(self):
diff --git a/holoviews/tests/operation/teststatsoperations.py b/holoviews/tests/operation/teststatsoperations.py
index 5966b32f5c..88baeb409f 100644
--- a/holoviews/tests/operation/teststatsoperations.py
+++ b/holoviews/tests/operation/teststatsoperations.py
@@ -56,18 +56,20 @@ def test_bivariate_kde(self):
         self.assertEqual(kde, img)
 
     def test_bivariate_kde_contours(self):
+        np.random.seed(1)
         bivariate = Bivariate(np.random.rand(100, 2))
         kde = bivariate_kde(bivariate, n_samples=100, x_range=(0, 1),
                             y_range=(0, 1), contours=True, levels=10)
         self.assertIsInstance(kde, Contours)
-        self.assertEqual(len(kde.data), 10)
+        self.assertEqual(len(kde.data), 6)
 
     def test_bivariate_kde_contours_filled(self):
+        np.random.seed(1)
         bivariate = Bivariate(np.random.rand(100, 2))
         kde = bivariate_kde(bivariate, n_samples=100, x_range=(0, 1),
                             y_range=(0, 1), contours=True, filled=True, levels=10)
         self.assertIsInstance(kde, Polygons)
-        self.assertEqual(len(kde.data), 10)
+        self.assertEqual(len(kde.data), 7)
 
     def test_bivariate_kde_nans(self):
         kde = bivariate_kde(self.bivariate_nans, n_samples=2, x_range=(0, 4),
diff --git a/holoviews/tests/plotting/bokeh/testpathplot.py b/holoviews/tests/plotting/bokeh/testpathplot.py
index c471b01468..7a9dba00b6 100644
--- a/holoviews/tests/plotting/bokeh/testpathplot.py
+++ b/holoviews/tests/plotting/bokeh/testpathplot.py
@@ -1,9 +1,11 @@
 import datetime as dt
+from unittest import SkipTest
 
 import numpy as np
 from holoviews.core import NdOverlay
 from holoviews.core.options import Cycle
 from holoviews.element import Path, Polygons, Contours
+from holoviews.plotting.bokeh.util import bokeh_version
 
 from .testplot import TestBokehPlot, bokeh_renderer
 
@@ -161,6 +163,36 @@ def test_empty_polygons_plot(self):
         self.assertEqual(len(source.data['ys']), 0)
         self.assertEqual(len(source.data['Intensity']), 0)
 
+    def test_polygon_with_hole_plot(self):
+        if bokeh_version < '1.0':
+            raise SkipTest('Plotting Polygons with holes requires bokeh >= 1.0')
+        xs = [1, 2, 3]
+        ys = [2, 0, 7]
+        holes = [[[(1.5, 2), (2, 3), (1.6, 1.6)], [(2.1, 4.5), (2.5, 5), (2.3, 3.5)]]]
+        poly = Polygons([{'x': xs, 'y': ys, 'holes': holes}])
+        plot = bokeh_renderer.get_plot(poly)
+        source = plot.handles['source']
+        self.assertEqual(source.data['xs'], [[[np.array([1, 2, 3]), np.array([1.5, 2, 1.6]),
+                                              np.array([2.1, 2.5, 2.3])]]])
+        self.assertEqual(source.data['ys'], [[[np.array([2, 0, 7]), np.array([2, 3, 1.6]),
+                                              np.array([4.5, 5, 3.5])]]])
+
+    def test_multi_polygon_hole_plot(self):
+        if bokeh_version < '1.0':
+            raise SkipTest('Plotting Polygons with holes requires bokeh >= 1.0')
+        xs = [1, 2, 3, np.nan, 6, 7, 3]
+        ys = [2, 0, 7, np.nan, 7, 5, 2]
+        holes = [
+            [[(1.5, 2), (2, 3), (1.6, 1.6)], [(2.1, 4.5), (2.5, 5), (2.3, 3.5)]],
+            []
+        ]
+        poly = Polygons([{'x': xs, 'y': ys, 'holes': holes}])
+        plot = bokeh_renderer.get_plot(poly)
+        source = plot.handles['source']
+        self.assertEqual(source.data['xs'], [[[np.array([1, 2, 3]), np.array([1.5, 2, 1.6]),
+                                               np.array([2.1, 2.5, 2.3])], [np.array([6, 7, 3])]]])
+        self.assertEqual(source.data['ys'], [[[np.array([2, 0, 7]), np.array([2, 3, 1.6]),
+                                               np.array([4.5, 5, 3.5])], [np.array([7, 5, 2])]]])
 
 
 class TestContoursPlot(TestBokehPlot):
diff --git a/holoviews/tests/plotting/matplotlib/testpathplot.py b/holoviews/tests/plotting/matplotlib/testpathplot.py
index be3fb01f55..3a8610ebd8 100644
--- a/holoviews/tests/plotting/matplotlib/testpathplot.py
+++ b/holoviews/tests/plotting/matplotlib/testpathplot.py
@@ -17,6 +17,45 @@ def test_polygons_colored(self):
             self.assertEqual(artist.get_array(), np.array([j]))
             self.assertEqual(artist.get_clim(), (0, 4))
 
+    def test_polygon_with_hole_plot(self):
+        xs = [1, 2, 3]
+        ys = [2, 0, 7]
+        holes = [[[(1.5, 2), (2, 3), (1.6, 1.6)], [(2.1, 4.5), (2.5, 5), (2.3, 3.5)]]]
+        poly = Polygons([{'x': xs, 'y': ys, 'holes': holes}])
+        plot = mpl_renderer.get_plot(poly)
+        artist = plot.handles['artist']
+        paths = artist.get_paths()
+        self.assertEqual(len(paths), 1)
+        path = paths[0]
+        self.assertEqual(path.vertices, np.array([
+            (1, 2), (2, 0), (3, 7), (1.5, 2), (2, 3), (1.6, 1.6),
+            (2.1, 4.5), (2.5, 5), (2.3, 3.5)])
+        )
+        self.assertEqual(path.codes, np.array([1, 2, 2, 1, 2, 2, 1, 2, 2]))
+
+    def test_multi_polygon_hole_plot(self):
+        xs = [1, 2, 3, np.nan, 6, 7, 3]
+        ys = [2, 0, 7, np.nan, 7, 5, 2]
+        holes = [
+            [[(1.5, 2), (2, 3), (1.6, 1.6)], [(2.1, 4.5), (2.5, 5), (2.3, 3.5)]],
+            []
+        ]
+        poly = Polygons([{'x': xs, 'y': ys, 'holes': holes, 'value': 1}], vdims=['value'])
+        plot = mpl_renderer.get_plot(poly)
+        artist = plot.handles['artist']
+        self.assertEqual(artist.get_array(), np.array([1, 1]))
+        paths = artist.get_paths()
+        self.assertEqual(len(paths), 2)
+        path = paths[0]
+        self.assertEqual(path.vertices, np.array([
+            (1, 2), (2, 0), (3, 7), (1.5, 2), (2, 3), (1.6, 1.6),
+            (2.1, 4.5), (2.5, 5), (2.3, 3.5)])
+        )
+        self.assertEqual(path.codes, np.array([1, 2, 2, 1, 2, 2, 1, 2, 2]))
+        path2 = paths[1]
+        self.assertEqual(path2.vertices, np.array([(6, 7), (7, 5), (3, 2)]))
+        self.assertEqual(path2.codes, np.array([1, 2, 2]))
+
 
 class TestContoursPlot(TestMPLPlot):