diff --git a/openmc/geometry.py b/openmc/geometry.py
index ad86cc56b..ee4a58a3b 100644
--- a/openmc/geometry.py
+++ b/openmc/geometry.py
@@ -1,11 +1,14 @@
+from __future__ import annotations
 import os
 import typing
 from collections import OrderedDict, defaultdict
-from collections.abc import Iterable
 from copy import deepcopy
+from collections.abc import Iterable
 from pathlib import Path
-from xml.etree import ElementTree as ET
 import warnings
+from xml.etree import ElementTree as ET
+
+import numpy as np
 
 import openmc
 import openmc._xml as xml
@@ -51,19 +54,19 @@ class Geometry:
                 self._root_universe = univ
 
     @property
-    def root_universe(self):
+    def root_universe(self) -> openmc.UniverseBase:
         return self._root_universe
 
     @property
-    def bounding_box(self):
+    def bounding_box(self) -> np.ndarray:
         return self.root_universe.bounding_box
 
     @property
-    def merge_surfaces(self):
+    def merge_surfaces(self) -> bool:
         return self._merge_surfaces
 
     @property
-    def surface_precision(self):
+    def surface_precision(self) -> int:
         return self._surface_precision
 
     @root_universe.setter
@@ -105,7 +108,7 @@ class Geometry:
                 if universe.id in volume_calc.volumes:
                     universe.add_volume_information(volume_calc)
 
-    def to_xml_element(self, remove_surfs=False):
+    def to_xml_element(self, remove_surfs=False) -> ET.Element:
         """Creates a 'geometry' element to be written to an XML file.
 
         Parameters
@@ -164,7 +167,7 @@ class Geometry:
         tree.write(str(p), xml_declaration=True, encoding='utf-8')
 
     @classmethod
-    def from_xml_element(cls, elem, materials=None):
+    def from_xml_element(cls, elem, materials=None) -> Geometry:
         """Generate geometry from an XML element
 
         Parameters
@@ -260,7 +263,7 @@ class Geometry:
         cls,
         path: PathLike = 'geometry.xml',
         materials: typing.Optional[typing.Union[PathLike, 'openmc.Materials']] = 'materials.xml'
-    ):
+    ) -> Geometry:
         """Generate geometry from XML file
 
         Parameters
@@ -290,7 +293,7 @@ class Geometry:
 
         return cls.from_xml_element(root, materials)
 
-    def find(self, point):
+    def find(self, point) -> list:
         """Find cells/universes/lattices which contain a given point
 
         Parameters
@@ -307,7 +310,7 @@ class Geometry:
         """
         return self.root_universe.find(point)
 
-    def get_instances(self, paths):
+    def get_instances(self, paths) -> typing.Union[int, typing.List[int]]:
         """Return the instance number(s) for a cell/material in a geometry path.
 
         The instance numbers are used as indices into distributed
@@ -354,7 +357,7 @@ class Geometry:
 
         return indices if return_list else indices[0]
 
-    def get_all_cells(self):
+    def get_all_cells(self) -> typing.OrderedDict[int, openmc.Cell]:
         """Return all cells in the geometry.
 
         Returns
@@ -368,7 +371,7 @@ class Geometry:
         else:
             return OrderedDict()
 
-    def get_all_universes(self):
+    def get_all_universes(self) -> typing.OrderedDict[int, openmc.Universe]:
         """Return all universes in the geometry.
 
         Returns
@@ -383,7 +386,7 @@ class Geometry:
         universes.update(self.root_universe.get_all_universes())
         return universes
 
-    def get_all_materials(self):
+    def get_all_materials(self) -> typing.OrderedDict[int, openmc.Material]:
         """Return all materials within the geometry.
 
         Returns
@@ -398,7 +401,7 @@ class Geometry:
         else:
             return OrderedDict()
 
-    def get_all_material_cells(self):
+    def get_all_material_cells(self) -> typing.OrderedDict[int, openmc.Cell]:
         """Return all cells filled by a material
 
         Returns
@@ -417,7 +420,7 @@ class Geometry:
 
         return material_cells
 
-    def get_all_material_universes(self):
+    def get_all_material_universes(self) -> typing.OrderedDict[int, openmc.Universe]:
         """Return all universes having at least one material-filled cell.
 
         This method can be used to find universes that have at least one cell
@@ -440,7 +443,7 @@ class Geometry:
 
         return material_universes
 
-    def get_all_lattices(self):
+    def get_all_lattices(self) -> typing.OrderedDict[int, openmc.Lattice]:
         """Return all lattices defined
 
         Returns
@@ -458,7 +461,7 @@ class Geometry:
 
         return lattices
 
-    def get_all_surfaces(self):
+    def get_all_surfaces(self) -> typing.OrderedDict[int, openmc.Surface]:
         """
         Return all surfaces used in the geometry
 
@@ -475,7 +478,7 @@ class Geometry:
                 surfaces = cell.region.get_surfaces(surfaces)
         return surfaces
 
-    def _get_domains_by_name(self, name, case_sensitive, matching, domain_type):
+    def _get_domains_by_name(self, name, case_sensitive, matching, domain_type) -> list:
         if not case_sensitive:
             name = name.lower()
 
@@ -492,7 +495,9 @@ class Geometry:
         domains.sort(key=lambda x: x.id)
         return domains
 
-    def get_materials_by_name(self, name, case_sensitive=False, matching=False):
+    def get_materials_by_name(
+        self, name, case_sensitive=False, matching=False
+    ) -> typing.List[openmc.Material]:
         """Return a list of materials with matching names.
 
         Parameters
@@ -513,7 +518,9 @@ class Geometry:
         """
         return self._get_domains_by_name(name, case_sensitive, matching, 'material')
 
-    def get_cells_by_name(self, name, case_sensitive=False, matching=False):
+    def get_cells_by_name(
+        self, name, case_sensitive=False, matching=False
+    ) -> typing.List[openmc.Cell]:
         """Return a list of cells with matching names.
 
         Parameters
@@ -534,7 +541,9 @@ class Geometry:
         """
         return self._get_domains_by_name(name, case_sensitive, matching, 'cell')
 
-    def get_surfaces_by_name(self, name, case_sensitive=False, matching=False):
+    def get_surfaces_by_name(
+        self, name, case_sensitive=False, matching=False
+    ) -> typing.List[openmc.Surface]:
         """Return a list of surfaces with matching names.
 
         .. versionadded:: 0.13.3
@@ -557,7 +566,9 @@ class Geometry:
         """
         return self._get_domains_by_name(name, case_sensitive, matching, 'surface')
 
-    def get_cells_by_fill_name(self, name, case_sensitive=False, matching=False):
+    def get_cells_by_fill_name(
+        self, name, case_sensitive=False, matching=False
+    ) -> typing.List[openmc.Cell]:
         """Return a list of cells with fills with matching names.
 
         Parameters
@@ -602,7 +613,9 @@ class Geometry:
 
         return sorted(cells, key=lambda x: x.id)
 
-    def get_universes_by_name(self, name, case_sensitive=False, matching=False):
+    def get_universes_by_name(
+        self, name, case_sensitive=False, matching=False
+    ) -> typing.List[openmc.Universe]:
         """Return a list of universes with matching names.
 
         Parameters
@@ -623,7 +636,9 @@ class Geometry:
         """
         return self._get_domains_by_name(name, case_sensitive, matching, 'universe')
 
-    def get_lattices_by_name(self, name, case_sensitive=False, matching=False):
+    def get_lattices_by_name(
+        self, name, case_sensitive=False, matching=False
+    ) -> typing.List[openmc.Lattice]:
         """Return a list of lattices with matching names.
 
         Parameters
@@ -644,7 +659,7 @@ class Geometry:
         """
         return self._get_domains_by_name(name, case_sensitive, matching, 'lattice')
 
-    def remove_redundant_surfaces(self):
+    def remove_redundant_surfaces(self) -> typing.Dict[int, openmc.Surface]:
         """Remove and return all of the redundant surfaces.
 
         Uses surface_precision attribute of Geometry instance for rounding and
@@ -707,7 +722,7 @@ class Geometry:
         # Recursively traverse the CSG tree to count all cell instances
         self.root_universe._determine_paths(instances_only=instances_only)
 
-    def clone(self):
+    def clone(self) -> Geometry:
         """Create a copy of this geometry with new unique IDs for all of its
         enclosed materials, surfaces, cells, universes and lattices."""
 
diff --git a/openmc/material.py b/openmc/material.py
index a2c6138d2..c84ff8c27 100644
--- a/openmc/material.py
+++ b/openmc/material.py
@@ -1,3 +1,4 @@
+from __future__ import annotations
 from collections import OrderedDict, defaultdict, namedtuple, Counter
 from collections.abc import Iterable
 from copy import deepcopy
@@ -6,11 +7,11 @@ from pathlib import Path
 import re
 import typing  # imported separately as py3.8 requires typing.Iterable
 import warnings
-from typing import Optional
+from typing import Optional, List, Union, Dict
 from xml.etree import ElementTree as ET
 
-import h5py
 import numpy as np
+import h5py
 
 import openmc
 import openmc.data
@@ -134,7 +135,7 @@ class Material(IDManagerMixin):
         # If specified, a list of table names
         self._sab = []
 
-    def __repr__(self):
+    def __repr__(self) -> str:
         string = 'Material\n'
         string += '{: <16}=\t{}\n'.format('\tID', self._id)
         string += '{: <16}=\t{}\n'.format('\tName', self._name)
@@ -166,34 +167,34 @@ class Material(IDManagerMixin):
         return string
 
     @property
-    def name(self):
+    def name(self) -> Optional[str]:
         return self._name
 
     @property
-    def temperature(self):
+    def temperature(self) -> Optional[float]:
         return self._temperature
 
     @property
-    def density(self):
+    def density(self) -> Optional[float]:
         return self._density
 
     @property
-    def density_units(self):
+    def density_units(self) -> str:
         return self._density_units
 
     @property
-    def depletable(self):
+    def depletable(self) -> bool:
         return self._depletable
 
     @property
-    def paths(self):
+    def paths(self) -> List[str]:
         if self._paths is None:
             raise ValueError('Material instance paths have not been determined. '
                              'Call the Geometry.determine_paths() method.')
         return self._paths
 
     @property
-    def num_instances(self):
+    def num_instances(self) -> int:
         if self._num_instances is None:
             raise ValueError(
                 'Number of material instances have not been determined. Call '
@@ -201,15 +202,15 @@ class Material(IDManagerMixin):
         return self._num_instances
 
     @property
-    def nuclides(self):
+    def nuclides(self) -> List[namedtuple]:
         return self._nuclides
 
     @property
-    def isotropic(self):
+    def isotropic(self) -> List[str]:
         return self._isotropic
 
     @property
-    def average_molar_mass(self):
+    def average_molar_mass(self) -> float:
         # Using the sum of specified atomic or weight amounts as a basis, sum
         # the mass and moles of the material
         mass = 0.
@@ -226,11 +227,11 @@ class Material(IDManagerMixin):
         return mass / moles
 
     @property
-    def volume(self):
+    def volume(self) -> Optional[float]:
         return self._volume
 
     @property
-    def ncrystal_cfg(self):
+    def ncrystal_cfg(self) -> Optional[str]:
         return self._ncrystal_cfg
 
     @name.setter
@@ -267,7 +268,7 @@ class Material(IDManagerMixin):
         self._isotropic = list(isotropic)
 
     @property
-    def fissionable_mass(self):
+    def fissionable_mass(self) -> float:
         if self.volume is None:
             raise ValueError("Volume must be set in order to determine mass.")
         density = 0.0
@@ -291,7 +292,7 @@ class Material(IDManagerMixin):
         return openmc.data.combine_distributions(dists, probs) if dists else None
 
     @classmethod
-    def from_hdf5(cls, group: h5py.Group):
+    def from_hdf5(cls, group: h5py.Group) -> Material:
         """Create material from HDF5 group
 
         Parameters
@@ -346,7 +347,7 @@ class Material(IDManagerMixin):
         return material
 
     @classmethod
-    def from_ncrystal(cls, cfg, **kwargs):
+    def from_ncrystal(cls, cfg, **kwargs) -> Material:
         """Create material from NCrystal configuration string.
 
         Density, temperature, and material composition, and (ultimately) thermal
@@ -881,7 +882,7 @@ class Material(IDManagerMixin):
     def make_isotropic_in_lab(self):
         self.isotropic = [x.name for x in self._nuclides]
 
-    def get_elements(self):
+    def get_elements(self) -> List[str]:
         """Returns all elements in the material
 
         .. versionadded:: 0.12
@@ -895,7 +896,7 @@ class Material(IDManagerMixin):
 
         return sorted({re.split(r'(\d+)', i)[0] for i in self.get_nuclides()})
 
-    def get_nuclides(self, element: Optional[str] = None):
+    def get_nuclides(self, element: Optional[str] = None) -> List[str]:
         """Returns a list of all nuclides in the material, if the element
         argument is specified then just nuclides of that element are returned.
 
@@ -925,7 +926,7 @@ class Material(IDManagerMixin):
 
         return matching_nuclides
 
-    def get_nuclide_densities(self):
+    def get_nuclide_densities(self) -> Dict[str, tuple]:
         """Returns all nuclides in the material and their densities
 
         Returns
@@ -944,7 +945,7 @@ class Material(IDManagerMixin):
 
         return nuclides
 
-    def get_nuclide_atom_densities(self, nuclide: Optional[str] = None):
+    def get_nuclide_atom_densities(self, nuclide: Optional[str] = None) -> Dict[str, float]:
         """Returns one or all nuclides in the material and their atomic
         densities in units of atom/b-cm
 
@@ -1030,7 +1031,7 @@ class Material(IDManagerMixin):
         return nuclides
 
     def get_activity(self, units: str = 'Bq/cm3', by_nuclide: bool = False,
-                     volume: Optional[float] = None):
+                     volume: Optional[float] = None) -> Union[Dict[str, float], float]:
         """Returns the activity of the material or for each nuclide in the
         material in units of [Bq], [Bq/g] or [Bq/cm3].
 
@@ -1077,7 +1078,7 @@ class Material(IDManagerMixin):
         return activity if by_nuclide else sum(activity.values())
 
     def get_decay_heat(self, units: str = 'W', by_nuclide: bool = False,
-                       volume: Optional[float] = None):
+                       volume: Optional[float] = None) -> Union[Dict[str, float], float]:
         """Returns the decay heat of the material or for each nuclide in the
         material in units of [W], [W/g] or [W/cm3].
 
@@ -1125,7 +1126,7 @@ class Material(IDManagerMixin):
 
         return decayheat if by_nuclide else sum(decayheat.values())
 
-    def get_nuclide_atoms(self, volume: Optional[float] = None):
+    def get_nuclide_atoms(self, volume: Optional[float] = None) -> Dict[str, float]:
         """Return number of atoms of each nuclide in the material
 
         .. versionadded:: 0.13.1
@@ -1154,7 +1155,7 @@ class Material(IDManagerMixin):
             atoms[nuclide] = 1.0e24 * atom_per_bcm * volume
         return atoms
 
-    def get_mass_density(self, nuclide: Optional[str] = None):
+    def get_mass_density(self, nuclide: Optional[str] = None) -> float:
         """Return mass density of one or all nuclides
 
         Parameters
@@ -1176,7 +1177,7 @@ class Material(IDManagerMixin):
             mass_density += density_i
         return mass_density
 
-    def get_mass(self, nuclide: Optional[str] = None, volume: Optional[float] = None):
+    def get_mass(self, nuclide: Optional[str] = None, volume: Optional[float] = None) -> float:
         """Return mass of one or all nuclides.
 
         Note that this method requires that the :attr:`Material.volume` has
@@ -1206,7 +1207,7 @@ class Material(IDManagerMixin):
             raise ValueError("Volume must be set in order to determine mass.")
         return volume*self.get_mass_density(nuclide)
 
-    def clone(self, memo: Optional[dict] = None):
+    def clone(self, memo: Optional[dict] = None) -> Material:
         """Create a copy of this material with a new unique ID.
 
         Parameters
@@ -1245,7 +1246,7 @@ class Material(IDManagerMixin):
 
         return memo[self]
 
-    def _get_nuclide_xml(self, nuclide: str):
+    def _get_nuclide_xml(self, nuclide: str) -> ET.Element:
         xml_element = ET.Element("nuclide")
         xml_element.set("name", nuclide.name)
 
@@ -1256,19 +1257,19 @@ class Material(IDManagerMixin):
 
         return xml_element
 
-    def _get_macroscopic_xml(self, macroscopic: str):
+    def _get_macroscopic_xml(self, macroscopic: str) -> ET.Element:
         xml_element = ET.Element("macroscopic")
         xml_element.set("name", macroscopic)
 
         return xml_element
 
-    def _get_nuclides_xml(self, nuclides: typing.Iterable[str]):
+    def _get_nuclides_xml(self, nuclides: typing.Iterable[str]) -> List[ET.Element]:
         xml_elements = []
         for nuclide in nuclides:
             xml_elements.append(self._get_nuclide_xml(nuclide))
         return xml_elements
 
-    def to_xml_element(self):
+    def to_xml_element(self) -> ET.Element:
         """Return XML representation of the material
 
         Returns
@@ -1338,7 +1339,7 @@ class Material(IDManagerMixin):
 
     @classmethod
     def mix_materials(cls, materials, fracs: typing.Iterable[float],
-                      percent_type: str = 'ao', name: Optional[str] = None):
+                      percent_type: str = 'ao', name: Optional[str] = None) -> Material:
         """Mix materials together based on atom, weight, or volume fractions
 
         .. versionadded:: 0.12
@@ -1436,7 +1437,7 @@ class Material(IDManagerMixin):
         return new_mat
 
     @classmethod
-    def from_xml_element(cls, elem: ET.Element):
+    def from_xml_element(cls, elem: ET.Element) -> Material:
         """Generate material from an XML element
 
         Parameters
@@ -1530,7 +1531,7 @@ class Materials(cv.CheckedList):
             self += materials
 
     @property
-    def cross_sections(self):
+    def cross_sections(self) -> Optional[Path]:
         return self._cross_sections
 
     @cross_sections.setter
@@ -1638,7 +1639,7 @@ class Materials(cv.CheckedList):
             self._write_xml(fh)
 
     @classmethod
-    def from_xml_element(cls, elem):
+    def from_xml_element(cls, elem) -> Material:
         """Generate materials collection from XML file
 
         Parameters
@@ -1665,7 +1666,7 @@ class Materials(cv.CheckedList):
         return materials
 
     @classmethod
-    def from_xml(cls, path: PathLike = 'materials.xml'):
+    def from_xml(cls, path: PathLike = 'materials.xml') -> Material:
         """Generate materials collection from XML file
 
         Parameters
diff --git a/openmc/model/model.py b/openmc/model/model.py
index 65aad6de9..933e1172e 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -1,3 +1,4 @@
+from __future__ import annotations
 from collections.abc import Iterable
 from contextlib import contextmanager
 from functools import lru_cache
@@ -7,6 +8,7 @@ from numbers import Integral
 from tempfile import NamedTemporaryFile
 import warnings
 from xml.etree import ElementTree as ET
+from typing import Optional, Dict
 
 import h5py
 
@@ -93,27 +95,27 @@ class Model:
             self.plots = plots
 
     @property
-    def geometry(self):
+    def geometry(self) -> Optional[openmc.Geometry]:
         return self._geometry
 
     @property
-    def materials(self):
+    def materials(self) -> Optional[openmc.Materials]:
         return self._materials
 
     @property
-    def settings(self):
+    def settings(self) -> Optional[openmc.Settings]:
         return self._settings
 
     @property
-    def tallies(self):
+    def tallies(self) -> Optional[openmc.Tallies]:
         return self._tallies
 
     @property
-    def plots(self):
+    def plots(self) -> Optional[openmc.Plots]:
         return self._plots
 
     @property
-    def is_initialized(self):
+    def is_initialized(self) -> bool:
         try:
             import openmc.lib
             return openmc.lib.is_initialized
@@ -122,7 +124,7 @@ class Model:
 
     @property
     @lru_cache(maxsize=None)
-    def _materials_by_id(self):
+    def _materials_by_id(self) -> dict:
         """Dictionary mapping material ID --> material"""
         if self.materials:
             mats = self.materials
@@ -132,14 +134,14 @@ class Model:
 
     @property
     @lru_cache(maxsize=None)
-    def _cells_by_id(self):
+    def _cells_by_id(self) -> dict:
         """Dictionary mapping cell ID --> cell"""
         cells = self.geometry.get_all_cells()
         return {cell.id: cell for cell in cells.values()}
 
     @property
     @lru_cache(maxsize=None)
-    def _cells_by_name(self):
+    def _cells_by_name(self) -> Dict[int, openmc.Cell]:
         # Get the names maps, but since names are not unique, store a set for
         # each name key. In this way when the user requests a change by a name,
         # the change will be applied to all of the same name.
@@ -152,7 +154,7 @@ class Model:
 
     @property
     @lru_cache(maxsize=None)
-    def _materials_by_name(self):
+    def _materials_by_name(self) -> Dict[int, openmc.Material]:
         if self.materials is None:
             mats = self.geometry.get_all_materials().values()
         else:
@@ -207,7 +209,7 @@ class Model:
     @classmethod
     def from_xml(cls, geometry='geometry.xml', materials='materials.xml',
                  settings='settings.xml', tallies='tallies.xml',
-                 plots='plots.xml'):
+                 plots='plots.xml') -> Model:
         """Create model from existing XML files
 
         Parameters
@@ -604,7 +606,6 @@ class Model:
         this method creates the XML files and runs OpenMC via a system call. In
         both cases this method returns the path to the last statepoint file
         generated.
-
         .. versionchanged:: 0.12
             Instead of returning the final k-effective value, this function now
             returns the path to the final statepoint written.
diff --git a/openmc/universe.py b/openmc/universe.py
index 4054485fd..67296dda8 100644
--- a/openmc/universe.py
+++ b/openmc/universe.py
@@ -7,6 +7,7 @@ from numbers import Integral, Real
 from pathlib import Path
 from tempfile import TemporaryDirectory
 from xml.etree import ElementTree as ET
+from warnings import warn
 
 import h5py
 import numpy as np
@@ -88,7 +89,8 @@ class UniverseBase(ABC, IDManagerMixin):
                 self._volume = volume_calc.volumes[self.id].n
                 self._atoms = volume_calc.atoms[self.id]
             else:
-                raise ValueError('No volume information found for this universe.')
+                raise ValueError(
+                    'No volume information found for this universe.')
         else:
             raise ValueError('No volume information found for this universe.')
 
@@ -169,7 +171,7 @@ class UniverseBase(ABC, IDManagerMixin):
             clone._cells = OrderedDict()
             for cell in self._cells.values():
                 clone.add_cell(cell.clone(clone_materials, clone_regions,
-                     memo))
+                                          memo))
 
             # Memoize the clone
             memo[self] = clone
@@ -294,9 +296,15 @@ class Universe(UniverseBase):
                     return [self, cell] + cell.fill.find(p)
         return []
 
+    # default kwargs that are passed to plt.legend in the plot method below.
+    _default_legend_kwargs = {'bbox_to_anchor': (
+        1.05, 1), 'loc': 2, 'borderaxespad': 0.0}
+
     def plot(self, origin=(0., 0., 0.), width=(1., 1.), pixels=(200, 200),
              basis='xy', color_by='cell', colors=None, seed=None,
-             openmc_exec='openmc', axes=None, outline=False, **kwargs):
+             openmc_exec='openmc', axes=None, legend=False,
+             legend_kwargs=_default_legend_kwargs, outline=False,
+             **kwargs):
         """Display a slice plot of the universe.
 
         Parameters
@@ -330,13 +338,20 @@ class Universe(UniverseBase):
             Axes to draw to
 
             .. versionadded:: 0.13.1
-        **kwargs
-            Keyword arguments passed to :func:`matplotlib.pyplot.imshow`
+        legend : bool
+            Whether a legend showing material or cell names should be drawn
 
+            .. versionadded:: 0.13.4
+        legend_kwargs : dict
+            Keyword arguments passed to :func:`matplotlib.pyplot.legend`.
+
+            .. versionadded:: 0.13.4
         outline : bool
             Whether outlines between color boundaries should be drawn
 
             .. versionadded:: 0.13.4
+        **kwargs
+            Keyword arguments passed to :func:`matplotlib.pyplot.imshow`
 
         Returns
         -------
@@ -345,15 +360,19 @@ class Universe(UniverseBase):
 
         """
         import matplotlib.image as mpimg
+        import matplotlib.patches as mpatches
         import matplotlib.pyplot as plt
 
         # Determine extents of plot
         if basis == 'xy':
             x, y = 0, 1
+            xlabel, ylabel = 'x [cm]', 'y [cm]'
         elif basis == 'yz':
             x, y = 1, 2
+            xlabel, ylabel = 'y [cm]', 'z [cm]'
         elif basis == 'xz':
             x, y = 0, 2
+            xlabel, ylabel = 'x [cm]', 'z [cm]'
         x_min = origin[x] - 0.5*width[0]
         x_max = origin[x] + 0.5*width[0]
         y_min = origin[y] - 0.5*width[1]
@@ -397,6 +416,8 @@ class Universe(UniverseBase):
             if axes is None:
                 px = 1/plt.rcParams['figure.dpi']
                 fig, axes = plt.subplots()
+                axes.set_xlabel(xlabel)
+                axes.set_ylabel(ylabel)
                 params = fig.subplotpars
                 width = pixels[0]*px/(params.right - params.left)
                 height = pixels[0]*px/(params.top - params.bottom)
@@ -417,6 +438,35 @@ class Universe(UniverseBase):
                     extent=(x_min, x_max, y_min, y_max),
                 )
 
+            # add legend showing which colors represent which material
+            # or cell if that was requested
+            if legend:
+                if plot.colors is None:
+                    raise ValueError("Must pass 'colors' dictionary if you "
+                                     "are adding a legend via legend=True.")
+
+                if color_by == "cell":
+                    expected_key_type = openmc.Cell
+                else:
+                    expected_key_type = openmc.Material
+
+                patches = []
+                for key, color in plot.colors.items():
+
+                    if isinstance(key, int):
+                        raise TypeError(
+                            "Cannot use IDs in colors dict for auto legend.")
+                    elif not isinstance(key, expected_key_type):
+                        raise TypeError(
+                            "Color dict key type does not match color_by")
+
+                    # this works whether we're doing cells or materials
+                    label = key.name if key.name != '' else key.id
+                    key_patch = mpatches.Patch(color=color, label=label)
+                    patches.append(key_patch)
+
+                axes.legend(handles=patches, **legend_kwargs)
+
             # Plot image and return the axes
             return axes.imshow(img, extent=(x_min, x_max, y_min, y_max), **kwargs)
 
@@ -759,8 +809,9 @@ class DAGMCUniverse(UniverseBase):
     @property
     def material_names(self):
         dagmc_file_contents = h5py.File(self.filename)
-        material_tags_hex=dagmc_file_contents['/tstt/tags/NAME'].get('values')
-        material_tags_ascii=[]
+        material_tags_hex = dagmc_file_contents['/tstt/tags/NAME'].get(
+            'values')
+        material_tags_ascii = []
         for tag in material_tags_hex:
             candidate_tag = tag.tobytes().decode().replace('\x00', '')
             # tags might be for temperature or reflective surfaces
@@ -926,7 +977,8 @@ class DAGMCUniverse(UniverseBase):
         openmc.Universe
             Universe instance
         """
-        bounding_cell = openmc.Cell(fill=self, cell_id=bounding_cell_id, region=self.bounding_region(**kwargs))
+        bounding_cell = openmc.Cell(
+            fill=self, cell_id=bounding_cell_id, region=self.bounding_region(**kwargs))
         return openmc.Universe(cells=[bounding_cell])
 
     @classmethod
diff --git a/openmc/weight_windows.py b/openmc/weight_windows.py
index 65b9b7bb0..2993db089 100644
--- a/openmc/weight_windows.py
+++ b/openmc/weight_windows.py
@@ -1,14 +1,18 @@
 from __future__ import annotations
 from collections.abc import Iterable
 from numbers import Real, Integral
-from typing import Iterable, List
+import pathlib
+import typing
+from typing import Iterable, List, Optional, Union, Dict
 
 from xml.etree import ElementTree as ET
 import numpy as np
+import h5py
 
 from openmc.filter import _PARTICLES
 from openmc.mesh import MeshBase, RectilinearMesh, UnstructuredMesh
 import openmc.checkvalue as cv
+from openmc.checkvalue import PathLike
 
 from ._xml import get_text
 from .mixin import IDManagerMixin
@@ -103,16 +107,20 @@ class WeightWindows(IDManagerMixin):
     next_id = 1
     used_ids = set()
 
-    def __init__(self, mesh, lower_ww_bounds,
-                 upper_ww_bounds=None,
-                 upper_bound_ratio=None,
-                 energy_bounds=None,
-                 particle_type='neutron',
-                 survival_ratio=3,
-                 max_lower_bound_ratio=None,
-                 max_split=10,
-                 weight_cutoff=1.e-38,
-                 id=None):
+    def __init__(
+        self, 
+        mesh: MeshBase, 
+        lower_ww_bounds: Iterable[float],
+        upper_ww_bounds: Optional[Iterable[float]] = None,
+        upper_bound_ratio: Optional[float] = None,
+        energy_bounds: Optional[Iterable[Real]] = None,
+        particle_type: str = 'neutron',
+        survival_ratio: float = 3,
+        max_lower_bound_ratio: Optional[float] = None,
+        max_split: int = 10,
+        weight_cutoff: float = 1.e-38,
+        id: Optional[int] = None
+    ):
         self.mesh = mesh
         self.id = id
         self.particle_type = particle_type
@@ -161,7 +169,7 @@ class WeightWindows(IDManagerMixin):
         string += '{: <16}=\t{}\n'.format('\tWeight Cutoff', self._weight_cutoff)
         return string
 
-    def __eq__(self, other) -> bool:
+    def __eq__(self, other: WeightWindows) -> bool:
         # ensure that `other` is a WeightWindows object
         if not isinstance(other, WeightWindows):
             return False
@@ -195,7 +203,7 @@ class WeightWindows(IDManagerMixin):
         return self._mesh
 
     @mesh.setter
-    def mesh(self, mesh):
+    def mesh(self, mesh: MeshBase):
         cv.check_type('Weight window mesh', mesh, MeshBase)
         self._mesh = mesh
 
@@ -204,7 +212,7 @@ class WeightWindows(IDManagerMixin):
         return self._particle_type
 
     @particle_type.setter
-    def particle_type(self, pt):
+    def particle_type(self, pt: str):
         cv.check_value('Particle type', pt, _PARTICLES)
         self._particle_type = pt
 
@@ -213,7 +221,7 @@ class WeightWindows(IDManagerMixin):
         return self._energy_bounds
 
     @energy_bounds.setter
-    def energy_bounds(self, bounds):
+    def energy_bounds(self, bounds: Iterable[float]):
         cv.check_type('Energy bounds', bounds, Iterable, Real)
         self._energy_bounds = np.asarray(bounds)
 
@@ -228,7 +236,7 @@ class WeightWindows(IDManagerMixin):
         return self._lower_ww_bounds
 
     @lower_ww_bounds.setter
-    def lower_ww_bounds(self, bounds):
+    def lower_ww_bounds(self, bounds: Iterable[float]):
         cv.check_iterable_type('Lower WW bounds',
                                bounds,
                                Real,
@@ -247,7 +255,7 @@ class WeightWindows(IDManagerMixin):
         return self._upper_ww_bounds
 
     @upper_ww_bounds.setter
-    def upper_ww_bounds(self, bounds):
+    def upper_ww_bounds(self, bounds: Iterable[float]):
         cv.check_iterable_type('Upper WW bounds',
                                bounds,
                                Real,
@@ -266,7 +274,7 @@ class WeightWindows(IDManagerMixin):
         return self._survival_ratio
 
     @survival_ratio.setter
-    def survival_ratio(self, val):
+    def survival_ratio(self, val: float):
         cv.check_type('Survival ratio', val, Real)
         cv.check_greater_than('Survival ratio', val, 1.0, True)
         self._survival_ratio = val
@@ -276,7 +284,7 @@ class WeightWindows(IDManagerMixin):
         return self._max_lower_bound_ratio
 
     @max_lower_bound_ratio.setter
-    def max_lower_bound_ratio(self, val):
+    def max_lower_bound_ratio(self, val: float):
         cv.check_type('Maximum lower bound ratio', val, Real)
         cv.check_greater_than('Maximum lower bound ratio', val, 1.0)
         self._max_lower_bound_ratio = val
@@ -286,7 +294,7 @@ class WeightWindows(IDManagerMixin):
         return self._max_split
 
     @max_split.setter
-    def max_split(self, val):
+    def max_split(self, val: int):
         cv.check_type('Max split', val, Integral)
         self._max_split = val
 
@@ -295,7 +303,7 @@ class WeightWindows(IDManagerMixin):
         return self._weight_cutoff
 
     @weight_cutoff.setter
-    def weight_cutoff(self, cutoff):
+    def weight_cutoff(self, cutoff: float):
         cv.check_type('Weight cutoff', cutoff, Real)
         cv.check_greater_than('Weight cutoff', cutoff, 0.0, True)
         self._weight_cutoff = cutoff
@@ -343,7 +351,7 @@ class WeightWindows(IDManagerMixin):
         return element
 
     @classmethod
-    def from_xml_element(cls, elem, root) -> WeightWindows:
+    def from_xml_element(cls, elem: ET.Element, root: ET.Element) -> WeightWindows:
         """Generate weight window settings from an XML element
 
         Parameters
@@ -398,7 +406,7 @@ class WeightWindows(IDManagerMixin):
         )
 
     @classmethod
-    def from_hdf5(cls, group, meshes) -> WeightWindows:
+    def from_hdf5(cls, group: h5py.Group, meshes: Dict[int, MeshBase]) -> WeightWindows:
         """Create weight windows from HDF5 group
 
         Parameters
@@ -443,7 +451,7 @@ class WeightWindows(IDManagerMixin):
         )
 
 
-def wwinp_to_wws(path) -> List[WeightWindows]:
+def wwinp_to_wws(path: PathLike) -> List[WeightWindows]:
     """Create WeightWindows instances from a wwinp file
 
     .. versionadded:: 0.13.1
diff --git a/src/tallies/trigger.cpp b/src/tallies/trigger.cpp
index 79af65879..87f298baa 100644
--- a/src/tallies/trigger.cpp
+++ b/src/tallies/trigger.cpp
@@ -39,7 +39,8 @@ std::pair<double, double> get_tally_uncertainty(
   auto mean = sum / n;
 
   // if the result has no contributions, return an invalid pair
-  if (mean == 0) return {-1 , -1};
+  if (mean == 0)
+    return {-1, -1};
 
   double std_dev = std::sqrt((sum_sq / n - mean * mean) / (n - 1));
   double rel_err = (mean != 0.) ? std_dev / std::abs(mean) : 0.;
@@ -193,8 +194,10 @@ void check_triggers()
       keff_ratio);
   } else {
     if (tally_ratio == INFINITY) {
-      msg = fmt::format("Triggers unsatisfied, no result tallied for score {} in tally {}", reaction_name(score), tally_id);
-    } else{
+      msg = fmt::format(
+        "Triggers unsatisfied, no result tallied for score {} in tally {}",
+        reaction_name(score), tally_id);
+    } else {
       msg = fmt::format(
         "Triggers unsatisfied, max unc./thresh. is {} for {} in tally {}",
         tally_ratio, reaction_name(score), tally_id);
@@ -211,13 +214,20 @@ void check_triggers()
     auto n_pred_batches = static_cast<int>(n_active * max_ratio * max_ratio) +
                           settings::n_inactive + 1;
 
-    std::string msg =
-      fmt::format("The estimated number of batches is {}", n_pred_batches);
-    if (n_pred_batches > settings::n_max_batches) {
-      msg.append(" --- greater than max batches");
-      warning(msg);
-    } else {
+    if (max_ratio == INFINITY) {
+      std::string msg =
+        fmt::format("One or more tallies with triggers have no scores. Unable "
+                    "to estimate the number of remaining batches.");
       write_message(msg, 7);
+    } else {
+      std::string msg =
+        fmt::format("The estimated number of batches is {}", n_pred_batches);
+      if (n_pred_batches > settings::n_max_batches) {
+        msg.append(" --- greater than max batches");
+        warning(msg);
+      } else {
+        write_message(msg, 7);
+      }
     }
   }
 }