diff --git a/holoviews/core/data/__init__.py b/holoviews/core/data/__init__.py
index c4e1764c75..1f70b2257a 100644
--- a/holoviews/core/data/__init__.py
+++ b/holoviews/core/data/__init__.py
@@ -14,6 +14,7 @@
 from .array import ArrayInterface
 from .dictionary import DictInterface
 from .grid import GridInterface
+from .image import ImageInterface
 from .ndelement import NdElementInterface
 
 datatypes = ['array', 'dictionary', 'grid', 'ndelement']
@@ -276,11 +277,12 @@ def select(self, selection_specs=None, **selection):
         if selection_specs and not any(self.matches(sp) for sp in selection_specs):
             return self
 
-        data = self.interface.select(self, **selection)
+        data, kwargs = self.interface.select(self, **selection)
+
         if np.isscalar(data):
             return data
         else:
-            return self.clone(data)
+            return self.clone(data, **kwargs)
 
 
     def reindex(self, kdims=None, vdims=None):
@@ -395,11 +397,19 @@ def aggregate(self, dimensions=None, function=None, spreadfn=None, **kwargs):
                 combined = combined.add_dimension(dim, ndims+i, dvals, True)
             return combined
 
+        ndims = len(dimensions)
+        min_d, max_d = self.params('kdims').bounds
         if np.isscalar(aggregated):
             return aggregated
         else:
-            return self.clone(aggregated, kdims=kdims, vdims=vdims)
-
+            new_type = None if ndims >= min_d and ndims <= max_d else Dataset
+            try:
+                return self.clone(aggregated, kdims=kdims, vdims=vdims,
+                                  new_type=new_type)
+            except:
+                datatype = self.params('datatype').default
+                return self.clone(aggregated, kdims=kdims, vdims=vdims,
+                                  new_type=new_type, datatype=datatype)
 
 
     def groupby(self, dimensions=[], container_type=HoloMap, group_type=None,
diff --git a/holoviews/core/data/array.py b/holoviews/core/data/array.py
index 6f8503db81..35701770a6 100644
--- a/holoviews/core/data/array.py
+++ b/holoviews/core/data/array.py
@@ -185,7 +185,7 @@ def select(cls, dataset, selection_mask=None, **selection):
         data = np.atleast_2d(dataset.data[selection_mask, :])
         if len(data) == 1 and indexed:
             data = data[0, dataset.ndims]
-        return data
+        return data, {}
 
 
     @classmethod
diff --git a/holoviews/core/data/dask.py b/holoviews/core/data/dask.py
index 65ccb40f59..188504ac20 100644
--- a/holoviews/core/data/dask.py
+++ b/holoviews/core/data/dask.py
@@ -126,7 +126,7 @@ def select(cls, columns, selection_mask=None, **selection):
         df = df if selection_mask is None else df[selection_mask]
         if indexed and len(df) == 1:
             return df[columns.vdims[0].alias].compute().iloc[0]
-        return df
+        return df, {}
     
     @classmethod
     def groupby(cls, columns, dimensions, container_type, group_type, **kwargs):
diff --git a/holoviews/core/data/dictionary.py b/holoviews/core/data/dictionary.py
index 031a807549..d4cf64e030 100644
--- a/holoviews/core/data/dictionary.py
+++ b/holoviews/core/data/dictionary.py
@@ -216,7 +216,7 @@ def select(cls, dataset, selection_mask=None, **selection):
                            for k, v in dataset.data.items())
         if indexed and len(list(data.values())[0]) == 1:
             return data[dataset.vdims[0].alias][0]
-        return data
+        return data, {}
 
 
     @classmethod
diff --git a/holoviews/core/data/grid.py b/holoviews/core/data/grid.py
index 250a7a6d90..dde86a2a39 100644
--- a/holoviews/core/data/grid.py
+++ b/holoviews/core/data/grid.py
@@ -273,9 +273,8 @@ def select(cls, dataset, selection_mask=None, **selection):
             data[vdim.alias] = dataset.data[vdim.alias][index]
 
         if indexed and len(data[dataset.vdims[0].alias]) == 1:
-            return data[dataset.vdims[0].alias][0]
-
-        return data
+            return data[dataset.vdims[0].name][0], {}
+        return data, {}
 
 
     @classmethod
diff --git a/holoviews/core/data/image.py b/holoviews/core/data/image.py
new file mode 100644
index 0000000000..02624b66cd
--- /dev/null
+++ b/holoviews/core/data/image.py
@@ -0,0 +1,148 @@
+import numpy as np
+
+from ..boundingregion import BoundingRegion, BoundingBox
+from ..dimension import Dimension 
+from ..ndmapping import OrderedDict
+from ..sheetcoords import SheetCoordinateSystem, Slice
+from .grid import GridInterface
+from .interface import Interface
+
+
+class ImageInterface(GridInterface):
+    """
+    Interface for 2 or 3D arrays representing images
+    of raw luminance values, RGB values or HSV values.
+    """
+
+    types = (np.ndarray,)
+
+    datatype = 'image'
+
+    @classmethod
+    def init(cls, eltype, data, kdims, vdims):
+        if kdims is None:
+            kdims = eltype.kdims
+        if vdims is None:
+            vdims = eltype.vdims
+
+        dimensions = [d.name if isinstance(d, Dimension) else
+                      d for d in kdims + vdims]
+        if not isinstance(data, np.ndarray) or data.ndim != 2:
+            raise ValueError('ImageIntereface expects a 2D array.')
+
+        if kdims is None:
+            kdims = eltype.kdims
+        if vdims is None:
+            vdims = eltype.vdims
+        return data, {'kdims':kdims, 'vdims':vdims}, {}
+
+
+    @classmethod
+    def validate(cls, dataset):
+        pass
+
+
+    @classmethod
+    def range(cls, obj, dim):
+        dim_idx = obj.get_dimension_index(dim)
+        if dim_idx in [0, 1] and obj.bounds:
+            l, b, r, t = obj.bounds.lbrt()
+            if dim_idx:
+                drange = (b, t)
+            else:
+                drange = (l, r)
+        elif 1 < dim_idx < len(obj.vdims) + 2:
+            dim_idx -= 2
+            data = np.atleast_3d(obj.data)[:, :, dim_idx]
+            drange = (np.nanmin(data), np.nanmax(data))
+        else:
+            drange = (None, None)
+        return drange
+
+    
+    @classmethod
+    def values(cls, dataset, dim, expanded=True, flat=True):
+        """
+        The set of samples available along a particular dimension.
+        """
+        dim_idx = dataset.get_dimension_index(dim)
+        if dim_idx in [0, 1]:
+            l, b, r, t = dataset.bounds.lbrt()
+            dim2, dim1 = dataset.data.shape[:2]
+            d1_half_unit = (r - l)/dim1/2.
+            d2_half_unit = (t - b)/dim2/2.
+            d1lin = np.linspace(l+d1_half_unit, r-d1_half_unit, dim1)
+            d2lin = np.linspace(b+d2_half_unit, t-d2_half_unit, dim2)
+            if expanded:
+                values = np.meshgrid(d2lin, d1lin)[abs(dim_idx-1)]
+                return values.flatten() if flat else values
+            else:
+                return d2lin if dim_idx else d1lin
+        elif dim_idx == 2:
+            # Raster arrays are stored with different orientation
+            # than expanded column format, reorient before expanding
+            data = np.flipud(dataset.data)
+            return data.flatten() if flat else data
+        else:
+            return None, None
+
+
+    @classmethod
+    def select(cls, dataset, selection_mask=None, **selection):
+        """
+        Slice the underlying numpy array in sheet coordinates.
+        """
+        selection = {k: slice(*sel) if isinstance(sel, tuple) else sel
+                     for k, sel in selection.items()}
+        coords = tuple(selection[kd.name] if kd.name in selection else slice(None)
+                       for kd in dataset.kdims)
+        if not any([isinstance(el, slice) for el in coords]):
+            data = dataset.data[dataset.sheet2matrixidx(*coords)]
+        xidx, yidx = coords
+        l, b, r, t = dataset.bounds.lbrt()
+        xunit = (1./dataset.xdensity)
+        yunit = (1./dataset.ydensity)    
+        if isinstance(xidx, slice):
+            l = l if xidx.start is None else max(l, xidx.start)
+            r = r if xidx.stop is None else min(r, xidx.stop)
+        if isinstance(yidx, slice):
+            b = b if yidx.start is None else max(b, yidx.start)
+            t = t if yidx.stop is None else min(t, yidx.stop)
+        bounds = BoundingBox(points=((l, b), (r, t)))
+        slc = Slice(bounds, dataset)
+        data = slc.submatrix(dataset.data)
+        l, b, r, t = slc.compute_bounds(dataset).lbrt()
+        if not isinstance(xidx, slice):
+            xc, _ = dataset.closest_cell_center(xidx, b)
+            l, r = xc-xunit/2, xc+xunit/2
+            _, x = dataset.sheet2matrixidx(xidx, b)
+            data = data[:, x][:, np.newaxis]
+        elif not isinstance(yidx, slice):
+            _, yc = dataset.closest_cell_center(l, yidx)
+            b, t = yc-yunit/2, yc+yunit/2
+            y, _ = dataset.sheet2matrixidx(l, yidx)
+            data = data[y, :][np.newaxis, :]
+        bounds = BoundingBox(points=((l, b), (r, t)))
+        return data, {'bounds': bounds}
+
+
+    @classmethod
+    def length(cls, dataset):
+        return np.product(dataset.data.shape)
+    
+
+    @classmethod
+    def aggregate(cls, dataset, kdims, function, **kwargs):
+        kdims = [kd.name if isinstance(kd, Dimension) else kd for kd in kdims]
+        axes = tuple(dataset.ndims-dataset.get_dimension_index(kdim)-1
+                     for kdim in dataset.kdims if kdim not in kdims)
+        
+        data = np.atleast_1d(function(dataset.data, axis=axes, **kwargs))
+        if np.isscalar(data):
+            return data
+        elif len(axes) == 1:
+            return {kdims[0]: cls.values(dataset, axes[0], expanded=False),
+                    dataset.vdims[0].name: data}
+
+
+Interface.register(ImageInterface)
diff --git a/holoviews/core/data/iris.py b/holoviews/core/data/iris.py
index fe9dd7d64d..e337e770c3 100644
--- a/holoviews/core/data/iris.py
+++ b/holoviews/core/data/iris.py
@@ -288,12 +288,12 @@ def select(cls, dataset, selection_mask=None, **selection):
         indexed = cls.indexed(dataset, selection)
         extracted = dataset.data.extract(constraint)
         if indexed and not extracted.dim_coords:
-            return extracted.data.item()
+            return extracted.data.item(), {}
         post_dim_coords = [c.name() for c in extracted.dim_coords]
         dropped = [c for c in pre_dim_coords if c not in post_dim_coords]
         for d in dropped:
             extracted = iris.util.new_axis(extracted, d)
-        return extracted
+        return extracted, {}
 
 
 Interface.register(CubeInterface)
diff --git a/holoviews/core/data/ndelement.py b/holoviews/core/data/ndelement.py
index c05573a1d8..48c898dbb1 100644
--- a/holoviews/core/data/ndelement.py
+++ b/holoviews/core/data/ndelement.py
@@ -118,9 +118,9 @@ def groupby(cls, columns, dimensions, container_type, group_type, **kwargs):
     @classmethod
     def select(cls, columns, selection_mask=None, **selection):
         if selection_mask is None:
-            return columns.data.select(**selection)
+            return columns.data.select(**selection), {}
         else:
-            return columns.data[selection_mask]
+            return columns.data[selection_mask], {}
 
     @classmethod
     def sample(cls, columns, samples=[]):
diff --git a/holoviews/core/data/pandas.py b/holoviews/core/data/pandas.py
index 2653f66f7e..5023d2943d 100644
--- a/holoviews/core/data/pandas.py
+++ b/holoviews/core/data/pandas.py
@@ -188,8 +188,8 @@ def select(cls, columns, selection_mask=None, **selection):
         indexed = cls.indexed(columns, selection)
         df = df.ix[selection_mask]
         if indexed and len(df) == 1:
-            return df[columns.vdims[0].alias].iloc[0]
-        return df
+            return df[columns.vdims[0].alias].iloc[0], {}
+        return df, {}
 
 
     @classmethod
diff --git a/holoviews/core/data/xarray.py b/holoviews/core/data/xarray.py
index 5f61bb6251..8adf9ec058 100644
--- a/holoviews/core/data/xarray.py
+++ b/holoviews/core/data/xarray.py
@@ -224,8 +224,8 @@ def select(cls, dataset, selection_mask=None, **selection):
         indexed = cls.indexed(dataset, selection)
         if (indexed and len(data.data_vars) == 1 and
             len(data[dataset.vdims[0].alias].shape) == 0):
-            return data[dataset.vdims[0].alias].item()
-        return data
+            return data[dataset.vdims[0].alias].item(), {}
+        return data, {}
 
     @classmethod
     def length(cls, dataset):
diff --git a/holoviews/element/raster.py b/holoviews/element/raster.py
index abb424838c..312d523f1d 100644
--- a/holoviews/element/raster.py
+++ b/holoviews/element/raster.py
@@ -393,7 +393,7 @@ def raster(self):
         return self.gridded.dimension_values(2, flat=False)
 
 
-class Image(SheetCoordinateSystem, Raster):
+class Image(Dataset, Element2D, SheetCoordinateSystem):
     """
     Image is the atomic unit as which 2D data is stored, along with
     its bounds object. The input data may be a numpy.matrix object or
@@ -403,45 +403,93 @@ class Image(SheetCoordinateSystem, Raster):
     access to the data, via the .data attribute.
     """
 
-    bounds = param.ClassSelector(class_=BoundingRegion, default=BoundingBox(), doc="""
+    datatype = param.List(default=['image', 'grid', 'xarray'])
+
+    bounds = param.ClassSelector(class_=BoundingRegion, default=None, doc="""
        The bounding region in sheet coordinates containing the data.""")
 
+    kdims = param.List(default=[Dimension('x'), Dimension('y')],
+                       bounds=(2, 2), constant=True, doc="""
+        The label of the x- and y-dimension of the Raster in form
+        of a string or dimension object.""")
+
     group = param.String(default='Image', constant=True)
 
     vdims = param.List(default=[Dimension('z')],
                        bounds=(1, 1), doc="""
         The dimension description of the data held in the matrix.""")
 
-
     def __init__(self, data, bounds=None, extents=None, xdensity=None, ydensity=None, **params):
-        bounds = bounds if bounds is not None else BoundingBox()
-        if np.isscalar(bounds):
+        extents = extents if extents else (None, None, None, None)
+        if data is None: data = np.array([[0]])
+        Dataset.__init__(self, data, extents=extents, bounds=None, **params)
+
+        (dim1, dim2) = self.shape[1], self.shape[0]
+        l, r = self.range(0)
+        b, t = self.range(1)
+        if bounds is None and None in (l, b, r, t):
+            bounds = BoundingBox()
+        if bounds is None:
+            TOL = 10e-9
+            if not xdensity:
+                xvals = self.dimension_values(0, False)
+                xdiff = np.diff(xvals)
+                xunit = np.unique(np.floor(xdiff/TOL).astype(int))*TOL
+                if len(xunit) > 1:
+                    raise Exception('')
+                xdensity = 1./xunit[0]
+            if not ydensity:
+                yvals = self.dimension_values(1, False)
+                ydiff = np.diff(yvals)
+                yunit = np.unique(np.floor(ydiff/TOL).astype(int))*TOL
+                if len(xunit) > 1:
+                    raise Exception('')
+                ydensity = 1./yunit[0]
+            halfx, halfy = 0.5/xdensity, 0.5/ydensity
+            l, r = l-halfx, r+halfx
+            b, t = b-halfy, t+halfy
+            bounds = BoundingBox(points=((l, b), (r, t)))
+        elif np.isscalar(bounds):
             bounds = BoundingBox(radius=bounds)
         elif isinstance(bounds, (tuple, list, np.ndarray)):
             l, b, r, t = bounds
             bounds = BoundingBox(points=((l, b), (r, t)))
-        if data is None: data = np.array([[0]])
-        l, b, r, t = bounds.lbrt()
-        extents = extents if extents else (None, None, None, None)
-        Element2D.__init__(self, data, extents=extents, bounds=bounds,
-                           **params)
 
-        (dim1, dim2) = self.data.shape[1], self.data.shape[0]
+        l, b, r, t = bounds.lbrt()
         xdensity = xdensity if xdensity else dim1/float(r-l)
         ydensity = ydensity if ydensity else dim2/float(t-b)
         SheetCoordinateSystem.__init__(self, bounds, xdensity, ydensity)
 
-        if len(self.data.shape) == 3:
-            if self.data.shape[2] != len(self.vdims):
+        if len(self.shape) == 3:
+            if self.shape[2] != len(self.vdims):
                 raise ValueError("Input array has shape %r but %d value dimensions defined"
-                                 % (self.data.shape, len(self.vdims)))
+                                 % (self.shape, len(self.vdims)))
 
 
-    def _convert_element(self, data):
-        if isinstance(data, (Raster, HeatMap)):
-            return data.data
+    def select(self, selection_specs=None, **selection):
+        """
+        Allows selecting data by the slices, sets and scalar values
+        along a particular dimension. The indices should be supplied as
+        keywords mapping between the selected dimension and
+        value. Additionally selection_specs (taking the form of a list
+        of type.group.label strings, types or functions) may be
+        supplied, which will ensure the selection is only applied if the
+        specs match the selected object.
+        """
+        if selection_specs and not any(self.matches(sp) for sp in selection_specs):
+            return self
+
+        data, kwargs = self.interface.select(self, **selection)
+
+        if np.isscalar(data):
+            return data
         else:
-            return super(Image, self)._convert_element(data)
+            return self.clone(data, xdensity=self.xdensity,
+                              ydensity=self.ydensity, bounds=None, **kwargs)
+
+
+    def _coord2matrix(self, coord):
+        return self.sheet2matrixidx(*coord)
 
 
     def closest(self, coords=[], **kwargs):
@@ -481,95 +529,9 @@ def closest(self, coords=[], **kwargs):
         else:
             return [getter(self.closest_cell_center(*el)) for el in coords]
 
-
-    def __getitem__(self, coords):
-        """
-        Slice the underlying numpy array in sheet coordinates.
-        """
-        if coords in self.dimensions(): return self.dimension_values(coords)
-        coords = util.process_ellipses(self,coords)
-        if coords is () or coords == slice(None, None):
-            return self
-
-        if not isinstance(coords, tuple):
-            coords = (coords, slice(None))
-        if len(coords) > (2 + self.depth):
-            raise KeyError("Can only slice %d dimensions" % 2 + self.depth)
-        elif len(coords) == 3 and coords[-1] not in [self.vdims[0].name, slice(None)]:
-            raise KeyError("%r is the only selectable value dimension" % self.vdims[0].name)
-
-        coords = coords[:2]
-        if not any([isinstance(el, slice) for el in coords]):
-            return self.data[self.sheet2matrixidx(*coords)]
-        if all([isinstance(c, slice) for c in coords]):
-            l, b, r, t = self.bounds.lbrt()
-            xcoords, ycoords = coords
-            xstart = l if xcoords.start is None else max(l, xcoords.start)
-            xend = r if xcoords.stop is None else min(r, xcoords.stop)
-            ystart = b if ycoords.start is None else max(b, ycoords.start)
-            yend = t if ycoords.stop is None else min(t, ycoords.stop)
-            bounds = BoundingBox(points=((xstart, ystart), (xend, yend)))
-        else:
-            raise KeyError('Indexing requires x- and y-slice ranges.')
-
-        return self.clone(Slice(bounds, self).submatrix(self.data),
-                          bounds=bounds)
-
-
-    def range(self, dim, data_range=True):
-        dim_idx = dim if isinstance(dim, int) else self.get_dimension_index(dim)
-        dim = self.get_dimension(dim_idx)
-        if dim.range != (None, None):
-            return dim.range
-        elif dim_idx in [0, 1]:
-            l, b, r, t = self.bounds.lbrt()
-            if dim_idx:
-                drange = (b, t)
-            else:
-                drange = (l, r)
-        elif dim_idx < len(self.vdims) + 2:
-            dim_idx -= 2
-            data = np.atleast_3d(self.data)[:, :, dim_idx]
-            drange = (np.nanmin(data), np.nanmax(data))
-        if data_range:
-            soft_range = [sr for sr in dim.soft_range if sr is not None]
-            if soft_range:
-                return util.max_range([drange, soft_range])
-            else:
-                return drange
-        else:
-            return dim.soft_range
-
-
-    def _coord2matrix(self, coord):
-        return self.sheet2matrixidx(*coord)
-
-
-    def dimension_values(self, dim, expanded=True, flat=True):
-        """
-        The set of samples available along a particular dimension.
-        """
-        dim_idx = self.get_dimension_index(dim)
-        if dim_idx in [0, 1]:
-            l, b, r, t = self.bounds.lbrt()
-            dim2, dim1 = self.data.shape[:2]
-            d1_half_unit = (r - l)/dim1/2.
-            d2_half_unit = (t - b)/dim2/2.
-            d1lin = np.linspace(l+d1_half_unit, r-d1_half_unit, dim1)
-            d2lin = np.linspace(b+d2_half_unit, t-d2_half_unit, dim2)
-            if expanded:
-                values = np.meshgrid(d2lin, d1lin)[abs(dim_idx-1)]
-                return values.flatten() if flat else values
-            else:
-                return d2lin if dim_idx else d1lin
-        elif dim_idx == 2:
-            # Raster arrays are stored with different orientation
-            # than expanded column format, reorient before expanding
-            data = np.flipud(self.data).T
-            return data.flatten() if flat else data
-        else:
-            super(Image, self).dimension_values(dim)
-
+    @property
+    def depth(self):
+        return len(self.vdims)
 
 
 class GridImage(Dataset, Element2D):