diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py
index db31e2f32..b4a14f8ce 100644
--- a/openmc/mgxs/mgxs.py
+++ b/openmc/mgxs/mgxs.py
@@ -62,6 +62,9 @@ _DOMAINS = (openmc.Cell,
 # Supported ScatterMatrixXS and NuScatterMatrixXS angular distribution types
 MU_TREATMENTS = ('legendre', 'histogram')
 
+# Maximum Legendre order supported by OpenMC
+MAX_LEGENDRE = 10
+
 
 @add_metaclass(ABCMeta)
 class MGXS(object):
@@ -3465,7 +3468,8 @@ class ScatterMatrixXS(MatrixMGXS):
         cv.check_type('legendre_order', legendre_order, Integral)
         cv.check_greater_than('legendre_order', legendre_order, 0,
                               equality=True)
-        cv.check_less_than('legendre_order', legendre_order, 10, equality=True)
+        cv.check_less_than('legendre_order', legendre_order, MAX_LEGENDRE,
+                           equality=True)
 
         if self.scatter_format == 'legendre':
             if self.correction == 'P0' and legendre_order > 0:
diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py
index ef1dd236c..6fdbb1de7 100644
--- a/openmc/mgxs_library.py
+++ b/openmc/mgxs_library.py
@@ -2132,3 +2132,43 @@ class MGXSLibrary(object):
             xsdata.to_hdf5(file)
 
         file.close()
+
+    @classmethod
+    def from_hdf5(cls, filename=None):
+        """Generate an MGXS Library from an HDF5 group or file
+        Parameters
+        ----------
+        filename : str, optional
+            Name of HDF5 file containing MGXS data. Default is None.
+            If not provided, the value of the OPENMC_MG_CROSS_SECTIONS
+            environmental variable will be used
+        Returns
+        -------
+        openmc.MGXSLibrary
+            Multi-group cross section data object.
+        """
+
+        # If filename is None, get the cross sections from the
+        # OPENMC_CROSS_SECTIONS environment variable
+        if filename is None:
+            filename = os.environ.get('OPENMC_MG_CROSS_SECTIONS')
+
+        # Check to make sure there was an environmental variable.
+        if filename is None:
+            raise ValueError("Either path or OPENMC_MG_CROSS_SECTIONS "
+                             "environmental variable must be set")
+
+        check_type('filename', filename, str)
+        file = h5py.File(filename, 'r')
+
+        group_structure = file.attrs['group structure']
+        num_delayed_groups = file.attrs['delayed_groups']
+        energy_groups = openmc.mgxs.EnergyGroups(group_structure)
+        data = cls(energy_groups, num_delayed_groups)
+
+        for group_name, group in file.items():
+            data.add_xsdata(openmc.XSdata.from_hdf5(group, group_name,
+                                                    energy_groups,
+                                                    num_delayed_groups))
+
+        return data
diff --git a/openmc/plotter.py b/openmc/plotter.py
index aa48ef944..334567d2a 100644
--- a/openmc/plotter.py
+++ b/openmc/plotter.py
@@ -10,10 +10,14 @@ import openmc.data
 PLOT_TYPES = ['total', 'scatter', 'elastic', 'inelastic', 'fission',
               'absorption', 'capture', 'nu-fission', 'nu-scatter', 'unity',
               'slowing-down power', 'damage']
+
 # Supported keywoards for multi-group cross section plotting
-PLOT_TYPES_MGXS = ['total', 'absorption', 'fission', 'kappa-fission',
-                   'chi', 'chi-prompt', 'nu-fission', 'prompt-nu-fission',
-                   'inverse-velocity', 'unity']
+PLOT_TYPES_MGXS = ['total', 'absorption', 'scatter', 'fission',
+                   'kappa-fission', 'chi', 'chi-prompt', 'nu-fission',
+                   'prompt-nu-fission', 'inverse-velocity', 'unity']
+# Add on values for scattering moments
+PLOT_TYPES_MGXS += ['scatter-' + str(i)
+                    for i in range(0, openmc.mgxs.MAX_LEGENDRE + 1)]
 
 # Special MT values
 UNITY_MT = -1
@@ -822,6 +826,21 @@ def _calculate_mgxs_nuc_macro(this, types, library, temperature=294.):
         for i, line in enumerate(types):
             if line == 'unity':
                 data[i, :] = 1.
+            elif line.startswith('scatter'):
+                # We have to remove the outgoing dependence
+                attr = line.replace(' ', '_').replace('-', '_')
+                matrix = xsdata.scatter_matrix[t]
+                # Sum over outgoing groups
+                vector = np.sum(matrix, axis=1)
+                # Now get the actual order of interest
+                if line == 'scatter':
+                    order = 0
+                else:
+                    order = int(line.split('-')[1])
+                if order < xsdata.xs_shapes["[G][G'][Order]"][-1]:
+                    data[i, :] = vector[:, order]
+                else:
+                    data[i, :] = 0.
             else:
                 attr = line.replace(' ', '_').replace('-', '_')
                 data[i, :] = getattr(xsdata, attr)[t]
@@ -883,10 +902,10 @@ def _calculate_mgxs_elem_mat(this, types, library, temperature=294.,
             nuclides = {this._macroscopic: (this._macroscopic, this.density)}
         else:
             # Expand elements in to nuclides with atomic densities
-            nuclides = this.get_nuclide_atom_densities(ce_cross_sections)
+            nuclides = this.get_nuclide_atom_densities()
 
         # For ease of processing split out nuc and nuc_density
-        nuc_multiplier = [nuclide[1][1] for nuclide in nuclides.items()]
+        nuc_fraction = [nuclide[1][1] for nuclide in nuclides.items()]
     else:
         T = temperature
         # Expand elements in to nuclides with atomic densities
@@ -894,7 +913,7 @@ def _calculate_mgxs_elem_mat(this, types, library, temperature=294.,
                                cross_sections=ce_cross_sections)
 
         # For ease of processing split out nuc and nuc_fractions
-        nuc_multiplier = [nuclide[1] for nuclide in nuclides]
+        nuc_fraction = [nuclide[1] for nuclide in nuclides]
 
     nuc_data = []
     for nuclide in nuclides.items():
@@ -908,6 +927,6 @@ def _calculate_mgxs_elem_mat(this, types, library, temperature=294.,
             data[line, :] = 1.
         else:
             for n in range(len(nuclides)):
-                data[line, :] += nuc_multiplier[n] * nuc_data[n][line, :]
+                data[line, :] += nuc_fraction[n] * nuc_data[n][line, :]
 
     return data