@@ -303,7 +303,7 @@ def find(self, point):
 
     def plot(self, origin=None, width=None, pixels=40000,
              basis='xy', color_by='cell', colors=None, seed=None,
-             openmc_exec='openmc', axes=None, legend=False,
+             openmc_exec='openmc', axes=None, legend=False, axis_units='cm',
              legend_kwargs=_default_legend_kwargs, outline=False,
              **kwargs):
         """Display a slice plot of the universe.
@@ -360,6 +360,10 @@ def plot(self, origin=None, width=None, pixels=40000,
         outline : bool
             Whether outlines between color boundaries should be drawn
 
+            .. versionadded:: 0.13.4
+        axis_units : {'km', 'm', 'cm', 'mm'}
+            Units used on the plot axis
+
             .. versionadded:: 0.13.4
         **kwargs
             Keyword arguments passed to :func:`matplotlib.pyplot.imshow`
@@ -377,13 +381,13 @@ def plot(self, origin=None, width=None, pixels=40000,
         # Determine extents of plot
         if basis == 'xy':
             x, y = 0, 1
-            xlabel, ylabel = 'x [cm]', 'y [cm]'
+            xlabel, ylabel = f'x [{axis_units}]', f'y [{axis_units}]'
         elif basis == 'yz':
             x, y = 1, 2
-            xlabel, ylabel = 'y [cm]', 'z [cm]'
+            xlabel, ylabel = f'y [{axis_units}]', f'z [{axis_units}]'
         elif basis == 'xz':
             x, y = 0, 2
-            xlabel, ylabel = 'x [cm]', 'z [cm]'
+            xlabel, ylabel = f'x [{axis_units}]', f'z [{axis_units}]'
 
         bb = self.bounding_box
         # checks to see if bounding box contains -inf or inf values
@@ -408,10 +412,12 @@ def plot(self, origin=None, width=None, pixels=40000,
             pixels_y = math.sqrt(pixels / aspect_ratio)
             pixels = (int(pixels / pixels_y), int(pixels_y))
 
-        x_min = origin[x] - 0.5*width[0]
-        x_max = origin[x] + 0.5*width[0]
-        y_min = origin[y] - 0.5*width[1]
-        y_max = origin[y] + 0.5*width[1]
+        axis_scaling_factor = {'km': 0.00001, 'm': 0.01, 'cm': 1, 'mm': 10}
+
+        x_min = (origin[x] - 0.5*width[0]) * axis_scaling_factor[axis_units]
+        x_max = (origin[x] + 0.5*width[0]) * axis_scaling_factor[axis_units]
+        y_min = (origin[y] - 0.5*width[1]) * axis_scaling_factor[axis_units]
+        y_max = (origin[y] + 0.5*width[1]) * axis_scaling_factor[axis_units]
 
         with TemporaryDirectory() as tmpdir:
             model = openmc.Model()