diff --git a/docs/source/pythonapi/examples/mgxs-part-iii.ipynb b/docs/source/pythonapi/examples/mgxs-part-iii.ipynb
index a38677945..5f0acde3f 100644
--- a/docs/source/pythonapi/examples/mgxs-part-iii.ipynb
+++ b/docs/source/pythonapi/examples/mgxs-part-iii.ipynb
@@ -260,7 +260,6 @@
    "source": [
     "# Create fuel assembly Lattice\n",
     "assembly = openmc.RectLattice(name='1.6% Fuel Assembly')\n",
-    "assembly.dimension = (17, 17)\n",
     "assembly.pitch = (1.26, 1.26)\n",
     "assembly.lower_left = [-1.26 * 17. / 2.0] * 2"
    ]
@@ -1597,21 +1596,21 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 2",
-   "language": "python2",
-   "name": "python2"
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
-    "version": 2
+    "version": 3
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",
    "name": "python",
    "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython2",
-   "version": "2.7.11"
+   "pygments_lexer": "ipython3",
+   "version": "3.5.1"
   }
  },
  "nbformat": 4,
diff --git a/docs/source/pythonapi/examples/mgxs-part-iv.ipynb b/docs/source/pythonapi/examples/mgxs-part-iv.ipynb
index b330e7ace..65b3f44dc 100644
--- a/docs/source/pythonapi/examples/mgxs-part-iv.ipynb
+++ b/docs/source/pythonapi/examples/mgxs-part-iv.ipynb
@@ -244,7 +244,6 @@
    "source": [
     "# Create fuel assembly Lattice\n",
     "assembly = openmc.RectLattice(name='1.6% Fuel Assembly')\n",
-    "assembly.dimension = (17, 17)\n",
     "assembly.pitch = (1.26, 1.26)\n",
     "assembly.lower_left = [-1.26 * 17. / 2.0] * 2"
    ]
diff --git a/docs/source/pythonapi/examples/pandas-dataframes.ipynb b/docs/source/pythonapi/examples/pandas-dataframes.ipynb
index b88cf9949..2c222ad6e 100644
--- a/docs/source/pythonapi/examples/pandas-dataframes.ipynb
+++ b/docs/source/pythonapi/examples/pandas-dataframes.ipynb
@@ -199,7 +199,6 @@
    "source": [
     "# Create fuel assembly Lattice\n",
     "assembly = openmc.RectLattice(name='1.6% Fuel - 0BA')\n",
-    "assembly.dimension = (17, 17)\n",
     "assembly.pitch = (1.26, 1.26)\n",
     "assembly.lower_left = [-1.26 * 17. / 2.0] * 2\n",
     "assembly.universes = [[pin_cell_universe] * 17] * 17"
@@ -2194,21 +2193,21 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 2",
+   "display_name": "Python 3",
    "language": "python",
-   "name": "python2"
+   "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
-    "version": 2
+    "version": 3
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",
    "name": "python",
    "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython2",
-   "version": "2.7.6"
+   "pygments_lexer": "ipython3",
+   "version": "3.5.1"
   }
  },
  "nbformat": 4,
diff --git a/examples/python/lattice/nested/build-xml.py b/examples/python/lattice/nested/build-xml.py
index b2d611d34..a964d882c 100644
--- a/examples/python/lattice/nested/build-xml.py
+++ b/examples/python/lattice/nested/build-xml.py
@@ -98,14 +98,12 @@ univ4.add_cell(cell2)
 
 # Instantiate nested Lattices
 lattice1 = openmc.RectLattice(lattice_id=4, name='4x4 assembly')
-lattice1.dimension = [2, 2]
 lattice1.lower_left = [-1., -1.]
 lattice1.pitch = [1., 1.]
 lattice1.universes = [[univ1, univ2],
                       [univ2, univ3]]
 
 lattice2 = openmc.RectLattice(lattice_id=6, name='4x4 core')
-lattice2.dimension = [2, 2]
 lattice2.lower_left = [-2., -2.]
 lattice2.pitch = [2., 2.]
 lattice2.universes = [[univ4, univ4],
diff --git a/examples/python/lattice/simple/build-xml.py b/examples/python/lattice/simple/build-xml.py
index 65c355479..4961b96b8 100644
--- a/examples/python/lattice/simple/build-xml.py
+++ b/examples/python/lattice/simple/build-xml.py
@@ -94,7 +94,6 @@ root.add_cell(cell1)
 
 # Instantiate a Lattice
 lattice = openmc.RectLattice(lattice_id=5)
-lattice.dimension = [4, 4]
 lattice.lower_left = [-2., -2.]
 lattice.pitch = [1., 1.]
 lattice.universes = [[univ1, univ2, univ1, univ2],
diff --git a/openmc/cell.py b/openmc/cell.py
index 8f8ebd936..29f85754a 100644
--- a/openmc/cell.py
+++ b/openmc/cell.py
@@ -97,7 +97,10 @@ class Cell(object):
             self.region = region
 
     def __contains__(self, point):
-        return point in self.region
+        if self.region is None:
+            return True
+        else:
+            return point in self.region
 
     def __eq__(self, other):
         if not isinstance(other, Cell):
diff --git a/openmc/geometry.py b/openmc/geometry.py
index ed437f6e1..7eddadfc7 100644
--- a/openmc/geometry.py
+++ b/openmc/geometry.py
@@ -62,6 +62,23 @@ class Geometry(object):
         tree.write("geometry.xml", xml_declaration=True, encoding='utf-8',
                    method="xml")
 
+    def find(self, point):
+        """Find cells/universes/lattices which contain a given point
+
+        Parameters
+        ----------
+        point : 3-tuple of float
+            Cartesian coordinatesof the point
+
+        Returns
+        -------
+        list
+            Sequence of universes, cells, and lattices which are traversed to
+            find the given point
+
+        """
+        return self.root_universe.find(point)
+
     def get_cell_instance(self, path):
         """Return the instance number for the final cell in a geometry path.
 
diff --git a/openmc/lattice.py b/openmc/lattice.py
index af6c14a6a..d8deddc42 100644
--- a/openmc/lattice.py
+++ b/openmc/lattice.py
@@ -1,8 +1,12 @@
+from __future__ import division
+
 import abc
 from collections import OrderedDict, Iterable
+from math import sqrt, floor
 from numbers import Real, Integral
 from xml.etree import ElementTree as ET
 import sys
+import warnings
 
 import numpy as np
 
@@ -113,12 +117,6 @@ class Lattice(object):
         cv.check_type('outer universe', outer, openmc.Universe)
         self._outer = outer
 
-    @universes.setter
-    def universes(self, universes):
-        cv.check_iterable_type('lattice universes', universes, openmc.Universe,
-                               min_depth=2, max_depth=3)
-        self._universes = np.asarray(universes)
-
     def get_unique_universes(self):
         """Determine all unique universes in the lattice
 
@@ -239,6 +237,11 @@ class Lattice(object):
 class RectLattice(Lattice):
     """A lattice consisting of rectangular prisms.
 
+    To completely define a rectangular lattice, the
+    :attr:`RectLattice.lower_left` :attr:`RectLattice.pitch`,
+    :attr:`RectLattice.outer`, and :attr:`RectLattice.universes` properties need
+    to be set.
+
     Parameters
     ----------
     lattice_id : int, optional
@@ -253,12 +256,6 @@ class RectLattice(Lattice):
         Unique identifier for the lattice
     name : str
         Name of the lattice
-    dimension : Iterable of int
-        An array of two or three integers representing the number of lattice
-        cells in the x- and y- (and z-) directions, respectively.
-    lower_left : Iterable of float
-        The coordinates of the lower-left corner of the lattice. If the lattice
-        is two-dimensional, only the x- and y-coordinates are specified.
     pitch : Iterable of float
         Pitch of the lattice in the x, y, and (if applicable) z directions in
         cm.
@@ -266,7 +263,25 @@ class RectLattice(Lattice):
         A universe to fill all space outside the lattice
     universes : Iterable of Iterable of openmc.Universe
         A two- or three-dimensional list/array of universes filling each element
-        of the lattice
+        of the lattice. The first dimension corresponds to the z-direction (if
+        applicable), the second dimension corresponds to the y-direction, and
+        the third dimension corresponds to the x-direction. Note that for the
+        y-direction, a higher index corresponds to a lower physical
+        y-value. Each z-slice in the array can be thought of as a top-down view
+        of the lattice.
+    lower_left : Iterable of float
+        The Cartesian coordinates of the lower-left corner of the lattice. If
+        the lattice is two-dimensional, only the x- and y-coordinates are
+        specified.
+    indices : list of tuple
+        A list of all possible (z,y,x) or (y,x) lattice element indices. These
+        indices correspond to indices in the :attr:`RectLattice.universes`
+        property.
+    ndim : int
+        The number of dimensions of the lattice
+    shape : Iterable of int
+        An array of two or three integers representing the number of lattice
+        cells in the x- and y- (and z-) directions, respectively.
 
     """
 
@@ -274,7 +289,6 @@ class RectLattice(Lattice):
         super(RectLattice, self).__init__(lattice_id, name)
 
         # Initialize Lattice class attributes
-        self._dimension = None
         self._lower_left = None
         self._offsets = None
 
@@ -283,7 +297,7 @@ class RectLattice(Lattice):
             return False
         elif not super(RectLattice, self).__eq__(other):
             return False
-        elif self.dimension != other.dimension:
+        elif self.shape != other.shape:
             return False
         elif self.lower_left != other.lower_left:
             return False
@@ -300,8 +314,8 @@ class RectLattice(Lattice):
         string = 'RectLattice\n'
         string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
         string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
-        string += '{0: <16}{1}{2}\n'.format('\tDimension', '=\t',
-                                            self._dimension)
+        string += '{0: <16}{1}{2}\n'.format('\tShape', '=\t',
+                                            self.shape)
         string += '{0: <16}{1}{2}\n'.format('\tLower Left', '=\t',
                                             self._lower_left)
         string += '{0: <16}{1}{2}\n'.format('\tPitch', '=\t', self._pitch)
@@ -320,7 +334,7 @@ class RectLattice(Lattice):
             string += '{0} '.format(universe._id)
 
             # Add a newline character every time we reach end of row of cells
-            if (i+1) % self._dimension[-1] == 0:
+            if (i+1) % self.shape[0] == 0:
                 string += '\n'
 
         string = string.rstrip('\n')
@@ -333,7 +347,7 @@ class RectLattice(Lattice):
                 string += '{0} '.format(offset)
 
                 # Add a newline character when we reach end of row of cells
-                if (i+1) % self._dimension[-1] == 0:
+                if (i+1) % self.shape[0] == 0:
                     string += '\n'
 
             string = string.rstrip('\n')
@@ -341,24 +355,29 @@ class RectLattice(Lattice):
         return string
 
     @property
-    def dimension(self):
-        return self._dimension
+    def indices(self):
+        if self.ndim == 2:
+            return list(np.broadcast(*np.ogrid[
+                :self.shape[1], :self.shape[0]]))
+        else:
+            return list(np.broadcast(*np.ogrid[
+                :self.shape[2], :self.shape[1], :self.shape[0]]))
 
     @property
     def lower_left(self):
         return self._lower_left
 
+    @property
+    def ndim(self):
+        return len(self.pitch)
+
     @property
     def offsets(self):
         return self._offsets
 
-    @dimension.setter
-    def dimension(self, dimension):
-        cv.check_type('lattice dimension', dimension, Iterable, Integral)
-        cv.check_length('lattice dimension', dimension, 2, 3)
-        for dim in dimension:
-            cv.check_greater_than('lattice dimension', dim, 0)
-        self._dimension = dimension
+    @property
+    def shape(self):
+        return self._universes.shape[::-1]
 
     @lower_left.setter
     def lower_left(self, lower_left):
@@ -379,8 +398,13 @@ class RectLattice(Lattice):
             cv.check_greater_than('lattice pitch', dim, 0.0)
         self._pitch = pitch
 
-    def get_cell_instance(self, path, distribcell_index):
+    @Lattice.universes.setter
+    def universes(self, universes):
+        cv.check_iterable_type('lattice universes', universes, openmc.Universe,
+                               min_depth=2, max_depth=3)
+        self._universes = np.asarray(universes)
 
+    def get_cell_instance(self, path, distribcell_index):
         # Extract the lattice element from the path
         next_index = path.index('-')
         lat_id_indices = path[:next_index]
@@ -395,7 +419,7 @@ class RectLattice(Lattice):
         lat_z = int(i.split(',')[2]) - 1
 
         # For 2D Lattices
-        if len(self._dimension) == 2:
+        if self.ndim == 2:
             offset = self._offsets[lat_z, lat_y, lat_x, distribcell_index-1]
             offset += self._universes[lat_x][lat_y].get_cell_instance(path,
                                                               distribcell_index)
@@ -408,6 +432,128 @@ class RectLattice(Lattice):
 
         return offset
 
+    def find_element(self, point):
+        """Determine index of lattice element and local coordinates for a point
+
+        Parameters
+        ----------
+        point : Iterable of float
+            Cartesian coordinates of point
+
+        Returns
+        -------
+        2- or 3-tuple of int
+            A tuple of the corresponding (x,y,z) lattice element indices
+        3-tuple of float
+            Carestian coordinates of the point in the corresponding lattice
+            element coordinate system
+
+        """
+        ix = floor((point[0] - self._lower_left[0])/self._pitch[0])
+        iy = floor((point[1] - self._lower_left[1])/self._pitch[1])
+        if self.ndim == 2:
+            idx = (ix, iy)
+        else:
+            iz = floor((point[2] - self._lower_left[2])/self._pitch[2])
+            idx = (ix, iy, iz)
+        return idx, self.get_local_coordinates(point, idx)
+
+    def get_local_coordinates(self, point, idx):
+        """Determine local coordinates of a point within a lattice element
+
+        Parameters
+        ----------
+        point : Iterable of float
+            Cartesian coordinates of point
+        idx : Iterable of int
+            (x,y,z) indices of lattice element. If the lattice is 2D, the z
+            index can be omitted.
+
+        Returns
+        -------
+        3-tuple of float
+            Cartesian coordinates of point in the lattice element coordinate
+            system
+
+        """
+        x = point[0] - (self._lower_left[0] + (idx[0] + 0.5)*self._pitch[0])
+        y = point[1] - (self._lower_left[1] + (idx[1] + 0.5)*self._pitch[1])
+        if self.ndim == 2:
+            z = point[2]
+        else:
+            z = point[2] - (self._lower_left[2] + (idx[2] + 0.5)*self._pitch[2])
+        return (x, y, z)
+
+    def get_universe_index(self, idx):
+        """Return index in the universes array corresponding to a lattice element index
+
+        Parameters
+        ----------
+        idx : Iterable of int
+            Lattice element indices in the :math:`(x,y,z)` coordinate system
+
+        Returns
+        -------
+        2- or 3-tuple of int
+            Indices used when setting the :attr:`RectLattice.universes` property
+
+        """
+        max_y = self.shape[1] - 1
+        if self.ndim == 2:
+            x, y = idx
+            return (max_y - y, x)
+        else:
+            x, y, z = idx
+            return (z, max_y - y, x)
+
+    def is_valid_index(self, idx):
+        """Determine whether lattice element index is within defined range
+
+        Parameters
+        ----------
+        idx : Iterable of int
+            Lattice element indices in the :math:`(x,y,z)` coordinate system
+
+        Returns
+        -------
+        bool
+            Whether index is valid
+
+        """
+        if self.ndim == 2:
+            return (0 <= idx[0] < self.shape[0] and
+                    0 <= idx[1] < self.shape[1])
+        else:
+            return (0 <= idx[0] < self.shape[0] and
+                    0 <= idx[1] < self.shape[1] and
+                    0 <= idx[2] < self.shape[2])
+
+    def find(self, point):
+        """Find cells/universes/lattices which contain a given point
+
+        Parameters
+        ----------
+        point : 3-tuple of float
+            Cartesian coordinatesof the point
+
+        Returns
+        -------
+        list
+            Sequence of universes, cells, and lattices which are traversed to
+            find the given point
+
+        """
+        idx, p = self.find_element(point)
+        if self.is_valid_index(idx):
+            idx_u = self.get_universe_index(idx)
+            u = self.universes[idx_u]
+        else:
+            if self.outer is not None:
+                u = self.outer
+            else:
+                return []
+        return [(self, idx)] + u.find(p)
+
     def create_xml_subelement(self, xml_element):
 
         # Determine if XML element already contains subelement for this Lattice
@@ -436,7 +582,7 @@ class RectLattice(Lattice):
 
         # Export Lattice cell dimensions
         dimension = ET.SubElement(lattice_subelement, "dimension")
-        dimension.text = ' '.join(map(str, self._dimension))
+        dimension.text = ' '.join(map(str, self.shape))
 
         # Export Lattice lower left
         lower_left = ET.SubElement(lattice_subelement, "lower_left")
@@ -446,10 +592,10 @@ class RectLattice(Lattice):
         universe_ids = '\n'
 
         # 3D Lattices
-        if len(self._dimension) == 3:
-            for z in range(self._dimension[2]):
-                for y in range(self._dimension[1]):
-                    for x in range(self._dimension[0]):
+        if self.ndim == 3:
+            for z in range(self.shape[2]):
+                for y in range(self.shape[1]):
+                    for x in range(self.shape[0]):
                         universe = self._universes[z][y][x]
 
                         # Append Universe ID to the Lattice XML subelement
@@ -466,8 +612,8 @@ class RectLattice(Lattice):
 
         # 2D Lattices
         else:
-            for y in range(self._dimension[1]):
-                for x in range(self._dimension[0]):
+            for y in range(self.shape[1]):
+                for x in range(self.shape[0]):
                     universe = self._universes[y][x]
 
                     # Append Universe ID to Lattice XML subelement
@@ -492,6 +638,10 @@ class RectLattice(Lattice):
 class HexLattice(Lattice):
     """A lattice consisting of hexagonal prisms.
 
+    To completely define a hexagonal lattice, the :attr:`HexLattice.center`,
+    :attr:`HexLattice.pitch`, :attr:`HexLattice.universes`, and
+    :attr:`HexLattice.outer` properties need to be set.
+
     Parameters
     ----------
     lattice_id : int, optional
@@ -506,26 +656,31 @@ class HexLattice(Lattice):
         Unique identifier for the lattice
     name : str
         Name of the lattice
-    num_rings : int
-        Number of radial ring positions in the xy-plane
-    num_axial : int
-        Number of positions along the z-axis.
-    center : Iterable of float
-        Coordinates of the center of the lattice. If the lattice does not have
-        axial sections then only the x- and y-coordinates are specified
     pitch : Iterable of float
         Pitch of the lattice in cm. The first item in the iterable specifies the
         pitch in the radial direction and, if the lattice is 3D, the second item
         in the iterable specifies the pitch in the axial direction.
     outer : openmc.Universe
         A universe to fill all space outside the lattice
-    universes : Iterable of Iterable of openmc.Universe
+    universes : Nested Iterable of openmc.Universe
         A two- or three-dimensional list/array of universes filling each element
         of the lattice. Each sub-list corresponds to one ring of universes and
         should be ordered from outermost ring to innermost ring. The universes
         within each sub-list are ordered from the "top" and proceed in a
         clockwise fashion. The :meth:`HexLattice.show_indices` method can be
         used to help figure out indices for this property.
+    center : Iterable of float
+        Coordinates of the center of the lattice. If the lattice does not have
+        axial sections then only the x- and y-coordinates are specified
+    indices : list of tuple
+        A list of all possible (z,r,i) or (r,i) lattice element indices that are
+        possible, where z is the axial index, r is in the ring index (starting
+        from the outermost ring), and i is the index with a ring starting from
+        the top and proceeding clockwise.
+    num_rings : int
+        Number of radial ring positions in the xy-plane
+    num_axial : int
+        Number of positions along the z-axis.
 
     """
 
@@ -597,17 +752,15 @@ class HexLattice(Lattice):
     def center(self):
         return self._center
 
-    @num_rings.setter
-    def num_rings(self, num_rings):
-        cv.check_type('number of rings', num_rings, Integral)
-        cv.check_greater_than('number of rings', num_rings, 0)
-        self._num_rings = num_rings
-
-    @num_axial.setter
-    def num_axial(self, num_axial):
-        cv.check_type('number of axial', num_axial, Integral)
-        cv.check_greater_than('number of axial', num_axial, 0)
-        self._num_axial = num_axial
+    @property
+    def indices(self):
+        if self.num_axial is None:
+            return [(r, i) for r in range(self._num_rings)
+                    for i in range(max(6*(self._num_rings - 1 - r), 1))]
+        else:
+            return [(z, r, i) for z in range(self._num_axial)
+                    for r in range(self._num_rings)
+                    for i in range(max(6*(self._num_rings - 1 - r), 1))]
 
     @center.setter
     def center(self, center):
@@ -625,8 +778,9 @@ class HexLattice(Lattice):
 
     @Lattice.universes.setter
     def universes(self, universes):
-        # Call Lattice.universes parent class setter property
-        Lattice.universes.fset(self, universes)
+        cv.check_iterable_type('lattice universes', universes, openmc.Universe,
+                               min_depth=2, max_depth=3)
+        self._universes = universes
 
         # NOTE: This routine assumes that the user creates a "ragged" list of
         # lists, where each sub-list corresponds to one ring of Universes.
@@ -649,14 +803,14 @@ class HexLattice(Lattice):
 
         # Set the number of axial positions.
         if n_dims == 3:
-            self.num_axial = len(self._universes)
+            self._num_axial = len(self._universes)
         else:
             self._num_axial = None
 
         # Set the number of rings and make sure this number is consistent for
         # all axial positions.
         if n_dims == 3:
-            self.num_rings = len(self._universes[0])
+            self._num_rings = len(self._universes[0])
             for rings in self._universes:
                 if len(rings) != self._num_rings:
                     msg = 'HexLattice ID={0:d} has an inconsistent number of ' \
@@ -664,7 +818,7 @@ class HexLattice(Lattice):
                     raise ValueError(msg)
 
         else:
-            self.num_rings = len(self._universes)
+            self._num_rings = len(self._universes)
 
         # Make sure there are the correct number of elements in each ring.
         if n_dims == 3:
@@ -705,6 +859,170 @@ class HexLattice(Lattice):
                                              6*(self._num_rings - 1 - r))
                     raise ValueError(msg)
 
+    def find_element(self, point):
+        """Determine index of lattice element and local coordinates for a point
+
+        Parameters
+        ----------
+        point : Iterable of float
+            Cartesian coordinates of point
+
+        Returns
+        -------
+        3-tuple of int
+            Indices of corresponding lattice element in (x,:math:`alpha`,z)
+            bases
+        numpy.ndarray
+            Carestian coordinates of the point in the corresponding lattice
+            element coordinate system
+
+        """
+        # Convert coordinates to skewed bases
+        x = point[0] - self._center[0]
+        y = point[1] - self._center[1]
+        if self._num_axial is None:
+            iz = 1
+        else:
+            z = point[2] - self._center[2]
+            iz = floor(z/self._pitch[1] + 0.5*self._num_axial)
+        alpha = y - x/sqrt(3.)
+        ix = floor(x/(sqrt(0.75) * self._pitch[0]))
+        ia = floor(alpha/self._pitch[0])
+
+        # Check four lattice elements to see which one is closest based on local
+        # coordinates
+        d_min = np.inf
+        for idx in [(ix, ia, iz), (ix + 1, ia, iz), (ix, ia + 1, iz),
+                  (ix + 1, ia + 1, iz)]:
+            p = self.get_local_coordinates(point, idx)
+            d = p[0]**2 + p[1]**2
+            if d < d_min:
+                d_min = d
+                idx_min = idx
+                p_min = p
+
+        return idx_min, p_min
+
+    def get_local_coordinates(self, point, idx):
+        """Determine local coordinates of a point within a lattice element
+
+        Parameters
+        ----------
+        point : Iterable of float
+            Cartesian coordinates of point
+        idx : Iterable of int
+            Indices of lattice element in (x,:math:`alpha`,z) bases
+
+        Returns
+        -------
+        3-tuple of float
+            Cartesian coordinates of point in the lattice element coordinate
+            system
+
+        """
+        x = point[0] - (self._center[0] + sqrt(0.75)*self._pitch[0]*idx[0])
+        y = point[1] - (self._center[1] + (0.5*idx[0] + idx[1])*self._pitch[0])
+        if self._num_axial is None:
+            z = point[2]
+        else:
+            z = point[2] - (self._center[2] + (idx[2] + 0.5 - 0.5*self._num_axial)*
+                            self._pitch[1])
+        return (x, y, z)
+
+    def get_universe_index(self, idx):
+        """Return index in the universes array corresponding to a lattice element index
+
+        Parameters
+        ----------
+        idx : Iterable of int
+            Lattice element indices in the :math:`(x,\alpha,z)` coordinate
+            system
+
+        Returns
+        -------
+        2- or 3-tuple of int
+            Indices used when setting the :attr:`HexLattice.universes` property
+
+        """
+
+        # First we determine which ring the index corresponds to.
+        x = idx[0]
+        a = idx[1]
+        z = -a - x
+        g = max(abs(x), abs(a), abs(z))
+
+        # Next we use a clever method to figure out where along the ring we are.
+        i_ring = self._num_rings - 1 - g
+        if x >= 0:
+            if a >= 0:
+                i_within = x
+            else:
+                i_within = 2*g + z
+        else:
+            if a <= 0:
+                i_within = 3*g - x
+            else:
+                i_within = 5*g - z
+
+        if self.num_axial is None:
+            return (i_ring, i_within)
+        else:
+            return (idx[2], i_ring, i_within)
+
+    def is_valid_index(self, idx):
+        """Determine whether lattice element index is within defined range
+
+        Parameters
+        ----------
+        idx : Iterable of int
+            Lattice element indices in the :math:`(x,\alpha,z)` coordinate
+            system
+
+        Returns
+        -------
+        bool
+            Whether index is valid
+
+        """
+        x = idx[0]
+        y = idx[1]
+        z = 0 - y - x
+        g = max(abs(x), abs(y), abs(z))
+        if self.num_axial is None:
+            return g < self.num_rings
+        else:
+            return g < self.num_rings and 0 <= idx[2] < self.num_axial
+
+    def find(self, point):
+        """Find cells/universes/lattices which contain a given point
+
+        Parameters
+        ----------
+        point : 3-tuple of float
+            Cartesian coordinatesof the point
+
+        Returns
+        -------
+        list
+            Sequence of universes, cells, and lattices which are traversed to
+            find the given point
+
+        """
+        idx, p = self.find_element(point)
+        if self.is_valid_index(idx):
+            idx_u = self.get_universe_index(idx)
+            if self.num_axial is None:
+                u = self.universes[idx_u[0]][idx_u[1]]
+            else:
+                u = self.universes[idx_u[0]][idx_u[1]][idx_u[2]]
+        else:
+            if self.outer is not None:
+                u = self.outer
+            else:
+                return []
+
+        return [(self, idx)] + u.find(p)
+
     def create_xml_subelement(self, xml_element):
         # Determine if XML element already contains subelement for this Lattice
         path = './hex_lattice[@id=\'{0}\']'.format(self._id)
@@ -736,8 +1054,8 @@ class HexLattice(Lattice):
             lattice_subelement.set("n_axial", str(self._num_axial))
 
         # Export Lattice cell center
-        dimension = ET.SubElement(lattice_subelement, "center")
-        dimension.text = ' '.join(map(str, self._center))
+        center = ET.SubElement(lattice_subelement, "center")
+        center.text = ' '.join(map(str, self._center))
 
         # Export the Lattice nested Universe IDs.
 
diff --git a/openmc/opencg_compatible.py b/openmc/opencg_compatible.py
index 562fe9cad..c6b364d7c 100644
--- a/openmc/opencg_compatible.py
+++ b/openmc/opencg_compatible.py
@@ -861,8 +861,8 @@ def get_opencg_lattice(openmc_lattice):
         universes = new_universes
 
     # Initialize an empty array for the OpenCG nested Universes in this Lattice
-    universe_array = np.ndarray(tuple(np.array(dimension)[::-1]),
-                                dtype=opencg.Universe)
+    universe_array = np.empty(tuple(np.array(dimension)[::-1]),
+                              dtype=opencg.Universe)
 
     # Create OpenCG Universes for each unique nested Universe in this Lattice
     unique_universes = openmc_lattice.get_unique_universes()
@@ -929,8 +929,8 @@ def get_openmc_lattice(opencg_lattice):
     outer = opencg_lattice.outside
 
     # Initialize an empty array for the OpenMC nested Universes in this Lattice
-    universe_array = np.ndarray(tuple(np.array(dimension)[::-1]),
-                                dtype=openmc.Universe)
+    universe_array = np.empty(tuple(np.array(dimension)[::-1]),
+                              dtype=openmc.Universe)
 
     # Create OpenMC Universes for each unique nested Universe in this Lattice
     unique_universes = opencg_lattice.get_unique_universes()
@@ -953,7 +953,6 @@ def get_openmc_lattice(opencg_lattice):
                    np.array(dimension, dtype=np.float64))) / -2.0
 
     openmc_lattice = openmc.RectLattice(lattice_id=lattice_id)
-    openmc_lattice.dimension = dimension
     openmc_lattice.pitch = width
     openmc_lattice.universes = universe_array
     openmc_lattice.lower_left = lower_left
diff --git a/openmc/summary.py b/openmc/summary.py
index 6fb6b7000..c3277809a 100644
--- a/openmc/summary.py
+++ b/openmc/summary.py
@@ -358,7 +358,6 @@ class Summary(object):
 
                 # Create the Lattice
                 lattice = openmc.RectLattice(lattice_id=lattice_id, name=name)
-                lattice.dimension = tuple(dimension)
                 lattice.lower_left = lower_left
                 lattice.pitch = pitch
 
@@ -368,7 +367,7 @@ class Summary(object):
 
                 # Build array of Universe pointers for the Lattice
                 universes = \
-                    np.ndarray(tuple(universe_ids.shape), dtype=openmc.Universe)
+                    np.empty(tuple(universe_ids.shape), dtype=openmc.Universe)
 
                 for z in range(universe_ids.shape[0]):
                     for y in range(universe_ids.shape[1]):
@@ -403,8 +402,6 @@ class Summary(object):
 
                 # Create the Lattice
                 lattice = openmc.HexLattice(lattice_id=lattice_id, name=name)
-                lattice.num_rings = n_rings
-                lattice.num_axial = n_axial
                 lattice.center = center
                 lattice.pitch = pitch
 
@@ -417,12 +414,12 @@ class Summary(object):
                 # (x, alpha, z) to the Python API's format of a ragged nested
                 # list of (z, ring, theta).
                 universes = []
-                for z in range(lattice.num_axial):
+                for z in range(n_axial):
                     # Add a list for this axial level.
                     universes.append([])
-                    x = lattice.num_rings - 1
-                    a = 2*lattice.num_rings - 2
-                    for r in range(lattice.num_rings - 1, 0, -1):
+                    x = n_rings - 1
+                    a = 2*n_rings - 2
+                    for r in range(n_rings - 1, 0, -1):
                         # Add a list for this ring.
                         universes[-1].append([])
 
diff --git a/openmc/universe.py b/openmc/universe.py
index 770e789da..d8e1c4dac 100644
--- a/openmc/universe.py
+++ b/openmc/universe.py
@@ -1,6 +1,7 @@
 from collections import OrderedDict, Iterable
 from numbers import Integral
 from xml.etree import ElementTree as ET
+import random
 import sys
 import warnings
 
@@ -124,6 +125,110 @@ class Universe(object):
         else:
             self._name = ''
 
+    def find(self, point):
+        """Find cells/universes/lattices which contain a given point
+
+        Parameters
+        ----------
+        point : 3-tuple of float
+            Cartesian coordinatesof the point
+
+        Returns
+        -------
+        list
+            Sequence of universes, cells, and lattices which are traversed to
+            find the given point
+
+        """
+        p = np.asarray(point)
+        for cell in self._cells.values():
+            if p in cell:
+                if cell._type in ('normal', 'void'):
+                    return [self, cell]
+                elif cell._type == 'fill':
+                    if cell.translation is not None:
+                        p -= cell.translation
+                    if cell.rotation is not None:
+                        p[:] = cell.rotation_matrix.dot(p)
+                    return [self, cell] + cell.fill.find(p)
+                else:
+                    return [self, cell] + cell.fill.find(p)
+        return []
+
+    def plot(self, center=(0., 0., 0.), width=(1., 1.), pixels=(200, 200),
+             basis='xy', color_by='cell'):
+        """Display a slice plot of the universe.
+
+        Parameters
+        ----------
+        center : Iterable of float
+            Coordinates at the center of the plot
+        width : Iterable of float
+            Width of the plot in each basis direction
+        pixels : Iterable of int
+            Number of pixels to use in each basis direction
+        basis : {'xy', 'xz', 'yz'}
+            The basis directions for the plot
+        color_by : {'cell', 'material'}
+            Indicate whether the plot should be colored by cell or by material
+
+        """
+        import matplotlib.pyplot as plt
+
+        if basis == 'xy':
+            x_min = center[0] - 0.5*width[0]
+            x_max = center[0] + 0.5*width[0]
+            y_min = center[1] - 0.5*width[1]
+            y_max = center[1] + 0.5*width[1]
+        elif basis == 'yz':
+            # The x-axis will correspond to physical y and the y-axis will correspond to physical z
+            x_min = center[1] - 0.5*width[0]
+            x_max = center[1] + 0.5*width[0]
+            y_min = center[2] - 0.5*width[1]
+            y_max = center[2] + 0.5*width[1]
+        elif basis == 'xz':
+            # The y-axis will correspond to physical z
+            x_min = center[0] - 0.5*width[0]
+            x_max = center[0] + 0.5*width[0]
+            y_min = center[2] - 0.5*width[1]
+            y_max = center[2] + 0.5*width[1]
+
+        # Determine locations to determine cells at
+        x_coords = np.linspace(x_min, x_max, pixels[0], endpoint=False) + \
+                   0.5*(x_max - x_min)/pixels[0]
+        y_coords = np.linspace(y_max, y_min, pixels[1], endpoint=False) - \
+                   0.5*(y_max - y_min)/pixels[1]
+
+        colors = {}
+        img = np.zeros(pixels + (4,))  # Use RGBA form
+        for i, x in enumerate(x_coords):
+            for j, y in enumerate(y_coords):
+                if basis == 'xy':
+                    path = self.find((x, y, center[2]))
+                elif basis == 'yz':
+                    path = self.find((center[0], x, y))
+                elif basis == 'xz':
+                    path = self.find((x, center[1], y))
+
+                if len(path) > 0:
+                    try:
+                        if color_by == 'cell':
+                            uid = path[-1].id
+                        elif color_by == 'material':
+                            if path[-1].fill_type == 'material':
+                                uid = path[-1].fill.id
+                            else:
+                                continue
+                    except AttributeError:
+                        continue
+                    if uid not in colors:
+                        colors[uid] = (random.random(), random.random(),
+                                       random.random(), 1.0)
+                    img[j,i,:] = colors[uid]
+
+        plt.imshow(img, extent=(x_min, x_max, y_min, y_max))
+        plt.show()
+
     def add_cell(self, cell):
         """Add a cell to the universe.
 
diff --git a/tests/input_set.py b/tests/input_set.py
index 2c6841e25..fefd2ca5f 100644
--- a/tests/input_set.py
+++ b/tests/input_set.py
@@ -350,7 +350,6 @@ class InputSet(object):
         # Define fuel lattices.
         l100 = openmc.RectLattice(name='Fuel assembly (lower half)',
                                   lattice_id=100)
-        l100.dimension = (17, 17)
         l100.lower_left = (-10.71, -10.71)
         l100.pitch = (1.26, 1.26)
         l100.universes = [
@@ -384,7 +383,6 @@ class InputSet(object):
 
         l101 = openmc.RectLattice(name='Fuel assembly (upper half)',
                                   lattice_id=101)
-        l101.dimension = (17, 17)
         l101.lower_left = (-10.71, -10.71)
         l101.pitch = (1.26, 1.26)
         l101.universes = [
@@ -444,7 +442,6 @@ class InputSet(object):
         # Define core lattices
         l200 = openmc.RectLattice(name='Core lattice (lower half)',
                                   lattice_id=200)
-        l200.dimension = (21, 21)
         l200.lower_left = (-224.91, -224.91)
         l200.pitch = (21.42, 21.42)
         l200.universes = [
@@ -472,7 +469,6 @@ class InputSet(object):
 
         l201 = openmc.RectLattice(name='Core lattice (lower half)',
                                   lattice_id=201)
-        l201.dimension = (21, 21)
         l201.lower_left = (-224.91, -224.91)
         l201.pitch = (21.42, 21.42)
         l201.universes = [
diff --git a/tests/test_asymmetric_lattice/test_asymmetric_lattice.py b/tests/test_asymmetric_lattice/test_asymmetric_lattice.py
index 504cc4746..40def1f62 100644
--- a/tests/test_asymmetric_lattice/test_asymmetric_lattice.py
+++ b/tests/test_asymmetric_lattice/test_asymmetric_lattice.py
@@ -24,7 +24,6 @@ class AsymmetricLatticeTestHarness(PyAPITestHarness):
 
         # Construct a 3x3 lattice of fuel assemblies
         core_lat = openmc.RectLattice(name='3x3 Core Lattice', lattice_id=202)
-        core_lat.dimension = (3, 3)
         core_lat.lower_left = (-32.13, -32.13)
         core_lat.pitch = (21.42, 21.42)
         core_lat.universes = [[fuel,  water, water],
diff --git a/tests/test_distribmat/test_distribmat.py b/tests/test_distribmat/test_distribmat.py
index d8f78c5cf..6700a96b6 100644
--- a/tests/test_distribmat/test_distribmat.py
+++ b/tests/test_distribmat/test_distribmat.py
@@ -53,7 +53,6 @@ class DistribmatTestHarness(PyAPITestHarness):
         fuel_univ.add_cells((c11, c12))
 
         lat = openmc.RectLattice(lattice_id=101)
-        lat.dimension = [2, 2]
         lat.lower_left = [-2.0, -2.0]
         lat.pitch = [2.0, 2.0]
         lat.universes = [[fuel_univ]*2]*2