diff --git a/.gitignore b/.gitignore
index c722179..5ee6cc2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -21,6 +21,8 @@ MANIFEST
 # Sphinx
 docs/api
 docs/_build
+docs/gallery
+docs/gen_modules
 
 # Eclipse editor project files
 .project
diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index 7796a44..713d601 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -58,9 +58,9 @@ stages:
         - linux: py38-oldestdeps
           libraries: {}
 
-        - linux: build_docs
-          posargs: " "
-          pytest: false
+#        - linux: build_docs
+#          posargs: " "
+#          pytest: false
 
   - stage: FigureTests
     displayName: Figure Tests
diff --git a/docs/conf.py b/docs/conf.py
index 11e673e..332ee1d 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -31,6 +31,7 @@
 # -- General configuration ---------------------------------------------------
 
 import sphinx_rtd_theme
+from sphinx_gallery.sorting import FileNameSortKey
 
 # Add any Sphinx extension module names here, as strings. They can be
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
@@ -47,6 +48,7 @@
     'sphinx.ext.napoleon',
     'sphinx.ext.todo',
     'sphinx.ext.viewcode',
+    'sphinx_gallery.gen_gallery',
 ]
 numpydoc_show_class_members = False
 nbsphinx_execute = 'never'
@@ -69,6 +71,22 @@
     "sklearn": ('https://scikit-learn.org/stable', None),
 }
 
+# sphinx-gallery configuration
+sphinx_gallery_conf = {
+    # path to your example scripts
+    'examples_dirs': ['../examples/gallery'],
+    # path to where to save gallery generated output
+    'gallery_dirs': ['gallery'],
+    # specify that examples should be ordered according to filename
+    'within_subsection_order': FileNameSortKey,
+    # directory where function/class granular galleries are stored
+    'backreferences_dir': 'gen_modules/backreferences',
+    # Modules for which function/class level galleries are created.
+    'doc_module': ('mcalf',),
+    # don't print e.g. sphinx_gallery_thumbnail_number = 2
+    'remove_config_comments': True,
+}
+
 # -- Options for HTML output -------------------------------------------------
 
 # The theme to use for HTML and HTML Help pages.  See the documentation for
diff --git a/docs/index.rst b/docs/index.rst
index ba63f0c..df6f951 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -20,6 +20,7 @@ a :any:`Module Index <modindex>` are available.
    :caption: Contents:
 
    guide/index
+   gallery/index
    code_ref/index
    code_of_conduct
    license
diff --git a/docs/requirements.txt b/docs/requirements.txt
index d112f87..7b6e023 100644
--- a/docs/requirements.txt
+++ b/docs/requirements.txt
@@ -4,3 +4,5 @@ sphinx==3.1.2
 ipykernel
 nbsphinx==0.7.1
 sphinx-rtd-theme==0.5.0
+sphinx-gallery
+pytest
diff --git a/examples/gallery/README.rst b/examples/gallery/README.rst
new file mode 100644
index 0000000..b6ef39c
--- /dev/null
+++ b/examples/gallery/README.rst
@@ -0,0 +1,4 @@
+Example Gallery
+***************
+
+Here are a collection of examples on how this package can be used.
diff --git a/examples/gallery/visualisation/README.rst b/examples/gallery/visualisation/README.rst
new file mode 100644
index 0000000..eb8328b
--- /dev/null
+++ b/examples/gallery/visualisation/README.rst
@@ -0,0 +1,4 @@
+Visualisation
+=============
+
+Below are examples of plots produced using functions within the visualisation module:
diff --git a/examples/gallery/visualisation/plot_bar.py b/examples/gallery/visualisation/plot_bar.py
new file mode 100644
index 0000000..03984d8
--- /dev/null
+++ b/examples/gallery/visualisation/plot_bar.py
@@ -0,0 +1,94 @@
+"""
+Plot a bar chart of classifications
+===================================
+This is an example showing how to produce a bar chart showing the percentage
+abundance of each classification in a 2D or 3D array of classifications.
+"""
+
+#%%
+# First we shall create a random 3D grid of classifications that can be plotted.
+# Usually you would use a method such as
+# :meth:`mcalf.models.ModelBase.classify_spectra`
+# to classify an array of spectra.
+
+from mcalf.tests.helpers import class_map as c
+
+t = 3  # Three images
+x = 50  # 50 coordinates along x-axis
+y = 20  # 20 coordinates along y-axis
+n = 5  # Possible classifications [0, 1, 2, 3, 4]
+
+class_map = c(t, x, y, n)  # 3D array of classifications (t, y, x)
+
+#%%
+# Next, we shall import :func:`mcalf.visualisation.bar`.
+
+from mcalf.visualisation import bar
+
+#%%
+# We can now simply plot the 3D array.
+# By default, the first dimension of a 3D array will be averaged to
+# produce a time average, selecting the most common classification
+# at each (x, y) coordinate.
+# This means the percentage abundances will correspond to the
+# most common classification at each coordinate.
+
+bar(class_map)
+
+#%%
+# Instead, the percentage abundances can be determined for the whole
+# 3D array of classifications by setting ``reduce=True``.
+# This skips the averaging process.
+
+bar(class_map, reduce=False)
+
+#%%
+# Alternatively, a 2D array can be passed to the function.
+# If a 2D array is passed, no averaging is needed, and
+# the ``reduce`` parameter is ignored.
+
+#%%
+# A narrower range of classifications to be plotted can be
+# requested with the ``vmin`` and ``vmax`` parameters.
+# To show bars for only classifcations 1, 2, and 3,
+
+bar(class_map, vmin=1, vmax=3)
+
+#%%
+# An alternative set of colours can be requested.
+# Passing a name of a matplotlib colormap to the
+# ``style`` parameter will produce a corresponding
+# list of colours for each of the bars.
+# For advanced use, explore the ``cmap`` parameter.
+
+bar(class_map, style='viridis')
+
+#%%
+# The bar function integrates well with matplotlib, allowing
+# extensive flexibility.
+
+import matplotlib.pyplot as plt
+
+fig, ax = plt.subplots(1, 2, constrained_layout=True)
+
+bar(class_map, vmax=2, style='viridis', ax=ax[0])
+bar(class_map, vmin=3, style='cividis', ax=ax[1])
+
+ax[0].set_title('first 3')
+ax[1].set_title('last 2')
+
+ax[1].set_ylabel('')
+
+plt.show()
+
+#%%
+# Note that ``vmin`` and ``vmax`` are applied before the
+# ``reduce`` value is applied. So setting these ranges
+# can change the calculated abundances for other
+# classifications if ``class_map`` is 3D and
+# ``reduce=True``.
+#
+# The bars do not add up to 100% as a bar for
+# negative, invalid classifications (and therefore
+# classifications out of the ``vmin`` and ``vmax``
+# range) is not shown.
diff --git a/examples/gallery/visualisation/plot_multiple_classifications.py b/examples/gallery/visualisation/plot_multiple_classifications.py
new file mode 100644
index 0000000..10b16b6
--- /dev/null
+++ b/examples/gallery/visualisation/plot_multiple_classifications.py
@@ -0,0 +1,83 @@
+"""
+Combine multiple classification plots
+=====================================
+This is an example showing how to use multiple classification
+plotting functions in a single figure.
+"""
+
+#%%
+# First we shall create a random 3D grid of classifications that can be plotted.
+# Usually you would use a method such as
+# :meth:`mcalf.models.ModelBase.classify_spectra`
+# to classify an array of spectra.
+
+from mcalf.tests.helpers import class_map as c
+
+t = 1  # One image
+x = 20  # 20 coordinates along x-axis
+y = 20  # 20 coordinates along y-axis
+n = 3  # Possible classifications [0, 1, 2]
+
+class_map = c(t, x, y, n)  # 3D array of classifications (t, y, x)
+
+#%%
+# Next we shall create a random array of spectra each labelled
+# with a random classifications.
+# Usually you would provide your own set of hand labelled spectra
+# taken from spectral imaging observations of the Sun.
+# Or you could provide a set of spectra labelled by the classifier.
+
+from mcalf.tests.visualisation.test_classifications import spectra as s
+
+n = 400  # 200 spectra
+w = 20  # 20 wavelength points for each spectrum
+low, high = 0, 3  # Possible classifications [0, 1, 2]
+
+# 2D array of spectra (n, w), 1D array of labels (n,)
+spectra, labels = s(n, w, low, high)
+
+#%%
+# If a GridSpec returned by the plot_classification function has
+# free space, a new axes can be added to the returned GridSpec.
+# We can then request plot_class_map to plot onto this
+# new axes.
+# The colorbar axes can be set to ``fig.axes`` such that
+# the colorbar takes the full height of the figure, as
+# in this case, its colours are the same as the line plots.
+
+import matplotlib.pyplot as plt
+from mcalf.visualisation import plot_classifications, plot_class_map
+
+fig = plt.figure(constrained_layout=True)
+
+gs = plot_classifications(spectra, labels, nrows=2, show_labels=False)
+
+ax = fig.add_subplot(gs[-1])
+plot_class_map(class_map, ax=ax, colorbar_settings={
+    'ax': fig.axes,
+    'label': 'classification',
+})
+
+#%%
+# The function :func:`mcalf.visualisation.init_class_data`` is
+# intended to be an internal function for generating data that
+# is common to multiple plotting functions. However, it may be
+# used externally if necessary.
+
+from mcalf.visualisation import init_class_data, bar
+
+fig, ax = plt.subplots(1, 2, constrained_layout=True)
+
+data = init_class_data(class_map, resolution=(0.25, 0.25), ax=ax[1])
+
+bar(data=data, ax=ax[0])
+plot_class_map(data=data, ax=ax[1], show_colorbar=False)
+
+#%%
+# The following example should be equivalent to the example above,
+
+fig, ax = plt.subplots(1, 2, constrained_layout=True)
+
+bar(class_map, ax=ax[0])
+plot_class_map(class_map, ax=ax[1], show_colorbar=False,
+               resolution=(0.25, 0.25))
diff --git a/examples/gallery/visualisation/plot_plot_class_map.py b/examples/gallery/visualisation/plot_plot_class_map.py
new file mode 100644
index 0000000..db5bcc9
--- /dev/null
+++ b/examples/gallery/visualisation/plot_plot_class_map.py
@@ -0,0 +1,80 @@
+"""
+Plot a map of classifications
+=============================
+This is an example showing how to produce a map showing the spatial
+distribution of spectral classifications in a 2D region of the Sun.
+"""
+
+#%%
+# First we shall create a random 3D grid of classifications that can be plotted.
+# Usually you would use a method such as
+# :meth:`mcalf.models.ModelBase.classify_spectra`
+# to classify an array of spectra.
+
+from mcalf.tests.helpers import class_map as c
+
+t = 3  # Three images
+x = 50  # 50 coordinates along x-axis
+y = 20  # 20 coordinates along y-axis
+n = 5  # Possible classifications [0, 1, 2, 3, 4]
+
+class_map = c(t, x, y, n)  # 3D array of classifications (t, y, x)
+
+#%%
+# Next, we shall import :func:`mcalf.visualisation.plot_class_map`.
+
+from mcalf.visualisation import plot_class_map
+
+#%%
+# We can now simply plot the 3D array.
+# By default, the first dimension of a 3D array will be averaged to
+# produce a time average, selecting the most common classification
+# at each (x, y) coordinate.
+
+plot_class_map(class_map)
+
+#%%
+# A spatial resolution with units can be specified for each axis.
+
+import astropy.units as u
+
+plot_class_map(class_map, resolution=(0.75 * u.km, 1.75 * u.Mm),
+               offset=(-25, -10),
+               dimension=('x distance', 'y distance'))
+
+#%%
+# A narrower range of classifications to be plotted can be
+# requested with the ``vmin`` and ``vmax`` parameters.
+# Classifications outside of the range will appear as grey,
+# the same as pixels with a negative, unassigned classification.
+
+plot_class_map(class_map, vmin=1, vmax=3)
+
+#%%
+# An alternative set of colours can be requested.
+# Passing a name of a matplotlib colormap to the
+# ``style`` parameter will produce a corresponding
+# list of colours for each of the classifications.
+# For advanced use, explore the ``cmap`` parameter.
+
+plot_class_map(class_map, style='viridis')
+
+#%%
+# The plot_class_map function integrates well with
+# matplotlib, allowing extensive flexibility.
+# This example also shows how you can plot a 2D
+# ``class_map`` and skip the averaging.
+
+import matplotlib.pyplot as plt
+
+fig, ax = plt.subplots(2, constrained_layout=True)
+
+plot_class_map(class_map[0], style='viridis', ax=ax[0],
+               show_colorbar=False)
+plot_class_map(class_map[1], style='viridis', ax=ax[1],
+               colorbar_settings={'ax': ax, 'label': 'classification'})
+
+ax[0].set_title('time $t=0$')
+ax[1].set_title('time $t=1$')
+
+plt.show()
diff --git a/examples/gallery/visualisation/plot_plot_classifications.py b/examples/gallery/visualisation/plot_plot_classifications.py
new file mode 100644
index 0000000..71649f2
--- /dev/null
+++ b/examples/gallery/visualisation/plot_plot_classifications.py
@@ -0,0 +1,54 @@
+"""
+Plot a grid of spectra grouped by classification
+================================================
+This is an example showing how to produce a grid of line plots of
+an array of spectra labelled with a classification.
+"""
+
+#%%
+# First we shall create a random array of spectra each labelled
+# with a random classifications.
+# Usually you would provide your own set of hand labelled spectra
+# taken from spectral imaging observations of the Sun.
+
+from mcalf.tests.visualisation.test_classifications import spectra as s
+
+n = 200  # 200 spectra
+w = 20  # 20 wavelength points for each spectrum
+low, high = 1, 5  # Possible classifications [1, 2, 3, 4]
+
+# 2D array of spectra (n, w), 1D array of labels (n,)
+spectra, labels = s(n, w, low, high)
+
+#%%
+# Next, we shall import :func:`mcalf.visualisation.plot_classifications`.
+
+from mcalf.visualisation import plot_classifications
+
+#%%
+# We can now plot a simple grid of the spectra grouped by their
+# classification. By default, a maximum of 20 spectra are plotted
+# for each classification. The first 20 spectra are selected.
+
+plot_classifications(spectra, labels)
+
+#%%
+# A specific number of rows or columns can be requested,
+
+plot_classifications(spectra, labels, ncols=1)
+
+#%%
+# The plot settings can be adjusted,
+
+plot_classifications(spectra, labels, show_labels=False, nlines=5,
+                     style='viridis', plot_settings={'ls': '--', 'marker': 'x'})
+
+#%%
+# By default, the y-axis goes from 0 to 1. This is because
+# labelled training data will typically be rescaled
+# between 0 and 1. However, if a particular classification
+# has spectra that are not within 0 and 1, the y-axis limits
+# are determined by matplotlib.
+
+spectra[labels == 2] += 0.5
+plot_classifications(spectra, labels)
diff --git a/examples/gallery/visualisation/plot_plot_ibis8542.py b/examples/gallery/visualisation/plot_plot_ibis8542.py
new file mode 100644
index 0000000..a060c42
--- /dev/null
+++ b/examples/gallery/visualisation/plot_plot_ibis8542.py
@@ -0,0 +1,129 @@
+"""
+Plot a fitted spectrum
+======================
+This is an example showing how to plot the result of fitting
+a spectrum using the :class:`~mcalf.models.IBIS8542Model` class.
+"""
+
+#%%
+# First we shall create a list of wavelengths, with a variable
+# wavelength spacing.
+# Next, we shall use the Voigt profile to generate spectral
+# intensities at each of the wavelength points.
+# Typically you would provide a spectrum obtained from observations.
+#
+# The data in this example are produced from randomly selected
+# parameters, so numerical values in this example should be ignored.
+
+import numpy as np
+from mcalf.models import IBIS8542Model
+from mcalf.profiles.voigt import double_voigt
+
+# Create the wavelength grid and intensity values
+wavelengths = np.linspace(8541, 8543, 20)
+wavelengths = np.delete(wavelengths, np.s_[1:6:2])
+wavelengths = np.delete(wavelengths, np.s_[-6::2])
+
+spectrum = double_voigt(wavelengths, -526, 8542, 0.1, 0.1,
+                        300, 8541.9, 0.2, 0.05, 1242)
+
+#%%
+
+from mcalf.visualisation import plot_spectrum
+
+plot_spectrum(wavelengths, spectrum, normalised=False)
+
+#%%
+# A basic model is created,
+
+model = IBIS8542Model(original_wavelengths=wavelengths)
+
+#%%
+# The spectrum is provided to the model's fit method.
+# A classifier has not been loaded so the classification must be
+# provided manually.
+# The fitting algorithm assumes that the intensity at the ends of the
+# spectrum is zero, so in this case we need to provide it with a
+# background value to subtract from the spectrum before fitting.
+fit = model.fit(spectrum=spectrum, classifications=4, background=1242)
+
+print(fit)
+
+#%%
+# The spectrum can now be plotted,
+
+# sphinx_gallery_thumbnail_number = 2
+model.plot(fit, spectrum=spectrum, background=1242)
+
+#%%
+# If an array of spectra and associated background values
+# had been loaded into the model with the
+# :meth:`~mcalf.models.ModelBase.load_array` and
+# :meth:`~mcalf.models.ModelBase.load_background`
+# methods respectively, the ``spectrum`` and ``background``
+# parameters would not have to be specified when plotting.
+# This is because the ``fit`` object would contain
+# indices that the ``model`` object would use to look
+# up the original loaded values.
+#
+# Equivalent to above, the plot method can be called on
+# the fit object directly. Remember to specify the ``model``
+# which is needed for additional information such as the
+# stationary line core value.
+
+fit.plot(model, spectrum=spectrum, background=1242)
+
+#%%
+# If the fit has multiple spectral components, such as
+# an active emission profile mixed with a quiescent
+# absorption profile, the follow method can be used
+# to plot the components separatly.
+#
+# If the fit only has a single component the ``plot``
+# method as shown above is used.
+
+model.plot_separate(fit, spectrum=spectrum, background=1242)
+
+#%%
+# If the fit has an emission component, it is subtracted
+# from the raw spectral data. Otherwise, the default
+# ``plot`` method is used.
+
+model.plot_subtraction(fit, spectrum=spectrum, background=1242)
+
+#%%
+# The same line on multiple plots is only labelled the first time it
+# plotted in the figure. This prevents duplicated entries in
+# legends.
+
+import matplotlib.pyplot as plt
+
+fig, ax = plt.subplots(1, 2)
+
+model.plot(fit, spectrum=spectrum, background=1242,
+           show_legend=False, ax=ax[0])
+model.plot_separate(fit, spectrum=spectrum, background=1242,
+                    show_legend=False, ax=ax[1])
+
+fig.subplots_adjust(top=0.75)  # Create space above for legend
+fig.legend(ncol=2, loc='upper center', bbox_to_anchor=(0.5, 0.97))
+
+plt.show()
+
+#%%
+# The underlying :func:`mcalf.visualisation.plot_ibis8542`
+# function can be used directly. However, it is recommended
+# to plot using the method detailed above as it will do
+# additional processing to the wavelengths and spectrum
+# and also pass additional parameters, such as sigma,
+# to this fitting function.
+
+from mcalf.visualisation import plot_ibis8542
+
+plot_ibis8542(wavelengths, spectrum, fit.parameters, 1242)
+
+#%%
+# The y-axis and legend can be easily hidden,
+
+model.plot(fit, spectrum=spectrum, background=1242,
+           show_intensity=False, show_legend=False)
diff --git a/examples/gallery/visualisation/plot_plot_map.py b/examples/gallery/visualisation/plot_plot_map.py
new file mode 100644
index 0000000..7b3b5fb
--- /dev/null
+++ b/examples/gallery/visualisation/plot_plot_map.py
@@ -0,0 +1,114 @@
+"""
+Plot a map of velocities
+========================
+This is an example showing how to produce a map showing the spatial
+distribution of velocities in a 2D region of the Sun.
+"""
+
+#%%
+# First we shall create a random 2D grid of velocities that can be plotted.
+# Usually you would use a method such as
+# :meth:`mcalf.models.FitResults.velocities`
+# to extract an array of velocities from fitted spectra.
+
+import numpy as np
+np.random.seed(0)
+
+x = 50  # 50 coordinates along x-axis
+y = 40  # 40 coordinates along y-axis
+low, high = -10, 10  # range of velocities (km/s)
+
+
+def a(x, y, low, high):
+    arr = np.random.normal(0, (high - low) / 2 * 0.3, (y, x))
+    arr[arr < low] = low
+    arr[arr > high] = high
+    i = np.random.randint(0, arr.size, arr.size // 100)
+    arr[np.unravel_index(i, arr.shape)] = np.nan
+    return arr
+
+
+arr = a(x, y, low, high)  # 2D array of velocities (y, x)
+
+#%%
+# Next, we shall import :func:`mcalf.visualisation.plot_map`.
+
+from mcalf.visualisation import plot_map
+
+#%%
+# We can now simply plot the 2D array.
+
+plot_map(arr)
+
+#%%
+# Notice that pixels with missing data (NaN) are shown in grey.
+
+#%%
+# By default, the velocity data are assumed to have units km/s.
+# If your data are not in km/s, you must either 1) rescale the
+# array such that it is in km/s, 2) attach an astropy unit
+# to the array to override the default, or 3) pass an
+# astropy unit to the ``unit`` parameter to override the
+# default. For example, we can change from km/s to m/s,
+
+import astropy.units as u
+
+plot_map(arr * 1000 * u.m / u.s)
+
+#%%
+# A spatial resolution with units can be specified for each axis.
+
+plot_map(arr, resolution=(0.5 * u.km, 0.6 * u.Mm), offset=(-25, -20))
+
+#%%
+# A narrower range of velocities to be plotted can be
+# requested with the ``vmin`` and ``vmax`` parameters.
+# Classifications outside of the range will appear saturated.
+# Providing only one of ``vmin`` and ``vmax`` with set the
+# other such that zero is the midpoint.
+
+plot_map(arr, vmax=4)
+
+#%%
+# A mask can be applied to the velocity array to isolate
+# a region of interest. This functionally is useful if,
+# for example, data only exist for a circular region
+# and you want to distinguish between the pixels that are
+# out of bounds and the data that were not successfully
+# fitted.
+
+from mcalf.utils.mask import genmask
+mask = genmask(50, 40, 18)
+
+plot_map(arr, ~mask)
+
+#%%
+# Notice how data out of bounds are grey, while data
+# which were not fitted successfully are now black.
+
+#%%
+# A region of interest, typically the umbra of a
+# sunspot, can be outlined by passing a different
+# mask.
+
+umbra_mask = genmask(50, 40, 5, 5, 5)
+
+plot_map(arr, ~mask, umbra_mask)
+
+#%%
+# The plot_map function integrates well with
+# matplotlib, allowing extensive flexibility.
+
+import matplotlib.pyplot as plt
+
+fig, ax = plt.subplots(1, 2, constrained_layout=True)
+
+plot_map(arr[:, :25], vmax=6, ax=ax[0], show_colorbar=False)
+im = plot_map(arr[:, 25:], vmax=6, ax=ax[1], show_colorbar=False)
+
+fig.colorbar(im, ax=[ax], location='bottom', label='velocity (km/s)')
+
+ax[0].set_title('first half')
+ax[1].set_title('second half')
+
+plt.show()
diff --git a/examples/gallery/visualisation/plot_plot_spectrum.py b/examples/gallery/visualisation/plot_plot_spectrum.py
new file mode 100644
index 0000000..5cc98e1
--- /dev/null
+++ b/examples/gallery/visualisation/plot_plot_spectrum.py
@@ -0,0 +1,47 @@
+"""
+Plot a spectrum
+===============
+This is an example showing how to plot a spectrum with the
+:func:`mcalf.visualisation.plot_spectrum` function.
+"""
+
+#%%
+# First we shall create a list of wavelengths, with a variable
+# wavelength spacing.
+# Next, we shall use the Voigt profile to generate spectral
+# intensities at each of the wavelength points.
+# Typically you would provide a spectrum obtained from observations.
+
+import numpy as np
+wavelengths = np.linspace(8541, 8543, 20)
+wavelengths = np.delete(wavelengths, np.s_[1:6:2])
+wavelengths = np.delete(wavelengths, np.s_[-6::2])
+
+from mcalf.profiles.voigt import voigt
+spectrum = voigt(wavelengths, -526, 8542, 0.1, 0.1, 1242)
+
+#%%
+# Next, we shall import :func:`mcalf.visualisation.plot_spectrum`.
+
+from mcalf.visualisation import plot_spectrum
+
+#%%
+# We can now simply plot the spectrum.
+
+plot_spectrum(wavelengths, spectrum)
+
+#%%
+# Notice how the spectrum above is normalised.
+# The normalisation is applied by dividing through
+# by the mean of the three rightmost points.
+# To plot the raw spectrum,
+
+plot_spectrum(wavelengths, spectrum, normalised=False)
+
+#%%
+# The line connecting the points provided in the ``spectrum``
+# array above is smooth. This is due to spline interpolation
+# being applied. Interpolation can be disabled, resulting
+# in a straight line between each of the points.
+
+plot_spectrum(wavelengths, spectrum, smooth=False)
diff --git a/setup.cfg b/setup.cfg
index fe6fe78..a6a3e62 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -51,6 +51,8 @@ docs =
   ipykernel
   nbsphinx==0.7.1
   sphinx-rtd-theme==0.5.0
+  sphinx-gallery
+  pytest
 
 [options.packages.find]
 where = src
diff --git a/src/mcalf/models/ibis.py b/src/mcalf/models/ibis.py
index e804ecf..2b598e8 100644
--- a/src/mcalf/models/ibis.py
+++ b/src/mcalf/models/ibis.py
@@ -356,6 +356,10 @@ def plot(self, fit=None, time=None, row=None, column=None, spectrum=None, classi
         plot_separate : Plot the fit parameters separately.
         plot_subtraction : Plot the spectrum with the emission fit subtracted from it.
         mcalf.models.FitResult.plot : Plotting method on the fit result.
+
+        Examples
+        --------
+        .. minigallery:: mcalf.visualisation.plot_ibis8542
         """
         if fit.__class__ == FitResult:  # If fit is a `FitResult`
 
diff --git a/src/mcalf/profiles/voigt.py b/src/mcalf/profiles/voigt.py
index 3b2ceb6..3cbe789 100644
--- a/src/mcalf/profiles/voigt.py
+++ b/src/mcalf/profiles/voigt.py
@@ -20,6 +20,11 @@
 except IndexError:  # File does not exist
     warnings.warn("Could not locate the external C library. Further use of `clib` will fail!")
 
+###
+# readthedocs.org does not support clib (force clib=False)
+import os
+not_on_rtd = os.environ.get('READTHEDOCS') != 'True'
+###
 
 # Parameters for `voigt_approx_nobg` and other approx. Voigt functions
 params = np.array([[-1.2150, -1.3509, -1.2150, -1.3509],
@@ -109,7 +114,7 @@ def voigt_nobg(x, a, b, s, g, clib=True):
     """
     warnings.filterwarnings("ignore", category=IntegrationWarning)
     u = x - b
-    if clib:
+    if clib and not_on_rtd:
         i = [quad(cvoigt, -np.inf, np.inf, args=(v, s, g), epsabs=1.49e-1, epsrel=1.49e-4)[0] for v in u]
     else:
         i = quad_vec(lambda y: np.exp(-y**2 / (2 * s**2)) / (g**2 + (u - y)**2), -np.inf, np.inf)[0]
diff --git a/src/mcalf/tests/baseline_mpl_334_ft_261/test_plot_spectrum.png b/src/mcalf/tests/baseline_mpl_334_ft_261/test_plot_spectrum.png
index 2006d9d..833f422 100644
Binary files a/src/mcalf/tests/baseline_mpl_334_ft_261/test_plot_spectrum.png and b/src/mcalf/tests/baseline_mpl_334_ft_261/test_plot_spectrum.png differ
diff --git a/src/mcalf/tests/baseline_mpl_334_ft_261/test_plot_spectrum_no_norm.png b/src/mcalf/tests/baseline_mpl_334_ft_261/test_plot_spectrum_no_norm.png
index 12485aa..b8550bd 100644
Binary files a/src/mcalf/tests/baseline_mpl_334_ft_261/test_plot_spectrum_no_norm.png and b/src/mcalf/tests/baseline_mpl_334_ft_261/test_plot_spectrum_no_norm.png differ
diff --git a/src/mcalf/tests/baseline_mpl_334_ft_261/test_plot_spectrum_no_smooth.png b/src/mcalf/tests/baseline_mpl_334_ft_261/test_plot_spectrum_no_smooth.png
index 4f843cc..0edea4a 100644
Binary files a/src/mcalf/tests/baseline_mpl_334_ft_261/test_plot_spectrum_no_smooth.png and b/src/mcalf/tests/baseline_mpl_334_ft_261/test_plot_spectrum_no_smooth.png differ
diff --git a/src/mcalf/tests/baseline_mpl_334_ft_261/test_plot_spectrum_no_smooth_norm.png b/src/mcalf/tests/baseline_mpl_334_ft_261/test_plot_spectrum_no_smooth_norm.png
index ee1facc..6ea3976 100644
Binary files a/src/mcalf/tests/baseline_mpl_334_ft_261/test_plot_spectrum_no_smooth_norm.png and b/src/mcalf/tests/baseline_mpl_334_ft_261/test_plot_spectrum_no_smooth_norm.png differ
diff --git a/src/mcalf/tests/figure_hashes_mpl_334_ft_261.json b/src/mcalf/tests/figure_hashes_mpl_334_ft_261.json
index ba91416..0536fe8 100644
--- a/src/mcalf/tests/figure_hashes_mpl_334_ft_261.json
+++ b/src/mcalf/tests/figure_hashes_mpl_334_ft_261.json
@@ -45,10 +45,10 @@
   "mcalf.tests.visualisation.test_classifications.test_plot_class_map_units": "3737733a86ca318cfa6bde65cf922934f269b72cba90c75361971221000fbd27",
   "mcalf.tests.visualisation.test_classifications.test_plot_class_map_allbad": "944e89a497b2e481fe6f61552964d1c183cd39e778113ac366a3ff55dce107a5",
   "mcalf.tests.visualisation.test_classifications.test_plot_bar_and_class_map": "af1291cf4c4a04efd98fc0a69fe92ba422742e332d8018fbe5ce5b18a748d6c2",
-  "mcalf.tests.visualisation.test_spec.test_plot_spectrum": "79afa94c229e75646304b14bb388df30f23d6fd06ac82cfa4443ab05fc4c4a4a",
-  "mcalf.tests.visualisation.test_spec.test_plot_spectrum_no_smooth": "c4bf62988715825cecaf75e03b66c05231854e5b84962fbd03f11e1c381e1273",
-  "mcalf.tests.visualisation.test_spec.test_plot_spectrum_no_norm": "f4a52acdbb2a614e0a775ff7833a95151bf6188126286d1eb87fedef7f324711",
-  "mcalf.tests.visualisation.test_spec.test_plot_spectrum_no_smooth_norm": "96ce04857fe6fdc479c1620e72a0d26c981dd82082da2627e92e9716c46b2182",
+  "mcalf.tests.visualisation.test_spec.test_plot_spectrum": "20633cadf2b0e975c174053eba72e379ee09c376c314ec69ed0f1ce6c9f74eeb",
+  "mcalf.tests.visualisation.test_spec.test_plot_spectrum_no_smooth": "e09508418236ce5a9906b5c9c63716c8c11075c8e38326fc3d271a54cb3b867d",
+  "mcalf.tests.visualisation.test_spec.test_plot_spectrum_no_norm": "be773c5aa34c160efa5caa12b9c0d46a6d5c58c859263f5fa3456affd1feb2df",
+  "mcalf.tests.visualisation.test_spec.test_plot_spectrum_no_smooth_norm": "317f5ba87d38eea86581f928100427546b2fddc0046318f896b3342ba6c9792e",
   "mcalf.tests.visualisation.test_velocity.test_plot_map_basic": "7aafc54666ffd269480caecd10e0606fd0f59b1f43ed180c230421671df8bd63",
   "mcalf.tests.visualisation.test_velocity.test_plot_map_mask": "a8d3af74f08a110ea6404aa782d5cf09820524b0cd9ec43ea22366757a922037",
   "mcalf.tests.visualisation.test_velocity.test_plot_map_umbra_mask": "c8e305a8ab5da653387a585ea21a74d1d53020bb7d43657aecf4b47429b7b813",
diff --git a/src/mcalf/utils/mask.py b/src/mcalf/utils/mask.py
index d5278fa..c10b9f3 100644
--- a/src/mcalf/utils/mask.py
+++ b/src/mcalf/utils/mask.py
@@ -24,6 +24,10 @@ def genmask(width, height, radius=np.inf, right_shift=0, up_shift=0):
     -------
     array : numpy.ndarray, shape=(height, width)
         The generated mask.
+
+    Examples
+    --------
+    .. minigallery:: mcalf.utils.mask.genmask
     """
     array = radial_distances(width, height) < radius  # Create mask
     array = np.roll(array, [up_shift, right_shift], [0, 1])  # Apply shift
diff --git a/src/mcalf/visualisation/classifications.py b/src/mcalf/visualisation/classifications.py
index 6c435f3..37845e6 100644
--- a/src/mcalf/visualisation/classifications.py
+++ b/src/mcalf/visualisation/classifications.py
@@ -55,6 +55,10 @@ def plot_classifications(spectra, labels, nrows=None, ncols=None, nlines=20, sty
     -------
     gs : matplotlib.gridspec.GridSpec
         The grid layout subplots are placed on within the figure.
+
+    Examples
+    --------
+    .. minigallery:: mcalf.visualisation.plot_classifications
     """
     if fig is None:
         fig = plt.gcf()
@@ -202,6 +206,10 @@ def bar(class_map=None, vmin=None, vmax=None, reduce=True, style='original', cma
     -----
     Visualisation assumes that all integers between `vmin` and `vmax` are valid
     classifications, even if they do not appear in `class_map`.
+
+    Examples
+    --------
+    .. minigallery:: mcalf.visualisation.bar
     """
     if ax is None:
         ax = plt.gca()
@@ -295,6 +303,10 @@ def plot_class_map(class_map=None, vmin=None, vmax=None, resolution=None, offset
     -----
     Visualisation assumes that all integers between `vmin` and `vmax` are valid
     classifications, even if they do not appear in `class_map`.
+
+    Examples
+    --------
+    .. minigallery:: mcalf.visualisation.plot_class_map
     """
     if ax is None:
         ax = plt.gca()
@@ -384,6 +396,10 @@ def init_class_data(class_map, vmin=None, vmax=None, reduce=True, resolution=Non
     mcalf.utils.smooth.mask_classifications : Mask 2D and 3D arrays of classifications.
     mcalf.utils.plot.calculate_extent : Calculate the extent from a particular data shape and resolution.
     mcalf.utils.plot.class_cmap : Create a listed colormap for a specific number of classifications.
+
+    Examples
+    --------
+    .. minigallery:: mcalf.visualisation.init_class_data
     """
     # Mask and average classification map according to the classification range and shape
     class_map, vmin, vmax = mask_classifications(class_map, vmin, vmax, reduce=reduce)
diff --git a/src/mcalf/visualisation/spec.py b/src/mcalf/visualisation/spec.py
index 1d4c4ef..7c7ce8d 100644
--- a/src/mcalf/visualisation/spec.py
+++ b/src/mcalf/visualisation/spec.py
@@ -61,6 +61,10 @@ class instead.
     mcalf.models.IBIS8542Model.plot_separate : Plot the fit parameters separately.
     mcalf.models.IBIS8542Model.plot_subtraction : Plot the spectrum with the emission fit subtracted from it.
     mcalf.models.FitResult.plot : Plotting method provided by the fit result.
+
+    Examples
+    --------
+    .. minigallery:: mcalf.visualisation.plot_ibis8542
     """
     if ax is None:
         ax = plt.gca()
@@ -176,6 +180,10 @@ def plot_spectrum(wavelengths, spectrum, normalised=True, smooth=True, ax=None):
     -------
     ax : matplotlib.axes.Axes
         Axes the lines are drawn on.
+
+    Examples
+    --------
+    .. minigallery:: mcalf.visualisation.plot_spectrum
     """
     # Inputted wavelengths and spectrum
     wavelengths_pts = wavelengths.copy()
@@ -208,5 +216,7 @@ def plot_spectrum(wavelengths, spectrum, normalised=True, smooth=True, ax=None):
     ylabel = 'normalised intensity ($I/I_c$)' if normalised else 'intensity ($I$)'
     ax.set_ylabel(ylabel)
     ax.minorticks_on()
+    for label in ax.xaxis.get_ticklabels()[::2]:
+        label.set_visible(False)
 
     return ax
diff --git a/src/mcalf/visualisation/velocity.py b/src/mcalf/visualisation/velocity.py
index 1c8c582..569efdf 100644
--- a/src/mcalf/visualisation/velocity.py
+++ b/src/mcalf/visualisation/velocity.py
@@ -61,6 +61,10 @@ def plot_map(arr, mask=None, umbra_mask=None, resolution=None, offset=(0, 0), vm
     See Also
     --------
     mcalf.models.FitResults.velocities : Calculate the Doppler velocities for an array of fits.
+
+    Examples
+    --------
+    .. minigallery:: mcalf.visualisation.plot_map
     """
     if ax is None:
         ax = plt.gca()