diff --git a/openmc/lib/plot.py b/openmc/lib/plot.py
index 471ae688b..196682394 100644
--- a/openmc/lib/plot.py
+++ b/openmc/lib/plot.py
@@ -87,7 +87,7 @@ _dll.openmc_slice_data.errcheck = _error_handler
 
 
 def slice_data(origin, width=None, basis='xy', u_span=None, v_span=None,
-                pixels=None, show_overlaps=False, level=-1, filter=None,
+                pixels=None, show_overlaps=False, level=None, filter=None,
                 include_properties=True):
     """Generate a 2D raster of geometry and property data for plotting.
 
@@ -111,7 +111,7 @@ def slice_data(origin, width=None, basis='xy', u_span=None, v_span=None,
     show_overlaps : bool, optional
         Whether to detect overlapping cells
     level : int, optional
-        Universe level (-1 for deepest)
+        Universe level (None for deepest)
     filter : openmc.lib.Filter, optional
         Filter for bin index lookup
     include_properties : bool, optional
@@ -127,6 +127,12 @@ def slice_data(origin, width=None, basis='xy', u_span=None, v_span=None,
         Array of shape (v_res, h_res, 2) with float64 dtype containing
         [temperature, density], or None if include_properties=False
     """
+    # Set deepest level as default
+    if level is None:
+        level = -1
+    if not isinstance(level, int):
+        raise TypeError("level must be an integer.")
+
     if pixels is None:
         raise ValueError("pixels must be specified.")
     if len(pixels) != 2:
diff --git a/src/plot.cpp b/src/plot.cpp
index 06c2b550a..96df1c5be 100644
--- a/src/plot.cpp
+++ b/src/plot.cpp
@@ -88,10 +88,7 @@ void PropertyData::set_value(size_t y, size_t x, const Particle& p, int level,
 {
   Cell* c = model::cells.at(p.lowest_coord().cell()).get();
   data_(y, x, 0) = (p.sqrtkT() * p.sqrtkT()) / K_BOLTZMANN;
-  if (c->type_ != Fill::UNIVERSE && p.material() != MATERIAL_VOID) {
-    Material* m = model::materials.at(p.material()).get();
-    data_(y, x, 1) = m->density_gpcc_;
-  }
+  data_(y, x, 1) = c->density(p.cell_instance());
 }
 
 void PropertyData::set_overlap(size_t y, size_t x, int /*overlap_idx*/)
@@ -150,7 +147,7 @@ void RasterData::set_value(size_t y, size_t x, const Particle& p, int level,
   // set density (g/cm³)
   if (c->type_ != Fill::UNIVERSE && p.material() != MATERIAL_VOID) {
     Material* m = model::materials.at(p.material()).get();
-    property_data_(y, x, 1) = m->density_gpcc_;
+    property_data_(y, x, 1) = c->density(p.cell_instance());
   }
 }
 
diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py
index f62306e83..4cfae0df2 100644
--- a/tests/unit_tests/test_lib.py
+++ b/tests/unit_tests/test_lib.py
@@ -926,6 +926,22 @@ def test_property_map(lib_init):
     assert np.allclose(expected_properties, properties, atol=1e-04)
 
 
+def test_slice_data(lib_init):
+    expected_properties = np.array(
+        [[(293.6, 0.740582), (293.6, 6.55), (293.6, 0.740582)],
+         [ (293.6, 6.55), (293.6, 10.29769),  (293.6, 6.55)],
+         [(293.6, 0.740582), (293.6, 6.55), (293.6, 0.740582)]], dtype='float')
+    origin = (0.0, 0.0, 0.0)
+    _, properties = openmc.lib.slice_data(
+        origin,
+        width=(1.26, 1.26),
+        basis='xy',
+        pixels=(3, 3),
+        include_properties=True
+    )
+    assert np.allclose(expected_properties, properties, atol=1e-04)
+
+
 def test_solid_raytrace_plot(lib_init, pincell_model):
     # Ensure plot mapping can be accessed and grows after allocation
     n0 = len(openmc.lib.plots)