diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py
index ef3e90fc41..ca9fa5d669 100644
--- a/openmc/mgxs/library.py
+++ b/openmc/mgxs/library.py
@@ -73,6 +73,9 @@ class Library(object):
     name : str, optional
         Name of the multi-group cross section library. Used as a label to
         identify tallies in OpenMC 'tallies.xml' file.
+    sparse : bool
+        Whether or not the Library's tallies use SciPy's LIL sparse matrix
+        format for compressed data storage
 
     """
 
@@ -92,6 +95,7 @@ class Library(object):
         self._all_mgxs = OrderedDict()
         self._sp_filename = None
         self._keff = None
+        self._sparse = False
 
         self.name = name
         self.openmc_geometry = openmc_geometry
@@ -119,6 +123,7 @@ class Library(object):
             clone._all_mgxs = self.all_mgxs
             clone._sp_filename = self._sp_filename
             clone._keff = self._keff
+            clone._sparse = self.sparse
 
             clone._all_mgxs = OrderedDict()
             for domain in self.domains:
@@ -208,6 +213,10 @@ class Library(object):
     def keff(self):
         return self._keff
 
+    @property
+    def sparse(self):
+        return self._sparse
+
     @openmc_geometry.setter
     def openmc_geometry(self, openmc_geometry):
         cv.check_type('openmc_geometry', openmc_geometry, openmc.Geometry)
@@ -286,6 +295,28 @@ class Library(object):
         cv.check_type('tally trigger', tally_trigger, openmc.Trigger)
         self._tally_trigger = tally_trigger
 
+    @sparse.setter
+    def sparse(self, sparse):
+        """Convert tally data from NumPy arrays to SciPy list of lists (LIL)
+        sparse matrices, and vice versa.
+
+        This property may be used to reduce the amount of data in memory during
+        tally data processing. The tally data will be stored as SciPy LIL
+        matrices internally within the Tally object. All tally data access
+        properties and methods will return data as a dense NumPy array.
+
+        """
+
+        cv.check_type('sparse', sparse, bool)
+        
+        # Sparsify or densify each MGXS in the Library
+        for domain in self.domains:
+            for mgxs_type in self.mgxs_types:
+                mgxs = self.get_mgxs(domain, mgxs_type)
+                mgxs.sparse = self.sparse
+
+        self._sparse = sparse
+
     def build_library(self):
         """Initialize MGXS objects in each domain and for each reaction type
         in the library.
@@ -381,6 +412,7 @@ class Library(object):
             for mgxs_type in self.mgxs_types:
                 mgxs = self.get_mgxs(domain, mgxs_type)
                 mgxs.load_from_statepoint(statepoint)
+                mgxs.sparse = self.sparse
 
     def get_mgxs(self, domain, mgxs_type):
         """Return the MGXS object for some domain and reaction rate type.
diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py
index 7f521a975c..968504a06d 100644
--- a/openmc/mgxs/mgxs.py
+++ b/openmc/mgxs/mgxs.py
@@ -103,6 +103,9 @@ class MGXS(object):
     nuclides : list of str or 'sum'
         A list of nuclide string names (e.g., 'U-238', 'O-16') when by_nuclide
         is True and 'sum' when by_nuclide is False.
+    sparse : bool
+        Whether or not the MGXS' tallies use SciPy's LIL sparse matrix format
+        for compressed data storage
 
     """
 
@@ -121,6 +124,7 @@ class MGXS(object):
         self._tally_trigger = None
         self._tallies = None
         self._xs_tally = None
+        self._sparse = False
 
         self.name = name
         self.by_nuclide = by_nuclide
@@ -146,6 +150,7 @@ class MGXS(object):
             clone._energy_groups = copy.deepcopy(self.energy_groups, memo)
             clone._tally_trigger = copy.deepcopy(self.tally_trigger, memo)
             clone._xs_tally = copy.deepcopy(self.xs_tally, memo)
+            clone._sparse = self.sparse
 
             clone._tallies = OrderedDict()
             for tally_type, tally in self.tallies.items():
@@ -199,6 +204,10 @@ class MGXS(object):
     def xs_tally(self):
         return self._xs_tally
 
+    @property
+    def sparse(self):
+        return self._sparse
+
     @property
     def num_subdomains(self):
         tally = list(self.tallies.values())[0]
@@ -249,6 +258,29 @@ class MGXS(object):
         cv.check_type('tally trigger', tally_trigger, openmc.Trigger)
         self._tally_trigger = tally_trigger
 
+    @sparse.setter
+    def sparse(self, sparse):
+        """Convert tally data from NumPy arrays to SciPy list of lists (LIL)
+        sparse matrices, and vice versa.
+
+        This property may be used to reduce the amount of data in memory during
+        tally data processing. The tally data will be stored as SciPy LIL
+        matrices internally within the Tally object. All tally data access
+        properties and methods will return data as a dense NumPy array.
+
+        """
+
+        cv.check_type('sparse', sparse, bool)
+
+        # Sparsify or densify the derived MGXS tally and its base tallies
+        if self.xs_tally:
+            self.xs_tally.sparse = sparse
+
+        for tally_name in self.tallies:
+                self.tallies[tally_name].sparse = sparse
+
+        self._sparse = sparse
+
     @staticmethod
     def get_mgxs(mgxs_type, domain=None, domain_type=None,
                  energy_groups=None, by_nuclide=False, name=''):
@@ -503,6 +535,7 @@ class MGXS(object):
         # Remove NaNs which may have resulted from divide-by-zero operations
         self._xs_tally._mean = np.nan_to_num(self.xs_tally.mean)
         self._xs_tally._std_dev = np.nan_to_num(self.xs_tally.std_dev)
+        self.xs_tally.sparse = self.sparse
 
     def load_from_statepoint(self, statepoint):
         """Extracts tallies in an OpenMC StatePoint with the data needed to
@@ -565,6 +598,7 @@ class MGXS(object):
                                             estimator=tally.estimator)
             sp_tally = sp_tally.get_slice(tally.scores, filters,
                                           filter_bins, tally.nuclides)
+            sp_tally.sparse = self.sparse
             self.tallies[tally_type] = sp_tally
 
         # Compute the cross section from the tallies
@@ -716,6 +750,7 @@ class MGXS(object):
 
         # Clone this MGXS to initialize the condensed version
         condensed_xs = copy.deepcopy(self)
+        condensed_xs.sparse = False
         condensed_xs.energy_groups = coarse_groups
 
         # Build energy indices to sum across
@@ -754,10 +789,8 @@ class MGXS(object):
                     std_dev = np.sqrt(std_dev)
 
             # Reshape condensed data arrays with one dimension for all filters
-            new_shape = \
-                (tally.num_filter_bins, tally.num_nuclides, tally.num_scores,)
-            mean = np.reshape(mean, new_shape)
-            std_dev = np.reshape(std_dev, new_shape)
+            mean = np.reshape(mean, tally.shape)
+            std_dev = np.reshape(std_dev, tally.shape)
 
             # Override tally's data with the new condensed data
             tally._mean = mean
@@ -765,6 +798,7 @@ class MGXS(object):
 
         # Compute the energy condensed multi-group cross section
         condensed_xs.compute_xs()
+        condensed_xs.sparse = self.sparse
         return condensed_xs
 
     def get_subdomain_avg_xs(self, subdomains='all'):
@@ -807,8 +841,9 @@ class MGXS(object):
 
         # Clone this MGXS to initialize the subdomain-averaged version
         avg_xs = copy.deepcopy(self)
+        avg_xs.sparse = False
 
-            # If domain is distribcell, make the new domain 'cell'
+        # If domain is distribcell, make the new domain 'cell'
         if self.domain_type == 'distribcell':
             avg_xs.domain_type = 'cell'
 
@@ -836,10 +871,8 @@ class MGXS(object):
                 domain_filter.num_bins = 1
 
             # Reshape averaged data arrays with one dimension for all filters
-            new_shape = \
-                (tally.num_filter_bins, tally.num_nuclides, tally.num_scores,)
-            mean = np.reshape(mean, new_shape)
-            std_dev = np.reshape(std_dev, new_shape)
+            mean = np.reshape(mean, tally.shape)
+            std_dev = np.reshape(std_dev, tally.shape)
 
             # Override tally's data with the new condensed data
             tally._mean = mean
@@ -847,7 +880,7 @@ class MGXS(object):
 
         # Compute the subdomain-averaged multi-group cross section
         avg_xs.compute_xs()
-
+        avg_xs.sparse = self.sparse
         return avg_xs
 
     def print_xs(self, subdomains='all', nuclides='all', xs_type='macro'):
diff --git a/openmc/opencg_compatible.py b/openmc/opencg_compatible.py
index 61485b5099..0bda48c160 100644
--- a/openmc/opencg_compatible.py
+++ b/openmc/opencg_compatible.py
@@ -901,7 +901,9 @@ def get_opencg_lattice(openmc_lattice):
 
     # Convert 2D universes array to 3D for OpenCG
     if len(universes.shape) == 2:
-        universes.shape = (1,) + universes.shape
+        new_universes = universes.copy()
+        new_universes.shape = (1,) + universes.shape
+        universes = new_universes
 
     # Initialize an empty array for the OpenCG nested Universes in this Lattice
     universe_array = np.ndarray(tuple(np.array(dimension)[::-1]),
diff --git a/openmc/statepoint.py b/openmc/statepoint.py
index 2b1e57dfa4..3c0759f124 100644
--- a/openmc/statepoint.py
+++ b/openmc/statepoint.py
@@ -75,6 +75,9 @@ class StatePoint(object):
         energy of the source site.
     source_present : bool
         Indicate whether source sites are present
+    sparse : bool
+        Whether or not the tallies uses SciPy's LIL sparse matrix format for
+        compressed data storage
     tallies : dict
         Dictionary whose keys are tally IDs and whose values are Tally objects
     tallies_present : bool
@@ -110,6 +113,7 @@ class StatePoint(object):
         self._tallies_read = False
         self._summary = False
         self._global_tallies = None
+        self._sparse = False
 
     def close(self):
         self._f.close()
@@ -318,6 +322,10 @@ class StatePoint(object):
     def source_present(self):
         return self._f['source_present'].value > 0
 
+    @property
+    def sparse(self):
+        return self._sparse
+
     @property
     def tallies(self):
         if not self._tallies_read:
@@ -340,7 +348,8 @@ class StatePoint(object):
             for tally_key in tally_keys:
 
                 # Read the Tally size specifications
-                n_realizations = self._f['{0}{1}/n_realizations'.format(base, tally_key)].value
+                n_realizations = \
+                    self._f['{0}{1}/n_realizations'.format(base, tally_key)].value
 
                 # Create Tally object and assign basic properties
                 tally = openmc.Tally(tally_id=tally_key)
@@ -350,7 +359,8 @@ class StatePoint(object):
                 tally.num_realizations = n_realizations
 
                 # Read the number of Filters
-                n_filters = self._f['{0}{1}/n_filters'.format(base, tally_key)].value
+                n_filters = \
+                    self._f['{0}{1}/n_filters'.format(base, tally_key)].value
 
                 subbase = '{0}{1}/filter '.format(base, tally_key)
 
@@ -358,7 +368,8 @@ class StatePoint(object):
                 for j in range(1, n_filters+1):
 
                     # Read the Filter type
-                    filter_type = self._f['{0}{1}/type'.format(subbase, j)].value.decode()
+                    filter_type = \
+                        self._f['{0}{1}/type'.format(subbase, j)].value.decode()
 
                     n_bins = self._f['{0}{1}/n_bins'.format(subbase, j)].value
 
@@ -380,7 +391,8 @@ class StatePoint(object):
                     tally.add_filter(filter)
 
                 # Read Nuclide bins
-                nuclide_names = self._f['{0}{1}/nuclides'.format(base, tally_key)].value
+                nuclide_names = \
+                    self._f['{0}{1}/nuclides'.format(base, tally_key)].value
 
                 # Add all Nuclides to the Tally
                 for name in nuclide_names:
@@ -415,6 +427,7 @@ class StatePoint(object):
                     tally.add_score(score)
 
                 # Add Tally to the global dictionary of all Tallies
+                tally.sparse = self.sparse
                 self._tallies[tally_key] = tally
 
             self._tallies_read = True
@@ -439,6 +452,26 @@ class StatePoint(object):
     def with_summary(self):
         return False if self.summary is None else True
 
+    @sparse.setter
+    def sparse(self, sparse):
+        """Convert tally data from NumPy arrays to SciPy list of lists (LIL)
+        sparse matrices, and vice versa.
+
+        This property may be used to reduce the amount of data in memory during
+        tally data processing. The tally data will be stored as SciPy LIL
+        matrices internally within each Tally object. All tally data access
+        properties and methods will return data as a dense NumPy array.
+
+        """
+
+        cv.check_type('sparse', sparse, bool)
+        self._sparse = sparse
+
+        # Update tally sparsities
+        if self._tallies_read:
+            for tally_id in self.tallies:
+                self.tallies[tally_id].sparse = self.sparse
+
     def get_tally(self, scores=[], filters=[], nuclides=[],
                   name=None, id=None, estimator=None):
         """Finds and returns a Tally object with certain properties.
diff --git a/openmc/tallies.py b/openmc/tallies.py
index 190f147279..7ec1ca50cb 100644
--- a/openmc/tallies.py
+++ b/openmc/tallies.py
@@ -73,6 +73,9 @@ class Tally(object):
         Total number of filter bins accounting for all filters
     num_bins : Integral
         Total number of bins for the tally
+    shape : 3-tuple of Integral
+        The shape of the tally data array ordered as the number of filter bins, 
+        nuclide bins and score bins
     num_realizations : Integral
         Total number of realizations
     with_summary : bool
@@ -86,6 +89,11 @@ class Tally(object):
         An array containing the sample mean for each bin
     std_dev : ndarray
         An array containing the sample standard deviation for each bin
+    derived : bool
+        Whether or not the tally is derived from one or more other tallies
+    sparse : bool
+        Whether or not the tally uses SciPy's LIL sparse matrix format for
+        compressed data storage
 
     """
 
@@ -108,6 +116,7 @@ class Tally(object):
         self._std_dev = None
         self._with_batch_statistics = False
         self._derived = False
+        self._sparse = False
 
         self._sp_filename = None
         self._results_read = False
@@ -129,6 +138,7 @@ class Tally(object):
             clone._with_summary = self.with_summary
             clone._with_batch_statistics = self.with_batch_statistics
             clone._derived = self.derived
+            clone._sparse = self.sparse
             clone._sp_filename = self._sp_filename
             clone._results_read = self._results_read
 
@@ -265,6 +275,10 @@ class Tally(object):
         num_bins *= self.num_scores
         return num_bins
 
+    @property
+    def shape(self):
+        return (self.num_filter_bins, self.num_nuclides, self.num_scores)
+
     @property
     def estimator(self):
         return self._estimator
@@ -298,29 +312,33 @@ class Tally(object):
             sum = data['sum']
             sum_sq = data['sum_sq']
 
-            # Define a routine to convert 0 to 1
-            def nonzero(val):
-                return 1 if not val else val
-
             # Reshape the results arrays
-            new_shape = (nonzero(self.num_filter_bins),
-                         nonzero(self.num_nuclides),
-                         nonzero(self.num_scores))
-
-            sum = np.reshape(sum, new_shape)
-            sum_sq = np.reshape(sum_sq, new_shape)
+            sum = np.reshape(sum, self.shape)
+            sum_sq = np.reshape(sum_sq, self.shape)
 
             # Set the data for this Tally
             self._sum = sum
             self._sum_sq = sum_sq
 
+            # Convert NumPy arrays to SciPy sparse LIL matrices
+            if self.sparse:
+                import scipy.sparse as sps
+
+                self._sum = \
+                    sps.lil_matrix(self._sum.flatten(), self._sum.shape)
+                self._sum_sq = \
+                    sps.lil_matrix(self._sum_sq.flatten(), self._sum_sq.shape)
+
             # Indicate that Tally results have been read
             self._results_read = True
 
             # Close the HDF5 statepoint file
             f.close()
 
-        return self._sum
+        if self.sparse:
+            return np.reshape(self._sum.toarray(), self.shape)
+        else:
+            return self._sum
 
     @property
     def sum_sq(self):
@@ -331,7 +349,10 @@ class Tally(object):
             # Force reading of sum and sum_sq
             self.sum
 
-        return self._sum_sq
+        if self.sparse:
+            return np.reshape(self._sum_sq.toarray(), self.shape)
+        else:
+            return self._sum_sq
 
     @property
     def mean(self):
@@ -340,7 +361,18 @@ class Tally(object):
                 return None
 
             self._mean = self.sum / self.num_realizations
-        return self._mean
+
+            # Convert NumPy array to SciPy sparse LIL matrix
+            if self.sparse:
+                import scipy.sparse as sps
+
+                self._mean = \
+                    sps.lil_matrix(self._mean.flatten(), self._mean.shape)
+
+        if self.sparse:
+            return np.reshape(self._mean.toarray(), self.shape)
+        else:
+            return self._mean
 
     @property
     def std_dev(self):
@@ -353,8 +385,20 @@ class Tally(object):
             self._std_dev = np.zeros_like(self.mean)
             self._std_dev[nonzero] = np.sqrt((self.sum_sq[nonzero]/n -
                                               self.mean[nonzero]**2)/(n - 1))
+
+            # Convert NumPy array to SciPy sparse LIL matrix
+            if self.sparse:
+                import scipy.sparse as sps
+
+                self._std_dev = \
+                    sps.lil_matrix(self._std_dev.flatten(), self._std_dev.shape)
+
             self.with_batch_statistics = True
-        return self._std_dev
+
+        if self.sparse:
+            return np.reshape(self._std_dev.toarray(), self.shape)
+        else:
+            return self._std_dev
 
     @property
     def with_batch_statistics(self):
@@ -364,6 +408,10 @@ class Tally(object):
     def derived(self):
         return self._derived
 
+    @property
+    def sparse(self):
+        return self._sparse
+
     @estimator.setter
     def estimator(self, estimator):
         cv.check_value('estimator', estimator,
@@ -491,6 +539,51 @@ class Tally(object):
         cv.check_type('sum_sq', sum_sq, Iterable)
         self._sum_sq = sum_sq
 
+    @sparse.setter
+    def sparse(self, sparse):
+        """Convert tally data from NumPy arrays to SciPy list of lists (LIL)
+        sparse matrices, and vice versa.
+
+        This property may be used to reduce the amount of data in memory during
+        tally data processing. The tally data will be stored as SciPy LIL
+        matrices internally within the Tally object. All tally data access
+        properties and methods will return data as a dense NumPy array.
+
+        """
+
+        cv.check_type('sparse', sparse, bool)
+
+        # Convert NumPy arrays to SciPy sparse LIL matrices
+        if sparse and not self.sparse:
+            import scipy.sparse as sps
+
+            if self._sum is not None:
+                self._sum = \
+                    sps.lil_matrix(self._sum.flatten(), self._sum.shape)
+            if self._sum_sq is not None:
+                self._sum_sq = \
+                    sps.lil_matrix(self._sum_sq.flatten(), self._sum_sq.shape)
+            if self._mean is not None:
+                self._mean = \
+                    sps.lil_matrix(self._mean.flatten(), self._mean.shape)
+            if self._std_dev is not None:
+                self._std_dev = \
+                    sps.lil_matrix(self._std_dev.flatten(), self._std_dev.shape)
+
+            self._sparse = True
+
+        # Convert SciPy sparse LIL matrices to NumPy arrays
+        elif not sparse and self.sparse:
+            if self._sum is not None:
+                self._sum = np.reshape(self._sum.toarray(), self.shape)
+            if self._sum_sq is not None:
+                self._sum_sq = np.reshape(self._sum_sq.toarray(), self.shape)
+            if self._mean is not None:
+                self._mean = np.reshape(self._mean.toarray(), self.shape)
+            if self._std_dev is not None:
+                self._std_dev = np.reshape(self._std_dev.toarray(), self.shape)
+            self._sparse = False
+
     def remove_score(self, score):
         """Remove a score from the tally
 
@@ -618,7 +711,8 @@ class Tally(object):
         """
 
         if not self.can_merge(tally):
-            msg = 'Unable to merge tally ID="{0}" with "{1}"'.format(tally.id, self.id)
+            msg = 'Unable to merge tally ID="{0}" with ' + \
+                   '"{1}"'.format(tally.id, self.id)
             raise ValueError(msg)
 
         # Create deep copy of tally to return as merged tally
@@ -1067,22 +1161,19 @@ class Tally(object):
         Raises
         ------
         ValueError
-            When this method is called before the Tally is populated with data
-            by the StatePoint.read_results() method. ValueError is also thrown
-            if the input parameters do not correspond to the Tally's attributes,
+            When this method is called before the Tally is populated with data,
+            or the input parameters do not correspond to the Tally's attributes,
             e.g., if the score(s) do not match those in the Tally.
 
         """
 
-        # Ensure that StatePoint.read_results() was called first
+        # Ensure that the tally has data
         if (value == 'mean' and self.mean is None) or \
            (value == 'std_dev' and self.std_dev is None) or \
            (value == 'rel_err' and self.mean is None) or \
            (value == 'sum' and self.sum is None) or \
            (value == 'sum_sq' and self.sum_sq is None):
-            msg = 'The Tally ID="{0}" has no data to return. Call the ' \
-                  'StatePoint.read_results() method before using ' \
-                  'Tally.get_values(...)'.format(self.id)
+            msg = 'The Tally ID="{0}" has no data to return'.format(self.id)
             raise ValueError(msg)
 
         # Get filter, nuclide and score indices
@@ -1149,17 +1240,14 @@ class Tally(object):
         ------
         KeyError
             When this method is called before the Tally is populated with data
-            by the StatePoint.read_results() method.
         ImportError
             When Pandas can not be found on the caller's system
 
         """
 
-        # Ensure that StatePoint.read_results() was called first
+        # Ensure that the tally has data
         if self.mean is None or self.std_dev is None:
-            msg = 'The Tally ID="{0}" has no data to return. Call the ' \
-                  'StatePoint.read_results() method before using ' \
-                  'Tally.get_pandas_dataframe(...)'.format(self.id)
+            msg = 'The Tally ID="{0}" has no data to return'.format(self.id)
             raise KeyError(msg)
 
         # If using Summary, ensure StatePoint.link_with_summary(...) was called
@@ -1305,16 +1393,13 @@ class Tally(object):
         Raises
         ------
         KeyError
-            When this method is called before the Tally is populated with data
-            by the StatePoint.read_results() method.
+            When this method is called before the Tally is populated with data.
 
         """
 
-        # Ensure that StatePoint.read_results() was called first
+        # Ensure that the tally has data
         if self._sum is None or self._sum_sq is None and not self.derived:
-            msg = 'The Tally ID="{0}" has no data to export. Call the ' \
-                  'StatePoint.read_results() method before using ' \
-                  'Tally.export_results(...)'.format(self.id)
+            msg = 'The Tally ID="{0}" has no data to export'.format(self.id)
             raise KeyError(msg)
 
         if not isinstance(filename, basestring):
@@ -1476,7 +1561,7 @@ class Tally(object):
         ------
         ValueError
             When this method is called before the other tally is populated
-            with data by the StatePoint.read_results() method.
+            with data.
 
         """
 
@@ -1531,6 +1616,9 @@ class Tally(object):
         self_copy = copy.deepcopy(self)
         other_copy = copy.deepcopy(other)
 
+        self_copy.sparse = False
+        other_copy.sparse = False
+
         # Align the tally data based on desired hybrid product
         data = self_copy._align_tally_data(other_copy, filter_product,
                                            nuclide_product, score_product)
@@ -1596,7 +1684,8 @@ class Tally(object):
         else:
             all_nuclides = [self_copy.nuclides, other_copy.nuclides]
             for self_nuclide, other_nuclide in itertools.product(*all_nuclides):
-                new_nuclide = CrossNuclide(self_nuclide, other_nuclide, binary_op)
+                new_nuclide = \
+                    CrossNuclide(self_nuclide, other_nuclide, binary_op)
                 new_tally.add_nuclide(new_nuclide)
 
         # Add scores to the new tally
@@ -1653,8 +1742,10 @@ class Tally(object):
         """
 
         # Get the set of filters that each tally is missing
-        other_missing_filters = set(self.filters).difference(set(other.filters))
-        self_missing_filters = set(other.filters).difference(set(self.filters))
+        other_missing_filters = \
+            set(self.filters).difference(set(other.filters))
+        self_missing_filters = \
+            set(other.filters).difference(set(self.filters))
 
         # Add filters present in self but not in other to other
         for filter in other_missing_filters:
@@ -1681,30 +1772,40 @@ class Tally(object):
         # Repeat and tile the data by nuclide in preparation for performing
         # the tensor product across nuclides.
         if nuclide_product == 'tensor':
-            self._mean = np.repeat(self.mean, other.num_nuclides, axis=1)
-            self._std_dev = np.repeat(self.std_dev, other.num_nuclides, axis=1)
-            other._mean = np.tile(other.mean, (1, self.num_nuclides, 1))
-            other._std_dev = np.tile(other.std_dev, (1, self.num_nuclides, 1))
+            self._mean = \
+                np.repeat(self.mean, other.num_nuclides, axis=1)
+            self._std_dev = \
+                np.repeat(self.std_dev, other.num_nuclides, axis=1)
+            other._mean = \
+                np.tile(other.mean, (1, self.num_nuclides, 1))
+            other._std_dev = \
+                np.tile(other.std_dev, (1, self.num_nuclides, 1))
 
         # Add nuclides to each tally such that each tally contains the complete
-        # set of nuclides necessary to perform an entrywise product. New nuclides
-        # added to a tally will have all their scores set to zero.
+        # set of nuclides necessary to perform an entrywise product. New 
+        # nuclides added to a tally will have all their scores set to zero.
         else:
 
             # Get the set of nuclides that each tally is missing
-            other_missing_nuclides = set(self.nuclides).difference(set(other.nuclides))
-            self_missing_nuclides = set(other.nuclides).difference(set(self.nuclides))
+            other_missing_nuclides = \
+                set(self.nuclides).difference(set(other.nuclides))
+            self_missing_nuclides = \
+                set(other.nuclides).difference(set(self.nuclides))
 
             # Add nuclides present in self but not in other to other
             for nuclide in other_missing_nuclides:
-                other._mean = np.insert(other.mean, other.num_nuclides, 0, axis=1)
-                other._std_dev = np.insert(other.std_dev, other.num_nuclides, 0, axis=1)
+                other._mean = \
+                    np.insert(other.mean, other.num_nuclides, 0, axis=1)
+                other._std_dev = \
+                    np.insert(other.std_dev, other.num_nuclides, 0, axis=1)
                 other.add_nuclide(nuclide)
 
             # Add nuclides present in other but not in self to self
             for nuclide in self_missing_nuclides:
-                self._mean = np.insert(self.mean, self.num_nuclides, 0, axis=1)
-                self._std_dev = np.insert(self.std_dev, self.num_nuclides, 0, axis=1)
+                self._mean = \
+                    np.insert(self.mean, self.num_nuclides, 0, axis=1)
+                self._std_dev = \
+                    np.insert(self.std_dev, self.num_nuclides, 0, axis=1)
                 self.add_nuclide(nuclide)
 
             # Align other nuclides with self nuclides
@@ -1729,8 +1830,10 @@ class Tally(object):
         else:
 
             # Get the set of scores that each tally is missing
-            other_missing_scores = set(self.scores).difference(set(other.scores))
-            self_missing_scores = set(other.scores).difference(set(self.scores))
+            other_missing_scores = \
+                set(self.scores).difference(set(other.scores))
+            self_missing_scores = \
+                set(other.scores).difference(set(self.scores))
 
             # Add scores present in self but not in other to other
             for score in other_missing_scores:
@@ -1792,7 +1895,7 @@ class Tally(object):
         ------
         ValueError
             If this is a derived tally or this method is called before the tally
-            is populated with data by the StatePoint.read_results() method.
+            is populated with data.
 
         """
 
@@ -1879,7 +1982,7 @@ class Tally(object):
         ------
         ValueError
             If this is a derived tally or this method is called before the tally
-            is populated with data by the StatePoint.read_results() method.
+            is populated with data.
 
         """
 
@@ -1946,7 +2049,7 @@ class Tally(object):
         ------
         ValueError
             If this is a derived tally or this method is called before the tally
-            is populated with data by the StatePoint.read_results() method.
+            is populated with data.
 
         """
 
@@ -2030,8 +2133,7 @@ class Tally(object):
         Raises
         ------
         ValueError
-            When this method is called before the Tally is populated with data
-            by the StatePoint.read_results() method.
+            When this method is called before the Tally is populated with data.
 
         """
 
@@ -2044,6 +2146,10 @@ class Tally(object):
         if isinstance(other, Tally):
             new_tally = self.hybrid_product(other, binary_op='+')
 
+            # If both tally operands were sparse, sparsify the new tally
+            if self.sparse and other.sparse:
+                new_tally.sparse = True
+
         elif isinstance(other, Real):
             new_tally = Tally(name='derived')
             new_tally._derived = True
@@ -2062,6 +2168,9 @@ class Tally(object):
             for score in self.scores:
                 new_tally.add_score(score)
 
+            # If this tally operand is sparse, sparsify the new tally
+            new_tally.sparse = self.sparse
+
         else:
             msg = 'Unable to add "{0}" to Tally ID="{1}"'.format(other, self.id)
             raise ValueError(msg)
@@ -2099,8 +2208,7 @@ class Tally(object):
         Raises
         ------
         ValueError
-            When this method is called before the Tally is populated with data
-            by the StatePoint.read_results() method.
+            When this method is called before the Tally is populated with data.
 
         """
 
@@ -2113,6 +2221,10 @@ class Tally(object):
         if isinstance(other, Tally):
             new_tally = self.hybrid_product(other, binary_op='-')
 
+            # If both tally operands were sparse, sparsify the new tally
+            if self.sparse and other.sparse:
+                new_tally.sparse = True
+
         elif isinstance(other, Real):
             new_tally = Tally(name='derived')
             new_tally._derived = True
@@ -2130,6 +2242,9 @@ class Tally(object):
             for score in self.scores:
                 new_tally.add_score(score)
 
+            # If this tally operand is sparse, sparsify the new tally
+            new_tally.sparse = self.sparse
+
         else:
             msg = 'Unable to subtract "{0}" from Tally ' \
                   'ID="{1}"'.format(other, self.id)
@@ -2168,8 +2283,7 @@ class Tally(object):
         Raises
         ------
         ValueError
-            When this method is called before the Tally is populated with data
-            by the StatePoint.read_results() method.
+            When this method is called before the Tally is populated with data.
 
         """
 
@@ -2182,6 +2296,10 @@ class Tally(object):
         if isinstance(other, Tally):
             new_tally = self.hybrid_product(other, binary_op='*')
 
+            # If original tally operands were sparse, sparsify the new tally
+            if self.sparse and other.sparse:
+                new_tally.sparse = True
+
         elif isinstance(other, Real):
             new_tally = Tally(name='derived')
             new_tally._derived = True
@@ -2199,6 +2317,9 @@ class Tally(object):
             for score in self.scores:
                 new_tally.add_score(score)
 
+            # If this tally operand is sparse, sparsify the new tally
+            new_tally.sparse = self.sparse
+
         else:
             msg = 'Unable to multiply Tally ID="{0}" ' \
                   'by "{1}"'.format(self.id, other)
@@ -2237,8 +2358,7 @@ class Tally(object):
         Raises
         ------
         ValueError
-            When this method is called before the Tally is populated with data
-            by the StatePoint.read_results() method.
+            When this method is called before the Tally is populated with data.
 
         """
 
@@ -2251,6 +2371,10 @@ class Tally(object):
         if isinstance(other, Tally):
             new_tally = self.hybrid_product(other, binary_op='/')
 
+            # If original tally operands were sparse, sparsify the new tally
+            if self.sparse and other.sparse:
+                new_tally.sparse = True
+
         elif isinstance(other, Real):
             new_tally = Tally(name='derived')
             new_tally._derived = True
@@ -2268,6 +2392,9 @@ class Tally(object):
             for score in self.scores:
                 new_tally.add_score(score)
 
+            # If this tally operand is sparse, sparsify the new tally
+            new_tally.sparse = self.sparse
+
         else:
             msg = 'Unable to divide Tally ID="{0}" ' \
                   'by "{1}"'.format(self.id, other)
@@ -2309,8 +2436,7 @@ class Tally(object):
         Raises
         ------
         ValueError
-            When this method is called before the Tally is populated with data
-            by the StatePoint.read_results() method.
+            When this method is called before the Tally is populated with data.
 
         """
 
@@ -2323,6 +2449,10 @@ class Tally(object):
         if isinstance(power, Tally):
             new_tally = self.hybrid_product(power, binary_op='^')
 
+            # If original tally operands were sparse, sparsify the new tally
+            if self.sparse and other.sparse:
+                new_tally.sparse = True
+
         elif isinstance(power, Real):
             new_tally = Tally(name='derived')
             new_tally._derived = True
@@ -2341,6 +2471,9 @@ class Tally(object):
             for score in self.scores:
                 new_tally.add_score(score)
 
+            # If original tally was sparse, sparsify the exponentiated tally
+            new_tally.sparse = self.sparse
+
         else:
             msg = 'Unable to raise Tally ID="{0}" to ' \
                   'power "{1}"'.format(self.id, power)
@@ -2491,18 +2624,18 @@ class Tally(object):
         Raises
         ------
         ValueError
-            When this method is called before the Tally is populated with data
-            by the StatePoint.read_results() method.
+            When this method is called before the Tally is populated with data.
 
         """
 
-        # Ensure that StatePoint.read_results() was called first
+        # Ensure that the tally has data
         if not self.derived and self.sum is None:
             msg = 'Unable to use tally arithmetic with Tally ID="{0}" ' \
                   'since it does not contain any results.'.format(self.id)
             raise ValueError(msg)
 
         new_tally = copy.deepcopy(self)
+        new_tally.sparse = False
 
         if self.sum is not None:
             new_sum = self.get_values(scores, filters, filter_bins,
@@ -2518,7 +2651,7 @@ class Tally(object):
             new_tally._mean = new_mean
         if self.std_dev is not None:
             new_std_dev = self.get_values(scores, filters, filter_bins,
-                                       nuclides, 'std_dev')
+                                          nuclides, 'std_dev')
             new_tally._std_dev = new_std_dev
 
         # SCORES
@@ -2578,6 +2711,8 @@ class Tally(object):
             filter.stride = stride
             stride *= filter.num_bins
 
+        # If original tally was sparse, sparsify the sliced tally
+        new_tally.sparse = self.sparse
         return new_tally
 
     def summation(self, scores=[], filter_type=None,
@@ -2680,6 +2815,8 @@ class Tally(object):
                     filters[i] = CrossFilter(filters[i-1], filters[i], '+')
                 tally_sum.add_filter(filters[-1])
 
+        # If original tally was sparse, sparsify the tally summation
+        tally_sum.sparse = self.sparse
         return tally_sum
 
     def diagonalize_filter(self, new_filter):
@@ -2716,12 +2853,6 @@ class Tally(object):
         new_tally = copy.deepcopy(self)
         new_tally.add_filter(new_filter)
 
-        # Determine the shape of data in the new diagonalized Tally
-        num_filter_bins = new_tally.num_filter_bins
-        num_nuclides = new_tally.num_nuclides
-        num_scores = new_tally.num_scores
-        new_shape = (num_filter_bins, num_nuclides, num_scores)
-
         # Determine "base" indices along the new "diagonal", and the factor
         # by which the "base" indices should be repeated to account for all
         # other filter bins in the diagonalized tally
@@ -2737,16 +2868,16 @@ class Tally(object):
 
         # Inject this Tally's data along the diagonal of the diagonalized Tally
         if self.sum is not None:
-            new_tally._sum = np.zeros(new_shape, dtype=np.float64)
+            new_tally._sum = np.zeros(new_tally.shape, dtype=np.float64)
             new_tally._sum[diag_indices, :, :] = self.sum
         if self.sum_sq is not None:
-            new_tally._sum_sq = np.zeros(new_shape, dtype=np.float64)
+            new_tally._sum_sq = np.zeros(new_tally.shape, dtype=np.float64)
             new_tally._sum_sq[diag_indices, :, :] = self.sum_sq
         if self.mean is not None:
-            new_tally._mean = np.zeros(new_shape, dtype=np.float64)
+            new_tally._mean = np.zeros(new_tally.shape, dtype=np.float64)
             new_tally._mean[diag_indices, :, :] = self.mean
         if self.std_dev is not None:
-            new_tally._std_dev = np.zeros(new_shape, dtype=np.float64)
+            new_tally._std_dev = np.zeros(new_tally.shape, dtype=np.float64)
             new_tally._std_dev[diag_indices, :, :] = self.std_dev
 
         # Correct each Filter's stride
@@ -2755,6 +2886,8 @@ class Tally(object):
             filter.stride = stride
             stride *= filter.num_bins
 
+        # If original tally was sparse, sparsify the diagonalized tally
+        new_tally.sparse = self.sparse
         return new_tally
 
 
diff --git a/openmc/universe.py b/openmc/universe.py
index b997e8845c..8417f27c71 100644
--- a/openmc/universe.py
+++ b/openmc/universe.py
@@ -1081,7 +1081,7 @@ class RectLattice(Lattice):
         # For 2D Lattices
         if len(self._dimension) == 2:
             offset = self._offsets[i[1]-1, i[2]-1, 0, distribcell_index-1]
-            offset += self._universes[i[1]][i[2]].get_cell_instance(path,
+            offset += self._universes[i[1]-1][i[2]-1].get_cell_instance(path,
                                                               distribcell_index)
 
         # For 3D Lattices
diff --git a/setup.py b/setup.py
index 907d80e31b..d401f5fbaf 100644
--- a/setup.py
+++ b/setup.py
@@ -37,6 +37,7 @@ if have_setuptools:
         # Optional dependencies
         'extras_require': {
             'pandas': ['pandas'],
+            'sparse' : ['scipy'],
             'vtk': ['vtk', 'silomesh'],
             'validate': ['lxml']
         }})