diff --git a/Project.toml b/Project.toml
index e4bf085..aa4219f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,6 +5,7 @@ version = "0.2.1-pre"
 
 [deps]
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
+EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
 GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
@@ -13,6 +14,7 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [compat]
 ArgParse = "1.1"
+EllipsisNotation = "1"
 GR = "0.51, 0.52"
 Glob = "1.3"
 HDF5 = "0.13"
diff --git a/src/Trixi2Img.jl b/src/Trixi2Img.jl
index 7da82fc..f909c3a 100644
--- a/src/Trixi2Img.jl
+++ b/src/Trixi2Img.jl
@@ -1,6 +1,7 @@
 module Trixi2Img
 
 # Include other packages
+using EllipsisNotation
 using Glob: glob
 using HDF5: h5open, attrs, exists
 using Plots: plot, plot!, gr, savefig, contourf!
diff --git a/src/convert.jl b/src/convert.jl
index 0c7980a..0341ebb 100644
--- a/src/convert.jl
+++ b/src/convert.jl
@@ -18,16 +18,22 @@ Convert two-dimensional Trixi-generated output files to image files (PNG or PDF)
 - `nvisnodes`: Number of visualization nodes per element (default: twice the number of DG nodes).
                A value of `0` (zero) uses the number of nodes in the DG elements.
 - `max_supported_level`: Maximum cell refinement level supported for plotting.
+- `slice_axis`: Axis orthogonal to the slice plane. Will be ignored in 2D. May be :x, :y or :z.
+- `slice_axis_intercept`: Point on the slice axis where it intersects with the slice plane.
 
 # Examples
 ```julia
 julia> trixi2img("out/solution_000*.h5")
 [...]
+
+julia> trixi2img("out/solution_000*.h5", slice_axis=:z, slice_axis_intercept=0.5) # Slice plane will be z=0.5
+[...]
 ```
 """
 function trixi2img(filename::AbstractString...;
                    format=:png, variables=[], verbose=false, grid_lines=false,
-                   output_directory=".", nvisnodes=nothing, max_supported_level=11)
+                   output_directory=".", nvisnodes=nothing, max_supported_level=11,
+                   slice_axis=:z, slice_axis_intercept=0)
   # Reset timer
   reset_timer!()
 
@@ -65,7 +71,7 @@ function trixi2img(filename::AbstractString...;
     # Read mesh
     verbose && println("| Reading mesh file...")
     @timeit "read mesh" (center_level_0, length_level_0,
-                         leaf_cells, coordinates, levels) = read_meshfile(meshfile)
+                         leaf_cells, coordinates, levels, ndims) = read_meshfile(meshfile)
 
     # Read data
     verbose && println("| Reading data file...")
@@ -78,7 +84,7 @@ function trixi2img(filename::AbstractString...;
             "maximum supported level $max_supported_level")
     end
     max_available_nodes_per_finest_element = 2^(max_supported_level - max_level)
-    if nvisnodes == nothing
+    if nvisnodes === nothing
       max_nvisnodes = 2 * n_nodes
     elseif nvisnodes == 0
       max_nvisnodes = n_nodes
@@ -92,6 +98,14 @@ function trixi2img(filename::AbstractString...;
     # level-0-cell) that contains the number of visualization nodes for any
     # refinement level to visualize on an equidistant grid
 
+    if ndims == 3
+      verbose && println("| Extracting 2D slice...")
+      # convert 3d unstructured data to 2d slice
+      @timeit "extract 2D slice" (unstructured_data, coordinates, levels,
+          center_level_0) = unstructured_2d_to_3d(unstructured_data, coordinates,
+          levels, length_level_0, center_level_0, slice_axis, slice_axis_intercept)
+    end
+
     # Normalize element coordinates: move center to (0, 0) and domain size to [-1, 1]²
     n_elements = length(levels)
     normalized_coordinates = similar(coordinates)
@@ -180,4 +194,3 @@ function trixi2img(filename::AbstractString...;
   print_timer()
   println()
 end
-
diff --git a/src/interpolate.jl b/src/interpolate.jl
index 52d7bde..e3f6804 100644
--- a/src/interpolate.jl
+++ b/src/interpolate.jl
@@ -40,26 +40,137 @@ function cell2node(cell_centered_data::AbstractArray{Float64})
 end
 
 
+# Convert 3d unstructured data to 2d slice.
+# Additional to the new unstructured data updated coordinates, levels and
+# center coordinates are returned.
+function unstructured_2d_to_3d(unstructured_data, coordinates, levels,
+                               length_level_0, center_level_0, slice_axis,
+                               slice_axis_intercept)
+
+  if slice_axis === :x
+    slice_axis_dimension = 1
+    other_dimensions = [2, 3]
+  elseif slice_axis === :y
+    slice_axis_dimension = 2
+    other_dimensions = [1, 3]
+  elseif slice_axis === :z
+    slice_axis_dimension = 3
+    other_dimensions = [1, 2]
+  else
+    error("illegal dimension '$slice_axis', supported dimensions are :x, :y, and :z")
+  end
+
+  # Limits of domain in slice_axis dimension
+  lower_limit = center_level_0[slice_axis_dimension] - length_level_0 / 2
+  upper_limit = center_level_0[slice_axis_dimension] + length_level_0 / 2
+
+  if slice_axis_intercept < lower_limit || slice_axis_intercept > upper_limit
+    error(string("Slice plane $slice_axis = $slice_axis_intercept outside of domain. ",
+        "$slice_axis must be between $lower_limit and $upper_limit"))
+  end
+
+  # Extract data shape information
+  n_nodes_in, _, _, n_elements, n_variables = size(unstructured_data)
+
+  # Get node coordinates for DG locations on reference element
+  nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in)
+
+  # New unstructured data has one dimension less.
+  # The redundant element ids are removed later.
+  @views new_unstructured_data = similar(unstructured_data[1, ..])
+
+  # Declare new empty arrays to fill in new coordinates and levels
+  new_coordinates = Array{Float64}(undef, 2, n_elements)
+  new_levels = Array{eltype(levels)}(undef, n_elements)
+
+  # Counter for new element ids
+  new_id = 0
+
+  # Save vandermonde matrices in a Dict to prevent redundant generation
+  vandermonde_to_2d = Dict()
+
+  # Permute dimensions such that the slice axis dimension is always the
+  # third dimension of the array. Below we can always interpolate in the 
+  # third dimension.
+  if slice_axis === :x
+    unstructured_data = permutedims(unstructured_data, [2, 3, 1, 4, 5])
+  elseif slice_axis === :y
+    unstructured_data = permutedims(unstructured_data, [1, 3, 2, 4, 5])
+  end
+
+  for element_id in 1:n_elements
+    # Distance from center to border of this element (half the length)
+    element_length = length_level_0 / 2^levels[element_id]
+    min_coordinate = coordinates[:, element_id] .- element_length / 2
+    max_coordinate = coordinates[:, element_id] .+ element_length / 2
+
+    # Check if slice plane and current element intersect.
+    # The first check uses a "greater but not equal" to only match one cell if the
+    # slice plane lies between two cells.
+    # The second check is needed if the slice plane is at the upper border of
+    # the domain due to this.
+    if !((min_coordinate[slice_axis_dimension] <= slice_axis_intercept &&
+          max_coordinate[slice_axis_dimension] > slice_axis_intercept) ||
+        (slice_axis_intercept == upper_limit &&
+          max_coordinate[slice_axis_dimension] == upper_limit))
+      # Continue for loop if they don't intersect
+      continue
+    end
+
+    # This element is of interest
+    new_id += 1
+
+    # Add element to new coordinates and levels
+    new_coordinates[:, new_id] = coordinates[other_dimensions, element_id]
+    new_levels[new_id] = levels[element_id]
+
+    # Construct vandermonde matrix (or load from Dict if possible)
+    normalized_intercept =
+        (slice_axis_intercept - min_coordinate[slice_axis_dimension]) /
+        element_length * 2 - 1
+
+    if haskey(vandermonde_to_2d, normalized_intercept)
+      vandermonde = vandermonde_to_2d[normalized_intercept]
+    else
+      # Generate vandermonde matrix to interpolate values at nodes_in to one value
+      vandermonde = polynomial_interpolation_matrix(nodes_in, [normalized_intercept])
+      vandermonde_to_2d[normalized_intercept] = vandermonde
+    end
+
+    # 1D interpolation to specified slice plane
+    # We permuted the dimensions above such that now the dimension in which
+    # we will interpolate is always the third one.
+    for i in 1:n_nodes_in
+      for ii in 1:n_nodes_in
+        # Interpolate in the third dimension
+        data = unstructured_data[i, ii, :, element_id, :]
+    
+        value = interpolate_nodes(permutedims(data), vandermonde, n_variables)
+        new_unstructured_data[i, ii, new_id, :] = value[:, 1]
+      end
+    end
+  end
+
+  # Remove redundant element ids
+  unstructured_data = new_unstructured_data[:, :, 1:new_id, :]
+  new_coordinates = new_coordinates[:, 1:new_id]
+  new_levels = new_levels[1:new_id]
+
+  center_level_0 = center_level_0[other_dimensions]
+
+  return unstructured_data, new_coordinates, new_levels, center_level_0
+end
+
+
 # Interpolate unstructured DG data to structured data (cell-centered)
-function unstructured2structured(unstructured_data::AbstractArray{Float64},
-                                 normalized_coordinates::AbstractArray{Float64},
-                                 levels::AbstractArray{Int}, resolution::Int,
-                                 nvisnodes_per_level::AbstractArray{Int})
+function unstructured2structured(unstructured_data, normalized_coordinates,
+                                 levels, resolution, nvisnodes_per_level)
   # Extract data shape information
   n_nodes_in, _, n_elements, n_variables = size(unstructured_data)
 
   # Get node coordinates for DG locations on reference element
   nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in)
 
-  #=# Calculate node coordinates for structured locations on reference element=#
-  #=max_level = length(nvisnodes_per_level) - 1=#
-  #=visnodes_per_level = []=#
-  #=for l in 0:max_level=#
-  #=  n_nodes_out = nvisnodes_per_level[l + 1]=#
-  #=  dx = 2 / n_nodes_out=#
-  #=  push!(visnodes_per_level, collect(range(-1 + dx/2, 1 - dx/2, length=n_nodes_out)))=#
-  #=end=#
-
   # Calculate interpolation vandermonde matrices for each level
   max_level = length(nvisnodes_per_level) - 1
   vandermonde_per_level = []
@@ -241,4 +352,3 @@ function calc_vertices(coordinates::AbstractArray{Float64, 2},
 
   return x, y
 end
-
diff --git a/src/io.jl b/src/io.jl
index 7447e0d..5d8f184 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -12,22 +12,22 @@ end
 
 # Read in mesh file and return relevant data
 function read_meshfile(filename::String)
-  n_children_per_cell = 2^ndim
   # Open file for reading
   h5open(filename, "r") do file
     # Extract basic information
     if exists(attrs(file), "ndims")
-      ndims_ = read(attrs(file)["ndims"])
+      ndims = read(attrs(file)["ndims"])
     else
-      ndims_ = read(attrs(file)["ndim"]) # FIXME once Trixi's 3D branch is merged & released
+      ndims = read(attrs(file)["ndim"]) # FIXME once Trixi's 3D branch is merged & released
     end
+    n_children_per_cell = 2^ndims
     n_cells = read(attrs(file)["n_cells"])
     n_leaf_cells = read(attrs(file)["n_leaf_cells"])
     center_level_0 = read(attrs(file)["center_level_0"])
     length_level_0 = read(attrs(file)["length_level_0"])
 
     # Extract coordinates, levels, child cells
-    coordinates = Array{Float64}(undef, ndim, n_cells)
+    coordinates = Array{Float64}(undef, ndims, n_cells)
     coordinates .= read(file["coordinates"])
     levels = Array{Int}(undef, n_cells)
     levels .= read(file["levels"])
@@ -50,7 +50,7 @@ function read_meshfile(filename::String)
     coordinates = coordinates[:, leaf_cells]
     levels = levels[leaf_cells]
 
-    return center_level_0, length_level_0, leaf_cells, coordinates, levels
+    return center_level_0, length_level_0, leaf_cells, coordinates, levels, ndims
   end
 end
 
@@ -58,6 +58,7 @@ end
 function read_datafile(filename::String)
   # Open file for reading
   h5open(filename, "r") do file
+    ndims = read(attrs(file)["ndims"])
     # Extract basic information
     if exists(attrs(file), "polydeg")
       polydeg = read(attrs(file)["polydeg"])
@@ -76,14 +77,21 @@ function read_datafile(filename::String)
 
     # Extract data arrays
     n_nodes = polydeg + 1
-    data = Array{Float64}(undef, n_nodes, n_nodes, n_elements, n_variables)
+
+    if ndims == 3
+      # Read 3d data
+      data = Array{Float64}(undef, n_nodes, n_nodes, n_nodes, n_elements, n_variables)
+    elseif ndims == 2
+      # Read 2d data
+      data = Array{Float64}(undef, n_nodes, n_nodes, n_elements, n_variables)
+    else
+      error("unsupported number of dimensions: $ndims")
+    end
     for v = 1:n_variables
       vardata = read(file["variables_$v"])
-      #=@show minimum(vardata), maximum(vardata)=#
-      @views data[:, :, :, v][:] .= vardata
+      @views data[.., :, v][:] .= vardata
     end
 
     return labels, data, n_nodes, time
   end
 end
-
diff --git a/test/runtests.jl b/test/runtests.jl
index 82e02aa..49a2d05 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1 +1,2 @@
 include("test_2d.jl")
+include("test_3d.jl")
diff --git a/test/test_3d.jl b/test/test_3d.jl
new file mode 100644
index 0000000..ba493de
--- /dev/null
+++ b/test/test_3d.jl
@@ -0,0 +1,80 @@
+module Test3D
+
+using Test
+using Trixi2Img
+
+include("test_trixi2img.jl")
+
+# Start with a clean environment: remove Trixi output directory if it exists
+outdir = "out"
+isdir(outdir) && rm(outdir, recursive=true)
+
+
+@testset "3D" begin
+  run_trixi(joinpath("3d", "parameters.toml"), n_steps_max=1)
+
+  @testset "uniform mesh as PNG" begin
+    @testset "default slice" begin
+      test_trixi2img("solution_000000.h5", outdir)
+        # To be able to compare hashes, we would need to make sure that exactly the same versions
+        # of all required libraries are installed everywhere, which does not seem to be a good option
+        # right now.
+    end
+
+    @testset "x-axis slice" begin
+      test_trixi2img("solution_000000.h5", outdir, slice_axis=:x, slice_axis_intercept=0.6)
+    end
+
+    @testset "y-axis slice" begin
+      test_trixi2img("solution_000000.h5", outdir, slice_axis=:y, slice_axis_intercept=-0.35)
+    end
+
+    @testset "z-axis slice" begin
+      test_trixi2img("solution_000000.h5", outdir, slice_axis=:z, slice_axis_intercept=-0.89)
+    end
+
+    @testset "negative border slice" begin
+      test_trixi2img("solution_000000.h5", outdir, slice_axis=:y, slice_axis_intercept=-1)
+    end
+
+    @testset "positive border slice" begin
+      test_trixi2img("solution_000000.h5", outdir, slice_axis=:z, slice_axis_intercept=1)
+    end
+
+    @testset "outside negative border slice" begin
+      slice_axis_intercept=-1.01
+      slice_axis=:x
+      @test_throws ErrorException(string(
+          "Slice plane $slice_axis = $slice_axis_intercept outside of domain. ",
+          "$slice_axis must be between -1.0 and 1.0")) trixi2img(
+          joinpath(outdir, "solution_000000.h5"); output_directory=outdir,
+          slice_axis=slice_axis, slice_axis_intercept=slice_axis_intercept)
+    end
+
+    @testset "outside positive border slice" begin
+      slice_axis_intercept=1.005
+      slice_axis=:y
+      @test_throws ErrorException(string(
+          "Slice plane $slice_axis = $slice_axis_intercept outside of domain. ",
+          "$slice_axis must be between -1.0 and 1.0")) trixi2img(
+          joinpath(outdir, "solution_000000.h5"); output_directory=outdir,
+          slice_axis=slice_axis, slice_axis_intercept=slice_axis_intercept)
+    end
+  end
+
+  # PDF tests do not work on Windows
+  if !Sys.iswindows()
+    @testset "uniform mesh as PDF with grid lines" begin
+      test_trixi2img("solution_000000.h5", outdir,
+          format=:pdf, grid_lines=true)
+          # To be able to compare hashes, we would need to make sure that exactly the same versions
+          # of all required libraries are installed everywhere, which does not seem to be a good option
+          # right now.
+    end
+  end
+end
+
+# Clean up afterwards: delete Trixi output directory
+@test_nowarn rm(outdir, recursive=true)
+
+end