From ee44095c1516f1d8ce5a6b52a24c3d4df91380ab Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 13:24:55 -0300
Subject: [PATCH 01/46] Draft function for generating pyvista objects out of
 prisms

Create a new visualization module to host it.
---
 harmonica/visualization/__init__.py |   0
 harmonica/visualization/prism.py    | 180 ++++++++++++++++++++++++++++
 2 files changed, 180 insertions(+)
 create mode 100644 harmonica/visualization/__init__.py
 create mode 100644 harmonica/visualization/prism.py

diff --git a/harmonica/visualization/__init__.py b/harmonica/visualization/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/harmonica/visualization/prism.py b/harmonica/visualization/prism.py
new file mode 100644
index 000000000..5c0f9f179
--- /dev/null
+++ b/harmonica/visualization/prism.py
@@ -0,0 +1,180 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Functions for visualizing prisms through pyvista
+"""
+import numpy as np
+
+try:
+    import vtk
+except ImportError:
+    vtk = None
+
+try:
+    import pyvista
+except ImportError:
+    pyvista = None
+
+
+def prisms_to_pyvista(prisms, properties=None):
+    """
+    Create a ``pyvista.UnstructuredGrid`` from a set of prisms
+
+    Parameters
+    ----------
+    prisms : list or 2d-array
+        List or 2d-array with the boundaries of the prisms.
+        Each row contains the boundaries of each prism in the following order:
+        ``west``, ``east``, ``south``, ``north``, ``bottom``, ``top``.
+    properties : dict or None (optional)
+        Dictionary with the physical properties of the prisms.
+        Each key should be a string and their value an array.
+        If None, no property will be added to the
+        :class:`pyvista.UnstructuredGrid`.
+        Default to None.
+
+    Returns
+    -------
+    pv_grid : :class:`pyvista.UnstructuredGrid`
+        :class:`pyvista.UnstructuredGrid` that represents the prisms with their
+        properties (if any).
+
+    Examples
+    --------
+
+    >>> pv_grid = prisms_to_pyvista(
+    ...     [0, 1e3, -2e3, 2e3, -5e3, -3e3],
+    ...     properties={"density": 2670},
+    ... )
+    >>> pv_grid.n_cells
+    1
+    >>> pv_grid.n_points
+    8
+    >>> pv_grid.get_array("density")
+    array([2670.])
+    >>> pv_grid.cell_bounds(0)
+    [0.0, 1000.0, -2000.0, 2000.0, -5000.0, -3000.0]
+
+
+    """
+    # Check if vkt and pyvista are installed
+    msg = "Missing optional dependency '{}' required for generating pyvista grids."
+    if vtk is None:
+        raise ValueError(msg.format("vtk"))
+    if pyvista is None:
+        raise ValueError(msg.format("pyvista"))
+    # Get prisms and number of prisms
+    prisms = np.atleast_2d(prisms)
+    n_prisms = prisms.shape[0]
+    # Get the vertices of the prisms
+    vertices = np.vstack([_get_prism_vertices(prism) for prism in prisms])
+    # Get the cells for these prisms
+    cells = _build_cells(n_prisms)
+    # Define an array with the types of the cells
+    cell_type = np.full(n_prisms, vtk.VTK_HEXAHEDRON)
+    # Build the UnstructuredGrid
+    pv_grid = pyvista.UnstructuredGrid(cells, cell_type, vertices)
+    # Add properties to the grid
+    if properties is not None:
+        for name, prop in properties.items():
+            # Convert the property to an array and ravel it
+            prop = np.atleast_1d(prop).ravel()
+            # Assign the property to the cell_data
+            pv_grid.cell_data[name] = prop
+    return pv_grid
+
+
+def _build_cells(n_prisms):
+    """
+    Build the VTK cells for a given number of prisms
+
+    We will represent each prism a single cell of ``VTK_HEXAHEDRON`` type.
+    Each cell is an array with 9 elements:
+      - the first one indicates the number of vertices that the cell will have
+        (8 for prisms), and
+      - the following 8 elements are the indices of the vertices that form the
+        prism.
+
+    The cells is a concatenation of every 9 elements cell array, resulting in
+    a 1d-array of 9 times ``n_prisms`` elements.
+
+    Parameters
+    ----------
+    n_prisms : int
+        Number of prisms in the set
+
+    Returns
+    -------
+    cells : 1d-array
+        Array representing the cells for each one of the prisms.
+
+    Examples
+    --------
+
+    >>> _build_cells(4)
+    array([ 8,  0,  1,  2,  3,  4,  5,  6,  7,  8,  8,  9, 10, 11, 12, 13, 14,
+           15,  8, 16, 17, 18, 19, 20, 21, 22, 23,  8, 24, 25, 26, 27, 28, 29,
+           30, 31])
+    """
+    # Get total number of vertices
+    n_vertices = n_prisms * 8
+    # Build the indices array as a 2D array: each row contains the indices for
+    # each prism
+    indices = np.arange(n_vertices).reshape(n_prisms, 8)
+    # Create a vertical array full of 8
+    n_vertices_per_cell = np.full(fill_value=8, shape=(n_prisms, 1))
+    # Stack the indices and n_vertices_per_cell arrays
+    cells = np.hstack((n_vertices_per_cell, indices)).ravel()
+    return cells
+
+
+def _get_prism_vertices(prism):
+    """
+    Return the vertices of the given prism
+
+    Parameters
+    ----------
+    prism : list or 1d-array
+        List or 1d-array with the boundaries of a single prism in the following
+        order: ``west``, ``east``, ``south``, ``north``, ``bottom``, ``top``.
+
+    Returns
+    -------
+    vertices : 2d-array
+        2D array with the coordinates of the vertices of the prism. Each row of
+        the array corresponds to the coordinate of a single vertex in the
+        following order: ``easting``, ``northing``, ``upward``.
+        The order of the vertices is fixed to be compatible with VTK.
+
+    Examples
+    --------
+
+    >>> _get_prism_vertices([-1, 1, -2, 2, -3, 3])
+    array([[-1, -2, -3],
+           [ 1, -2, -3],
+           [ 1,  2, -3],
+           [-1,  2, -3],
+           [-1, -2,  3],
+           [ 1, -2,  3],
+           [ 1,  2,  3],
+           [-1,  2,  3]])
+
+    """
+    w, e, s, n, bottom, top = prism[:]
+    vertices = np.array(
+        [
+            [w, s, bottom],
+            [e, s, bottom],
+            [e, n, bottom],
+            [w, n, bottom],
+            [w, s, top],
+            [e, s, top],
+            [e, n, top],
+            [w, n, top],
+        ]
+    )
+    return vertices

From b1cca3521344fbd152bd52a64b7f547dff12cd91 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 13:26:11 -0300
Subject: [PATCH 02/46] Add a to_pyvista method to PrismLayer

The method returns a pyvista.Unstructured object that represents the
layer.
---
 harmonica/forward/prism_layer.py | 19 +++++++++++++++++++
 1 file changed, 19 insertions(+)

diff --git a/harmonica/forward/prism_layer.py b/harmonica/forward/prism_layer.py
index 8063f67ef..ed3cc67ec 100644
--- a/harmonica/forward/prism_layer.py
+++ b/harmonica/forward/prism_layer.py
@@ -14,6 +14,7 @@
 import xarray as xr
 
 from .prism import prism_gravity
+from ..visualization.prism import prisms_to_pyvista
 
 
 def prism_layer(
@@ -446,3 +447,21 @@ def get_prism(self, indices):
         bottom = self._obj.bottom.values[indices]
         top = self._obj.top.values[indices]
         return west, east, south, north, bottom, top
+
+    def to_pyvista(self):
+        """
+        Return a pyvista UnstructuredGrid to plot the PrismLayer
+
+        Returns
+        -------
+        pv_grid : :class:`pyvista.UnstructuredGrid`
+            :class:`pyvista.UnstructuredGrid` containing each prism of the
+            layer as a hexahedron along with their properties.
+        """
+        prisms = self._to_prisms()
+        properties = None
+        if self._obj.data_vars:
+            properties = {
+                data_var: self._obj[data_var] for data_var in self._obj.data_vars
+            }
+        return prisms_to_pyvista(prisms, properties=properties)

From b92a2170ae74613be0dc3daf06e136878bf273b4 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 13:27:11 -0300
Subject: [PATCH 03/46] Add vtk and pyvista to environment.yml as optional
 dependencies

---
 environment.yml | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/environment.yml b/environment.yml
index d8f8104e0..5b32fdf6f 100644
--- a/environment.yml
+++ b/environment.yml
@@ -14,6 +14,9 @@ dependencies:
   - pooch>=0.7.0
   - verde>=1.5.0
   - xarray
+  # Optional requirements
+  - pyvista
+  - vtk>9
   # Development requirements
   - boule
   - pyproj

From 3833a77a785b2bc384437ddc502f989366fb68ab Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 13:27:34 -0300
Subject: [PATCH 04/46] Add link to pyvista docs in doc/conf.py

---
 doc/conf.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/doc/conf.py b/doc/conf.py
index a0c05deeb..4bd02692e 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -49,6 +49,7 @@
     "verde": ("https://www.fatiando.org/verde/latest/", None),
     "boule": ("https://www.fatiando.org/boule/latest/", None),
     "matplotlib": ("https://matplotlib.org/", None),
+    "pyvista": ("https://docs.pyvista.org", None),
 }
 
 # Autosummary pages will be generated by sphinx-autogen instead of sphinx-build

From 6ae56c4b7d599058442d34ccbf7335be9a21d7e5 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 13:27:47 -0300
Subject: [PATCH 05/46] Add a new example that plots a prism layer in 3D

---
 examples/visualization/prism_layer.py | 45 +++++++++++++++++++++++++++
 1 file changed, 45 insertions(+)
 create mode 100644 examples/visualization/prism_layer.py

diff --git a/examples/visualization/prism_layer.py b/examples/visualization/prism_layer.py
new file mode 100644
index 000000000..c314f7fc5
--- /dev/null
+++ b/examples/visualization/prism_layer.py
@@ -0,0 +1,45 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+3D plot of layer of prisms
+==========================
+
+"""
+import pyproj
+import verde as vd
+import harmonica as hm
+import pyvista as pv
+
+# Read South Africa topography
+south_africa_topo = hm.datasets.fetch_south_africa_topography()
+
+# Project the grid
+projection = pyproj.Proj(proj="merc", lat_ts=south_africa_topo.latitude.values.mean())
+south_africa_topo = vd.project_grid(south_africa_topo.topography, projection=projection)
+
+# Create a 2d array with the density of the prisms Points above the geoid will
+# have a density of 2670 kg/m^3 Points below the geoid will have a density
+# contrast equal to the difference between the density of the ocean and the
+# density of the upper crust: # 1000 kg/m^3 - 2900 kg/m^3
+density = south_africa_topo.copy()  # copy topography to a new xr.DataArray
+density.values[:] = 2670.0  # replace every value for the density of the topography
+# Change density values of ocean points
+density = density.where(south_africa_topo >= 0, 1000 - 2900)
+
+# Create layer of prisms
+prisms = hm.prism_layer(
+    (south_africa_topo.easting, south_africa_topo.northing),
+    surface=south_africa_topo,
+    reference=0,
+    properties={"density": density},
+)
+
+# Plot with pyvista
+plotter = pv.Plotter()
+plotter.add_mesh(prisms.prism_layer.to_pyvista())
+plotter.set_scale(zscale=100)  # exaggerate the vertical coordinate
+plotter.show()

From 0415eb1e65fad4a5d88d1762d4b2614009748bbd Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 13:32:26 -0300
Subject: [PATCH 06/46] Simplify the generation of cells through the dictionary
 interface

---
 harmonica/visualization/prism.py | 57 +++-----------------------------
 1 file changed, 5 insertions(+), 52 deletions(-)

diff --git a/harmonica/visualization/prism.py b/harmonica/visualization/prism.py
index 5c0f9f179..3a8710337 100644
--- a/harmonica/visualization/prism.py
+++ b/harmonica/visualization/prism.py
@@ -67,71 +67,24 @@ def prisms_to_pyvista(prisms, properties=None):
         raise ValueError(msg.format("vtk"))
     if pyvista is None:
         raise ValueError(msg.format("pyvista"))
+
     # Get prisms and number of prisms
     prisms = np.atleast_2d(prisms)
     n_prisms = prisms.shape[0]
     # Get the vertices of the prisms
     vertices = np.vstack([_get_prism_vertices(prism) for prism in prisms])
-    # Get the cells for these prisms
-    cells = _build_cells(n_prisms)
-    # Define an array with the types of the cells
-    cell_type = np.full(n_prisms, vtk.VTK_HEXAHEDRON)
+    # Generate the cells for the hexahedrons
+    cells = np.arange(n_prisms * 8).reshape([n_prisms, 8])
     # Build the UnstructuredGrid
-    pv_grid = pyvista.UnstructuredGrid(cells, cell_type, vertices)
+    pv_grid = pyvista.UnstructuredGrid({vtk.VTK_HEXAHEDRON: cells}, vertices)
     # Add properties to the grid
     if properties is not None:
         for name, prop in properties.items():
-            # Convert the property to an array and ravel it
-            prop = np.atleast_1d(prop).ravel()
             # Assign the property to the cell_data
-            pv_grid.cell_data[name] = prop
+            pv_grid.cell_data[name] = np.atleast_1d(prop).ravel()
     return pv_grid
 
 
-def _build_cells(n_prisms):
-    """
-    Build the VTK cells for a given number of prisms
-
-    We will represent each prism a single cell of ``VTK_HEXAHEDRON`` type.
-    Each cell is an array with 9 elements:
-      - the first one indicates the number of vertices that the cell will have
-        (8 for prisms), and
-      - the following 8 elements are the indices of the vertices that form the
-        prism.
-
-    The cells is a concatenation of every 9 elements cell array, resulting in
-    a 1d-array of 9 times ``n_prisms`` elements.
-
-    Parameters
-    ----------
-    n_prisms : int
-        Number of prisms in the set
-
-    Returns
-    -------
-    cells : 1d-array
-        Array representing the cells for each one of the prisms.
-
-    Examples
-    --------
-
-    >>> _build_cells(4)
-    array([ 8,  0,  1,  2,  3,  4,  5,  6,  7,  8,  8,  9, 10, 11, 12, 13, 14,
-           15,  8, 16, 17, 18, 19, 20, 21, 22, 23,  8, 24, 25, 26, 27, 28, 29,
-           30, 31])
-    """
-    # Get total number of vertices
-    n_vertices = n_prisms * 8
-    # Build the indices array as a 2D array: each row contains the indices for
-    # each prism
-    indices = np.arange(n_vertices).reshape(n_prisms, 8)
-    # Create a vertical array full of 8
-    n_vertices_per_cell = np.full(fill_value=8, shape=(n_prisms, 1))
-    # Stack the indices and n_vertices_per_cell arrays
-    cells = np.hstack((n_vertices_per_cell, indices)).ravel()
-    return cells
-
-
 def _get_prism_vertices(prism):
     """
     Return the vertices of the given prism

From 58367523b8caf8c310e62fe4db968e732bcc554d Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 14:09:50 -0300
Subject: [PATCH 07/46] Configure sphinx to show pyvista plots in gallery

Also improve the example
---
 doc/conf.py                           | 14 ++++++++++++++
 examples/visualization/README.txt     |  2 ++
 examples/visualization/prism_layer.py | 10 ++++++++--
 3 files changed, 24 insertions(+), 2 deletions(-)
 create mode 100644 examples/visualization/README.txt

diff --git a/doc/conf.py b/doc/conf.py
index 4bd02692e..eeccd8713 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -5,6 +5,7 @@
 # This code is part of the Fatiando a Terra project (https://www.fatiando.org)
 #
 import datetime
+import pyvista
 
 from sphinx_gallery.sorting import ExampleTitleSortKey
 
@@ -95,9 +96,22 @@
     "doc_module": "harmonica",
     # Insert links to documentation of objects in the examples
     "reference_url": {"harmonica": None},
+    # Add pyvista to the image scrapers
+    "image_scrapers": ("pyvista"),
 }
 
 
+# Pyvista configurations
+# -----------------------------------------------------------------------------
+# necessary when building the sphinx gallery
+pyvista.BUILDING_GALLERY = True
+pyvista.OFF_SCREEN = True
+
+# Optional - set parameters like theme or window size
+pyvista.set_plot_theme("document")
+pyvista.global_theme.window_size = (1024 * 2, 768 * 2)
+
+
 # HTML output configuration
 # -----------------------------------------------------------------------------
 html_title = f'{project} <span class="project-version">{version}</span>'
diff --git a/examples/visualization/README.txt b/examples/visualization/README.txt
new file mode 100644
index 000000000..e68465716
--- /dev/null
+++ b/examples/visualization/README.txt
@@ -0,0 +1,2 @@
+3D visualizations
+-----------------
diff --git a/examples/visualization/prism_layer.py b/examples/visualization/prism_layer.py
index c314f7fc5..bbd18650f 100644
--- a/examples/visualization/prism_layer.py
+++ b/examples/visualization/prism_layer.py
@@ -38,8 +38,14 @@
     properties={"density": density},
 )
 
+# Create a pyvista UnstructuredGrid from the prism layer
+pv_grid = prisms.prism_layer.to_pyvista()
+print(pv_grid)
+
 # Plot with pyvista
 plotter = pv.Plotter()
-plotter.add_mesh(prisms.prism_layer.to_pyvista())
-plotter.set_scale(zscale=100)  # exaggerate the vertical coordinate
+plotter.add_mesh(pv_grid, scalars="density")
+plotter.set_scale(zscale=75)  # exaggerate the vertical coordinate
+plotter.camera_position = "xz"
+plotter.camera.elevation = 15
 plotter.show()

From 6ed8fc6556c5d73e548d7b1486e8eff7dfbff8c9 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 14:23:36 -0300
Subject: [PATCH 08/46] Minor change on the prisms_to_pyvista examples

---
 harmonica/visualization/prism.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/harmonica/visualization/prism.py b/harmonica/visualization/prism.py
index 3a8710337..2da21e610 100644
--- a/harmonica/visualization/prism.py
+++ b/harmonica/visualization/prism.py
@@ -54,8 +54,8 @@ def prisms_to_pyvista(prisms, properties=None):
     1
     >>> pv_grid.n_points
     8
-    >>> pv_grid.get_array("density")
-    array([2670.])
+    >>> pv_grid.cell_data["density"]
+    pyvista_ndarray([2670])
     >>> pv_grid.cell_bounds(0)
     [0.0, 1000.0, -2000.0, 2000.0, -5000.0, -3000.0]
 

From 814640030dd6715fcec2e62c531ce34447772a9b Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 16:12:49 -0300
Subject: [PATCH 09/46] Add tests for prisms_to_pyvista function

---
 harmonica/tests/test_visualizations.py | 139 +++++++++++++++++++++++++
 1 file changed, 139 insertions(+)
 create mode 100644 harmonica/tests/test_visualizations.py

diff --git a/harmonica/tests/test_visualizations.py b/harmonica/tests/test_visualizations.py
new file mode 100644
index 000000000..bd31991ee
--- /dev/null
+++ b/harmonica/tests/test_visualizations.py
@@ -0,0 +1,139 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Test functions from the visualization module.
+"""
+import pytest
+import numpy as np
+import numpy.testing as npt
+import xarray as xr
+
+from ..visualization.prism import prisms_to_pyvista
+
+try:
+    import vtk
+except ImportError:
+    vtk = None
+
+try:
+    import pyvista
+except ImportError:
+    pyvista = None
+
+
+@pytest.mark.skipif(vtk is not None, reason="vtk must be missing")
+def test_prisms_to_pyvista_missing_vtk():
+    """
+    Check error raise after calling prisms_to_pyvista when vtk is missing
+    """
+    prism = [0, 1, 0, 1, 0, 1]
+    with pytest.raises(ValueError) as exception:
+        prisms_to_pyvista(prism)
+    assert "'vtk'" in str(exception.value)
+
+
+@pytest.mark.skipif(pyvista is not None, reason="pyvista must be missing")
+def test_prisms_to_pyvista_missing_pyvista():
+    """
+    Check error raise after calling prisms_to_pyvista when pyvista is missing
+    """
+    prism = [0, 1, 0, 1, 0, 1]
+    with pytest.raises(ValueError) as exception:
+        prisms_to_pyvista(prism)
+    assert "'pyvista'" in str(exception.value)
+
+
+@pytest.fixture(name="prisms")
+def fixture_prisms():
+    """
+    Define a set of sample prisms
+    """
+    return [
+        [1e3, 5e3, -2e3, 4e3, -10e3, -7e3],
+        [200, 500, 600, 900, 100, 300],
+        [-2e3, 0, 5e3, 15e3, -4e3, 2e3],
+        [-10e3, -8e3, -10e3, -8e3, -10e3, -8e3],
+    ]
+
+
+@pytest.fixture(name="density", params=("numpy", "xarray"))
+def fixture_density(request):
+    """
+    Return a density array for the sample prisms
+    """
+    density = [2670.0, 2900.0, 3300.0, 3200.0]
+    if request.param == "xarray":
+        density = xr.DataArray(density)
+    return density
+
+
+@pytest.fixture(name="susceptibility", params=("numpy", "xarray"))
+def fixture_susceptibility(request):
+    """
+    Return a susceptibility array for the sample prisms
+    """
+    susceptibility = [2e-4, 2.5e-4, 5e-4, 1e-3]
+    if request.param == "xarray":
+        susceptibility = xr.DataArray(susceptibility)
+    return susceptibility
+
+
+def test_prisms_to_pyvista(prisms, density):
+    """
+    Test output of prisms_to_pyvista
+    """
+    pv_grid = prisms_to_pyvista(prisms, properties={"density": density})
+    assert pv_grid.n_cells == 4
+    assert pv_grid.n_points == 32
+    # Check coordinates of prisms
+    easting = pv_grid.points[:, 0]
+    northing = pv_grid.points[:, 1]
+    upward = pv_grid.points[:, 2]
+    npt.assert_allclose(easting.min(), -10e3)
+    npt.assert_allclose(easting.max(), 5e3)
+    npt.assert_allclose(northing.min(), -10e3)
+    npt.assert_allclose(northing.max(), 15e3)
+    npt.assert_allclose(upward.min(), -10e3)
+    npt.assert_allclose(upward.max(), 2e3)
+    # Check properties of the prisms
+    assert pv_grid.n_arrays == 1
+    assert pv_grid.array_names == ["density"]
+    npt.assert_allclose(pv_grid.get_array("density"), density)
+
+
+@pytest.mark.parametrize("n_props", [0, 1, 2])
+def test_prisms_to_pyvista_properties(n_props, prisms, density, susceptibility):
+    """
+    Test prisms_to_pyvista for different number of properties
+    """
+    properties = None
+    if n_props == 1:
+        properties = {"density": density}
+    elif n_props == 2:
+        properties = {"density": density, "susceptibility": susceptibility}
+    pv_grid = prisms_to_pyvista(prisms, properties=properties)
+    # Check if the grid has the right properties
+    if properties is None:
+        assert pv_grid.n_arrays == 0
+        assert pv_grid.array_names == []
+    else:
+        assert pv_grid.n_arrays == len(properties)
+        assert set(pv_grid.array_names) == set(properties.keys())
+        for prop in properties:
+            npt.assert_allclose(pv_grid.get_array(prop), properties[prop])
+
+
+def test_prisms_to_pyvista_2d_property(prisms, density):
+    """
+    Test if prisms_to_pyvista handles a 2d property array
+    """
+    density_2d = np.array(density).reshape((2, 2))
+    pv_grid = prisms_to_pyvista(prisms, properties={"density": density_2d})
+    assert pv_grid.n_arrays == 1
+    assert pv_grid.array_names == ["density"]
+    assert pv_grid.get_array("density").ndim == 1
+    npt.assert_allclose(pv_grid.get_array("density"), density)

From 2a55f8665de75a109d531863709d0899e688f7f3 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 16:12:59 -0300
Subject: [PATCH 10/46] Add tests for PrismLayer to_pyvista method

---
 harmonica/tests/test_prism_layer.py | 40 +++++++++++++++++++++++++++++
 1 file changed, 40 insertions(+)

diff --git a/harmonica/tests/test_prism_layer.py b/harmonica/tests/test_prism_layer.py
index e531d5e20..56f183a06 100644
--- a/harmonica/tests/test_prism_layer.py
+++ b/harmonica/tests/test_prism_layer.py
@@ -17,6 +17,16 @@
 
 from .. import prism_gravity, prism_layer
 
+try:
+    import vtk
+except ImportError:
+    vtk = None
+
+try:
+    import pyvista
+except ImportError:
+    pyvista = None
+
 
 @pytest.fixture(params=("numpy", "xarray"))
 def dummy_layer(request):
@@ -385,3 +395,33 @@ def test_prism_layer_gravity_density_nans(field, dummy_layer, prism_layer_with_h
         result,
         prism_gravity(coordinates, prisms, rho, field=field),
     )
+
+
+@pytest.mark.skipif(vtk is None or pyvista is None, reason="requires vtk and pyvista")
+def test_to_pyvista(dummy_layer):
+    """
+    Test the conversion of the prism layer to pyvista.UnstructuredGrid
+    """
+    (easting, northing), surface, reference, density = dummy_layer
+    layer = prism_layer(
+        (easting, northing), surface, reference, properties={"density": density}
+    )
+    pv_grid = layer.prism_layer.to_pyvista()
+    assert pv_grid.n_cells == 20
+    assert pv_grid.n_points == 20 * 8
+    # Check coordinates of prisms
+    pv_easting = pv_grid.points[:, 0]
+    pv_northing = pv_grid.points[:, 1]
+    pv_upward = pv_grid.points[:, 2]
+    d_easting, d_northing = 1, 2
+    npt.assert_allclose(layer.easting.min() - d_easting / 2, pv_easting.min())
+    npt.assert_allclose(layer.easting.max() + d_easting / 2, pv_easting.max())
+    npt.assert_allclose(layer.northing.min() - d_northing / 2, pv_northing.min())
+    npt.assert_allclose(layer.northing.max() + d_northing / 2, pv_northing.max())
+    npt.assert_allclose(layer.bottom.min(), pv_upward.min())
+    npt.assert_allclose(layer.top.max(), pv_upward.max())
+    # Check properties of the prisms
+    assert pv_grid.n_arrays == 1
+    assert pv_grid.array_names == ["density"]
+    assert pv_grid.get_array("density").ndim == 1
+    npt.assert_allclose(pv_grid.get_array("density"), layer.density.values.ravel())

From eca7c28da04158d411f9eb5eb9bcf0489ca4c205 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 16:14:03 -0300
Subject: [PATCH 11/46] Add skipif decorations

---
 harmonica/tests/test_visualizations.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/harmonica/tests/test_visualizations.py b/harmonica/tests/test_visualizations.py
index bd31991ee..427d3d7ca 100644
--- a/harmonica/tests/test_visualizations.py
+++ b/harmonica/tests/test_visualizations.py
@@ -82,6 +82,7 @@ def fixture_susceptibility(request):
     return susceptibility
 
 
+@pytest.mark.skipif(vtk is None or pyvista is None, reason="requires vtk and pyvista")
 def test_prisms_to_pyvista(prisms, density):
     """
     Test output of prisms_to_pyvista
@@ -105,6 +106,7 @@ def test_prisms_to_pyvista(prisms, density):
     npt.assert_allclose(pv_grid.get_array("density"), density)
 
 
+@pytest.mark.skipif(vtk is None or pyvista is None, reason="requires vtk and pyvista")
 @pytest.mark.parametrize("n_props", [0, 1, 2])
 def test_prisms_to_pyvista_properties(n_props, prisms, density, susceptibility):
     """
@@ -127,6 +129,7 @@ def test_prisms_to_pyvista_properties(n_props, prisms, density, susceptibility):
             npt.assert_allclose(pv_grid.get_array(prop), properties[prop])
 
 
+@pytest.mark.skipif(vtk is None or pyvista is None, reason="requires vtk and pyvista")
 def test_prisms_to_pyvista_2d_property(prisms, density):
     """
     Test if prisms_to_pyvista handles a 2d property array

From 6c155bad8e0e21857d0f36e1c10b6ce594e67cc4 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 16:17:44 -0300
Subject: [PATCH 12/46] Extend to_pyvista test without properties

---
 harmonica/tests/test_prism_layer.py | 26 ++++++++++++++++++--------
 1 file changed, 18 insertions(+), 8 deletions(-)

diff --git a/harmonica/tests/test_prism_layer.py b/harmonica/tests/test_prism_layer.py
index 56f183a06..237aa76a5 100644
--- a/harmonica/tests/test_prism_layer.py
+++ b/harmonica/tests/test_prism_layer.py
@@ -398,15 +398,21 @@ def test_prism_layer_gravity_density_nans(field, dummy_layer, prism_layer_with_h
 
 
 @pytest.mark.skipif(vtk is None or pyvista is None, reason="requires vtk and pyvista")
-def test_to_pyvista(dummy_layer):
+@pytest.mark.parametrize("properties", (False, True))
+def test_to_pyvista(dummy_layer, properties):
     """
     Test the conversion of the prism layer to pyvista.UnstructuredGrid
     """
     (easting, northing), surface, reference, density = dummy_layer
-    layer = prism_layer(
-        (easting, northing), surface, reference, properties={"density": density}
-    )
+    # Build the layer with or without properties
+    if properties:
+        properties = {"density": density}
+    else:
+        properties = None
+    layer = prism_layer((easting, northing), surface, reference, properties=properties)
+    # Convert the layer to pyvista UnstructuredGrid
     pv_grid = layer.prism_layer.to_pyvista()
+    # Check properties of the pyvista grid
     assert pv_grid.n_cells == 20
     assert pv_grid.n_points == 20 * 8
     # Check coordinates of prisms
@@ -421,7 +427,11 @@ def test_to_pyvista(dummy_layer):
     npt.assert_allclose(layer.bottom.min(), pv_upward.min())
     npt.assert_allclose(layer.top.max(), pv_upward.max())
     # Check properties of the prisms
-    assert pv_grid.n_arrays == 1
-    assert pv_grid.array_names == ["density"]
-    assert pv_grid.get_array("density").ndim == 1
-    npt.assert_allclose(pv_grid.get_array("density"), layer.density.values.ravel())
+    if properties is None:
+        assert pv_grid.n_arrays == 0
+        assert pv_grid.array_names == []
+    else:
+        assert pv_grid.n_arrays == 1
+        assert pv_grid.array_names == ["density"]
+        assert pv_grid.get_array("density").ndim == 1
+        npt.assert_allclose(pv_grid.get_array("density"), layer.density.values.ravel())

From 08a3eb8f9bc6dad08d11b5c044708faf50206ed3 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 16:29:19 -0300
Subject: [PATCH 13/46] Don't check if vtk is installed, it's a dependency of
 pyvista

---
 harmonica/tests/test_prism_layer.py    |  7 +------
 harmonica/tests/test_visualizations.py | 22 +++-------------------
 harmonica/visualization/prism.py       | 17 ++++++-----------
 3 files changed, 10 insertions(+), 36 deletions(-)

diff --git a/harmonica/tests/test_prism_layer.py b/harmonica/tests/test_prism_layer.py
index 237aa76a5..ec7f2bc81 100644
--- a/harmonica/tests/test_prism_layer.py
+++ b/harmonica/tests/test_prism_layer.py
@@ -17,11 +17,6 @@
 
 from .. import prism_gravity, prism_layer
 
-try:
-    import vtk
-except ImportError:
-    vtk = None
-
 try:
     import pyvista
 except ImportError:
@@ -397,7 +392,7 @@ def test_prism_layer_gravity_density_nans(field, dummy_layer, prism_layer_with_h
     )
 
 
-@pytest.mark.skipif(vtk is None or pyvista is None, reason="requires vtk and pyvista")
+@pytest.mark.skipif(pyvista is None, reason="requires pyvista")
 @pytest.mark.parametrize("properties", (False, True))
 def test_to_pyvista(dummy_layer, properties):
     """
diff --git a/harmonica/tests/test_visualizations.py b/harmonica/tests/test_visualizations.py
index 427d3d7ca..64e2853f3 100644
--- a/harmonica/tests/test_visualizations.py
+++ b/harmonica/tests/test_visualizations.py
@@ -14,28 +14,12 @@
 
 from ..visualization.prism import prisms_to_pyvista
 
-try:
-    import vtk
-except ImportError:
-    vtk = None
-
 try:
     import pyvista
 except ImportError:
     pyvista = None
 
 
-@pytest.mark.skipif(vtk is not None, reason="vtk must be missing")
-def test_prisms_to_pyvista_missing_vtk():
-    """
-    Check error raise after calling prisms_to_pyvista when vtk is missing
-    """
-    prism = [0, 1, 0, 1, 0, 1]
-    with pytest.raises(ValueError) as exception:
-        prisms_to_pyvista(prism)
-    assert "'vtk'" in str(exception.value)
-
-
 @pytest.mark.skipif(pyvista is not None, reason="pyvista must be missing")
 def test_prisms_to_pyvista_missing_pyvista():
     """
@@ -82,7 +66,7 @@ def fixture_susceptibility(request):
     return susceptibility
 
 
-@pytest.mark.skipif(vtk is None or pyvista is None, reason="requires vtk and pyvista")
+@pytest.mark.skipif(pyvista is None, reason="requires pyvista")
 def test_prisms_to_pyvista(prisms, density):
     """
     Test output of prisms_to_pyvista
@@ -106,7 +90,7 @@ def test_prisms_to_pyvista(prisms, density):
     npt.assert_allclose(pv_grid.get_array("density"), density)
 
 
-@pytest.mark.skipif(vtk is None or pyvista is None, reason="requires vtk and pyvista")
+@pytest.mark.skipif(pyvista is None, reason="requires pyvista")
 @pytest.mark.parametrize("n_props", [0, 1, 2])
 def test_prisms_to_pyvista_properties(n_props, prisms, density, susceptibility):
     """
@@ -129,7 +113,7 @@ def test_prisms_to_pyvista_properties(n_props, prisms, density, susceptibility):
             npt.assert_allclose(pv_grid.get_array(prop), properties[prop])
 
 
-@pytest.mark.skipif(vtk is None or pyvista is None, reason="requires vtk and pyvista")
+@pytest.mark.skipif(pyvista is None, reason="requires pyvista")
 def test_prisms_to_pyvista_2d_property(prisms, density):
     """
     Test if prisms_to_pyvista handles a 2d property array
diff --git a/harmonica/visualization/prism.py b/harmonica/visualization/prism.py
index 2da21e610..002c264ff 100644
--- a/harmonica/visualization/prism.py
+++ b/harmonica/visualization/prism.py
@@ -9,15 +9,12 @@
 """
 import numpy as np
 
-try:
-    import vtk
-except ImportError:
-    vtk = None
-
 try:
     import pyvista
 except ImportError:
     pyvista = None
+else:
+    import vtk
 
 
 def prisms_to_pyvista(prisms, properties=None):
@@ -61,13 +58,11 @@ def prisms_to_pyvista(prisms, properties=None):
 
 
     """
-    # Check if vkt and pyvista are installed
-    msg = "Missing optional dependency '{}' required for generating pyvista grids."
-    if vtk is None:
-        raise ValueError(msg.format("vtk"))
+    # Check if pyvista are installed
     if pyvista is None:
-        raise ValueError(msg.format("pyvista"))
-
+        raise ValueError(
+            "Missing optional dependency 'pyvista' required for building pyvista grids."
+        )
     # Get prisms and number of prisms
     prisms = np.atleast_2d(prisms)
     n_prisms = prisms.shape[0]

From 81db0129f1539f2d896be156c7d503b5b646035b Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 16:30:54 -0300
Subject: [PATCH 14/46] Run make format

---
 doc/conf.py                            | 2 +-
 examples/visualization/prism_layer.py  | 3 ++-
 harmonica/forward/prism_layer.py       | 2 +-
 harmonica/tests/test_visualizations.py | 2 +-
 harmonica/visualization/__init__.py    | 6 ++++++
 5 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/doc/conf.py b/doc/conf.py
index eeccd8713..2fdb388ca 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -5,8 +5,8 @@
 # This code is part of the Fatiando a Terra project (https://www.fatiando.org)
 #
 import datetime
-import pyvista
 
+import pyvista
 from sphinx_gallery.sorting import ExampleTitleSortKey
 
 import harmonica
diff --git a/examples/visualization/prism_layer.py b/examples/visualization/prism_layer.py
index bbd18650f..1a274c680 100644
--- a/examples/visualization/prism_layer.py
+++ b/examples/visualization/prism_layer.py
@@ -10,9 +10,10 @@
 
 """
 import pyproj
+import pyvista as pv
 import verde as vd
+
 import harmonica as hm
-import pyvista as pv
 
 # Read South Africa topography
 south_africa_topo = hm.datasets.fetch_south_africa_topography()
diff --git a/harmonica/forward/prism_layer.py b/harmonica/forward/prism_layer.py
index ed3cc67ec..dd1514f5c 100644
--- a/harmonica/forward/prism_layer.py
+++ b/harmonica/forward/prism_layer.py
@@ -13,8 +13,8 @@
 import verde as vd
 import xarray as xr
 
-from .prism import prism_gravity
 from ..visualization.prism import prisms_to_pyvista
+from .prism import prism_gravity
 
 
 def prism_layer(
diff --git a/harmonica/tests/test_visualizations.py b/harmonica/tests/test_visualizations.py
index 64e2853f3..3c7d46d73 100644
--- a/harmonica/tests/test_visualizations.py
+++ b/harmonica/tests/test_visualizations.py
@@ -7,9 +7,9 @@
 """
 Test functions from the visualization module.
 """
-import pytest
 import numpy as np
 import numpy.testing as npt
+import pytest
 import xarray as xr
 
 from ..visualization.prism import prisms_to_pyvista
diff --git a/harmonica/visualization/__init__.py b/harmonica/visualization/__init__.py
index e69de29bb..b348c0fa7 100644
--- a/harmonica/visualization/__init__.py
+++ b/harmonica/visualization/__init__.py
@@ -0,0 +1,6 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#

From 98fbe3dc0b55161d0e0a8b7deea8ba691b3c4c1f Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 16:46:16 -0300
Subject: [PATCH 15/46] Skip testing prisms_to_pyvista example

It needs pyvista and it's not always installed on CIs
---
 harmonica/visualization/prism.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/harmonica/visualization/prism.py b/harmonica/visualization/prism.py
index 002c264ff..fa01a4d7c 100644
--- a/harmonica/visualization/prism.py
+++ b/harmonica/visualization/prism.py
@@ -43,7 +43,7 @@ def prisms_to_pyvista(prisms, properties=None):
     Examples
     --------
 
-    >>> pv_grid = prisms_to_pyvista(
+    >>> pv_grid = prisms_to_pyvista( # doctest: +SKIP
     ...     [0, 1e3, -2e3, 2e3, -5e3, -3e3],
     ...     properties={"density": 2670},
     ... )

From d9640ea2ed40c195e997b35bafa394be66316a6d Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 16:50:57 -0300
Subject: [PATCH 16/46] Remove the example from the docstring

There's no need for that example to exists, it's not a public function,
so it's purpose was for testing. It's better to leave the testing
functions inside tests.
---
 harmonica/visualization/prism.py | 18 ------------------
 1 file changed, 18 deletions(-)

diff --git a/harmonica/visualization/prism.py b/harmonica/visualization/prism.py
index fa01a4d7c..7d8d9a896 100644
--- a/harmonica/visualization/prism.py
+++ b/harmonica/visualization/prism.py
@@ -39,24 +39,6 @@ def prisms_to_pyvista(prisms, properties=None):
     pv_grid : :class:`pyvista.UnstructuredGrid`
         :class:`pyvista.UnstructuredGrid` that represents the prisms with their
         properties (if any).
-
-    Examples
-    --------
-
-    >>> pv_grid = prisms_to_pyvista( # doctest: +SKIP
-    ...     [0, 1e3, -2e3, 2e3, -5e3, -3e3],
-    ...     properties={"density": 2670},
-    ... )
-    >>> pv_grid.n_cells
-    1
-    >>> pv_grid.n_points
-    8
-    >>> pv_grid.cell_data["density"]
-    pyvista_ndarray([2670])
-    >>> pv_grid.cell_bounds(0)
-    [0.0, 1000.0, -2000.0, 2000.0, -5000.0, -3000.0]
-
-
     """
     # Check if pyvista are installed
     if pyvista is None:

From 6d216347a2ba550ed761be002661274574059b89 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 16:56:48 -0300
Subject: [PATCH 17/46] Test the boundaries of each prism instead min-max
 values

---
 harmonica/tests/test_prism_layer.py    | 12 ++----------
 harmonica/tests/test_visualizations.py | 11 ++---------
 2 files changed, 4 insertions(+), 19 deletions(-)

diff --git a/harmonica/tests/test_prism_layer.py b/harmonica/tests/test_prism_layer.py
index ec7f2bc81..ddcbe57c4 100644
--- a/harmonica/tests/test_prism_layer.py
+++ b/harmonica/tests/test_prism_layer.py
@@ -411,16 +411,8 @@ def test_to_pyvista(dummy_layer, properties):
     assert pv_grid.n_cells == 20
     assert pv_grid.n_points == 20 * 8
     # Check coordinates of prisms
-    pv_easting = pv_grid.points[:, 0]
-    pv_northing = pv_grid.points[:, 1]
-    pv_upward = pv_grid.points[:, 2]
-    d_easting, d_northing = 1, 2
-    npt.assert_allclose(layer.easting.min() - d_easting / 2, pv_easting.min())
-    npt.assert_allclose(layer.easting.max() + d_easting / 2, pv_easting.max())
-    npt.assert_allclose(layer.northing.min() - d_northing / 2, pv_northing.min())
-    npt.assert_allclose(layer.northing.max() + d_northing / 2, pv_northing.max())
-    npt.assert_allclose(layer.bottom.min(), pv_upward.min())
-    npt.assert_allclose(layer.top.max(), pv_upward.max())
+    for i, prism in enumerate(layer.prism_layer._to_prisms()):
+        npt.assert_allclose(prism, pv_grid.cell_bounds(i))
     # Check properties of the prisms
     if properties is None:
         assert pv_grid.n_arrays == 0
diff --git a/harmonica/tests/test_visualizations.py b/harmonica/tests/test_visualizations.py
index 3c7d46d73..8b211b7e8 100644
--- a/harmonica/tests/test_visualizations.py
+++ b/harmonica/tests/test_visualizations.py
@@ -75,15 +75,8 @@ def test_prisms_to_pyvista(prisms, density):
     assert pv_grid.n_cells == 4
     assert pv_grid.n_points == 32
     # Check coordinates of prisms
-    easting = pv_grid.points[:, 0]
-    northing = pv_grid.points[:, 1]
-    upward = pv_grid.points[:, 2]
-    npt.assert_allclose(easting.min(), -10e3)
-    npt.assert_allclose(easting.max(), 5e3)
-    npt.assert_allclose(northing.min(), -10e3)
-    npt.assert_allclose(northing.max(), 15e3)
-    npt.assert_allclose(upward.min(), -10e3)
-    npt.assert_allclose(upward.max(), 2e3)
+    for i, prism in enumerate(prisms):
+        npt.assert_allclose(prism, pv_grid.cell_bounds(i))
     # Check properties of the prisms
     assert pv_grid.n_arrays == 1
     assert pv_grid.array_names == ["density"]

From 3c88f935686d1147f7836ec4b4adea7af4473334 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 17:12:56 -0300
Subject: [PATCH 18/46] Require vtk >= 9 in environment.yml

---
 environment.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/environment.yml b/environment.yml
index 5b32fdf6f..3d994c29c 100644
--- a/environment.yml
+++ b/environment.yml
@@ -16,7 +16,7 @@ dependencies:
   - xarray
   # Optional requirements
   - pyvista
-  - vtk>9
+  - vtk>=9
   # Development requirements
   - boule
   - pyproj

From 10657382c8ac9da3c0cf79b46bc699827e4c9740 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 17:13:11 -0300
Subject: [PATCH 19/46] Add optional dependencies to installation instructions

---
 doc/install.rst | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/doc/install.rst b/doc/install.rst
index 75e015429..503fbece8 100644
--- a/doc/install.rst
+++ b/doc/install.rst
@@ -19,6 +19,12 @@ doesn't interfere with any other Python installations in your system.
 Dependencies
 ------------
 
+The required dependencies should be installed automatically when you install
+Harmonica using ``conda`` or ``pip``. Optional dependencies have to be
+installed manually.
+
+Required:
+
 * `numpy <http://www.numpy.org/>`__
 * `pandas <http://pandas.pydata.org/>`__
 * `numba <https://numba.pydata.org/>`__
@@ -28,6 +34,12 @@ Dependencies
 * `pooch <http://www.fatiando.org/pooch/>`__
 * `verde <http://www.fatiando.org/verde/>`__
 
+Optional:
+
+* `pyvista <https://www.pyvista.org/>`__ and
+  `vtk <https://vtk.org/>`__ (>= 9): for 3D visualizations.
+  See :meth:`harmonica.DatasetAccessorPrismLayer.to_pyvista`.
+
 The examples in the :ref:`gallery` also use:
 
 * `boule <http://www.fatiando.org/boule/>`__

From 31a64c7da8f956e94becee01fb2f334abb9fcb89 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 17:13:36 -0300
Subject: [PATCH 20/46] Add pyvista and vtk to env/requirements-docs.txt

---
 env/requirements-docs.txt | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/env/requirements-docs.txt b/env/requirements-docs.txt
index 506e2be6a..0aad2e042 100644
--- a/env/requirements-docs.txt
+++ b/env/requirements-docs.txt
@@ -6,3 +6,5 @@ boule
 pyproj
 matplotlib
 cartopy>=0.18
+pyvista
+vtk>=9

From 8bff006030ffce8623cb9efb3816155ce5bb8443 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 17:13:50 -0300
Subject: [PATCH 21/46] Add a requirements-optional.txt file

---
 requirements-optional.txt | 2 ++
 1 file changed, 2 insertions(+)
 create mode 100644 requirements-optional.txt

diff --git a/requirements-optional.txt b/requirements-optional.txt
new file mode 100644
index 000000000..0e8bb24cf
--- /dev/null
+++ b/requirements-optional.txt
@@ -0,0 +1,2 @@
+pyvista
+vtk>=9

From d72c1c56028e0d18198714dc976df4c1333738a5 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Wed, 1 Dec 2021 17:19:03 -0300
Subject: [PATCH 22/46] Add a dependencies option to the test matrix in Actions

Now the tests are run with and without optional dependencies.
---
 .github/workflows/test.yml | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index ec613e8d5..d5e54d102 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -30,7 +30,7 @@ jobs:
   #############################################################################
   # Run tests and upload to codecov
   test:
-    name: ${{ matrix.os }} py${{ matrix.python }}
+    name: ${{ matrix.os }} py${{ matrix.python }} ${{ matrix.dependencies }}
     runs-on: ${{ matrix.os }}-latest
     strategy:
       # Otherwise, the workflow would stop if a single job fails. We want to
@@ -41,9 +41,11 @@ jobs:
         python: ["3.6", "3.9"]
         # If "optional", will install non-required dependencies in the build
         # environment. Otherwise, only required dependencies are installed.
-        dependencies: [""]
+        dependencies: ["", "optional"]
     env:
       REQUIREMENTS: requirements.txt env/requirements-tests.txt
+      REQUIREMENTS_OPTIONAL: requirements-optional.txt
+      DEPENDENCIES: ${{ matrix.dependencies }}
       # Used to tag codecov submissions
       OS: ${{ matrix.os }}
       PYTHON: ${{ matrix.python }}
@@ -102,6 +104,10 @@ jobs:
                 cat $requirement >> $requirements_file
               done
           fi
+          if [ "$DEPENDENCIES" == "optional" ]; then
+              echo "Capturing optional dependencies from $REQUIREMENTS_OPTIONAL"
+              cat $REQUIREMENTS_OPTIONAL >> $requirements_file
+          fi
           if [ -f $requirements_file ]; then
               echo "Collected dependencies:"
               cat $requirements_file
@@ -159,7 +165,7 @@ jobs:
         with:
           token: ${{ secrets.CODECOV_TOKEN }}
           file: ./coverage.xml
-          env_vars: OS,PYTHON
+          env_vars: OS,PYTHON,DEPENDENCIES
           # Don't mark the job as failed if the upload fails for some reason.
           # It does sometimes but shouldn't be the reason for running
           # everything again unless something else is broken.

From 678f4dd87d07b9dcd63a60dab89eb2c2f773a8c4 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Thu, 2 Dec 2021 11:03:18 -0300
Subject: [PATCH 23/46] Set DISPLAY before building docs on CIs

---
 .github/workflows/docs.yml | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index c8f1e7e14..cca311e55 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -122,7 +122,8 @@ jobs:
           HARMONICA_DATA_DIR: ${{ runner.temp }}/cache/harmonica
 
       - name: Build the documentation
-        run: make -C doc clean all
+        # Configure DISPLAY env variable tu run pyvista on headless system
+        run: DISPLAY=:99.0; make -C doc clean all
 
       # Store the docs as a build artifact so we can deploy it later
       - name: Upload HTML documentation as an artifact

From d65181f509a847b9c7def2749d6c893ec92288b5 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Thu, 2 Dec 2021 11:14:38 -0300
Subject: [PATCH 24/46] Fix typo

---
 .github/workflows/docs.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index cca311e55..700088125 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -123,7 +123,7 @@ jobs:
 
       - name: Build the documentation
         # Configure DISPLAY env variable tu run pyvista on headless system
-        run: DISPLAY=:99.0; make -C doc clean all
+        run: DISPLAY=:99.0 make -C doc clean all
 
       # Store the docs as a build artifact so we can deploy it later
       - name: Upload HTML documentation as an artifact

From 72c6a2abcfeb5bb104e9a3970df9afede336ccb1 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Thu, 2 Dec 2021 11:30:21 -0300
Subject: [PATCH 25/46] Install xvfb and follow pyvista guide to run in on
 headless systems

---
 .github/workflows/docs.yml | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 700088125..31f3d7f7d 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -122,8 +122,17 @@ jobs:
           HARMONICA_DATA_DIR: ${{ runner.temp }}/cache/harmonica
 
       - name: Build the documentation
-        # Configure DISPLAY env variable tu run pyvista on headless system
-        run: DISPLAY=:99.0 make -C doc clean all
+        run: |
+          # Install xvfb and run some commands to allow pyvista to run on
+          # a headless system.
+          sudo apt-get install xvfb
+          export DISPLAY=:99.0
+          export PYVISTA_OFF_SCREEN=true
+          export PYVISTA_USE_IPYVTK=true
+          Xvfb :99 -screen 0 1024x768x24 > /dev/null 2>&1 &
+          sleep 3
+          # Build the docs
+          make -C doc clean all
 
       # Store the docs as a build artifact so we can deploy it later
       - name: Upload HTML documentation as an artifact

From 603657e6edfde807c8c7ba7fddb910808b671abf Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Thu, 2 Dec 2021 15:45:04 -0300
Subject: [PATCH 26/46] Add introductory text to the example

---
 examples/visualization/prism_layer.py | 23 +++++++++++++++++++++--
 1 file changed, 21 insertions(+), 2 deletions(-)

diff --git a/examples/visualization/prism_layer.py b/examples/visualization/prism_layer.py
index 1a274c680..e9572024b 100644
--- a/examples/visualization/prism_layer.py
+++ b/examples/visualization/prism_layer.py
@@ -5,8 +5,27 @@
 # This code is part of the Fatiando a Terra project (https://www.fatiando.org)
 #
 """
-3D plot of layer of prisms
-==========================
+Plot a prism layer in 3D
+========================
+
+The :func:`harmonica.prism_layer` allows to easily create a layer of prisms
+whose top and bottom boundaries might drape a certain surface, like topography,
+bathymetry or the Moho discontinuity. It returns a :class:`xarray.Dataset`
+object with the horizontal coordinates of the center of each prism, their top
+and bottom boundaries and their physical properties, like their density.
+Through the ``prism_layer`` accessor (see
+:class:`harmonica.DatasetAccessorPrismLayer`) we can call some methods for our
+prism layer. For example, the
+:meth:`harmonica.DatasetAccessorPrismLayer.gravity` method lets us compute the
+gravitational fields of the layer on any set of observation points.
+Another interesting method is the
+:meth:`harmonica.DatasetAccessorPrismLayer.to_pyvista`: it converts the prism
+layer into a :class:`pyvista.UnstructuredGrid` that could be easily plotted
+through ``pyvista``.
+
+In this example we will show how we can build a prism layer out of a topography
+and bathymetry grid for South Africa and how we can visualize the layer as a 3D
+plot using ``pyvista``.
 
 """
 import pyproj

From 9ce454c45d15e9cd2e740c70dbdbc371f6e4533f Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Thu, 2 Dec 2021 15:45:29 -0300
Subject: [PATCH 27/46] Add a ceiling light to the pyvista plot

---
 examples/visualization/prism_layer.py | 18 +++++++++++++++++-
 1 file changed, 17 insertions(+), 1 deletion(-)

diff --git a/examples/visualization/prism_layer.py b/examples/visualization/prism_layer.py
index e9572024b..6e6cb3bdc 100644
--- a/examples/visualization/prism_layer.py
+++ b/examples/visualization/prism_layer.py
@@ -63,9 +63,25 @@
 print(pv_grid)
 
 # Plot with pyvista
-plotter = pv.Plotter()
+plotter = pv.Plotter(lighting="three_lights")
 plotter.add_mesh(pv_grid, scalars="density")
 plotter.set_scale(zscale=75)  # exaggerate the vertical coordinate
 plotter.camera_position = "xz"
 plotter.camera.elevation = 15
+plotter.camera.azimuth = 35
+plotter.camera.zoom(1.2)
+
+# Add a ceiling light
+west, east, south, north = vd.get_region((prisms.easting, prisms.northing))
+easting_center, northing_center = (east + west) / 2, (north + south) / 2
+light = pv.Light(
+    position=(easting_center, northing_center, 10e3),
+    focal_point=(easting_center, northing_center, 0),
+    intensity=0.3,
+    light_type="scene light",  # the light doesn't move with the camera
+    positional=False,  # the light comes from infinity
+)
+plotter.add_light(light)
+
+plotter.show_axes()
 plotter.show()

From 4eec55c14baef3870c2cfc08c6ebf555af759066 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Thu, 2 Dec 2021 15:48:28 -0300
Subject: [PATCH 28/46] Fix window_size of the plot

---
 examples/visualization/prism_layer.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/examples/visualization/prism_layer.py b/examples/visualization/prism_layer.py
index 6e6cb3bdc..50c395331 100644
--- a/examples/visualization/prism_layer.py
+++ b/examples/visualization/prism_layer.py
@@ -63,11 +63,11 @@
 print(pv_grid)
 
 # Plot with pyvista
-plotter = pv.Plotter(lighting="three_lights")
+plotter = pv.Plotter(lighting="three_lights", window_size=(1000, 800))
 plotter.add_mesh(pv_grid, scalars="density")
 plotter.set_scale(zscale=75)  # exaggerate the vertical coordinate
 plotter.camera_position = "xz"
-plotter.camera.elevation = 15
+plotter.camera.elevation = 20
 plotter.camera.azimuth = 35
 plotter.camera.zoom(1.2)
 

From d4a1f0013274d4c458150aa753c2585b1db34d12 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Tue, 7 Dec 2021 09:32:44 -0300
Subject: [PATCH 29/46] Use repr capturing instead of printing the pv_grid

Co-authored-by: Bane Sullivan <banesullivan@gmail.com>
---
 examples/visualization/prism_layer.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/examples/visualization/prism_layer.py b/examples/visualization/prism_layer.py
index 50c395331..740609515 100644
--- a/examples/visualization/prism_layer.py
+++ b/examples/visualization/prism_layer.py
@@ -60,7 +60,9 @@
 
 # Create a pyvista UnstructuredGrid from the prism layer
 pv_grid = prisms.prism_layer.to_pyvista()
-print(pv_grid)
+pv_grid
+
+###############################################################################
 
 # Plot with pyvista
 plotter = pv.Plotter(lighting="three_lights", window_size=(1000, 800))

From e6ddfa8f33d36ff5acf56048450c9d0e632374d7 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Tue, 7 Dec 2021 09:33:15 -0300
Subject: [PATCH 30/46] Raise ImportError if pyvista is missing

Replace ValueError for ImportError.

Co-authored-by: Bane Sullivan <banesullivan@gmail.com>
---
 harmonica/visualization/prism.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/harmonica/visualization/prism.py b/harmonica/visualization/prism.py
index 7d8d9a896..3cdb2efff 100644
--- a/harmonica/visualization/prism.py
+++ b/harmonica/visualization/prism.py
@@ -42,7 +42,7 @@ def prisms_to_pyvista(prisms, properties=None):
     """
     # Check if pyvista are installed
     if pyvista is None:
-        raise ValueError(
+        raise ImportError(
             "Missing optional dependency 'pyvista' required for building pyvista grids."
         )
     # Get prisms and number of prisms

From 172ad1558abe3207704d48b29d3f6eae0097ba8c Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Tue, 7 Dec 2021 09:41:49 -0300
Subject: [PATCH 31/46] Add matplotlib to image_scrapers on Sphinx configs

Co-authored-by: Bane Sullivan <banesullivan@gmail.com>
---
 doc/conf.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/conf.py b/doc/conf.py
index 2fdb388ca..00117c152 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -97,7 +97,7 @@
     # Insert links to documentation of objects in the examples
     "reference_url": {"harmonica": None},
     # Add pyvista to the image scrapers
-    "image_scrapers": ("pyvista"),
+    "image_scrapers": ("pyvista", "matplotlib"),
 }
 
 

From 16cb554225c1597116c4baa1874aac9656fbac0e Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Tue, 7 Dec 2021 10:23:36 -0300
Subject: [PATCH 32/46] Test if ImportError is raised when pyvista is missing

---
 harmonica/tests/test_visualizations.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/harmonica/tests/test_visualizations.py b/harmonica/tests/test_visualizations.py
index 8b211b7e8..3873ffcf7 100644
--- a/harmonica/tests/test_visualizations.py
+++ b/harmonica/tests/test_visualizations.py
@@ -26,7 +26,7 @@ def test_prisms_to_pyvista_missing_pyvista():
     Check error raise after calling prisms_to_pyvista when pyvista is missing
     """
     prism = [0, 1, 0, 1, 0, 1]
-    with pytest.raises(ValueError) as exception:
+    with pytest.raises(ImportError) as exception:
         prisms_to_pyvista(prism)
     assert "'pyvista'" in str(exception.value)
 

From 93986528051c6a188cba453c5f214fe712416e3a Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Fri, 4 Mar 2022 12:46:50 -0300
Subject: [PATCH 33/46] Make prism_to_pyvista to get only 1d arrays as
 properties

Raise ValueError if any property is passed as multidimensional array.
Update test function to check if error is raised.
Make the DataAccessorPrismLayer.to_pyvista() method to ravel the
properties arrays before passing them to prism_to_pyvista.
---
 harmonica/forward/prism_layer.py       |  3 ++-
 harmonica/tests/test_visualizations.py | 11 ++++-------
 harmonica/visualization/prism.py       | 11 +++++++++--
 3 files changed, 15 insertions(+), 10 deletions(-)

diff --git a/harmonica/forward/prism_layer.py b/harmonica/forward/prism_layer.py
index dd1514f5c..660648e82 100644
--- a/harmonica/forward/prism_layer.py
+++ b/harmonica/forward/prism_layer.py
@@ -462,6 +462,7 @@ def to_pyvista(self):
         properties = None
         if self._obj.data_vars:
             properties = {
-                data_var: self._obj[data_var] for data_var in self._obj.data_vars
+                data_var: np.asarray(self._obj[data_var]).ravel()
+                for data_var in self._obj.data_vars
             }
         return prisms_to_pyvista(prisms, properties=properties)
diff --git a/harmonica/tests/test_visualizations.py b/harmonica/tests/test_visualizations.py
index 3873ffcf7..06bdd3f4b 100644
--- a/harmonica/tests/test_visualizations.py
+++ b/harmonica/tests/test_visualizations.py
@@ -107,13 +107,10 @@ def test_prisms_to_pyvista_properties(n_props, prisms, density, susceptibility):
 
 
 @pytest.mark.skipif(pyvista is None, reason="requires pyvista")
-def test_prisms_to_pyvista_2d_property(prisms, density):
+def test_prisms_to_pyvista_error_2d_property(prisms, density):
     """
-    Test if prisms_to_pyvista handles a 2d property array
+    Test if prisms_to_pyvista raises error on property as 2D array
     """
     density_2d = np.array(density).reshape((2, 2))
-    pv_grid = prisms_to_pyvista(prisms, properties={"density": density_2d})
-    assert pv_grid.n_arrays == 1
-    assert pv_grid.array_names == ["density"]
-    assert pv_grid.get_array("density").ndim == 1
-    npt.assert_allclose(pv_grid.get_array("density"), density)
+    with pytest.raises(ValueError):
+        prisms_to_pyvista(prisms, properties={"density": density_2d})
diff --git a/harmonica/visualization/prism.py b/harmonica/visualization/prism.py
index 3cdb2efff..abc59d4fd 100644
--- a/harmonica/visualization/prism.py
+++ b/harmonica/visualization/prism.py
@@ -29,7 +29,7 @@ def prisms_to_pyvista(prisms, properties=None):
         ``west``, ``east``, ``south``, ``north``, ``bottom``, ``top``.
     properties : dict or None (optional)
         Dictionary with the physical properties of the prisms.
-        Each key should be a string and their value an array.
+        Each key should be a string and their value as a 1D array.
         If None, no property will be added to the
         :class:`pyvista.UnstructuredGrid`.
         Default to None.
@@ -57,8 +57,15 @@ def prisms_to_pyvista(prisms, properties=None):
     # Add properties to the grid
     if properties is not None:
         for name, prop in properties.items():
+            # Check if the property is given as 1d array
+            prop = np.atleast_1d(prop)
+            if prop.ndim > 1:
+                raise ValueError(
+                    f"Multidimensional array found in '{name}' property. "
+                    + "Please, pass prism properties as 1d arrays."
+                )
             # Assign the property to the cell_data
-            pv_grid.cell_data[name] = np.atleast_1d(prop).ravel()
+            pv_grid.cell_data[name] = prop
     return pv_grid
 
 

From e3a88559e7b56760040882fe6c4e1b6b5e031051 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Fri, 4 Mar 2022 14:29:55 -0300
Subject: [PATCH 34/46] Optimize generation of prisms vertices

Instead of generating each set of vertices per prism using a Python for
loop, we have defined a new function that runs this conversion for
a full set of prisms boundaries using Numpy for optimization.
Only 24 for loops iterations are used, independently of the size of the
prisms list.
---
 harmonica/visualization/prism.py | 123 +++++++++++++++++++++++--------
 1 file changed, 94 insertions(+), 29 deletions(-)

diff --git a/harmonica/visualization/prism.py b/harmonica/visualization/prism.py
index abc59d4fd..c0bf961aa 100644
--- a/harmonica/visualization/prism.py
+++ b/harmonica/visualization/prism.py
@@ -49,7 +49,7 @@ def prisms_to_pyvista(prisms, properties=None):
     prisms = np.atleast_2d(prisms)
     n_prisms = prisms.shape[0]
     # Get the vertices of the prisms
-    vertices = np.vstack([_get_prism_vertices(prism) for prism in prisms])
+    vertices = _prisms_boundaries_to_vertices(prisms)
     # Generate the cells for the hexahedrons
     cells = np.arange(n_prisms * 8).reshape([n_prisms, 8])
     # Build the UnstructuredGrid
@@ -69,14 +69,42 @@ def prisms_to_pyvista(prisms, properties=None):
     return pv_grid
 
 
-def _get_prism_vertices(prism):
+def _prisms_boundaries_to_vertices(prisms):
     """
-    Return the vertices of the given prism
+    Converts prisms boundaries to set of vertices for each prism
+
+    The vertices for each prism will be in the following order
+
+    .. code-block::
+
+                                7-------6
+                               /|      /|
+                              / |     / |
+                             4-------5  |
+        up   northing        |  |    |  |
+        ^   ~                |  3----|--2
+        |  /                 | /     | /
+        | /                  |/      |/
+         ------> easting     0-------1
+
+    So the vertices of a single prism should be like:
+
+    .. code-block:: python
+
+        [w, s, bottom],
+        [e, s, bottom],
+        [e, n, bottom],
+        [w, n, bottom],
+        [w, s, top],
+        [e, s, top],
+        [e, n, top],
+        [w, n, top],
+
 
     Parameters
     ----------
-    prism : list or 1d-array
-        List or 1d-array with the boundaries of a single prism in the following
+    prism : 2d-array
+        2d-array with the boundaries of a set of prisms in the following
         order: ``west``, ``east``, ``south``, ``north``, ``bottom``, ``top``.
 
     Returns
@@ -89,29 +117,66 @@ def _get_prism_vertices(prism):
 
     Examples
     --------
+    >>> _prisms_boundaries_to_vertices(np.array([[-1, 1, -2, 2, -3, 3]]))
+    array([[-1., -2., -3.],
+           [ 1., -2., -3.],
+           [ 1.,  2., -3.],
+           [-1.,  2., -3.],
+           [-1., -2.,  3.],
+           [ 1., -2.,  3.],
+           [ 1.,  2.,  3.],
+           [-1.,  2.,  3.]])
+    >>> _prisms_boundaries_to_vertices(
+    ...     np.array([[-1, 1, -2, 2, -3, 3], [-4, 4, -5, 5, -6, 6]])
+    ... )
+    array([[-1., -2., -3.],
+           [ 1., -2., -3.],
+           [ 1.,  2., -3.],
+           [-1.,  2., -3.],
+           [-1., -2.,  3.],
+           [ 1., -2.,  3.],
+           [ 1.,  2.,  3.],
+           [-1.,  2.,  3.],
+           [-4., -5., -6.],
+           [ 4., -5., -6.],
+           [ 4.,  5., -6.],
+           [-4.,  5., -6.],
+           [-4., -5.,  6.],
+           [ 4., -5.,  6.],
+           [ 4.,  5.,  6.],
+           [-4.,  5.,  6.]])
+    """
+    # Get number of prisms
+    n_prisms = prisms.shape[0]
 
-    >>> _get_prism_vertices([-1, 1, -2, 2, -3, 3])
-    array([[-1, -2, -3],
-           [ 1, -2, -3],
-           [ 1,  2, -3],
-           [-1,  2, -3],
-           [-1, -2,  3],
-           [ 1, -2,  3],
-           [ 1,  2,  3],
-           [-1,  2,  3]])
+    # Allocate vertices array
+    vertices = np.empty((n_prisms, 8, 3))
 
-    """
-    w, e, s, n, bottom, top = prism[:]
-    vertices = np.array(
-        [
-            [w, s, bottom],
-            [e, s, bottom],
-            [e, n, bottom],
-            [w, n, bottom],
-            [w, s, top],
-            [e, s, top],
-            [e, n, top],
-            [w, n, top],
-        ]
-    )
-    return vertices
+    # Define a dictionary with the indices of the vertices that contain each
+    # boundary of the prism.
+    # For example, the west boundary is present only in the vertices
+    # number 0, 3, 4 and 7.
+    indices = {
+        "west": (0, 3, 4, 7),
+        "east": (1, 2, 5, 6),
+        "south": (0, 1, 4, 5),
+        "north": (2, 3, 6, 7),
+        "bottom": (0, 1, 2, 3),
+        "top": (4, 5, 6, 7),
+    }
+
+    # Assign the values to each vertex
+    for i, boundary in enumerate(indices):
+        # Determine at which component of the vertices should the current
+        # boundary be assigned to.
+        #   The west and east (i = 0 and i = 1) should go to the component 0.
+        #   The south and north (i = 2 and i = 3) should go to the component 1.
+        #   The bottom and top (i = 4 and i = 5) should go to the component 2.
+        component = i // 2
+        # Assign vertices components
+        for vertex in indices[boundary]:
+            vertices[:, vertex, component] = prisms[:, i]
+
+    # Reshape the vertices array so it has (M, 3) shape, where M is the total
+    # number of vertices.
+    return vertices.reshape((n_prisms * 8, 3))

From 557f4212fed63801b431efa69ae0d9bf33d9f001 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Fri, 4 Mar 2022 14:38:46 -0300
Subject: [PATCH 35/46] Improve docstring of the
 _prisms_boundaries_to_vertices() function

---
 harmonica/visualization/prism.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/harmonica/visualization/prism.py b/harmonica/visualization/prism.py
index c0bf961aa..99eb6430f 100644
--- a/harmonica/visualization/prism.py
+++ b/harmonica/visualization/prism.py
@@ -103,8 +103,9 @@ def _prisms_boundaries_to_vertices(prisms):
 
     Parameters
     ----------
-    prism : 2d-array
-        2d-array with the boundaries of a set of prisms in the following
+    prisms : 2d-array
+        2d-array with the boundaries of a set of prisms. Each row of the array
+        should contain the boundaries of a single prism in the following
         order: ``west``, ``east``, ``south``, ``north``, ``bottom``, ``top``.
 
     Returns
@@ -113,6 +114,8 @@ def _prisms_boundaries_to_vertices(prisms):
         2D array with the coordinates of the vertices of the prism. Each row of
         the array corresponds to the coordinate of a single vertex in the
         following order: ``easting``, ``northing``, ``upward``.
+        The shape of this array is ``(M, 3)``, where ``M`` is the total number
+        of vertices in the whole set of prisms (number of prisms times 8).
         The order of the vertices is fixed to be compatible with VTK.
 
     Examples

From 408ecbfdeaddfd181383ef8e6f3fd7e49c6a0570 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Fri, 18 Mar 2022 13:27:39 -0300
Subject: [PATCH 36/46] Add PyVista and VTK to requirements-test.txt

---
 env/requirements-tests.txt | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/env/requirements-tests.txt b/env/requirements-tests.txt
index a968055da..605da369d 100644
--- a/env/requirements-tests.txt
+++ b/env/requirements-tests.txt
@@ -3,3 +3,5 @@ boule
 pytest
 pytest-cov
 coverage
+pyvista
+vtk>=9

From f91f410eb64dad0a2dea1b6dd388b03044f78dd6 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Mon, 4 Apr 2022 11:40:09 -0300
Subject: [PATCH 37/46] Add pyvista and vtk to setup.cfg

Ditch the old requirements-optional.txt file.
---
 requirements-optional.txt | 2 --
 setup.cfg                 | 5 +++++
 2 files changed, 5 insertions(+), 2 deletions(-)
 delete mode 100644 requirements-optional.txt

diff --git a/requirements-optional.txt b/requirements-optional.txt
deleted file mode 100644
index 0e8bb24cf..000000000
--- a/requirements-optional.txt
+++ /dev/null
@@ -1,2 +0,0 @@
-pyvista
-vtk>=9
diff --git a/setup.cfg b/setup.cfg
index 8c7b7eb43..e3132b3da 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -49,6 +49,11 @@ install_requires =
     xarray
     verde>=1.5.0
 
+[options.extras_require]
+visualizations =
+    pyvista
+    vtk>=9
+
 [options.package_data]
 harmonica.tests = data/*, baseline/*
 harmonica.datasets = registry.txt

From aa22a17ba006b6472dc2442157ec00bf3a952779 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Mon, 4 Apr 2022 11:50:46 -0300
Subject: [PATCH 38/46] Add test for error raising when pyvista is missing

---
 harmonica/tests/test_visualizations.py | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/harmonica/tests/test_visualizations.py b/harmonica/tests/test_visualizations.py
index 06bdd3f4b..2d62b5897 100644
--- a/harmonica/tests/test_visualizations.py
+++ b/harmonica/tests/test_visualizations.py
@@ -10,6 +10,7 @@
 import numpy as np
 import numpy.testing as npt
 import pytest
+from unittest.mock import patch
 import xarray as xr
 
 from ..visualization.prism import prisms_to_pyvista
@@ -114,3 +115,13 @@ def test_prisms_to_pyvista_error_2d_property(prisms, density):
     density_2d = np.array(density).reshape((2, 2))
     with pytest.raises(ValueError):
         prisms_to_pyvista(prisms, properties={"density": density_2d})
+
+
+@patch("harmonica.visualization.prism.pyvista", None)
+def test_prisms_pyvista_missing_error(prisms, density):
+    """
+    Check if prisms_to_pyvista raises error if pyvista is missing
+    """
+    # Check if error is raised
+    with pytest.raises(ImportError):
+        prisms_to_pyvista(prisms, properties={"density": density})

From fc7c2e77c355f4864bdccf0a179354b7b75cf62f Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Mon, 4 Apr 2022 11:53:23 -0300
Subject: [PATCH 39/46] Run isort

---
 harmonica/tests/test_visualizations.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/harmonica/tests/test_visualizations.py b/harmonica/tests/test_visualizations.py
index 2d62b5897..6bc3b071d 100644
--- a/harmonica/tests/test_visualizations.py
+++ b/harmonica/tests/test_visualizations.py
@@ -7,10 +7,11 @@
 """
 Test functions from the visualization module.
 """
+from unittest.mock import patch
+
 import numpy as np
 import numpy.testing as npt
 import pytest
-from unittest.mock import patch
 import xarray as xr
 
 from ..visualization.prism import prisms_to_pyvista

From 4d2d877f6e48144a62d56ebe128b12efef4bcc18 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Mon, 4 Apr 2022 15:57:21 -0300
Subject: [PATCH 40/46] Rename example file to avoid Sphinx issues with links

---
 examples/visualization/{prism_layer.py => prism_layer_pyvista.py} | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename examples/visualization/{prism_layer.py => prism_layer_pyvista.py} (100%)

diff --git a/examples/visualization/prism_layer.py b/examples/visualization/prism_layer_pyvista.py
similarity index 100%
rename from examples/visualization/prism_layer.py
rename to examples/visualization/prism_layer_pyvista.py

From 63b5b15e731fd708c0b245d1b16b2dfb65b7646e Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Fri, 22 Apr 2022 15:10:39 -0300
Subject: [PATCH 41/46] Make visualization a public module

Make prisms_to_pyvista public and available through
harmonica.visualization module.
Add visualization module to the API Reference.
---
 doc/api/index.rst                   | 8 ++++++++
 harmonica/visualization/__init__.py | 1 +
 harmonica/visualization/prism.py    | 5 ++++-
 3 files changed, 13 insertions(+), 1 deletion(-)

diff --git a/doc/api/index.rst b/doc/api/index.rst
index 40365b84f..9a2374511 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -55,6 +55,14 @@ Input and Output
 
     load_icgem_gdf
 
+Visualization
+-------------
+
+.. autosummary::
+   :toctree: generated/
+
+    visualization.prisms_to_pyvista
+
 Synthetic models and surveys
 ----------------------------
 .. autosummary::
diff --git a/harmonica/visualization/__init__.py b/harmonica/visualization/__init__.py
index b348c0fa7..82e740392 100644
--- a/harmonica/visualization/__init__.py
+++ b/harmonica/visualization/__init__.py
@@ -4,3 +4,4 @@
 #
 # This code is part of the Fatiando a Terra project (https://www.fatiando.org)
 #
+from .prism import prisms_to_pyvista
diff --git a/harmonica/visualization/prism.py b/harmonica/visualization/prism.py
index 99eb6430f..dc446daaf 100644
--- a/harmonica/visualization/prism.py
+++ b/harmonica/visualization/prism.py
@@ -19,7 +19,10 @@
 
 def prisms_to_pyvista(prisms, properties=None):
     """
-    Create a ``pyvista.UnstructuredGrid`` from a set of prisms
+    Create a :class:`pyvista.UnstructuredGrid` from a set of prisms
+
+    Builds a :class:`pyvista.UnstructuredGrid` out of a set of prisms that
+    could be used to plot a 3D representation through :mod:`pyvista`.
 
     Parameters
     ----------

From 9a6f7c00a483ea3d14c2d3b6b2eb5f17b720a736 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Fri, 22 Apr 2022 15:14:37 -0300
Subject: [PATCH 42/46] Rename prisms_to_pyvista to prism_to_pyvista

---
 doc/api/index.rst                      |  2 +-
 harmonica/forward/prism_layer.py       |  4 ++--
 harmonica/tests/test_visualizations.py | 30 +++++++++++++-------------
 harmonica/visualization/__init__.py    |  2 +-
 harmonica/visualization/prism.py       |  4 ++--
 5 files changed, 21 insertions(+), 21 deletions(-)

diff --git a/doc/api/index.rst b/doc/api/index.rst
index 9a2374511..44f5ab97c 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -61,7 +61,7 @@ Visualization
 .. autosummary::
    :toctree: generated/
 
-    visualization.prisms_to_pyvista
+    visualization.prism_to_pyvista
 
 Synthetic models and surveys
 ----------------------------
diff --git a/harmonica/forward/prism_layer.py b/harmonica/forward/prism_layer.py
index 660648e82..a23a21cb1 100644
--- a/harmonica/forward/prism_layer.py
+++ b/harmonica/forward/prism_layer.py
@@ -13,7 +13,7 @@
 import verde as vd
 import xarray as xr
 
-from ..visualization.prism import prisms_to_pyvista
+from ..visualization.prism import prism_to_pyvista
 from .prism import prism_gravity
 
 
@@ -465,4 +465,4 @@ def to_pyvista(self):
                 data_var: np.asarray(self._obj[data_var]).ravel()
                 for data_var in self._obj.data_vars
             }
-        return prisms_to_pyvista(prisms, properties=properties)
+        return prism_to_pyvista(prisms, properties=properties)
diff --git a/harmonica/tests/test_visualizations.py b/harmonica/tests/test_visualizations.py
index 6bc3b071d..cad20a320 100644
--- a/harmonica/tests/test_visualizations.py
+++ b/harmonica/tests/test_visualizations.py
@@ -14,7 +14,7 @@
 import pytest
 import xarray as xr
 
-from ..visualization.prism import prisms_to_pyvista
+from ..visualization.prism import prism_to_pyvista
 
 try:
     import pyvista
@@ -23,13 +23,13 @@
 
 
 @pytest.mark.skipif(pyvista is not None, reason="pyvista must be missing")
-def test_prisms_to_pyvista_missing_pyvista():
+def test_prism_to_pyvista_missing_pyvista():
     """
-    Check error raise after calling prisms_to_pyvista when pyvista is missing
+    Check error raise after calling prism_to_pyvista when pyvista is missing
     """
     prism = [0, 1, 0, 1, 0, 1]
     with pytest.raises(ImportError) as exception:
-        prisms_to_pyvista(prism)
+        prism_to_pyvista(prism)
     assert "'pyvista'" in str(exception.value)
 
 
@@ -69,11 +69,11 @@ def fixture_susceptibility(request):
 
 
 @pytest.mark.skipif(pyvista is None, reason="requires pyvista")
-def test_prisms_to_pyvista(prisms, density):
+def test_prism_to_pyvista(prisms, density):
     """
-    Test output of prisms_to_pyvista
+    Test output of prism_to_pyvista
     """
-    pv_grid = prisms_to_pyvista(prisms, properties={"density": density})
+    pv_grid = prism_to_pyvista(prisms, properties={"density": density})
     assert pv_grid.n_cells == 4
     assert pv_grid.n_points == 32
     # Check coordinates of prisms
@@ -87,16 +87,16 @@ def test_prisms_to_pyvista(prisms, density):
 
 @pytest.mark.skipif(pyvista is None, reason="requires pyvista")
 @pytest.mark.parametrize("n_props", [0, 1, 2])
-def test_prisms_to_pyvista_properties(n_props, prisms, density, susceptibility):
+def test_prism_to_pyvista_properties(n_props, prisms, density, susceptibility):
     """
-    Test prisms_to_pyvista for different number of properties
+    Test prism_to_pyvista for different number of properties
     """
     properties = None
     if n_props == 1:
         properties = {"density": density}
     elif n_props == 2:
         properties = {"density": density, "susceptibility": susceptibility}
-    pv_grid = prisms_to_pyvista(prisms, properties=properties)
+    pv_grid = prism_to_pyvista(prisms, properties=properties)
     # Check if the grid has the right properties
     if properties is None:
         assert pv_grid.n_arrays == 0
@@ -109,20 +109,20 @@ def test_prisms_to_pyvista_properties(n_props, prisms, density, susceptibility):
 
 
 @pytest.mark.skipif(pyvista is None, reason="requires pyvista")
-def test_prisms_to_pyvista_error_2d_property(prisms, density):
+def test_prism_to_pyvista_error_2d_property(prisms, density):
     """
-    Test if prisms_to_pyvista raises error on property as 2D array
+    Test if prism_to_pyvista raises error on property as 2D array
     """
     density_2d = np.array(density).reshape((2, 2))
     with pytest.raises(ValueError):
-        prisms_to_pyvista(prisms, properties={"density": density_2d})
+        prism_to_pyvista(prisms, properties={"density": density_2d})
 
 
 @patch("harmonica.visualization.prism.pyvista", None)
 def test_prisms_pyvista_missing_error(prisms, density):
     """
-    Check if prisms_to_pyvista raises error if pyvista is missing
+    Check if prism_to_pyvista raises error if pyvista is missing
     """
     # Check if error is raised
     with pytest.raises(ImportError):
-        prisms_to_pyvista(prisms, properties={"density": density})
+        prism_to_pyvista(prisms, properties={"density": density})
diff --git a/harmonica/visualization/__init__.py b/harmonica/visualization/__init__.py
index 82e740392..d08adbdde 100644
--- a/harmonica/visualization/__init__.py
+++ b/harmonica/visualization/__init__.py
@@ -4,4 +4,4 @@
 #
 # This code is part of the Fatiando a Terra project (https://www.fatiando.org)
 #
-from .prism import prisms_to_pyvista
+from .prism import prism_to_pyvista
diff --git a/harmonica/visualization/prism.py b/harmonica/visualization/prism.py
index dc446daaf..3683f4efa 100644
--- a/harmonica/visualization/prism.py
+++ b/harmonica/visualization/prism.py
@@ -17,9 +17,9 @@
     import vtk
 
 
-def prisms_to_pyvista(prisms, properties=None):
+def prism_to_pyvista(prisms, properties=None):
     """
-    Create a :class:`pyvista.UnstructuredGrid` from a set of prisms
+    Create a :class:`pyvista.UnstructuredGrid` out of prisms
 
     Builds a :class:`pyvista.UnstructuredGrid` out of a set of prisms that
     could be used to plot a 3D representation through :mod:`pyvista`.

From 4260a2553e4cea190af3d05a931d9a4a76adb221 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Fri, 22 Apr 2022 15:56:36 -0300
Subject: [PATCH 43/46] Add example in prism_to_pyvista docstring

The example uses the plot-pyvista directive to run the code when
documentation builds and show the generated plot.
---
 doc/conf.py                      |  2 +-
 harmonica/visualization/prism.py | 27 +++++++++++++++++++++++++++
 2 files changed, 28 insertions(+), 1 deletion(-)

diff --git a/doc/conf.py b/doc/conf.py
index 00117c152..47a82c92a 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -31,9 +31,9 @@
     "sphinx.ext.viewcode",
     "sphinx.ext.extlinks",
     "sphinx.ext.intersphinx",
-    "matplotlib.sphinxext.plot_directive",
     "sphinx.ext.napoleon",
     "sphinx_gallery.gen_gallery",
+    "pyvista.ext.plot_directive",
 ]
 
 # Configuration to include links to other project docs when referencing
diff --git a/harmonica/visualization/prism.py b/harmonica/visualization/prism.py
index 3683f4efa..17b6150f7 100644
--- a/harmonica/visualization/prism.py
+++ b/harmonica/visualization/prism.py
@@ -42,6 +42,33 @@ def prism_to_pyvista(prisms, properties=None):
     pv_grid : :class:`pyvista.UnstructuredGrid`
         :class:`pyvista.UnstructuredGrid` that represents the prisms with their
         properties (if any).
+
+    Examples
+    --------
+
+    .. pyvista-plot::
+
+       Define a set of prisms and their densities:
+
+       >>> prisms = [
+       ...     [0, 4, 0, 5, -10, 0],
+       ...     [0, 4, 7, 9, -12, -3],
+       ...     [6, 9, 2, 6, -7, 3],
+       ... ]
+       >>> densities = [2900, 3000, 2670]
+
+       Generate a :class:`pyvista.UnstructuredGrid` out of them:
+
+       >>> import harmonica as hm
+       >>> pv_grid = hm.visualization.prism_to_pyvista(
+       ...     prisms, properties={"density": densities}
+       ... )
+       >>> pv_grid # doctest: +SKIP
+
+       Plot it using :mod:`pyvista`:
+
+       >>> pv_grid.plot() # doctest: +SKIP
+
     """
     # Check if pyvista are installed
     if pyvista is None:

From 58214e1079ccd4a841520281eadc5b0a5e585fe6 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Fri, 22 Apr 2022 16:00:38 -0300
Subject: [PATCH 44/46] Add pyvista-plot directive to flake8 configuration

---
 setup.cfg | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/setup.cfg b/setup.cfg
index e3132b3da..299de10b4 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -87,6 +87,9 @@ rst-roles =
     mod,
     meth,
     ref,
+# Add some directives used in our docstrings
+rst-directives =
+    pyvista-plot,
 # Ignore "Unknown target name" raised on citations
 extend-ignore = RST306
 

From b6dbb9c764c6ab3915e3620e725771c708fd1bba Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Fri, 22 Apr 2022 16:44:25 -0300
Subject: [PATCH 45/46] Fix typo

---
 harmonica/visualization/prism.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/harmonica/visualization/prism.py b/harmonica/visualization/prism.py
index 17b6150f7..79644ca30 100644
--- a/harmonica/visualization/prism.py
+++ b/harmonica/visualization/prism.py
@@ -101,7 +101,7 @@ def prism_to_pyvista(prisms, properties=None):
 
 def _prisms_boundaries_to_vertices(prisms):
     """
-    Converts prisms boundaries to set of vertices for each prism
+    Converts prisms boundaries to sets of vertices for each prism
 
     The vertices for each prism will be in the following order
 

From dbcf92c123a0fe1817c9021ba9995c9fa5a96ea9 Mon Sep 17 00:00:00 2001
From: Santiago Soler <santiago.r.soler@gmail.com>
Date: Fri, 22 Apr 2022 16:46:39 -0300
Subject: [PATCH 46/46] Improve docstring

---
 harmonica/visualization/prism.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/harmonica/visualization/prism.py b/harmonica/visualization/prism.py
index 79644ca30..1553bdd8a 100644
--- a/harmonica/visualization/prism.py
+++ b/harmonica/visualization/prism.py
@@ -32,7 +32,7 @@ def prism_to_pyvista(prisms, properties=None):
         ``west``, ``east``, ``south``, ``north``, ``bottom``, ``top``.
     properties : dict or None (optional)
         Dictionary with the physical properties of the prisms.
-        Each key should be a string and their value as a 1D array.
+        Each key should be a string and its corresponding value a 1D array.
         If None, no property will be added to the
         :class:`pyvista.UnstructuredGrid`.
         Default to None.