diff --git a/examples/forward/tesseroid_layer.py b/examples/forward/tesseroid_layer.py
new file mode 100644
index 000000000..ad9977cf1
--- /dev/null
+++ b/examples/forward/tesseroid_layer.py
@@ -0,0 +1,67 @@
+# 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)
+#
+"""
+Layer of tesseroids
+===================
+"""
+import boule as bl
+import ensaio
+import numpy as np
+import pygmt
+import verde as vd
+import xarray as xr
+
+import harmonica as hm
+
+fname = ensaio.fetch_earth_topography(version=1)
+topo = xr.load_dataarray(fname)
+
+region = (-78, -53, -57, -20)
+topo = topo.sel(latitude=slice(*region[2:]), longitude=slice(*region[:2]))
+
+ellipsoid = bl.WGS84
+
+longitude, latitude = np.meshgrid(topo.longitude, topo.latitude)
+reference = ellipsoid.geocentric_radius(latitude)
+surface = topo + reference
+density = xr.where(topo > 0, 2670.0, 1040.0 - 2670.0)
+
+tesseroids = hm.tesseroid_layer(
+    coordinates=(topo.longitude, topo.latitude),
+    surface=surface,
+    reference=reference,
+    properties={"density": density},
+)
+
+# Create a regular grid of computation points located at 10km above reference
+grid_longitude, grid_latitude = vd.grid_coordinates(region=region, spacing=0.5)
+grid_radius = ellipsoid.geocentric_radius(grid_latitude) + 10e3
+grid_coords = (grid_longitude, grid_latitude, grid_radius)
+
+# Compute gravity field of tesseroids on a regular grid of observation points
+gravity = tesseroids.tesseroid_layer.gravity(grid_coords, field="g_z")
+gravity = vd.make_xarray_grid(
+    grid_coords,
+    gravity,
+    data_names="g_z",
+    dims=("latitude", "longitude"),
+    extra_coords_names="radius",
+)
+
+# Plot gravity field
+fig = pygmt.Figure()
+maxabs = vd.maxabs(gravity.g_z)
+pygmt.makecpt(cmap="polar", series=(-maxabs, maxabs))
+fig.grdimage(
+    gravity.g_z,
+    projection="M15c",
+    nan_transparent=True,
+)
+fig.basemap(frame=True)
+fig.colorbar(frame='af+l"Gravity [mGal]"', position="JCR")
+fig.coast(shorelines="0.5p,black", borders=["1/0.5p,black"])
+fig.show()
diff --git a/harmonica/__init__.py b/harmonica/__init__.py
index 7c2d4eb04..f745da1c7 100644
--- a/harmonica/__init__.py
+++ b/harmonica/__init__.py
@@ -15,6 +15,7 @@
 from .forward.prism import prism_gravity
 from .forward.prism_layer import DatasetAccessorPrismLayer, prism_layer
 from .forward.tesseroid import tesseroid_gravity
+from .forward.tesseroid_layer import DatasetAccessorTesseroidLayer, tesseroid_layer
 from .gravity_corrections import bouguer_correction
 from .io import load_icgem_gdf
 from .isostasy import isostasy_airy, isostatic_moho_airy
diff --git a/harmonica/forward/tesseroid_layer.py b/harmonica/forward/tesseroid_layer.py
new file mode 100644
index 000000000..6c121f07b
--- /dev/null
+++ b/harmonica/forward/tesseroid_layer.py
@@ -0,0 +1,411 @@
+# 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)
+#
+"""
+Define a layer of tesseroids
+"""
+import warnings
+
+import numpy as np
+import verde as vd
+import xarray as xr
+
+from .tesseroid import tesseroid_gravity
+
+
+def tesseroid_layer(coordinates, surface, reference, properties=None):
+    """
+    Create a layer of tesseroids of equal size
+
+    Parameters
+    ----------
+    coordinates : tuple
+        List containing the coordinates of the centers of the tesseroids in
+        spherical coordinates in the following order ``longitude`` and
+        ``latitude``.
+    surface : 2d-array
+        Array used to create the uppermost boundary of the tesserois layer. All
+        radii should be in meters. On every point where ``surface`` is below
+        ``reference``, the ``surface`` value will be used to set the ``bottom``
+        boundary of that tesseroid, while the ``reference`` value will be used
+        to set the ``top`` boundary of the tesseroid.
+    reference : 2d-array or float
+        Reference surface used to create the lowermost boundary of the
+        tesseroids layer. It can be either a plane or an irregular surface
+        passed as 2d array. Radii must be in meters.
+    properties : dict or None
+        Dictionary containing the physical properties of the tesseroids. The
+        keys must be strings that will be used to name the corresponding
+        ``data_var`` inside the :class:`xarray.Dataset`, while the values must
+        be 2d-arrays. All physical properties must be passed in SI units. If
+        None, no ``data_var`` will be added to the :class:`xarray.Dataset`.
+        Default is None.
+
+    Returns
+    -------
+    dataset : :class:`xarray.Dataset`
+        Dataset containing the coordinates of the center of each tesseroid, the
+        height of its top and bottom boundaries ans its corresponding physical
+        properties.
+
+    See also
+    --------
+    harmonica.DatasetAccessorsTesseroidLayer
+
+    """
+    dims = ("latitude", "longitude")
+    # Initialize data and data_names as None
+    data, data_names = None, None
+    # If properties were passed, then replace data_names and data for its keys
+    # and values, respectively
+    if properties:
+        data_names = tuple(p for p in properties.keys())
+        data = tuple(np.asarray(p) for p in properties.values())
+    # Create xr.Dataset for tesseroids
+    tesseroids = vd.make_xarray_grid(
+        coordinates, data=data, data_names=data_names, dims=dims
+    )
+    _check_regular_grid(tesseroids.longitude.values, tesseroids.latitude.values)
+    # Check if tesseroid boundaries are overlapped
+    _check_overlap(tesseroids.longitude.values)
+    # Append some attributes to the xr.Dataset
+    attrs = {
+        "longitude_units": "degrees",
+        "latitude_units": "degrees",
+        "radius_units": "meters",
+        "properties_units": "SI",
+    }
+    tesseroids.attrs = attrs
+    # Create the top and bottom coordinates of the prisms
+    tesseroids.tesseroid_layer.update_top_bottom(surface, reference)
+    return tesseroids
+
+
+def _check_regular_grid(longitude, latitude):
+    """
+    Check if the longitude and latitude coordinates define a regular grid
+
+    .. note:
+
+        This function should live inside Verde in the future
+
+    """
+    if not np.allclose(longitude[1] - longitude[0], longitude[1:] - longitude[:-1]):
+        raise ValueError("Passed longitude coordinates are note evenly spaced.")
+    if not np.allclose(latitude[1] - latitude[0], latitude[1:] - latitude[:-1]):
+        raise ValueError("Passed latitude coordinates are note evenly spaced.")
+
+
+def _check_overlap(longitude):
+    """
+    Check if the prisms boundaries are overlapped
+    """
+    spacing = longitude[1] - longitude[0]
+    if longitude.max() - longitude.min() >= 360 - spacing:
+        raise ValueError(
+            "Found invalid longitude coordinates that would create overlapping tesseroids around the globe."
+        )
+
+
+@xr.register_dataset_accessor("tesseroid_layer")
+class DatasetAccessorTesseroidLayer:
+    """
+    Define dataset accessor for layer of tesseroids
+
+    .. warning::
+
+        This class in not intended to be initialized.
+        Use the `tesseroid_layer` accessor for accessing the methods and
+        attributes of this class.
+
+    See also
+    --------
+    harmonica.tesseroid_layer
+    """
+
+    def __init__(self, xarray_obj):
+        self._obj = xarray_obj
+
+    @property
+    def dims(self):
+        """
+        Return the dims tuple of the prism layer
+
+        The tuple follows the xarray order: ``"latitude"``, ``"longitude"``.
+        """
+        return ("latitude", "longitude")
+
+    @property
+    def spacing(self):
+        """
+        Spacing between center of tesseroids
+
+        Returns
+        -------
+        s_latitude : float
+            Spacing between center of the tesseroids on the latitude direction.
+        s_longitude : float
+            Spacing between center of the tesseroids on the longitude
+            direction.
+        """
+        latitude, longitude = self._obj.latitude.values, self._obj.longitude.values
+        _check_regular_grid(longitude, latitude)
+        s_latitude, s_longitude = latitude[1] - latitude[0], longitude[1] - longitude[0]
+        return s_latitude, s_longitude
+
+    @property
+    def size(self):
+        """
+        Return the total number of tesseroids on the layer
+
+        Returns
+        -------
+        size :  int
+            Total number of tesseroids in the layer.
+        """
+        return self._obj.latitude.size * self._obj.longitude.size
+
+    @property
+    def shape(self):
+        """
+        Return the number of tesseroids on each directions
+
+        Returns
+        -------
+        n_latitude : int
+            Number of tesseroids on the latitude direction.
+        n_longitude : int
+            Number of tesserods on the longitude direction.
+        """
+        return (self._obj.latitude.size, self._obj.longitude.size)
+
+    @property
+    def boundaries(self):
+        """
+        Boundaries of the layer
+
+        Returns
+        -------
+        boundaries : tuple
+        Boundaries of the layer of tesseroids in the following order:
+        ``longitude_w``, ``longitude_e``, ``latitude_s``, ``latitude_n``
+        """
+        s_latitude, s_longitude = self.spacing
+        longitude_w = self._obj.longitude.values.min() - s_longitude / 2
+        longitude_e = self._obj.longitude.values.max() + s_longitude / 2
+        latitude_s = self._obj.latitude.values.min() - s_latitude / 2
+        latitude_n = self._obj.latitude.values.max() + s_latitude / 2
+        return longitude_w, longitude_e, latitude_s, latitude_n
+
+    def update_top_bottom(self, surface, reference):
+        """
+        Update top and bottom boundaries of the layer
+
+        Change the values of the ``top`` and ``bottom`` coordinates based on
+        the passed ``surface`` and ``reference``. The ``top`` and ``bottom``
+        boundaries of every tesseroid will be equal to the corresponding
+        ``surface`` and ``reference`` values, respectively, if ``surface`` is
+        above the ``reference`` on that point. Otherwise the ``top`` and
+        ``bottom`` boundaries of the tesseroid will be equal to its
+        corresponding ``reference`` and ``surface``, respectively.
+
+        Parameters
+        ----------
+        surface : 2d-array
+            Array used to create the uppermost boundary of the tesseroid layer.
+            All heights should be in meters. On every point where ``surface``
+            is below ``reference``, the ``surface`` value will be used to set
+            the ``bottom`` boundary of that tesseroid, while the ``reference``
+            value will be used to set the ``top`` boundary of the tesseroid.
+
+        reference : 2d-array or float
+            Reference surface used to create the lowermost boundary of the
+            tesseroid layer. It can be either a plane or an irregular surface
+            passed as 2d array. Height(s) must be in meters.
+        """
+        surface, reference = np.asarray(surface), np.asarray(reference)
+        if surface.shape != self.shape:
+            raise ValueError(
+                f"Invalid surface array with shape '{surface.shape}'. "
+                + "Its shape should be compatible with the coordinates "
+                + "of the layer of tesseroids."
+            )
+        if reference.ndim != 0:
+            if reference.shape != self.shape:
+                raise ValueError(
+                    f"Invalid reference array with shape '{reference.shape}'. "
+                    + "Its shape should be compatible with the coordinates "
+                    + "of the layer of tesseroids."
+                )
+        else:
+            reference = reference * np.ones(self.shape)
+        top = surface.copy()
+        bottom = reference.copy()
+        reverse = surface < reference
+        top[reverse] = reference[reverse]
+        bottom[reverse] = surface[reverse]
+        self._obj.coords["top"] = (self.dims, top)
+        self._obj.coords["bottom"] = (self.dims, bottom)
+
+    def gravity(self, coordinates, field, density_name="density", **kwargs):
+        """
+        Computes the gravity generated by the layer of tesseroids
+
+        Parameters
+        ----------
+        coordinates : list or 1d-array
+            List of array containing ``latitude``, ``longitude`` and ``upward``
+            of the computation points defined on a spherical coordinates
+            system. ``upward`` coordinate should be in meters.
+        field : str
+            Gravitational field that wants to be computed.
+            The variable fields are:
+            - Gravitational potential: ``potential``
+            - Downward acceleration: ``g_z``
+        density_name : str (optional)
+            Name of the property layer (or ``data_var`` of the
+            :class:`xarray.Dataset`) that will be used for the density of each
+            tesseroid in the layer. Default to ``"density"``
+
+        Returns
+        -------
+        result : array
+            Gravitational field generated by the tesseroid on the computation
+            point in mGal
+
+        See also
+        --------
+        harmonica.tesseroid_gravity
+        """
+        # Get boundaries and density of the tesseroids
+        boundaries = self._to_tesseroids()
+        density = self._obj[density_name].values
+        # Get the mask for selecting  only the tesseroid whose top boundary,
+        # bottom boundary and density have no nans
+        mask = self._get_nonans_mask(property_name=density_name)
+        # Select only the boundaries and density elements for masked tesseroid
+        boundaries = boundaries[mask.ravel()]
+        density = density[mask]
+        # Return gravity field of tesserids
+        return tesseroid_gravity(
+            coordinates,
+            tesseroids=boundaries,
+            density=density,
+            field=field,
+            **kwargs,
+        )
+
+    def _get_nonans_mask(self, property_name=None):
+        """
+        Build a mask for tesseroid with no nans on top, bottom or a property
+
+        Parameters
+        ----------
+        mask : 2d-array
+            Array of bools that can be used as a mask for selecting tesseroids
+            with no nans on top boundaries, bottom boundaries ans the passed
+            property.
+        """
+        # Mask the tesseroid that contains no nans on top and bottom boundaries
+        mask = np.logical_and(
+            np.logical_not(np.isnan(self._obj.top.values)),
+            np.logical_not(np.isnan(self._obj.bottom.values)),
+        )
+        # Mask the tesseroids that contains nans on the selected property
+        if property_name is not None:
+            mask_property = np.logical_not(np.isnan(self._obj[property_name].values))
+            # Warn if a nan is found within the masked property
+            if not mask_property[mask].all():
+                warnings.warn(
+                    'Found missing values in "{}" property '.format(property_name)
+                    + "of the tesseroid layer. "
+                    + "The tesseroids with nan as "
+                    + '"{}" will be ignored.'.format(property_name)
+                )
+            mask = np.logical_and(mask, mask_property)
+        return mask
+
+    def _to_tesseroids(self):
+        """
+        Return the boundaries of each tesseroid of the layer
+
+        Returns
+        -------
+        tesseroids : 2d-array
+            Array containing the boundaries of each tesseroid of the layer.
+            Each row contains the boundaries of each tesseroid in the following
+            order: ``longitude_w``, ``longitude_e``, ``latitude_s``,
+            ``latitude_n``, ``bottom``, ``top``.
+        """
+        longitude, latitude = np.meshgrid(
+            self._obj.longitude.values, self._obj.latitude.values
+        )
+        (
+            longitude_w,
+            longitude_e,
+            latitude_s,
+            latitude_n,
+        ) = self._get_tesseroid_horizontal_boundaries(
+            longitude.ravel(), latitude.ravel()
+        )
+        bottom = self._obj.bottom.values.ravel()
+        top = self._obj.top.values.ravel()
+        tesseroids = np.vstack(
+            (longitude_w, longitude_e, latitude_s, latitude_n, bottom, top)
+        ).T
+        return tesseroids
+
+    def _get_tesseroid_horizontal_boundaries(self, longitude, latitude):
+        """
+        Compute the horizontal boundaries of the tesseroid
+
+        Parameters
+        ----------
+        latitude: float or array
+            Longitude coordinates of the center of the tesseroid
+        longitude : float or array
+            Longitude coordinates of the center of the tesseroid
+        """
+        spacing = self.spacing
+        longitude_w = longitude - spacing[1] / 2
+        longitude_e = longitude + spacing[1] / 2
+        latitude_s = latitude - spacing[0] / 2
+        latitude_n = latitude + spacing[0] / 2
+        return longitude_w, longitude_e, latitude_s, latitude_n
+
+    def get_tesseroid(self, indices):
+        """
+        Return the boundaries of the chosen tesseroid
+
+        Parameters
+        ----------
+        indices : tuple
+            Indices of the desired tesseroid of the layer in the following
+            order: ``(index_northing, index_easting)``.
+
+        Returns
+        -------
+        tesseroid : tuple
+           Boundaries of the prisms in the following order:
+           ``longitude_w``, ``longitude_e``, ``latitude_s``, ``latitude_n``,
+           ``bottom``, ``top``.
+        """
+        # Get the center of the tesseroid
+        center_longitude = self._obj.longitude.values[indices[1]]
+        center_latitude = self._obj.latitude.values[indices[0]]
+        # Calculate the boundaries of the tesseroid
+        # (
+        #     longitude_w,
+        #     Longitude_e,
+        #     latitude_s,
+        #     latitude_n,
+        boundaries = self._get_tesseroid_horizontal_boundaries(
+            center_longitude, center_latitude
+        )
+        bottom = self._obj.bottom.values[indices]
+        top = self._obj.top.values[indices]
+        # return longitude_w, longitude_e, latitude_s, latitude_n, bottom, top
+        return boundaries[0], boundaries[1], boundaries[2], boundaries[3], bottom, top
diff --git a/harmonica/tests/test_tesseroid_layer.py b/harmonica/tests/test_tesseroid_layer.py
new file mode 100644
index 000000000..75ba11e38
--- /dev/null
+++ b/harmonica/tests/test_tesseroid_layer.py
@@ -0,0 +1,467 @@
+# 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 tesseroids layer
+"""
+#
+import warnings
+
+import boule
+import numpy as np
+import numpy.testing as npt
+import pytest
+import verde as vd
+import xarray as xr
+
+from .. import tesseroid_gravity, tesseroid_layer
+
+
+@pytest.fixture
+def mean_earth_radius():
+    """
+    Return mean earth radius given by WGS84
+    """
+    return boule.WGS84.mean_radius
+
+
+@pytest.fixture(params=("numpy", "reference-as-array", "xarray"))
+def dummy_layer(mean_earth_radius, request):
+    """
+    Generate dummy array for defining tesseroid layers
+    """
+    latitude = np.linspace(-10, 10, 6)
+    longitude = np.linspace(-10, 10, 5)
+    shape = (latitude.size, longitude.size)
+    surface = mean_earth_radius * np.ones(shape) + 1e3
+    reference = mean_earth_radius
+    if request.param == "reference-as-array":
+        reference *= np.ones(shape)
+    density = 2670 * np.ones(shape)
+    if request.param == "xarray":
+        latitude = xr.DataArray(latitude, dims=("latitude",))
+        longitude = xr.DataArray(longitude, dims=("longitude",))
+        reference, surface = xr.DataArray(reference), xr.DataArray(surface)
+        density = xr.DataArray(density)
+    return (longitude, latitude), surface, reference, density
+
+
+@pytest.fixture
+def tesseroid_layer_with_holes(dummy_layer):
+    """
+    Return a set of tesseroids with some missing elements
+
+    The tesseroids are returned as a tuple of boundaries, ready to be passed to
+    ``hm.tesseroid_gravity``.
+    They would represent the same tesseroids that the ``dummy_layer``
+    generated, but with two missing tesseroids: the ``(3, 3)`` and the
+    ``(2, 1)``.
+    """
+    (longitude, latitude), surface, reference, density = dummy_layer
+    layer = tesseroid_layer(
+        (longitude, latitude), surface, reference, properties={"density": density}
+    )
+    indices = [(3, 3), (2, 1)]
+    tesseroids = list(
+        layer.tesseroid_layer.get_tesseroid((i, j))
+        for i in range(6)
+        for j in range(5)
+        if (i, j) not in indices
+    )
+    density = list(
+        density[i, j] for i in range(6) for j in range(5) if (i, j) not in indices
+    )
+    return tesseroids, density
+
+
+@pytest.mark.parametrize(
+    "west, east",
+    [
+        (0, 480),
+        (0, 360),
+        (-180, 180),
+        (0, 360 - 18 / 2),
+    ],
+)
+def test_tesseroid_overlap_invalid_coords(west, east, mean_earth_radius):
+    """
+    Check if error is raised when longitude coordinates create overlapped
+    tesseroid
+    """
+    latitude = np.linspace(-10, 10, 6)
+    longitude = np.linspace(west, east, 21)
+    shape = (latitude.size, longitude.size)
+    surface = mean_earth_radius * np.ones(shape) + 1e3
+    reference = mean_earth_radius * np.ones(shape)
+    message = (
+        "Found invalid longitude coordinates that would create "
+        + "overlapping tesseroids around the globe."
+    )
+    with pytest.raises(
+        ValueError,
+        match=message,
+    ):
+        tesseroid_layer((longitude, latitude), surface, reference)
+
+
+@pytest.mark.parametrize(
+    "west, east",
+    [
+        (0, 360 - 18),
+        (-180, 180 - 18),
+    ],
+)
+def test_tesseroid_overlap_valid_coords(west, east, mean_earth_radius):
+    """
+    Check if tesseroid_layer works properly when tesseroids are not overlapped
+    """
+    latitude = np.linspace(-10, 10, 6)
+    longitude = np.linspace(west, east, 21)
+    shape = (latitude.size, longitude.size)
+    surface = mean_earth_radius * np.ones(shape) + 1e3
+    reference = mean_earth_radius * np.ones(shape)
+    tesseroid_layer((longitude, latitude), surface, reference)
+
+
+def test_tesseroid_layer(dummy_layer, mean_earth_radius):
+    """
+    Check if the layer of tesseroids is property constructed
+    """
+    (longitude, latitude), surface, reference, _ = dummy_layer
+    layer = tesseroid_layer((longitude, latitude), surface, reference)
+    assert "latitude" in layer.coords
+    assert "longitude" in layer.coords
+    assert "top" in layer.coords
+    assert "bottom" in layer.coords
+    npt.assert_allclose(layer.latitude, latitude)
+    npt.assert_allclose(layer.longitude, longitude)
+    npt.assert_allclose(layer.top, surface)
+    npt.assert_allclose(layer.bottom, reference)
+    # Surface below reference on a single point
+    surface[1, 1] = mean_earth_radius - 1e3  # reference is on mean_earth_radius
+    expected_top = surface.copy()
+    expected_bottom = mean_earth_radius * np.ones_like(surface)
+    expected_top[1, 1], expected_bottom[1, 1] = mean_earth_radius, surface[1, 1]
+    layer = tesseroid_layer((longitude, latitude), surface, reference)
+    assert "latitude" in layer.coords
+    assert "longitude" in layer.coords
+    assert "top" in layer.coords
+    assert "bottom" in layer.coords
+    npt.assert_allclose(layer.latitude, latitude)
+    npt.assert_allclose(layer.longitude, longitude)
+    npt.assert_allclose(layer.top, expected_top)
+    npt.assert_allclose(layer.bottom, expected_bottom)
+
+
+def test_tesseroid_layer_invalid_surface_reference(dummy_layer):
+    """
+    Check if invalid surface and/or reference are caught
+    """
+    coordinates, surface, reference, _ = dummy_layer
+    # Surface with wrong shape
+    surface_invalid = np.arange(20, dtype=float)
+    with pytest.raises(ValueError, match="Invalid surface array with shape"):
+        tesseroid_layer(coordinates, surface_invalid, reference)
+    # Reference with wrong shape
+    reference_invalid = np.zeros(20)
+    with pytest.raises(ValueError, match="Invalid reference array with shape"):
+        tesseroid_layer(coordinates, surface, reference_invalid)
+
+
+def test_tesseroid_leyer_properties(dummy_layer):
+    """
+    Check passing physical properties to the tesseroid layer
+    """
+    coordinates, surface, reference, density = dummy_layer
+    suceptibility = 0 * density + 1e-3
+    layer = tesseroid_layer(
+        coordinates,
+        surface,
+        reference,
+        properties={"density": density, "suceptibility": suceptibility},
+    )
+    npt.assert_allclose(layer.density, density)
+    npt.assert_allclose(layer.suceptibility, suceptibility)
+
+
+def test_tesseroid_layer_no_regular_grid(
+    dummy_layer,
+):
+    """
+    Check if error is raised if the latitude or longitude are not regular
+    """
+
+    (longitude, latitude), surface, reference, _ = dummy_layer
+    # Longitude as not evenly spaced set of coordinates
+    longitude_invalid = longitude.copy()
+    longitude_invalid[3] = -22
+    with pytest.raises(ValueError):
+        tesseroid_layer(
+            (longitude_invalid, latitude),
+            surface,
+            reference,
+        )
+    # Latitude as not evenly spaced set of coordinates
+    latitude_invalid = latitude.copy()
+    latitude_invalid[3] = -22
+    with pytest.raises(ValueError):
+        tesseroid_layer(
+            (longitude, latitude_invalid),
+            surface,
+            reference,
+        )
+
+
+def test_tesseroi_layer_attibutes():
+    """
+    Check attributes of the DatasetAccessorTesseroidLayer class
+    """
+    latitude = np.linspace(-10, 10, 6)
+    longitude = np.linspace(-10, 10, 5)
+    shape = (latitude.size, longitude.size)
+    ellipsoid = boule.WGS84
+    surface = ellipsoid.mean_radius * np.ones(shape)
+    reference = (surface - 1e3) * np.ones(shape)
+    layer = tesseroid_layer((longitude, latitude), surface, reference)
+    assert layer.tesseroid_layer.dims == ("latitude", "longitude")
+    assert layer.tesseroid_layer.spacing == (4, 5)
+    assert layer.tesseroid_layer.boundaries == (
+        longitude[0] - 2.5,
+        longitude[-1] + 2.5,
+        latitude[0] - 2,
+        latitude[-1] + 2,
+    )
+    assert layer.tesseroid_layer.size == 30
+    assert layer.tesseroid_layer.shape == (6, 5)
+
+
+def test_tesseroid_layer_to_tesseroid():
+    """
+    Check the _to_tesseroid() method
+    """
+    latitude = np.linspace(-1, 1, 2)
+    longitude = np.linspace(-2, 2, 2)
+    shape = (latitude.size, longitude.size)
+    ellipsoid = boule.WGS84
+    surface = ellipsoid.mean_radius * np.ones(shape)
+    reference = (surface - 1e3) * np.ones(shape)
+    layer = tesseroid_layer((longitude, latitude), surface, reference)
+    expected_tesseroids = [
+        [-4.0, 0.0, -2.0, 0.0, ellipsoid.mean_radius - 1e3, ellipsoid.mean_radius],
+        [0.0, 4.0, -2.0, 0.0, ellipsoid.mean_radius - 1e3, ellipsoid.mean_radius],
+        [-4.0, 0.0, 0.0, 2.0, ellipsoid.mean_radius - 1e3, ellipsoid.mean_radius],
+        [0.0, 4.0, 0.0, 2.0, ellipsoid.mean_radius - 1e3, ellipsoid.mean_radius],
+    ]
+    npt.assert_allclose(expected_tesseroids, layer.tesseroid_layer._to_tesseroids())
+
+
+def test_tesseroid_layer_get_tesseroid_by_index():
+    """
+    Check if the right tesseroid is returned after index
+    """
+    latitude = np.linspace(-1, 1, 2)
+    longitude = np.linspace(-2, 2, 2)
+    shape = (latitude.size, longitude.size)
+    ellipsoid = boule.WGS84
+    surface = ellipsoid.mean_radius * np.ones(shape)
+    reference = (surface - 1e3) * np.ones(shape)
+    layer = tesseroid_layer((longitude, latitude), surface, reference)
+    expected_tesseroids = [
+        [
+            [-4.0, 0.0, -2.0, 0.0, ellipsoid.mean_radius - 1e3, ellipsoid.mean_radius],
+            [0.0, 4.0, -2.0, 0.0, ellipsoid.mean_radius - 1e3, ellipsoid.mean_radius],
+        ],
+        [
+            [-4.0, 0.0, 0.0, 2.0, ellipsoid.mean_radius - 1e3, ellipsoid.mean_radius],
+            [0.0, 4.0, 0.0, 2.0, ellipsoid.mean_radius - 1e3, ellipsoid.mean_radius],
+        ],
+    ]
+    print(layer)
+    for i in range(2):
+        for j in range(2):
+            print(i, j)
+            print(
+                layer.tesseroid_layer.get_tesseroid((i, j)), expected_tesseroids[i][j]
+            )
+            npt.assert_allclose(
+                layer.tesseroid_layer.get_tesseroid((i, j)), expected_tesseroids[i][j]
+            )
+
+
+def test_nonans_tesseroid_mask(dummy_layer):
+    """
+    Check if the mask for nonans tesseroid is correctly created
+    """
+    (longitude, latitude), surface, reference, _ = dummy_layer
+    shape = (latitude.size, longitude.size)
+    # No nan in top or bottom
+    layer = tesseroid_layer((longitude, latitude), surface, reference)
+    expected_mask = np.ones(shape, dtype=bool)
+    mask = layer.tesseroid_layer._get_nonans_mask()
+    npt.assert_allclose(mask, expected_mask)
+    # Nans in the top only
+    layer = tesseroid_layer((longitude, latitude), surface, reference)
+    expected_mask = np.ones(shape, dtype=bool)
+    for index in ((2, 1), (3, 2)):
+        layer.top[index] = np.nan
+        expected_mask[index] = False
+    mask = layer.tesseroid_layer._get_nonans_mask()
+    npt.assert_allclose(mask, expected_mask)
+    # Nans in the bottom only
+    layer = tesseroid_layer((longitude, latitude), surface, reference)
+    expected_mask = np.ones(shape, dtype=bool)
+    for index in ((2, 1), (3, 2)):
+        layer.bottom[index] = np.nan
+        expected_mask[index] = False
+    mask = layer.tesseroid_layer._get_nonans_mask()
+    npt.assert_allclose(mask, expected_mask)
+    # Nans in the top and bottom
+    layer = tesseroid_layer((longitude, latitude), surface, reference)
+    expected_mask = np.ones(shape, dtype=bool)
+    for index in ((1, 2), (2, 3)):
+        layer.top[index] = np.nan
+        expected_mask[index] = False
+    for index in ((1, 2), (2, 1), (3, 2)):
+        layer.bottom[index] = np.nan
+        expected_mask[index] = False
+    mask = layer.tesseroid_layer._get_nonans_mask()
+    npt.assert_allclose(mask, expected_mask)
+
+
+def test_nonans_tesseroid_mask_property(
+    dummy_layer,
+):
+    """
+    Check if the method masks the property and raises a warning
+    """
+    coordinates, surface, reference, density = dummy_layer
+    shape = (coordinates[1].size, coordinates[0].size)
+    # Nans in top and property (on the same tesseroid)
+    expected_mask = np.ones_like(surface, dtype=bool)
+    indices = ((1, 2), (2, 3))
+    # Set some elements of surface and density as nans
+    for index in indices:
+        surface[index] = np.nan
+        density[index] = np.nan
+        expected_mask[index] = False
+    layer = tesseroid_layer(
+        coordinates, surface, reference, properties={"density": density}
+    )
+    # Check if no warning is raised
+    with warnings.catch_warnings(record=True) as warn:
+        mask = layer.tesseroid_layer._get_nonans_mask(property_name="density")
+        assert len(warn) == 0
+    npt.assert_allclose(mask, expected_mask)
+    # Nans in top and property (not precisely on the same tesseroid)
+    surface = np.arange(30, dtype=float).reshape(shape)
+    density = 2670 * np.ones_like(surface)
+    expected_mask = np.ones_like(surface, dtype=bool)
+    # Set some elements of surface as nans
+    indices = ((1, 2), (2, 3))
+    for index in indices:
+        surface[index] = np.nan
+        expected_mask[index] = False
+    # Set a different set of elements of density as nans
+    indices = ((2, 2), (0, 1))
+    for index in indices:
+        density[index] = np.nan
+        expected_mask[index] = False
+    layer = tesseroid_layer(
+        coordinates, surface, reference, properties={"density": density}
+    )
+    # Check if warning is raised
+    with warnings.catch_warnings(record=True) as warn:
+        mask = layer.tesseroid_layer._get_nonans_mask(property_name="density")
+        assert len(warn) == 1
+        assert issubclass(warn[-1].category, UserWarning)
+    npt.assert_allclose(mask, expected_mask)
+
+
+@pytest.mark.use_numba
+@pytest.mark.parametrize("field", ["potential", "g_z"])
+def test_tesseroid_layer_gravity(field, dummy_layer):
+    """
+    Check if gravity method works as expected
+    """
+    (longitude, latitude), surface, reference, density = dummy_layer
+    # Create a regular grid of computation points located at 10km above surface
+    grid_coords = vd.grid_coordinates(
+        (-10, 10, -10, 10), spacing=7, extra_coords=(surface[0] + 10e3)
+    )
+    layer = tesseroid_layer(
+        (longitude, latitude), surface, reference, properties={"density": density}
+    )
+    expected_result = tesseroid_gravity(
+        grid_coords,
+        tesseroids=layer.tesseroid_layer._to_tesseroids(),
+        density=density,
+        field=field,
+    )
+    npt.assert_allclose(
+        expected_result, layer.tesseroid_layer.gravity(grid_coords, field=field)
+    )
+
+
+@pytest.mark.use_numba
+@pytest.mark.parametrize("field", ["potential", "g_z"])
+def test_tesseroid_layer_gravity_surface_nans(
+    field, dummy_layer, tesseroid_layer_with_holes
+):
+    """
+    Check if gravity method works as expected when surface has nans
+    """
+    (longitude, latitude), surface, reference, density = dummy_layer
+    grid_coords = vd.grid_coordinates(
+        (-10, 10, -10, 10), spacing=7, extra_coords=(surface[0] + 10e3)
+    )
+    # Create one layer that has nans on the surface array
+    surface_w_nans = surface.copy()
+    indices = [(3, 3), (2, 1)]
+    for index in indices:
+        surface_w_nans[index] = np.nan
+    layer = tesseroid_layer(
+        (longitude, latitude),
+        surface_w_nans,
+        reference,
+        properties={"density": density},
+    )
+    # Check if it generates the expected gravity field
+    tesseroids, rho = tesseroid_layer_with_holes
+    npt.assert_allclose(
+        layer.tesseroid_layer.gravity(grid_coords, field=field),
+        tesseroid_gravity(grid_coords, tesseroids, rho, field=field),
+    )
+
+
+@pytest.mark.use_numba
+@pytest.mark.parametrize("field", ["potential", "g_z"])
+def test_tesseroid_layer_gravity_density_nans(
+    field, dummy_layer, tesseroid_layer_with_holes
+):
+    """
+    Check if tesseroid is ignored after a nan is found in density array
+    """
+    (longitude, latitude), surface, reference, density = dummy_layer
+    grid_coords = vd.grid_coordinates(
+        (-10, 10, -10, 10), spacing=7, extra_coords=(surface[0] + 10e3)
+    )
+    # Create one layer that has nans on the density array
+    indices = [(3, 3), (2, 1)]
+    for index in indices:
+        density[index] = np.nan
+    layer = tesseroid_layer(
+        (longitude, latitude), surface, reference, properties={"density": density}
+    )
+    # Check if warning is raised after passing density with nans
+    with warnings.catch_warnings(record=True) as warn:
+        result = layer.tesseroid_layer.gravity(grid_coords, field=field)
+        assert len(warn) == 1
+    # Check if it generates the expected gravity field
+    tesseroids, rho = tesseroid_layer_with_holes
+    npt.assert_allclose(
+        result,
+        tesseroid_gravity(grid_coords, tesseroids, rho, field=field),
+    )