Repository: fatiando/harmonica (python)

Issue:
Add function to create tesseroid layers
**Description of the desired feature**

A tesseroid layer could be defined as a collection of tesseroids sorted in a regular way within a region, whose longitudinal and latitudinal boundaries touch each other leaving no empty space in between nor overlapping each other.
Tesseroids layers are very useful for gravity inversion (e.g. Moho inversion). So the function should allow us to modify the top and bottom boundaries of the tesseroid layer according to a given topography data or to a predicted radii by an inversion algorithm.

The tesseroid layer should be returned as an xarray object, what will ease the process of computing the forward model of the layer.

**Are you willing to help implement and maintain this feature?** Yes


Pull request:
Add function to create a tesseroid layer
Add function that creates a tesseroid layer, similar to the prism layer. Create
an Xarray DataAccesor that makes available some methods to deal with the prism
layer, included one to compute the gravitational field generated by the whole
layer in a set of computation points. Add checks to avoid generating layers
with tesseroids that overlap around the earth. Add tests and a new gallery
example.

Fixes #83 

Patch:
diff --git a/examples/forward/tesseroid_layer.py b/examples/forward/tesseroid_layer.py
new file mode 100644
index 0000000..ad9977c
--- /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 7c2d4eb..f745da1 100644
--- a/harmonica/__init__.py
+++ b/harmonica/__init__.py
@@ -15,6 +15,7 @@ from .forward.point import point_gravity
 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 0000000..6c121f0
--- /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


Test patch:
diff --git a/harmonica/tests/test_tesseroid_layer.py b/harmonica/tests/test_tesseroid_layer.py
new file mode 100644
index 0000000..75ba11e
--- /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),
+    )