diff --git a/.github/environment.yml b/.github/environment.yml
index e53af9e88..f8878cee6 100644
--- a/.github/environment.yml
+++ b/.github/environment.yml
@@ -16,6 +16,7 @@ dependencies:
   - xrft>=1.0
   - choclo>=0.1
   - boule>=0.6.0
+  - bordado>=0.4.0
   # Build
   - python-build
   - twine
@@ -30,7 +31,6 @@ dependencies:
   - ipykernel
   - boule
   - pyproj
-  - bordado>=0.4.0
   - ensaio>=0.5.*
   - netcdf4
   - matplotlib
diff --git a/doc/conf.py b/doc/conf.py
index 68fc5486e..ad4d4f42c 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -53,6 +53,7 @@
     "pooch": ("https://www.fatiando.org/pooch/latest/", None),
     "ensaio": ("https://www.fatiando.org/ensaio/latest/", None),
     "verde": ("https://www.fatiando.org/verde/latest/", None),
+    "bordado": ("https://www.fatiando.org/bordado/latest/", None),
     "boule": ("https://www.fatiando.org/boule/latest/", None),
     "matplotlib": ("https://matplotlib.org/", None),
     "pyproj": ("https://pyproj4.github.io/pyproj/stable/", None),
diff --git a/doc/install.rst b/doc/install.rst
index c91c9ff8f..8874051f6 100644
--- a/doc/install.rst
+++ b/doc/install.rst
@@ -74,6 +74,8 @@ Required:
 * `scikit-learn <https://scikit-learn.org>`__
 * `pooch <http://www.fatiando.org/pooch/>`__
 * `verde <http://www.fatiando.org/verde/>`__
+* `bordado <http://www.fatiando.org/bordado/>`__
+* `boule <http://www.fatiando.org/boule/>`__
 * `xrft <https://xrft.readthedocs.io/>`__
 
 Optional:
@@ -92,7 +94,6 @@ Optional:
 
 The examples in the :ref:`gallery` also use:
 
-* `boule <http://www.fatiando.org/boule/>`__
 * `ensaio <http://www.fatiando.org/ensaio/>`__ for downloading sample datasets
 * `pygmt <https://www.pygmt.org/>`__ for plotting maps
 * `pyproj <https://jswhit.github.io/pyproj/>`__ for cartographic projections
diff --git a/doc/user_guide/coordinate_systems.rst b/doc/user_guide/coordinate_systems.rst
index 5162cd666..83721e061 100644
--- a/doc/user_guide/coordinate_systems.rst
+++ b/doc/user_guide/coordinate_systems.rst
@@ -33,7 +33,7 @@ height:
 
 .. jupyter-execute::
 
-    import verde as vd
+    import bordado as bd
 
     # Define boundaries of the rectangular prism (in meters)
     west, east, south, north = -20, 20, -20, 20
@@ -45,8 +45,8 @@ Or a regular grid of points at 100 meters above the zeroth plane:
 .. jupyter-execute::
 
     # Define a regular grid of observation points (coordinates in meters)
-    coordinates = vd.grid_coordinates(
-        region=(-40, 40, -40, 40), shape=(5, 5), extra_coords=100
+    coordinates = bd.grid_coordinates(
+        region=(-40, 40, -40, 40), shape=(5, 5), non_dimensional_coords=100
     )
     easting, northing, upward = coordinates[:]
 
@@ -103,12 +103,12 @@ Spatial data are usually given in geodetic coordinates, along with the
 reference ellipsoid on which they are defined.
 
 For example, let's define a regular grid of points (separated by equal
-angles) at 2km above the ellipsoid using :mod:`verde`.
+angles) at 2km above the ellipsoid using :mod:`bordado`.
 
 .. jupyter-execute::
 
-    coordinates = vd.grid_coordinates(
-        region=(-70, -65, -35, -30), shape=(6, 6), extra_coords=2e3
+    coordinates = bd.grid_coordinates(
+        region=(-70, -65, -35, -30), shape=(6, 6), non_dimensional_coords=2e3
     )
     longitude, latitude, height = coordinates[:]
     print("longitude:", longitude)
@@ -199,10 +199,10 @@ the same radius equal to the *mean radius of the Earth*.
 
 .. jupyter-execute::
 
-    coordinates = vd.grid_coordinates(
+    coordinates = bd.grid_coordinates(
         region=(-70, -65, -35, -30),
         shape=(6, 6),
-        extra_coords=ellipsoid.mean_radius,
+        non_dimensional_coords=ellipsoid.mean_radius,
     )
     longitude, sph_latitude, radius = coordinates[:]
     print("longitude:", longitude)
diff --git a/doc/user_guide/equivalent_sources/block-averaged-eqs.rst b/doc/user_guide/equivalent_sources/block-averaged-eqs.rst
index c27e08dbb..59db0f641 100644
--- a/doc/user_guide/equivalent_sources/block-averaged-eqs.rst
+++ b/doc/user_guide/equivalent_sources/block-averaged-eqs.rst
@@ -44,10 +44,10 @@ the data:
 
 .. jupyter-execute::
 
-    import verde as vd
+    import bordado as bd
 
     region = (-5.5, -4.7, 57.8, 58.5)
-    inside = vd.inside((data.longitude, data.latitude), region)
+    inside = bd.inside((data.longitude, data.latitude), region)
     data = data[inside]
     data
 
@@ -60,7 +60,7 @@ And project the geographic coordinates to plain Cartesian ones:
     projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean())
     easting, northing = projection(data.longitude.values, data.latitude.values)
     coordinates = (easting, northing, data.height_m)
-    xy_region=vd.get_region(coordinates)
+    xy_region = bd.get_region((easting, northing))
 
 
 .. jupyter-execute::
@@ -77,9 +77,10 @@ And project the geographic coordinates to plain Cartesian ones:
 
 .. jupyter-execute::
 
+    import verde as vd
     import pygmt
 
-    maxabs = vd.maxabs(data.total_field_anomaly_nt)*.8
+    maxabs = vd.maxabs(data.total_field_anomaly_nt) * .8
 
     # Set figure properties
     w, e, s, n = xy_region
@@ -167,10 +168,10 @@ we are efectivelly upward continuing the data.
 
 .. jupyter-execute::
 
-    grid_coords = vd.grid_coordinates(
-        region=vd.get_region(coordinates),
+    grid_coords = bd.grid_coordinates(
+        region=bd.get_region((easting, northing)),
         spacing=500,
-        extra_coords=1500,
+        non_dimensional_coords=1500,
     )
     grid = eqs.grid(grid_coords, data_names=["magnetic_anomaly"])
     grid
diff --git a/doc/user_guide/equivalent_sources/eq-sources-spherical.rst b/doc/user_guide/equivalent_sources/eq-sources-spherical.rst
index a27d7f529..2166659c8 100644
--- a/doc/user_guide/equivalent_sources/eq-sources-spherical.rst
+++ b/doc/user_guide/equivalent_sources/eq-sources-spherical.rst
@@ -98,15 +98,17 @@ of 6 arcminutes.
 
 .. jupyter-execute::
 
+    import bordado as bd
+
     # Get the bounding region of the data in geodetic coordinates
-    region = vd.get_region((data.longitude, data.latitude))
+    region = bd.get_region((data.longitude, data.latitude))
 
     # Get the maximum height of the data coordinates
     max_height = data.height_sea_level_m.max()
 
     # Define a regular grid of points in geodetic coordinates
-    grid_coords = vd.grid_coordinates(
-        region=region, spacing=6 / 60, extra_coords=max_height
+    grid_coords = bd.grid_coordinates(
+        region=region, spacing=6 / 60, non_dimensional_coords=max_height
     )
 
 But before we can tell the equivalent sources to predict the
diff --git a/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst b/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst
index 22b5c2c57..e091a4606 100644
--- a/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst
+++ b/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst
@@ -175,12 +175,14 @@ And grid the data using the two equivalent sources:
 
 .. jupyter-execute::
 
+   import bordado as bd
+
    # Define grid coordinates
-   region = vd.get_region(coordinates)
-   grid_coords = vd.grid_coordinates(
+   region = bd.get_region((easting, northing))
+   grid_coords = bd.grid_coordinates(
        region=region,
        spacing=2e3,
-       extra_coords=2.5e3,
+       non_dimensional_coords=2.5e3,
    )
 
    grid_first_guess = eqs_first_guess.grid(grid_coords)
diff --git a/doc/user_guide/equivalent_sources/gradient-boosted-eqs.rst b/doc/user_guide/equivalent_sources/gradient-boosted-eqs.rst
index e44a3d218..c60055144 100644
--- a/doc/user_guide/equivalent_sources/gradient-boosted-eqs.rst
+++ b/doc/user_guide/equivalent_sources/gradient-boosted-eqs.rst
@@ -125,12 +125,13 @@ And then predict the field on a regular grid of computation points:
 
 .. jupyter-execute::
 
-    import verde as vd
-    region = vd.get_region(coordinates)
-    grid_coords = vd.grid_coordinates(
+    import bordado as bd
+
+    region = bd.get_region((easting, northing))
+    grid_coords = bd.grid_coordinates(
         region=region,
         spacing=5e3,
-        extra_coords=2.5e3,
+        non_dimensional_coords=2.5e3,
     )
     grid = eqs.grid(grid_coords, data_names=["gravity_disturbance"])
     grid
@@ -141,6 +142,8 @@ point. We can do so through the :func:`verde.distance_mask` function.
 
 .. jupyter-execute::
 
+    import verde as vd
+
     grid_masked = vd.distance_mask(coordinates, maxdist=50e3, grid=grid)
 
 And plot it:
diff --git a/doc/user_guide/equivalent_sources/index.rst b/doc/user_guide/equivalent_sources/index.rst
index 266a9e3e7..07b07a553 100644
--- a/doc/user_guide/equivalent_sources/index.rst
+++ b/doc/user_guide/equivalent_sources/index.rst
@@ -59,11 +59,11 @@ coordinates, so we need to project these gravity observations:
 .. jupyter-execute::
 
    import pyproj
-   import verde as vd
+   import bordado as bd
 
    projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.values.mean())
    easting, northing = projection(data.longitude.values, data.latitude.values)
-   region = vd.get_region((easting, northing))
+   region = bd.get_region((easting, northing))
 
 Now we can initialize the :class:`harmonica.EquivalentSources` class.
 
@@ -131,6 +131,7 @@ And plot it:
 .. jupyter-execute::
 
    import pygmt
+   import verde as vd
 
    # Get max absolute value for the observed gravity disturbance
    maxabs = vd.maxabs(data.gravity_disturbance_mgal)
@@ -176,7 +177,7 @@ And plot it:
 
 We can also *grid* and *upper continue* the field by predicting its values on
 a regular grid at a constant height higher than the observations. To do so we
-can use the :func:`verde.grid_coordinates` function to create the coordinates
+can use the :func:`bordado.grid_coordinates` function to create the coordinates
 of the grid and then use the :meth:`harmonica.EquivalentSources.grid` method.
 
 First, lets get the maximum height of the observations:
@@ -191,7 +192,9 @@ and use the equivalent sources to generate a gravity disturbance grid.
 .. jupyter-execute::
 
    # Build the grid coordinates
-   grid_coords = vd.grid_coordinates(region=region, spacing=2e3, extra_coords=2.2e3)
+   grid_coords = bd.grid_coordinates(
+       region=region, spacing=2e3, non_dimensional_coords=2.2e3
+   )
 
    # Grid the gravity disturbances
    grid = equivalent_sources.grid(grid_coords, data_names=["gravity_disturbance"])
diff --git a/doc/user_guide/forward_modelling/dipole.rst b/doc/user_guide/forward_modelling/dipole.rst
index e7a5e9a6d..eada12bb4 100644
--- a/doc/user_guide/forward_modelling/dipole.rst
+++ b/doc/user_guide/forward_modelling/dipole.rst
@@ -60,13 +60,13 @@ Let's define a regular grid of observation points:
 
 .. jupyter-execute::
 
-   import verde as vd
+   import bordado as bd
 
    region = (-100, 100, -100, 100)
    spacing = 1
    height = 0
-   coordinates = vd.grid_coordinates(
-      region=region, spacing=spacing, extra_coords=height
+   coordinates = bd.grid_coordinates(
+      region=region, spacing=spacing, non_dimensional_coords=height
    )
 
 And a set of dipoles with their magnetic moments:
diff --git a/doc/user_guide/forward_modelling/point.rst b/doc/user_guide/forward_modelling/point.rst
index 1b4942c9e..d1f9574f5 100644
--- a/doc/user_guide/forward_modelling/point.rst
+++ b/doc/user_guide/forward_modelling/point.rst
@@ -70,10 +70,10 @@ height:
 
 .. jupyter-execute::
 
-   import verde as vd
+   import bordado as bd
 
-   coordinates = vd.grid_coordinates(
-       region=(-250, 1250, -250, 1250), shape=(40, 40), extra_coords=0
+   coordinates = bd.grid_coordinates(
+       region=(-250, 1250, -250, 1250), shape=(40, 40), non_dimensional_coords=0
    )
 
 And finally calculate the vertical component of the gravitational acceleration
@@ -106,6 +106,7 @@ Lets plot this gravitational field:
 .. jupyter-execute::
 
    import pygmt
+   import verde as vd
 
    grid = vd.make_xarray_grid(
       coordinates, g_z, data_names="g_z", extra_coords_names="extra")
@@ -200,10 +201,10 @@ coordinates and located at 1km above the ellipsoid:
 
 .. jupyter-execute::
 
-   coordinates = vd.grid_coordinates(
+   coordinates = bd.grid_coordinates(
        region=(-72, -68, -46, -42),
        shape=(101, 101),
-       extra_coords=20e3,
+       non_dimensional_coords=20e3,
    )
 
 Before we can start forward modelling these sources, we need to convert them to
diff --git a/doc/user_guide/forward_modelling/prism.rst b/doc/user_guide/forward_modelling/prism.rst
index de520f41b..b6c7383e3 100644
--- a/doc/user_guide/forward_modelling/prism.rst
+++ b/doc/user_guide/forward_modelling/prism.rst
@@ -61,12 +61,14 @@ Build a regular grid of computation points located 10m above the zero height:
 
 .. jupyter-execute::
 
-   import verde as vd
+   import bordado as bd
 
    region = (-10e3, 10e3, -10e3, 10e3)
    shape = (51, 51)
    height = 10
-   coordinates = vd.grid_coordinates(region, shape=shape, extra_coords=height)
+   coordinates = bd.grid_coordinates(
+       region, shape=shape, non_dimensional_coords=height
+   )
 
 Define a single prism:
 
@@ -164,10 +166,10 @@ height:
 
 .. jupyter-execute::
 
-   import verde as vd
+   import bordado as bd
 
-   coordinates = vd.grid_coordinates(
-       region=(0, 10e3, 0, 10e3), shape=(40, 40), extra_coords=0
+   coordinates = bd.grid_coordinates(
+       region=(0, 10e3, 0, 10e3), shape=(40, 40), non_dimensional_coords=0
    )
 
 And finally calculate the vertical component of the gravitational acceleration
@@ -200,6 +202,8 @@ Lets plot this gravitational field:
 .. jupyter-execute::
 
    import pygmt
+   import verde as vd
+
    grid = vd.make_xarray_grid(
       coordinates, g_z, data_names="g_z", extra_coords_names="extra")
    fig = pygmt.Figure()
@@ -253,7 +257,9 @@ points by choosing ``field="b"``:
    region = (-10e3, 10e3, -10e3, 10e3)
    shape = (51, 51)
    height = 10
-   coordinates = vd.grid_coordinates(region, shape=shape, extra_coords=height)
+   coordinates = bd.grid_coordinates(
+       region, shape=shape, non_dimensional_coords=height
+   )
 
 .. jupyter-execute::
 
@@ -364,7 +370,7 @@ Then we can define a regular grid where the centers of the prisms will fall:
 
 .. jupyter-execute::
 
-   easting, northing = vd.grid_coordinates(region=region, spacing=spacing)
+   easting, northing = bd.grid_coordinates(region=region, spacing=spacing)
 
 We need to define a 2D array for the uppermost *surface* of the layer. We will
 sample a trigonometric function for this simple example:
@@ -396,9 +402,9 @@ Let's define a grid of observation points at 1 km above the zeroth height:
 
 .. jupyter-execute::
 
-   region_pad = vd.pad_region(region, 10e3)
-   coordinates = vd.grid_coordinates(
-       region_pad, spacing=spacing, extra_coords=1e3
+   region_pad = bd.pad_region(region, 10e3)
+   coordinates = bd.grid_coordinates(
+       region_pad, spacing=spacing, non_dimensional_coords=1e3
    )
 
 
diff --git a/doc/user_guide/forward_modelling/tesseroid.rst b/doc/user_guide/forward_modelling/tesseroid.rst
index 292ab2683..69374ceb0 100644
--- a/doc/user_guide/forward_modelling/tesseroid.rst
+++ b/doc/user_guide/forward_modelling/tesseroid.rst
@@ -75,12 +75,12 @@ the *top* surface of the tesseroid:
 
 .. jupyter-execute::
 
-   import verde as vd
+   import bordado as bd
 
-   coordinates = vd.grid_coordinates(
+   coordinates = bd.grid_coordinates(
        region=[-80, -40, -50, -10],
        shape=(80, 80),
-       extra_coords=100e3 + mean_radius,
+       non_dimensional_coords=100e3 + mean_radius,
    )
 
 Lets compute the *downward* component of the gravitational acceleration it
@@ -116,6 +116,8 @@ And finally plot the computed gravitational field
 .. jupyter-execute::
 
    import pygmt
+   import verde as vd
+
    grid = vd.make_xarray_grid(
       coordinates, gravity, data_names="gravity", extra_coords_names="extra")
 
@@ -132,7 +134,7 @@ And finally plot the computed gravitational field
 
    fig.colorbar(cmap=True, frame=["a200f50", "x+lmGal"])
    fig.coast(shorelines="1p,black")
-   
+
    # Plot edges of tesseroid
    fig.plot(
       x=[tesseroid[0], tesseroid[1], tesseroid[1], tesseroid[0], tesseroid[0]],
@@ -168,10 +170,10 @@ Compute their gravitational effect on a grid of computation points:
 
 .. jupyter-execute::
 
-   coordinates = vd.grid_coordinates(
+   coordinates = bd.grid_coordinates(
        region=[-80, -40, -50, -10],
        shape=(80, 80),
-       extra_coords=100e3 + mean_radius,
+       non_dimensional_coords=100e3 + mean_radius,
    )
    gravity = hm.tesseroid_gravity(coordinates, tesseroids, densities, field="g_z")
 
@@ -263,10 +265,10 @@ above the mean Earth radius:
 
 .. jupyter-execute::
 
-   coordinates = vd.grid_coordinates(
+   coordinates = bd.grid_coordinates(
        region=[-80, -40, -50, -10],
        shape=(80, 80),
-       extra_coords=100e3 + ellipsoid.mean_radius,
+       non_dimensional_coords=100e3 + ellipsoid.mean_radius,
    )
 
 And compute the gravitational fields the tesseroids generate:
diff --git a/doc/user_guide/topographic_correction.rst b/doc/user_guide/topographic_correction.rst
index ec7b60c07..fb5f8893a 100644
--- a/doc/user_guide/topographic_correction.rst
+++ b/doc/user_guide/topographic_correction.rst
@@ -165,8 +165,10 @@ And then crop it to a slightly larger region than the gravity observations:
 
 .. jupyter-execute::
 
-   region = vd.get_region((data.longitude, data.latitude))
-   region_pad = vd.pad_region(region, pad=1)
+   import bordado as bd
+
+   region = bd.get_region((data.longitude, data.latitude))
+   region_pad = bd.pad_region(region, pad=1)
 
    topography = topography.sel(
        longitude=slice(region_pad[0], region_pad[1]),
diff --git a/environment.yml b/environment.yml
index 956cf3305..5a7272946 100644
--- a/environment.yml
+++ b/environment.yml
@@ -18,6 +18,7 @@ dependencies:
   - xrft>=1.0
   - choclo>=0.1
   - boule>=0.6
+  - bordado>=0.4.0
   # Optional requirements
   - pyvista
   - vtk
@@ -38,7 +39,6 @@ dependencies:
   - ipykernel
   - pyproj
   - matplotlib
-  - bordado>=0.4.0
   - ensaio>=0.5.0
   - netcdf4
   - pygmt
diff --git a/pyproject.toml b/pyproject.toml
index cad85d5de..9d727e537 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -39,7 +39,8 @@ dependencies = [
     "verde >= 1.8.1",
     "xrft >= 1.0",
     "choclo >= 0.1",
-    "boule >= 0.6.0"
+    "boule >= 0.6.0",
+    "bordado >= 0.4.0",
 ]
 
 [project.urls]
diff --git a/src/harmonica/_equivalent_sources/cartesian.py b/src/harmonica/_equivalent_sources/cartesian.py
index 0b1f895ca..1e9b982d5 100644
--- a/src/harmonica/_equivalent_sources/cartesian.py
+++ b/src/harmonica/_equivalent_sources/cartesian.py
@@ -12,6 +12,7 @@
 
 import warnings
 
+import bordado as bd
 import numpy as np
 import verde as vd
 import verde.base as vdb
@@ -221,7 +222,7 @@ def fit(self, coordinates, data, weights=None):
             coordinates, data, weights, self.dtype
         )
         # Capture the data region to use as a default when gridding.
-        self.region_ = vd.get_region(coordinates[:2])
+        self.region_ = bd.get_region(coordinates[:2])
         coordinates = vdb.n_1d_arrays(coordinates, 3)
         if self.points is None:
             self.points_ = tuple(
@@ -268,7 +269,9 @@ def _build_points(self, coordinates):
         if self.block_size is not None:
             coordinates = self._block_average_coordinates(coordinates)
         if self.depth == "default":
-            self.depth_ = 4.5 * np.mean(vd.median_distance(coordinates, k_nearest=1))
+            self.depth_ = 4.5 * np.mean(
+                bd.neighbor_distance_statistics(coordinates[:2], "median", k=1)
+            )
         else:
             self.depth_ = self.depth
         return (
@@ -382,9 +385,9 @@ def grid(
         The coordinates of the regular grid must be passed through the
         ``coordinates`` argument as a tuple containing three arrays in the
         following order: ``(easting, nothing, upward)``. They can be easily
-        created through the :func:`verde.grid_coordinates` function. If the
+        created through the :func:`bordado.grid_coordinates` function. If the
         grid points must be all at the same height, it can be specified in the
-        ``extra_coords`` argument of :func:`verde.grid_coordinates`.
+        ``extra_coords`` argument of :func:`bordado.grid_coordinates`.
 
         Use the *dims* and *data_names* arguments to set custom names for the
         dimensions and the data field(s) in the output :class:`xarray.Dataset`.
@@ -427,7 +430,7 @@ def grid(
 
         See Also
         --------
-        :func:`verde.grid_coordinates`
+        :func:`bordado.grid_coordinates`
 
         """
         # We override the grid method from BaseGridder to change the docstring
diff --git a/src/harmonica/_equivalent_sources/gradient_boosted.py b/src/harmonica/_equivalent_sources/gradient_boosted.py
index 76b37de1e..ab522e394 100644
--- a/src/harmonica/_equivalent_sources/gradient_boosted.py
+++ b/src/harmonica/_equivalent_sources/gradient_boosted.py
@@ -12,10 +12,10 @@
 
 import warnings
 
+import bordado as bd
 import numpy as np
 import verde.base as vdb
 from sklearn import utils
-from verde import get_region, rolling_window
 
 from .cartesian import EquivalentSources
 from .utils import cast_fit_input, predict_numba_parallel
@@ -159,17 +159,17 @@ def estimate_required_memory(self, coordinates):
 
         Examples
         --------
-        >>> import verde as vd
-        >>> coordinates = vd.scatter_points(
+        >>> import bordado as bd
+        >>> coordinates = bd.random_coordinates(
         ...     region=(-1e3, 3e3, 2e3, 5e3),
         ...     size=100,
-        ...     extra_coords=100,
-        ...     random_state=42,
+        ...     non_dimensional_coords=100,
+        ...     random_seed=42,
         ... )
         >>> eqs = EquivalentSourcesGB(window_size=2e3)
         >>> n_bytes = eqs.estimate_required_memory(coordinates)
         >>> int(n_bytes)
-        9800
+        10952
         """
         # Build the sources and assign the points_ attribute
         coordinates = vdb.n_1d_arrays(coordinates, 3)
@@ -219,7 +219,7 @@ def fit(self, coordinates, data, weights=None):
             coordinates, data, weights, self.dtype
         )
         # Capture the data region to use as a default when gridding.
-        self.region_ = get_region(coordinates[:2])
+        self.region_ = bd.get_region(coordinates[:2])
         # Ravel coordinates, data and weights to 1d-arrays
         coordinates = vdb.n_1d_arrays(coordinates, 3)
         data = np.ravel(data)
@@ -342,21 +342,19 @@ def _create_windows(self, coordinates, shuffle=True):
             self.window_size_ = np.sqrt(window_area)
         else:
             self.window_size_ = self.window_size
-        # Compute window spacing based on overlapping
-        window_spacing = self.window_size_ * (1 - self.overlapping)
         # The windows for sources and data points are the same, but the
-        # verde.rolling_window function creates indices for the given
+        # bordado.rolling_window function creates indices for the given
         # coordinates. That's why we need to create two set of window indices:
         # one for the sources and one for the data points.
         # We pass the same region, size and spacing to be sure that both set of
         # windows are the same.
         kwargs = {
             "region": region,
-            "size": self.window_size_,
-            "spacing": window_spacing,
+            "window_size": self.window_size_,
+            "overlap": self.overlapping,
         }
-        _, source_windows = rolling_window(self.points_, **kwargs)
-        _, data_windows = rolling_window(coordinates, **kwargs)
+        _, source_windows = bd.rolling_window(self.points_[:2], **kwargs)
+        _, data_windows = bd.rolling_window(coordinates[:2], **kwargs)
         # Ravel the indices
         source_windows = [i[0] for i in source_windows.ravel()]
         data_windows = [i[0] for i in data_windows.ravel()]
@@ -388,8 +386,8 @@ def _get_region_data_sources(coordinates, points):
     -------
     region : tuple
     """
-    data_region = get_region(coordinates)
-    sources_region = get_region(points)
+    data_region = bd.get_region(coordinates)
+    sources_region = bd.get_region(points)
     region = (
         min(data_region[0], sources_region[0]),
         max(data_region[1], sources_region[1]),
diff --git a/src/harmonica/_equivalent_sources/spherical.py b/src/harmonica/_equivalent_sources/spherical.py
index 8c205f68e..52368daff 100644
--- a/src/harmonica/_equivalent_sources/spherical.py
+++ b/src/harmonica/_equivalent_sources/spherical.py
@@ -10,8 +10,8 @@
 
 import warnings
 
+import bordado as bd
 import numpy as np
-import verde as vd
 import verde.base as vdb
 from numba import jit
 from sklearn.utils.validation import check_is_fitted
@@ -152,7 +152,7 @@ def fit(self, coordinates, data, weights=None):
         """
         coordinates, data, weights = vdb.check_fit_input(coordinates, data, weights)
         # Capture the data region to use as a default when gridding.
-        self.region_ = vd.get_region(coordinates[:2])
+        self.region_ = bd.get_region(coordinates[:2])
         coordinates = vdb.n_1d_arrays(coordinates, 3)
         if self.points is None:
             self.points_ = (
@@ -244,9 +244,9 @@ def grid(
         The coordinates of the regular grid must be passed through the
         ``coordinates`` argument as a tuple containing three arrays in the
         following order: ``(longitude, latitude, radius)``. They can be easily
-        created through the :func:`verde.grid_coordinates` function. If the
+        created through the :func:`bordado.grid_coordinates` function. If the
         grid points must be all at the same radius, it can be specified in the
-        ``extra_coords`` argument of :func:`verde.grid_coordinates`.
+        ``extra_coords`` argument of :func:`bordado.grid_coordinates`.
 
         Use the *dims* and *data_names* arguments to set custom names for the
         dimensions and the data field(s) in the output :class:`xarray.Dataset`.
diff --git a/src/harmonica/_spherical_harmonics/igrf.py b/src/harmonica/_spherical_harmonics/igrf.py
index 697980022..42b13311c 100644
--- a/src/harmonica/_spherical_harmonics/igrf.py
+++ b/src/harmonica/_spherical_harmonics/igrf.py
@@ -11,6 +11,7 @@
 import datetime
 import pathlib
 
+import bordado as bd
 import boule
 import numba
 import numpy as np
@@ -437,8 +438,12 @@ def grid(self, region, height, shape=None, spacing=None, adjust="spacing"):
         """
         if spacing is None and shape is None:
             spacing = calculate_ideal_spacing(self.max_degree)
-        longitude, latitude, height = vd.grid_coordinates(
-            region, spacing=spacing, shape=shape, adjust=adjust, extra_coords=height
+        longitude, latitude, height = bd.grid_coordinates(
+            region,
+            spacing=spacing,
+            shape=shape,
+            adjust=adjust,
+            non_dimensional_coords=height,
         )
         longitude, latitude_sph, radius = self.ellipsoid.geodetic_to_spherical(
             (longitude, latitude, height)
diff --git a/test/ellipsoids/test_gravity.py b/test/ellipsoids/test_gravity.py
index 11b01b275..999060fcd 100644
--- a/test/ellipsoids/test_gravity.py
+++ b/test/ellipsoids/test_gravity.py
@@ -17,6 +17,7 @@
 import re
 from copy import copy
 
+import bordado as bd
 import numpy as np
 import pytest
 import verde as vd
@@ -72,11 +73,11 @@ def test_opposite_planes():
 
     # define observation points (2D grid) at surface height (z axis,
     # 'Upward') = 5
-    coordinates1 = vd.grid_coordinates(
-        region=(-20, 20, -20, 20), spacing=0.5, extra_coords=5
+    coordinates1 = bd.grid_coordinates(
+        region=(-20, 20, -20, 20), spacing=0.5, non_dimensional_coords=5
     )
-    coordinates2 = vd.grid_coordinates(
-        region=(-20, 20, -20, 20), spacing=0.5, extra_coords=-5
+    coordinates2 = bd.grid_coordinates(
+        region=(-20, 20, -20, 20), spacing=0.5, non_dimensional_coords=-5
     )
 
     _, _, gz1 = ellipsoid_gravity(coordinates1, triaxial_example)
@@ -452,8 +453,8 @@ def test_symmetry_when_flipping(self, ellipsoid):
         by the ellipsoid should be the same as before the rotation.
         """
         # Define observation points
-        coordinates = vd.grid_coordinates(
-            region=(-20, 20, -20, 20), spacing=0.5, extra_coords=5
+        coordinates = bd.grid_coordinates(
+            region=(-20, 20, -20, 20), spacing=0.5, non_dimensional_coords=5
         )
 
         # Generate a flipped ellipsoid
@@ -479,8 +480,8 @@ class TestMultipleEllipsoids:
     def coordinates(self):
         """Sample grid coordinates."""
         region = (-30, 30, -30, 30)
-        coordinates = vd.grid_coordinates(
-            region=region, shape=(21, 21), extra_coords=10
+        coordinates = bd.grid_coordinates(
+            region=region, shape=(21, 21), non_dimensional_coords=10
         )
         return coordinates
 
diff --git a/test/ellipsoids/test_magnetic.py b/test/ellipsoids/test_magnetic.py
index a8bb974a0..eafa9a5e9 100644
--- a/test/ellipsoids/test_magnetic.py
+++ b/test/ellipsoids/test_magnetic.py
@@ -18,6 +18,7 @@
 import re
 from copy import copy
 
+import bordado as bd
 import numpy as np
 import pytest
 import verde as vd
@@ -54,11 +55,11 @@ def test_magnetic_symmetry():
     susceptibility = 0.1
     ellipsoid = Ellipsoid(4, 3, 2, susceptibility=susceptibility)
 
-    coordinates = vd.grid_coordinates(
-        region=(-20, 20, -20, 20), spacing=0.5, extra_coords=5
+    coordinates = bd.grid_coordinates(
+        region=(-20, 20, -20, 20), spacing=0.5, non_dimensional_coords=5
     )
-    coordinates2 = vd.grid_coordinates(
-        region=(-20, 20, -20, 20), spacing=0.5, extra_coords=-5
+    coordinates2 = bd.grid_coordinates(
+        region=(-20, 20, -20, 20), spacing=0.5, non_dimensional_coords=-5
     )
 
     inducing_field = hm.magnetic_angles_to_vec(10_000, 0, 0)
@@ -79,8 +80,8 @@ def test_flipped_h0():
     susceptibility = 0.1
     oblate = Ellipsoid(a, a, c, susceptibility=susceptibility)
 
-    coordinates = vd.grid_coordinates(
-        region=(-20, 20, -20, 20), spacing=0.5, extra_coords=5
+    coordinates = bd.grid_coordinates(
+        region=(-20, 20, -20, 20), spacing=0.5, non_dimensional_coords=5
     )
 
     inducing_field = np.asarray(
@@ -100,8 +101,8 @@ def test_zero_susceptibility():
     """
     susceptibility = 0
     ellipsoid = Ellipsoid(2, 2, 1, susceptibility=susceptibility)
-    coordinates = vd.grid_coordinates(
-        region=(-10, 10, -10, 10), spacing=1.0, extra_coords=5
+    coordinates = bd.grid_coordinates(
+        region=(-10, 10, -10, 10), spacing=1.0, non_dimensional_coords=5
     )
     inducing_field = hm.magnetic_angles_to_vec(55_000, 0.0, 90.0)
 
@@ -118,8 +119,8 @@ def test_zero_field():
     """
     susceptibility = 0.01
     ellipsoid = Ellipsoid(2, 1, 1, susceptibility=susceptibility)
-    coordinates = vd.grid_coordinates(
-        region=(-10, 10, -10, 10), spacing=1.0, extra_coords=5
+    coordinates = bd.grid_coordinates(
+        region=(-10, 10, -10, 10), spacing=1.0, non_dimensional_coords=5
     )
 
     inducing_field = (0, 0, 0)
@@ -177,8 +178,8 @@ def test_mag_flipped_ellipsoid():
 
     # define observation points (2D grid) at surface height (z axis,
     # 'Upward') = 5
-    x, y, z = vd.grid_coordinates(
-        region=(-20, 20, -20, 20), spacing=0.5, extra_coords=5
+    x, y, z = bd.grid_coordinates(
+        region=(-20, 20, -20, 20), spacing=0.5, non_dimensional_coords=5
     )
 
     # ignore internal field as this won't be 'flipped' in the same natr
@@ -202,8 +203,8 @@ def test_euler_rotation_symmetry_mag():
     a, b, c = 5, 4, 3
     inducing_field = hm.magnetic_angles_to_vec(55_000, 0.0, 90.0)
     susceptibility = 0.01
-    coordinates = x, y, z = vd.grid_coordinates(
-        region=(-5, 5, -5, 5), spacing=1.0, extra_coords=5
+    x, y, z = bd.grid_coordinates(
+        region=(-5, 5, -5, 5), spacing=1.0, non_dimensional_coords=5
     )
     internal_mask = ((x**2) / (a**2) + (y**2) / (b**2) + (z**2) / (c**2)) < 1
     coordinates = tuple(c[internal_mask] for c in (x, y, z))
@@ -432,7 +433,9 @@ def coordinates(self, request):
         region = (-200, 200, -200, 200)
         shape = (151, 151)
         height = request.param
-        coordinates = vd.grid_coordinates(region, shape=shape, extra_coords=height)
+        coordinates = bd.grid_coordinates(
+            region, shape=shape, non_dimensional_coords=height
+        )
         return coordinates
 
     def get_ellipsoid(self, ellipsoid_type: str):
@@ -630,8 +633,8 @@ def test_symmetry_when_flipping(self, ellipsoid, magnetization_type):
         system, it won't rotate with the ellipsoid.
         """
         # Define observation points
-        coordinates = vd.grid_coordinates(
-            region=(-20, 20, -20, 20), spacing=0.5, extra_coords=5
+        coordinates = bd.grid_coordinates(
+            region=(-20, 20, -20, 20), spacing=0.5, non_dimensional_coords=5
         )
 
         # Define physical properties
@@ -668,8 +671,8 @@ class TestMultipleEllipsoids:
     def coordinates(self):
         """Sample grid coordinates."""
         region = (-30, 30, -30, 30)
-        coordinates = vd.grid_coordinates(
-            region=region, shape=(21, 21), extra_coords=10
+        coordinates = bd.grid_coordinates(
+            region=region, shape=(21, 21), non_dimensional_coords=10
         )
         return coordinates
 
@@ -1132,7 +1135,9 @@ def test_invariance_semiaxes_order(self):
         rotation matrix.
         """
         region = (-50, 50, -50, 50)
-        coordinates = vd.grid_coordinates(region, shape=(151, 151), extra_coords=0)
+        coordinates = bd.grid_coordinates(
+            region, shape=(151, 151), non_dimensional_coords=0
+        )
         inducing_field = np.array([0, 0, -20_000])  # pointing on z
         # Define susceptibility as a tensor with only component on the vertical axes
         # of the local coordinate system.
diff --git a/test/test_dipole.py b/test/test_dipole.py
index 17f05c625..028b4c001 100644
--- a/test/test_dipole.py
+++ b/test/test_dipole.py
@@ -8,11 +8,11 @@
 Test magnetic forward functions for dipoles.
 """
 
+import bordado as bd
 import choclo
 import numpy as np
 import numpy.testing as npt
 import pytest
-import verde as vd
 from choclo.dipole import magnetic_field
 
 try:
@@ -61,8 +61,8 @@ def test_progress_bar(field):
     """
     dipoles = ([1, 1], [-1, 0], [-10, -2])
     magnetic_moments = ([1.0, 1.0], [1.0, -1.0], [1.0, 5.0])
-    coordinates = vd.grid_coordinates(
-        region=(-100, 100, -100, 100), spacing=20, extra_coords=10
+    coordinates = bd.grid_coordinates(
+        region=(-100, 100, -100, 100), spacing=20, non_dimensional_coords=10
     )
     result_progress_true = dipole_magnetic(
         coordinates, dipoles, magnetic_moments, field=field, progressbar=True
@@ -140,8 +140,8 @@ def test_dipoles_parallel_vs_serial(
         Check results of parallelized and serials runs
         Run a large problem only with Numba enabled.
         """
-        coordinates = vd.grid_coordinates(
-            region=(-100, 100, -100, 100), spacing=20, extra_coords=10
+        coordinates = bd.grid_coordinates(
+            region=(-100, 100, -100, 100), spacing=20, non_dimensional_coords=10
         )
         dipoles = ([-100, 0], [0, 100], [-20, -50])
         magnetic_moments = [
diff --git a/test/test_eq_sources_cartesian.py b/test/test_eq_sources_cartesian.py
index 284732f86..d7267d44e 100644
--- a/test/test_eq_sources_cartesian.py
+++ b/test/test_eq_sources_cartesian.py
@@ -10,6 +10,7 @@
 
 from collections.abc import Iterable
 
+import bordado as bd
 import numpy as np
 import numpy.testing as npt
 import pytest
@@ -37,7 +38,7 @@ def fixture_coordinates(region):
     Return a set of sample coordinates at zero height.
     """
     shape = (40, 40)
-    return vd.grid_coordinates(region=region, shape=shape, extra_coords=0)
+    return bd.grid_coordinates(region=region, shape=shape, non_dimensional_coords=0)
 
 
 @pytest.fixture(name="points")
@@ -45,7 +46,9 @@ def fixture_points(region):
     """
     Return the coordinates of some sample point masses.
     """
-    points = vd.grid_coordinates(region=region, shape=(6, 6), extra_coords=-1e3)
+    points = bd.grid_coordinates(
+        region=region, shape=(6, 6), non_dimensional_coords=-1e3
+    )
     return points
 
 
@@ -79,7 +82,7 @@ def fixture_coordinates_small(region):
     Return a small set of 25 coordinates and variable elevation.
     """
     shape = (5, 5)
-    easting, northing = vd.grid_coordinates(region=region, shape=shape)
+    easting, northing = bd.grid_coordinates(region=region, shape=shape)
     upward = np.arange(25, dtype=float).reshape(shape)
     coordinates = (easting, northing, upward)
     return coordinates
@@ -99,7 +102,7 @@ def fixture_coordinates_9x9(region):
     Return a small set of 81 coordinates and variable elevation.
     """
     shape = (9, 9)
-    easting, northing = vd.grid_coordinates(region, shape=shape)
+    easting, northing = bd.grid_coordinates(region, shape=shape)
     upward = np.arange(shape[0] * shape[1], dtype=float).reshape(shape)
     coordinates = (easting, northing, upward)
     return coordinates
@@ -146,7 +149,9 @@ def test_equivalent_sources_cartesian(
     # to synthetic values
     upward = 0
     shape = (60, 60)
-    grid_coords = vd.grid_coordinates(region=region, shape=shape, extra_coords=upward)
+    grid_coords = bd.grid_coordinates(
+        region=region, shape=shape, non_dimensional_coords=upward
+    )
     true = point_gravity(grid_coords, points, masses, field="g_z")
     npt.assert_allclose(true, eqs.predict(grid_coords), atol=atol)
 
@@ -170,7 +175,9 @@ def test_equivalent_sources_small_data_cartesian(region, points, masses):
     Use Cartesian coordinates.
     """
     # Define a small set of observation points
-    coordinates = vd.grid_coordinates(region=region, shape=(8, 8), extra_coords=0)
+    coordinates = bd.grid_coordinates(
+        region=region, shape=(8, 8), non_dimensional_coords=0
+    )
     # Get synthetic data
     data = point_gravity(coordinates, points, masses, field="g_z")
 
@@ -188,7 +195,9 @@ def test_equivalent_sources_small_data_cartesian(region, points, masses):
     # to synthetic values
     upward = 20
     shape = (8, 8)
-    grid_coords = vd.grid_coordinates(region=region, shape=shape, extra_coords=upward)
+    grid_coords = bd.grid_coordinates(
+        region=region, shape=shape, non_dimensional_coords=upward
+    )
     true = point_gravity(grid_coords, points, masses, field="g_z")
     npt.assert_allclose(true, eqs.predict(grid_coords), rtol=0.08)
 
@@ -279,7 +288,8 @@ def test_equivalent_sources_custom_points_cartesian(region, coordinates, data):
     """
     # Pass a custom set of point sources
     points_custom = tuple(
-        i.ravel() for i in vd.grid_coordinates(region, shape=(3, 3), extra_coords=-550)
+        i.ravel()
+        for i in bd.grid_coordinates(region, shape=(3, 3), non_dimensional_coords=-550)
     )
     eqs = EquivalentSources(points=points_custom)
     eqs.fit(coordinates, data)
@@ -303,8 +313,8 @@ def test_equivalent_sources_jacobian_cartesian():
     Test Jacobian matrix under symmetric system of point sources.
     Use Cartesian coordinates.
     """
-    easting, northing, upward = vd.grid_coordinates(
-        region=[-100, 100, -100, 100], shape=(2, 2), extra_coords=0
+    easting, northing, upward = bd.grid_coordinates(
+        region=[-100, 100, -100, 100], shape=(2, 2), non_dimensional_coords=0
     )
     points = vdb.n_1d_arrays((easting, northing, upward + 100), n=3)
     coordinates = vdb.n_1d_arrays((easting, northing, upward), n=3)
@@ -366,7 +376,7 @@ def test_dtype(
     # Define the points argument for EquivalentSources
     points = None
     if custom_points:
-        points = vd.grid_coordinates(region, spacing=300, extra_coords=-2e3)
+        points = bd.grid_coordinates(region, spacing=300, non_dimensional_coords=-2e3)
     # Define the points argument for EquivalentSources.fit()
     if weights_none:
         weights = None
@@ -393,11 +403,13 @@ def test_jacobian_dtype(region, dtype):
     """
     # Build a set of custom coordinates
     coordinates = tuple(
-        c.ravel() for c in vd.grid_coordinates(region, shape=(10, 10), extra_coords=0)
+        c.ravel()
+        for c in bd.grid_coordinates(region, shape=(10, 10), non_dimensional_coords=0)
     )
     # Create custom set of point sources
     points = tuple(
-        p.ravel() for p in vd.grid_coordinates(region, shape=(6, 6), extra_coords=-2e3)
+        p.ravel()
+        for p in bd.grid_coordinates(region, shape=(6, 6), non_dimensional_coords=-2e3)
     )
     # Ravel the coordinates
     coordinates = tuple(c.ravel() for c in coordinates)
@@ -425,7 +437,9 @@ def test_error_deprecated_args(coordinates_small, data_small, region, deprecated
     # Define sample equivalent sources and fit against synthetic data
     eqs = EquivalentSources().fit(coordinates_small, data_small)
     # Build a target grid
-    grid_coords = vd.grid_coordinates(region=region, shape=(4, 4), extra_coords=2e3)
+    grid_coords = bd.grid_coordinates(
+        region=region, shape=(4, 4), non_dimensional_coords=2e3
+    )
     # Try to grid passing deprecated arguments
     msg = "The 'upward', 'region', 'shape' and 'spacing' arguments have been"
     with pytest.raises(ValueError, match=msg):
@@ -439,7 +453,9 @@ def test_error_ignored_args(coordinates_small, data_small, region):
     # Define sample equivalent sources and fit against synthetic data
     eqs = EquivalentSources().fit(coordinates_small, data_small)
     # Build a target grid
-    grid_coords = vd.grid_coordinates(region=region, shape=(4, 4), extra_coords=2e3)
+    grid_coords = bd.grid_coordinates(
+        region=region, shape=(4, 4), non_dimensional_coords=2e3
+    )
     # Try to grid passing kwarg arguments that will be ignored
     msg = "The 'bla' arguments are being ignored."
     with pytest.warns(FutureWarning, match=msg):
diff --git a/test/test_eq_sources_spherical.py b/test/test_eq_sources_spherical.py
index a8c783f95..b12509a1c 100644
--- a/test/test_eq_sources_spherical.py
+++ b/test/test_eq_sources_spherical.py
@@ -8,6 +8,7 @@
 Test the EquivalentSourcesSph gridder.
 """
 
+import bordado as bd
 import numpy as np
 import numpy.testing as npt
 import pytest
@@ -32,7 +33,9 @@ def fixture_points(region):
     Return the coordinates of some sample point masses.
     """
     radius = 6400e3
-    points = vd.grid_coordinates(region=region, shape=(6, 6), extra_coords=radius - 1e3)
+    points = bd.grid_coordinates(
+        region=region, shape=(6, 6), non_dimensional_coords=radius - 1e3
+    )
     return points
 
 
@@ -50,7 +53,7 @@ def fixture_coordinates_small(region):
     Return a small set of 25 coordinates and variable elevation.
     """
     shape = (5, 5)
-    longitude, latitude = vd.grid_coordinates(region=region, shape=shape)
+    longitude, latitude = bd.grid_coordinates(region=region, shape=shape)
     radius = 6400e3 + np.arange(25, dtype=float).reshape(shape)
     coordinates = (longitude, latitude, radius)
     return coordinates
@@ -75,13 +78,13 @@ def test_equivalent_sources_spherical():  # pragma: no cover
     region = (-70, -60, -40, -30)
     radius = 6400e3
     # Build synthetic point masses
-    points = vd.grid_coordinates(
-        region=region, shape=(6, 6), extra_coords=radius - 500e3
+    points = bd.grid_coordinates(
+        region=region, shape=(6, 6), non_dimensional_coords=radius - 500e3
     )
     masses = vd.synthetic.CheckerBoard(amplitude=1e13, region=region).predict(points)
     # Define a set of observation points
-    coordinates = vd.grid_coordinates(
-        region=region, shape=(40, 40), extra_coords=radius
+    coordinates = bd.grid_coordinates(
+        region=region, shape=(40, 40), non_dimensional_coords=radius
     )
     # Get synthetic data
     data = point_gravity(
@@ -97,7 +100,9 @@ def test_equivalent_sources_spherical():  # pragma: no cover
     # to synthetic values
     upward = radius
     shape = (60, 60)
-    grid_coords = vd.grid_coordinates(region=region, shape=shape, extra_coords=upward)
+    grid_coords = bd.grid_coordinates(
+        region=region, shape=shape, non_dimensional_coords=upward
+    )
     true = point_gravity(
         grid_coords, points, masses, field="g_z", coordinate_system="spherical"
     )
@@ -116,12 +121,14 @@ def test_equivalent_sources_small_data_spherical():
     region = (-70, -60, -40, -30)
     radius = 6400e3
     # Build synthetic point masses
-    points = vd.grid_coordinates(
-        region=region, shape=(6, 6), extra_coords=radius - 500e3
+    points = bd.grid_coordinates(
+        region=region, shape=(6, 6), non_dimensional_coords=radius - 500e3
     )
     masses = vd.synthetic.CheckerBoard(amplitude=1e13, region=region).predict(points)
     # Define a set of observation points
-    coordinates = vd.grid_coordinates(region=region, shape=(8, 8), extra_coords=radius)
+    coordinates = bd.grid_coordinates(
+        region=region, shape=(8, 8), non_dimensional_coords=radius
+    )
     # Get synthetic data
     data = point_gravity(
         coordinates, points, masses, field="g_z", coordinate_system="spherical"
@@ -141,7 +148,9 @@ def test_equivalent_sources_small_data_spherical():
     # to synthetic values
     upward = radius + 2e3
     shape = (8, 8)
-    grid_coords = vd.grid_coordinates(region=region, shape=shape, extra_coords=upward)
+    grid_coords = bd.grid_coordinates(
+        region=region, shape=shape, non_dimensional_coords=upward
+    )
     true = point_gravity(
         grid_coords, points, masses, field="g_z", coordinate_system="spherical"
     )
@@ -160,12 +169,14 @@ def test_equivalent_sources_custom_points_spherical():
     region = (-70, -60, -40, -30)
     radius = 6400e3
     # Build synthetic point masses
-    points = vd.grid_coordinates(
-        region=region, shape=(6, 6), extra_coords=radius - 500e3
+    points = bd.grid_coordinates(
+        region=region, shape=(6, 6), non_dimensional_coords=radius - 500e3
     )
     masses = vd.synthetic.CheckerBoard(amplitude=1e13, region=region).predict(points)
     # Define a set of observation points
-    coordinates = vd.grid_coordinates(region=region, shape=(5, 5), extra_coords=radius)
+    coordinates = bd.grid_coordinates(
+        region=region, shape=(5, 5), non_dimensional_coords=radius
+    )
     # Get synthetic data
     data = point_gravity(
         coordinates, points, masses, field="g_z", coordinate_system="spherical"
@@ -174,8 +185,8 @@ def test_equivalent_sources_custom_points_spherical():
     # Pass a custom set of point sources
     points_custom = tuple(
         i.ravel()
-        for i in vd.grid_coordinates(
-            region=region, shape=(3, 3), extra_coords=radius - 500e3
+        for i in bd.grid_coordinates(
+            region=region, shape=(3, 3), non_dimensional_coords=radius - 500e3
         )
     )
     eqs = EquivalentSourcesSph(points=points_custom)
@@ -220,13 +231,13 @@ def test_equivalent_sources_spherical_parallel():  # pragma: no cover
     region = (-70, -60, -40, -30)
     radius = 6400e3
     # Build synthetic point masses
-    points = vd.grid_coordinates(
-        region=region, shape=(6, 6), extra_coords=radius - 500e3
+    points = bd.grid_coordinates(
+        region=region, shape=(6, 6), non_dimensional_coords=radius - 500e3
     )
     masses = vd.synthetic.CheckerBoard(amplitude=1e13, region=region).predict(points)
     # Define a set of observation points
-    coordinates = vd.grid_coordinates(
-        region=region, shape=(40, 40), extra_coords=radius
+    coordinates = bd.grid_coordinates(
+        region=region, shape=(40, 40), non_dimensional_coords=radius
     )
     # Get synthetic data
     data = point_gravity(
@@ -242,7 +253,9 @@ def test_equivalent_sources_spherical_parallel():  # pragma: no cover
 
     upward = radius
     shape = (60, 60)
-    grid_coords = vd.grid_coordinates(region=region, shape=shape, extra_coords=upward)
+    grid_coords = bd.grid_coordinates(
+        region=region, shape=shape, non_dimensional_coords=upward
+    )
     grid_serial = eqs_serial.grid(grid_coords)
     grid_parallel = eqs_parallel.grid(grid_coords)
     npt.assert_allclose(grid_serial.scalars, grid_parallel.scalars, rtol=1e-7)
@@ -264,7 +277,9 @@ def test_error_deprecated_args(coordinates_small, data_small, region, deprecated
     # Define sample equivalent sources and fit against synthetic data
     eqs = EquivalentSourcesSph().fit(coordinates_small, data_small)
     # Build a target grid
-    grid_coords = vd.grid_coordinates(region=region, shape=(4, 4), extra_coords=2e3)
+    grid_coords = bd.grid_coordinates(
+        region=region, shape=(4, 4), non_dimensional_coords=2e3
+    )
     # Try to grid passing deprecated arguments
     msg = "The 'upward', 'region', 'shape' and 'spacing' arguments have been"
     with pytest.raises(ValueError, match=msg):
@@ -278,7 +293,9 @@ def test_error_ignored_args(coordinates_small, data_small, region):
     # Define sample equivalent sources and fit against synthetic data
     eqs = EquivalentSourcesSph().fit(coordinates_small, data_small)
     # Build a target grid
-    grid_coords = vd.grid_coordinates(region=region, shape=(4, 4), extra_coords=2e3)
+    grid_coords = bd.grid_coordinates(
+        region=region, shape=(4, 4), non_dimensional_coords=2e3
+    )
     # Try to grid passing kwarg arguments that will be ignored
     msg = "The 'bla' arguments are being ignored."
     with pytest.warns(FutureWarning, match=msg):
diff --git a/test/test_eq_sources_utils.py b/test/test_eq_sources_utils.py
index a86a5020d..b3b1682b3 100644
--- a/test/test_eq_sources_utils.py
+++ b/test/test_eq_sources_utils.py
@@ -12,7 +12,7 @@
 
 import numpy as np
 import pytest
-from verde import scatter_points
+import verde as vd
 
 from harmonica._equivalent_sources.utils import cast_fit_input, pop_extra_coords
 
@@ -42,7 +42,7 @@ def test_cast_fit_input(weights_none, dtype):
     Test cast_fit_input function.
     """
     region = (-7e3, 4e3, 10e3, 25e3)
-    coordinates = scatter_points(region=region, size=100, random_state=42)
+    coordinates = vd.scatter_points(region=region, size=100, random_state=42)
     data = np.arange(coordinates[0].size, dtype="float64")
     weights = None if weights_none else np.ones_like(data)
     coordinates, data, weights = cast_fit_input(coordinates, data, weights, dtype)
diff --git a/test/test_euler_methods.py b/test/test_euler_methods.py
index 6e098ae4c..8dfd2d771 100644
--- a/test/test_euler_methods.py
+++ b/test/test_euler_methods.py
@@ -4,6 +4,7 @@
 #
 # This code is part of the Fatiando a Terra project (https://www.fatiando.org)
 #
+import bordado as bd
 import numpy as np
 import numpy.testing as npt
 import pytest
@@ -49,7 +50,7 @@ def test_euler_inversion_with_numeric_derivatives(euler):
     inc, dec = -40, 15
     dipole_moments = magnetic_angles_to_vec(1.0e14, inc, dec)
     region = [-100e3, 100e3, -80e3, 80e3]
-    coordinates = vd.grid_coordinates(region, spacing=500, extra_coords=500)
+    coordinates = bd.grid_coordinates(region, spacing=500, non_dimensional_coords=500)
     b = dipole_magnetic(coordinates, dipole_coordinates, dipole_moments, field="b")
     # Add a fixed base level
     true_base_level = 200  # nT
@@ -88,7 +89,7 @@ def test_euler_inversion_with_analytic_derivatives(euler):
     masses_coordinates = (10e3, 15e3, -7e3)
     masses = 1.0e14
     region = [-100e3, 100e3, -80e3, 80e3]
-    coordinates = vd.grid_coordinates(region, spacing=500, extra_coords=500)
+    coordinates = bd.grid_coordinates(region, spacing=500, non_dimensional_coords=500)
     gz = point_gravity(coordinates, masses_coordinates, masses, field="g_z")
     random = np.random.default_rng(42)
     gz += random.normal(0, 0.01 * np.abs(gz).max(), gz.shape)
@@ -136,7 +137,7 @@ def test_euler_inversion_convergence_warning(euler):
     inc, dec = -40, 15
     dipole_moments = magnetic_angles_to_vec(1.0e14, inc, dec)
     region = [-100e3, 100e3, -80e3, 80e3]
-    coordinates = vd.grid_coordinates(region, spacing=500, extra_coords=500)
+    coordinates = bd.grid_coordinates(region, spacing=500, non_dimensional_coords=500)
     b = dipole_magnetic(coordinates, dipole_coordinates, dipole_moments, field="b")
     # Add a fixed base level
     true_base_level = 200  # nT
@@ -164,7 +165,7 @@ def test_euler_deconvolution_with_numeric_derivatives():
     inc, dec = -40, 15
     dipole_moments = magnetic_angles_to_vec(1.0e14, inc, dec)
     region = [-100e3, 100e3, -80e3, 80e3]
-    coordinates = vd.grid_coordinates(region, spacing=500, extra_coords=500)
+    coordinates = bd.grid_coordinates(region, spacing=500, non_dimensional_coords=500)
     b = dipole_magnetic(coordinates, dipole_coordinates, dipole_moments, field="b")
     # Add a fixed base level
     true_base_level = 200  # nT
@@ -193,7 +194,7 @@ def test_euler_deconvolution_with_analytic_derivatives():
     masses_coordinates = (10e3, 15e3, -10e3)
     masses = 1.0e12
     region = [-100e3, 100e3, -80e3, 80e3]
-    coordinates = vd.grid_coordinates(region, spacing=500, extra_coords=500)
+    coordinates = bd.grid_coordinates(region, spacing=500, non_dimensional_coords=500)
     gz = point_gravity(coordinates, masses_coordinates, masses, field="g_z")
     # Convert Eötvös to mGal because derivatives must be in mGal/m
     eotvos2mgal = 1.0e-4
diff --git a/test/test_filters.py b/test/test_filters.py
index 07ff4b9ce..132db4652 100644
--- a/test/test_filters.py
+++ b/test/test_filters.py
@@ -10,12 +10,13 @@
 
 import re
 
+import bordado as bd
 import numpy as np
 import numpy.testing as npt
 import pytest
+import verde as vd
 import xarray as xr
 import xarray.testing as xrt
-from verde import grid_coordinates, make_xarray_grid
 
 from harmonica.filters._fft import fft, ifft
 from harmonica.filters._filters import (
@@ -47,9 +48,9 @@ def fixture_sample_grid(region):
     """
     Return a sample grid as an :class:`xarray.DataArray`.
     """
-    easting, northing = grid_coordinates(region, spacing=500)
+    easting, northing = bd.grid_coordinates(region, spacing=500)
     data = np.sin(easting / 1e3) + np.cos(northing / 1e3)
-    return make_xarray_grid((easting, northing), data, data_names="sample").sample
+    return vd.make_xarray_grid((easting, northing), data, data_names="sample").sample
 
 
 @pytest.fixture(name="sample_grid_multiple_coords")
@@ -57,7 +58,9 @@ def fixture_sample_grid_multiple_coords(region):
     """
     Return a sample grid as an :class:`xarray.DataArray` with multiple coords.
     """
-    easting, northing, upward = grid_coordinates(region, spacing=500, extra_coords=100)
+    easting, northing, upward = bd.grid_coordinates(
+        region, spacing=500, non_dimensional_coords=100
+    )
     data = np.sin(easting / 1e3) + np.cos(northing / 1e3)
     easting, northing = np.unique(easting), np.unique(northing)
     easting_km, northing_km = easting * 1e-3, northing * 1e-3
@@ -112,9 +115,9 @@ def fixture_sample_fft_grid():
     Return a sample fft_grid to be used in test functions.
     """
     domain = (-9e-4, 9e-4, -8e-4, 8e-4)
-    freq_easting, freq_northing = grid_coordinates(region=domain, spacing=8e-4)
+    freq_easting, freq_northing = bd.grid_coordinates(region=domain, spacing=8e-4)
     dummy_fft = np.ones_like(freq_easting)
-    fft_grid = make_xarray_grid(
+    fft_grid = vd.make_xarray_grid(
         (freq_easting, freq_northing),
         dummy_fft,
         data_names=["sample_fft"],
diff --git a/test/test_gradient_boosted_eqs.py b/test/test_gradient_boosted_eqs.py
index fd72f4ea2..5038bd6eb 100644
--- a/test/test_gradient_boosted_eqs.py
+++ b/test/test_gradient_boosted_eqs.py
@@ -8,6 +8,7 @@
 Test functions for gradient-boosted equivalent sources.
 """
 
+import bordado as bd
 import numpy as np
 import numpy.testing as npt
 import pytest
@@ -33,7 +34,7 @@ def fixture_coordinates(region):
     Return a set of sample coordinates at zero height.
     """
     shape = (40, 40)
-    return vd.grid_coordinates(region=region, shape=shape, extra_coords=0)
+    return bd.grid_coordinates(region=region, shape=shape, non_dimensional_coords=0)
 
 
 @pytest.fixture(name="points")
@@ -41,7 +42,9 @@ def fixture_points(region):
     """
     Return the coordinates of some sample point masses.
     """
-    points = vd.grid_coordinates(region=region, shape=(6, 6), extra_coords=-1e3)
+    points = bd.grid_coordinates(
+        region=region, shape=(6, 6), non_dimensional_coords=-1e3
+    )
     return points
 
 
@@ -75,7 +78,7 @@ def fixture_coordinates_small(region):
     Return a small set of 25 coordinates and variable elevation.
     """
     shape = (8, 8)
-    return vd.grid_coordinates(region=region, shape=shape, extra_coords=0)
+    return bd.grid_coordinates(region=region, shape=shape, non_dimensional_coords=0)
 
 
 @pytest.fixture(name="data_small")
@@ -102,8 +105,8 @@ def test_get_region_data_sources(data_region, sources_region, expected_region):
     Test the EquivalentSourcesGB._get_region_data_sources method.
     """
     shape = (8, 10)
-    coordinates = vd.grid_coordinates(data_region, shape=shape)
-    points = vd.grid_coordinates(sources_region, shape=shape)
+    coordinates = bd.grid_coordinates(data_region, shape=shape)
+    points = bd.grid_coordinates(sources_region, shape=shape)
     region = _get_region_data_sources(coordinates, points)
     npt.assert_allclose(region, expected_region)
 
@@ -115,7 +118,9 @@ def test_custom_points(region, coordinates_small, data_small):
     # Pass a custom set of point sources
     points_custom = tuple(
         i.ravel()
-        for i in vd.grid_coordinates(region=region, shape=(3, 3), extra_coords=-550)
+        for i in bd.grid_coordinates(
+            region=region, shape=(3, 3), non_dimensional_coords=-550
+        )
     )
     eqs = EquivalentSourcesGB(points=points_custom, window_size=500, depth=500)
     eqs.fit(coordinates_small, data_small)
@@ -131,7 +136,9 @@ def test_memory_estimation(spacing, window_size, dtype, itemsize):
     Test the estimate_required_memory class method.
     """
     region = (-1e3, 5e3, 2e3, 8e3)
-    coordinates = vd.grid_coordinates(region=region, spacing=spacing, extra_coords=0)
+    coordinates = bd.grid_coordinates(
+        region=region, spacing=spacing, non_dimensional_coords=0
+    )
     # Compute expected required memory
     sources_p_window = (int(window_size / spacing) + 1) ** 2
     data_p_window = (int(window_size / spacing) + 1) ** 2
@@ -181,7 +188,7 @@ def test_gradient_boosted_eqs_single_window(
     npt.assert_allclose(data, eqs.predict(coordinates), rtol=1e-5, atol=atol)
     # Gridding onto a denser grid should be reasonably accurate when compared
     # to synthetic values
-    grid = vd.grid_coordinates(region, shape=(60, 60), extra_coords=0)
+    grid = bd.grid_coordinates(region, shape=(60, 60), non_dimensional_coords=0)
     true = point_gravity(grid, points, masses, field="g_z")
     npt.assert_allclose(true, eqs.predict(grid), rtol=1e-3, atol=atol)
 
@@ -203,7 +210,7 @@ def test_gradient_boosted_eqs_predictions(
 
     # Gridding onto a denser grid should be reasonably accurate when compared
     # to synthetic values
-    grid = vd.grid_coordinates(region, shape=(60, 60), extra_coords=0)
+    grid = bd.grid_coordinates(region, shape=(60, 60), non_dimensional_coords=0)
     true = point_gravity(grid, points, masses, field="g_z")
     # Error tolerance is 2% of the maximum data.
     npt.assert_allclose(true, eqs.predict(grid), rtol=0, atol=0.02 * vd.maxabs(true))
@@ -231,11 +238,13 @@ def test_same_number_of_windows_data_and_sources():
     spacing = 1
     # Create data points on a large region
     region = (1, 3, 1, 3)
-    coordinates = vd.grid_coordinates(region=region, spacing=spacing, extra_coords=0)
+    coordinates = bd.grid_coordinates(
+        region=region, spacing=spacing, non_dimensional_coords=0
+    )
     # Create source points on a smaller region
     sources_region = (1.5, 2.5, 1.5, 2.5)
-    points = vd.grid_coordinates(
-        region=sources_region, spacing=spacing, extra_coords=-10
+    points = bd.grid_coordinates(
+        region=sources_region, spacing=spacing, non_dimensional_coords=-10
     )
     # Create EquivalentSourcesGB
     eqs = EquivalentSourcesGB(window_size=spacing)
@@ -254,11 +263,13 @@ def test_same_windows_data_and_sources():
     spacing = 1
     # Create data points on a large region
     region = (1, 3, 1, 3)
-    coordinates = vd.grid_coordinates(region=region, spacing=spacing, extra_coords=0)
+    coordinates = bd.grid_coordinates(
+        region=region, spacing=spacing, non_dimensional_coords=0
+    )
     # Create source points on a subregion
     sources_region = (1, 2, 1, 3)
-    points = vd.grid_coordinates(
-        region=sources_region, spacing=spacing, extra_coords=-10
+    points = bd.grid_coordinates(
+        region=sources_region, spacing=spacing, non_dimensional_coords=-10
     )
     # Create EquivalentSourcesGB
     eqs = EquivalentSourcesGB(window_size=spacing)
@@ -308,7 +319,7 @@ def test_dtype(
     # Define the points argument for EquivalentSources
     points = None
     if custom_points:
-        points = vd.grid_coordinates(region, spacing=300, extra_coords=-2e3)
+        points = bd.grid_coordinates(region, spacing=300, non_dimensional_coords=-2e3)
     # Define the points argument for EquivalentSources.fit()
     if weights_none:
         weights = None
@@ -361,7 +372,9 @@ def test_error_deprecated_args(coordinates_small, data_small, region, deprecated
     # Define sample equivalent sources and fit against synthetic data
     eqs = EquivalentSourcesGB(window_size=500).fit(coordinates_small, data_small)
     # Build a target grid
-    grid_coords = vd.grid_coordinates(region=region, shape=(4, 4), extra_coords=2e3)
+    grid_coords = bd.grid_coordinates(
+        region=region, shape=(4, 4), non_dimensional_coords=2e3
+    )
     # Try to grid passing deprecated arguments
     msg = "The 'upward', 'region', 'shape' and 'spacing' arguments have been"
     with pytest.raises(ValueError, match=msg):
@@ -375,7 +388,9 @@ def test_error_ignored_args(coordinates_small, data_small, region):
     # Define sample equivalent sources and fit against synthetic data
     eqs = EquivalentSourcesGB(window_size=500).fit(coordinates_small, data_small)
     # Build a target grid
-    grid_coords = vd.grid_coordinates(region=region, shape=(4, 4), extra_coords=2e3)
+    grid_coords = bd.grid_coordinates(
+        region=region, shape=(4, 4), non_dimensional_coords=2e3
+    )
     # Try to grid passing kwarg arguments that will be ignored
     msg = "The 'bla' arguments are being ignored."
     with pytest.warns(FutureWarning, match=msg):
@@ -384,7 +399,9 @@ def test_error_ignored_args(coordinates_small, data_small, region):
 
 def test_window_size_less_than_5000():
     region = (0, 10e3, -5e3, 5e3)
-    grid_coords = vd.grid_coordinates(region=region, shape=(64, 64), extra_coords=0)
+    grid_coords = bd.grid_coordinates(
+        region=region, shape=(64, 64), non_dimensional_coords=0
+    )
     grid_coords = [c.ravel() for c in grid_coords]
     eqs = EquivalentSourcesGB()
     eqs.points_ = eqs._build_points(
@@ -404,7 +421,9 @@ def test_window_size_less_than_5000():
 
 def test_window_size():
     region = (0, 10e3, -5e3, 5e3)
-    grid_coords = vd.grid_coordinates(region=region, shape=(100, 100), extra_coords=0)
+    grid_coords = bd.grid_coordinates(
+        region=region, shape=(100, 100), non_dimensional_coords=0
+    )
     eqs = EquivalentSourcesGB()
     eqs.points_ = eqs._build_points(
         grid_coords
diff --git a/test/test_igrf.py b/test/test_igrf.py
index bf124db52..a11c3596e 100644
--- a/test/test_igrf.py
+++ b/test/test_igrf.py
@@ -9,10 +9,10 @@
 import datetime
 import pathlib
 
+import bordado as bd
 import numpy as np
 import numpy.testing as npt
 import pytest
-import verde as vd
 
 from harmonica._spherical_harmonics.igrf import (
     IGRF14,
@@ -247,7 +247,9 @@ def test_igrf_grid(date):
     region = (0, 360, -90, 90)
     spacing = 20
     height = 2043
-    coordinates = vd.grid_coordinates(region, spacing=spacing, extra_coords=height)
+    coordinates = bd.grid_coordinates(
+        region, spacing=spacing, non_dimensional_coords=height
+    )
     igrf = IGRF14(date)
     predicted = igrf.predict(coordinates)
     grid = igrf.grid(region, height, spacing=spacing)
diff --git a/test/test_point_gravity.py b/test/test_point_gravity.py
index 82516d76d..80b5a014f 100644
--- a/test/test_point_gravity.py
+++ b/test/test_point_gravity.py
@@ -11,10 +11,10 @@
 import re
 from pathlib import Path
 
+import bordado as bd
 import numpy as np
 import numpy.testing as npt
 import pytest
-import verde as vd
 from choclo.point import (
     gravity_e,
     gravity_ee,
@@ -335,9 +335,13 @@ def test_point_mass_cartesian_parallel(
     Check if parallel and serial runs return the same result.
     """
     region = (2e3, 10e3, -3e3, 5e3)
-    points = vd.scatter_points(region, size=30, extra_coords=-1e3, random_state=0)
+    points = bd.random_coordinates(
+        region, size=30, non_dimensional_coords=-1e3, random_seed=0
+    )
     masses = np.arange(points[0].size)
-    coordinates = vd.grid_coordinates(region=region, spacing=1e3, extra_coords=0)
+    coordinates = bd.grid_coordinates(
+        region=region, spacing=1e3, non_dimensional_coords=0
+    )
     result_serial = point_gravity(
         coordinates, points, masses, field=field, parallel=False
     )
@@ -355,9 +359,13 @@ def test_laplace_equation_cartesian():
     Use Cartesian coordinates.
     """
     region = (2e3, 10e3, -3e3, 5e3)
-    points = vd.scatter_points(region, size=30, extra_coords=-1e3, random_state=0)
+    points = bd.random_coordinates(
+        region, size=30, non_dimensional_coords=-1e3, random_seed=0
+    )
     masses = np.arange(points[0].size)
-    coordinates = vd.grid_coordinates(region=region, spacing=1e3, extra_coords=0)
+    coordinates = bd.grid_coordinates(
+        region=region, spacing=1e3, non_dimensional_coords=0
+    )
     g_ee = point_gravity(coordinates, points, masses, field="g_ee")
     g_nn = point_gravity(coordinates, points, masses, field="g_nn")
     g_zz = point_gravity(coordinates, points, masses, field="g_zz")
@@ -374,9 +382,13 @@ def test_tensor_non_diagonal_components(field, flipped_field):
     Check if function computes g_xy as the same as g_yx.
     """
     region = (2e3, 10e3, -3e3, 5e3)
-    points = vd.scatter_points(region, size=30, extra_coords=-1e3, random_state=0)
+    points = bd.random_coordinates(
+        region, size=30, non_dimensional_coords=-1e3, random_seed=0
+    )
     masses = np.arange(points[0].size)
-    coordinates = vd.grid_coordinates(region=region, spacing=1e3, extra_coords=0)
+    coordinates = bd.grid_coordinates(
+        region=region, spacing=1e3, non_dimensional_coords=0
+    )
     npt.assert_allclose(
         point_gravity(coordinates, points, masses, field=field),
         point_gravity(coordinates, points, masses, field=flipped_field),
@@ -763,11 +775,13 @@ def test_point_mass_spherical_parallel():  # pragma: no cover
     """
     region = (2, 10, -3, 5)
     radius = 6400e3
-    points = vd.scatter_points(
-        region, size=30, extra_coords=radius - 10e3, random_state=0
+    points = bd.random_coordinates(
+        region, size=30, non_dimensional_coords=radius - 10e3, random_seed=0
     )
     masses = np.arange(points[0].size)
-    coordinates = vd.grid_coordinates(region=region, spacing=1, extra_coords=radius)
+    coordinates = bd.grid_coordinates(
+        region=region, spacing=1, non_dimensional_coords=radius
+    )
     for field in ("potential", "g_z"):
         result_serial = point_gravity(
             coordinates,
diff --git a/test/test_prism.py b/test/test_prism.py
index 077e1470f..1abe20c87 100644
--- a/test/test_prism.py
+++ b/test/test_prism.py
@@ -11,10 +11,10 @@
 import re
 from unittest.mock import patch
 
+import bordado as bd
 import numpy as np
 import numpy.testing as npt
 import pytest
-import verde as vd
 from choclo.prism import (
     gravity_e,
     gravity_ee,
@@ -122,8 +122,8 @@ def test_forward_with_null_prisms():
     Test if the forward model with null prisms gives sensible results.
     """
     # Create a set of observation points
-    coordinates = vd.grid_coordinates(
-        region=(-50, 50, -50, 50), shape=(3, 3), extra_coords=0
+    coordinates = bd.grid_coordinates(
+        region=(-50, 50, -50, 50), shape=(3, 3), non_dimensional_coords=0
     )
     # Build a set of prisms that includes null ones (no volume or zero density)
     prisms = [
@@ -254,7 +254,7 @@ def test_laplace():  # pragma: no cover
     Test if the diagonal components satisfy Laplace equation.
     """
     region = (-10e3, 10e3, -10e3, 10e3)
-    coords = vd.grid_coordinates(region, shape=(10, 10), extra_coords=300)
+    coords = bd.grid_coordinates(region, shape=(10, 10), non_dimensional_coords=300)
     prisms = [
         [1e3, 7e3, -5e3, 2e3, -1e3, -500],
         [-4e3, 1e3, 4e3, 10e3, -2e3, 200],
@@ -313,8 +313,8 @@ def test_prisms_parallel_vs_serial():
         [0, 100, 0, 100, -10, 0],
     ]
     densities = [2000, 3000, 4000, 5000]
-    coordinates = vd.grid_coordinates(
-        region=(-100, 100, -100, 100), shape=(3, 3), extra_coords=10
+    coordinates = bd.grid_coordinates(
+        region=(-100, 100, -100, 100), shape=(3, 3), non_dimensional_coords=10
     )
     for field in ("potential", "g_z"):
         result_parallel = prism_gravity(
@@ -346,8 +346,8 @@ def densities(self):
     @pytest.fixture
     def coordinates(self):
         """Sample coordinates."""
-        coordinates = vd.grid_coordinates(
-            region=(-100, 100, -100, 100), spacing=20, extra_coords=10
+        coordinates = bd.grid_coordinates(
+            region=(-100, 100, -100, 100), spacing=20, non_dimensional_coords=10
         )
         return coordinates
 
diff --git a/test/test_prism_layer.py b/test/test_prism_layer.py
index acd560dd0..7625db419 100644
--- a/test/test_prism_layer.py
+++ b/test/test_prism_layer.py
@@ -11,10 +11,10 @@
 import warnings
 from unittest.mock import patch
 
+import bordado as bd
 import numpy as np
 import numpy.testing as npt
 import pytest
-import verde as vd
 import xarray as xr
 
 from harmonica import prism_gravity, prism_layer
@@ -345,7 +345,9 @@ def test_prism_layer_gravity(field, dummy_layer):
     """
     Check if gravity method works as expected.
     """
-    coordinates = vd.grid_coordinates((1, 3, 7, 10), spacing=1, extra_coords=30.0)
+    coordinates = bd.grid_coordinates(
+        (1, 3, 7, 10), spacing=1, non_dimensional_coords=30.0
+    )
     (easting, northing), surface, reference, density = dummy_layer
     layer = prism_layer(
         (easting, northing), surface, reference, properties={"density": density}
@@ -367,7 +369,9 @@ def test_prism_layer_gravity_surface_nans(field, dummy_layer, prism_layer_with_h
     """
     Check if gravity method works as expected when surface has nans.
     """
-    coordinates = vd.grid_coordinates((1, 3, 7, 10), spacing=1, extra_coords=30.0)
+    coordinates = bd.grid_coordinates(
+        (1, 3, 7, 10), spacing=1, non_dimensional_coords=30.0
+    )
     (easting, northing), surface, reference, density = dummy_layer
     # Create one layer that has nans on the surface array
     surface_w_nans = surface.copy()
@@ -391,7 +395,9 @@ def test_prism_layer_gravity_density_nans(field, dummy_layer, prism_layer_with_h
     """
     Check if prisms is ignored after a nan is found in density array.
     """
-    coordinates = vd.grid_coordinates((1, 3, 7, 10), spacing=1, extra_coords=30.0)
+    coordinates = bd.grid_coordinates(
+        (1, 3, 7, 10), spacing=1, non_dimensional_coords=30.0
+    )
     prisms_coords, surface, reference, density = dummy_layer
     # Create one layer that has nans on the density array
     indices = [(3, 3), (2, 1)]
@@ -480,7 +486,9 @@ def test_progress_bar(dummy_layer):
     """
     Check if forward gravity results with and without progress bar match.
     """
-    coordinates = vd.grid_coordinates((1, 3, 7, 10), spacing=1, extra_coords=30.0)
+    coordinates = bd.grid_coordinates(
+        (1, 3, 7, 10), spacing=1, non_dimensional_coords=30.0
+    )
     (easting, northing), surface, reference, density = dummy_layer
     layer = prism_layer(
         (easting, northing), surface, reference, properties={"density": density}
@@ -501,7 +509,9 @@ def test_numba_progress_missing_error(dummy_layer):
     Check if error is raised when progressbar=True and numba_progress package
     is not installed.
     """
-    coordinates = vd.grid_coordinates((1, 3, 7, 10), spacing=1, extra_coords=30.0)
+    coordinates = bd.grid_coordinates(
+        (1, 3, 7, 10), spacing=1, non_dimensional_coords=30.0
+    )
     (easting, northing), surface, reference, density = dummy_layer
     layer = prism_layer(
         (easting, northing), surface, reference, properties={"density": density}
@@ -515,7 +525,9 @@ def test_gravity_discarded_thin_prisms(dummy_layer):
     """
     Check if gravity of prism layer after discarding thin prisms is correct.
     """
-    coordinates = vd.grid_coordinates((1, 3, 7, 10), spacing=1, extra_coords=30.0)
+    coordinates = bd.grid_coordinates(
+        (1, 3, 7, 10), spacing=1, non_dimensional_coords=30.0
+    )
     prism_coords, surface, reference, density = dummy_layer
 
     layer = prism_layer(
diff --git a/test/test_prism_magnetic.py b/test/test_prism_magnetic.py
index 89f847ead..d8db86297 100644
--- a/test/test_prism_magnetic.py
+++ b/test/test_prism_magnetic.py
@@ -8,11 +8,11 @@
 Test forward functions for magnetic field of prisms.
 """
 
+import bordado as bd
 import choclo
 import numpy as np
 import numpy.testing as npt
 import pytest
-import verde as vd
 from choclo.prism import magnetic_field
 
 try:
@@ -47,8 +47,8 @@ def test_progress_bar(field):
         [0, 100, -100, 0, -10, 0],
     ]
     magnetizations = ([1.0, 1.0], [1.0, -1.0], [1.0, 5.0])
-    coordinates = vd.grid_coordinates(
-        region=(-100, 100, -100, 100), spacing=20, extra_coords=10
+    coordinates = bd.grid_coordinates(
+        region=(-100, 100, -100, 100), spacing=20, non_dimensional_coords=10
     )
     result_progress_true = prism_magnetic(
         coordinates, prisms, magnetizations, field=field, progressbar=True
@@ -104,8 +104,8 @@ def test_prisms_parallel_vs_serial(self, field):
             [1.0, -1.0, 1.0, 4.0],
             [1.0, 5.0, 3.0, 1.0],
         )
-        coordinates = vd.grid_coordinates(
-            region=(-100, 100, -100, 100), spacing=20, extra_coords=10
+        coordinates = bd.grid_coordinates(
+            region=(-100, 100, -100, 100), spacing=20, non_dimensional_coords=10
         )
         parallel = prism_magnetic(
             coordinates, prisms, magnetizations, field=field, parallel=True
diff --git a/test/test_tesseroid.py b/test/test_tesseroid.py
index daa9f66a2..9f4e2865e 100644
--- a/test/test_tesseroid.py
+++ b/test/test_tesseroid.py
@@ -11,12 +11,11 @@
 import re
 from unittest.mock import patch
 
+import bordado as bd
 import boule
 import numpy as np
 import numpy.testing as npt
 import pytest
-import verde as vd
-from verde import grid_coordinates
 
 from harmonica._forward._tesseroid_utils import (
     _discard_null_tesseroids,
@@ -691,7 +690,9 @@ def test_spherical_shell_two_dim_adaptive_discret(field):  # pragma: no cover
     # Define computation point located on the equator at the mean Earth radius
     ellipsoid = boule.WGS84
     radius = ellipsoid.mean_radius
-    coordinates = grid_coordinates([0, 350, -90, 90], spacing=10, extra_coords=radius)
+    coordinates = bd.grid_coordinates(
+        [0, 350, -90, 90], spacing=10, non_dimensional_coords=radius
+    )
     # Define lon and lat coordinates of spherical shell model made of
     # tesseroids
     shape = (12, 6)
@@ -788,8 +789,8 @@ def densities(self):
     @pytest.fixture
     def coordinates(self):
         """Sample coordinates."""
-        coordinates = vd.grid_coordinates(
-            region=(-15, 55, -80, 40), spacing=10, extra_coords=6.5e4
+        coordinates = bd.grid_coordinates(
+            region=(-15, 55, -80, 40), spacing=10, non_dimensional_coords=6.5e4
         )
         return coordinates
 
diff --git a/test/test_tesseroid_layer.py b/test/test_tesseroid_layer.py
index 04a35accb..6504aed27 100644
--- a/test/test_tesseroid_layer.py
+++ b/test/test_tesseroid_layer.py
@@ -11,11 +11,11 @@
 import warnings
 from unittest.mock import patch
 
+import bordado as bd
 import boule
 import numpy as np
 import numpy.testing as npt
 import pytest
-import verde as vd
 import xarray as xr
 
 from harmonica import tesseroid_gravity, tesseroid_layer
@@ -391,8 +391,8 @@ def test_tesseroid_layer_gravity(field, dummy_layer):
     """
     (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)
+    grid_coords = bd.grid_coordinates(
+        (-10, 10, -10, 10), spacing=7, non_dimensional_coords=(surface[0] + 10e3)
     )
     layer = tesseroid_layer(
         (longitude, latitude), surface, reference, properties={"density": density}
@@ -417,8 +417,8 @@ def test_tesseroid_layer_gravity_surface_nans(
     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)
+    grid_coords = bd.grid_coordinates(
+        (-10, 10, -10, 10), spacing=7, non_dimensional_coords=(surface[0] + 10e3)
     )
     # Create one layer that has nans on the surface array
     surface_w_nans = surface.copy()
@@ -448,8 +448,8 @@ def test_tesseroid_layer_gravity_density_nans(
     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)
+    grid_coords = bd.grid_coordinates(
+        (-10, 10, -10, 10), spacing=7, non_dimensional_coords=(surface[0] + 10e3)
     )
     # Create one layer that has nans on the density array
     indices = [(3, 3), (2, 1)]
@@ -477,8 +477,8 @@ def test_progress_bar(dummy_layer):
     Check if forward gravity results with and without progress bar match.
     """
     (longitude, latitude), surface, reference, density = dummy_layer
-    coordinates = vd.grid_coordinates(
-        (-10, 10, -10, 10), spacing=7, extra_coords=(surface[0] + 10e3)
+    coordinates = bd.grid_coordinates(
+        (-10, 10, -10, 10), spacing=7, non_dimensional_coords=(surface[0] + 10e3)
     )
     layer = tesseroid_layer(
         (longitude, latitude), surface, reference, properties={"density": density}
@@ -500,8 +500,8 @@ def test_numba_progress_missing_error(dummy_layer):
     is not installed.
     """
     (longitude, latitude), surface, reference, density = dummy_layer
-    coordinates = vd.grid_coordinates(
-        (-10, 10, -10, 10), spacing=7, extra_coords=(surface[0] + 10e3)
+    coordinates = bd.grid_coordinates(
+        (-10, 10, -10, 10), spacing=7, non_dimensional_coords=(surface[0] + 10e3)
     )
     layer = tesseroid_layer(
         (longitude, latitude), surface, reference, properties={"density": density}
@@ -516,8 +516,8 @@ def test_gravity_discarded_thin_tesseroids(dummy_layer):
     Check if gravity of tesseroid layer after discarding thin tesseroids is correct.
     """
     (longitude, latitude), surface, reference, density = dummy_layer
-    coordinates = vd.grid_coordinates(
-        (-10, 10, -10, 10), spacing=7, extra_coords=(surface[0] + 10e3)
+    coordinates = bd.grid_coordinates(
+        (-10, 10, -10, 10), spacing=7, non_dimensional_coords=(surface[0] + 10e3)
     )
     layer = tesseroid_layer(
         (longitude, latitude), surface, reference, properties={"density": density}
diff --git a/test/test_tesseroid_variable_density.py b/test/test_tesseroid_variable_density.py
index 9ec4de904..4707028b9 100644
--- a/test/test_tesseroid_variable_density.py
+++ b/test/test_tesseroid_variable_density.py
@@ -10,12 +10,11 @@
 
 from unittest.mock import patch
 
+import bordado as bd
 import numpy as np
 import numpy.testing as npt
 import pytest
-import verde as vd
 from numba import jit
-from verde import grid_coordinates
 
 import harmonica
 from harmonica import tesseroid_gravity
@@ -340,7 +339,9 @@ def constant_density(
         return density
 
     # Define a set of observation points
-    coordinates = grid_coordinates(region=(-5, 5, -5, 5), spacing=1, extra_coords=top)
+    coordinates = bd.grid_coordinates(
+        region=(-5, 5, -5, 5), spacing=1, non_dimensional_coords=top
+    )
 
     # Compare effects against the constant density implementation
     npt.assert_allclose(
@@ -464,7 +465,9 @@ def linear_density(radius):
     # Create a set of observation points located on the shell outer radius
     region = (-180, 180, -90, 90)
     shape = (19, 13)
-    coordinates = grid_coordinates(region=region, shape=shape, extra_coords=top)
+    coordinates = bd.grid_coordinates(
+        region=region, shape=shape, non_dimensional_coords=top
+    )
     # Get analytic solution
     analytic = analytical_spherical_shell_linear(
         coordinates[-1], bottom, top, slope, constant_term
@@ -508,7 +511,9 @@ def exponential_density(radius):
     # Create a set of observation points located on the shell outer radius
     region = (-180, 180, -90, 90)
     shape = (19, 13)
-    coordinates = grid_coordinates(region=region, shape=shape, extra_coords=top)
+    coordinates = bd.grid_coordinates(
+        region=region, shape=shape, non_dimensional_coords=top
+    )
     # Get analytic solution
     analytic = analytical_spherical_shell_exponential(
         coordinates[-1], bottom, top, a_factor, b_factor, constant_term
@@ -545,8 +550,8 @@ def density(r):
     @pytest.fixture
     def coordinates(self):
         """Sample coordinates."""
-        coordinates = vd.grid_coordinates(
-            region=(-15, 55, -80, 40), spacing=10, extra_coords=6.5e4
+        coordinates = bd.grid_coordinates(
+            region=(-15, 55, -80, 40), spacing=10, non_dimensional_coords=6.5e4
         )
         return coordinates
 
diff --git a/test/test_transformations.py b/test/test_transformations.py
index 84b01defb..19209f6f3 100644
--- a/test/test_transformations.py
+++ b/test/test_transformations.py
@@ -11,6 +11,7 @@
 import re
 from pathlib import Path
 
+import bordado as bd
 import numpy as np
 import numpy.testing as npt
 import pytest
@@ -63,8 +64,8 @@ def fixture_sample_grid_coords():
     """
     Define sample grid coordinates.
     """
-    grid_coords = vd.grid_coordinates(
-        region=(-300e3, 300e3, -300e3, 300e3), spacing=5000, extra_coords=0
+    grid_coords = bd.grid_coordinates(
+        region=(-300e3, 300e3, -300e3, 300e3), spacing=5000, non_dimensional_coords=0
     )
     return grid_coords
 
@@ -74,8 +75,8 @@ def fixture_upward_grid_coords():
     """
     Define upward grid coordinates.
     """
-    grid_coords = vd.grid_coordinates(
-        region=(-300e3, 300e3, -300e3, 300e3), spacing=5000, extra_coords=10e3
+    grid_coords = bd.grid_coordinates(
+        region=(-300e3, 300e3, -300e3, 300e3), spacing=5000, non_dimensional_coords=10e3
     )
     return grid_coords
 
@@ -473,8 +474,8 @@ def test_reduction_to_pole_remanent():
     """
     Test reduction_to_pole against an analytical solution with remanent magnetization.
     """
-    coordinates = vd.grid_coordinates(
-        (-70e3, 20e3, -20e3, 60e3), spacing=0.5e3, extra_coords=500
+    coordinates = bd.grid_coordinates(
+        (-70e3, 20e3, -20e3, 60e3), spacing=0.5e3, non_dimensional_coords=500
     )
     finc, fdec = -45, 13
     minc, mdec = -14, -24
@@ -510,8 +511,8 @@ def test_reduction_to_pole_induced():
     """
     Test reduction_to_pole against an analytical solution with induced magnetization.
     """
-    coordinates = vd.grid_coordinates(
-        (-70e3, 20e3, -20e3, 60e3), spacing=0.5e3, extra_coords=500
+    coordinates = bd.grid_coordinates(
+        (-70e3, 20e3, -20e3, 60e3), spacing=0.5e3, non_dimensional_coords=500
     )
     finc, fdec = -45, 13
     dipole = [-25e3, 20e3, -5000]
@@ -675,7 +676,7 @@ def test_gaussian_lowpass(self):
         Test gaussian_lowpass against a synthetic.
         """
         # Make a synthetic with only 2 known wavelengths
-        coordinates = vd.grid_coordinates((0, 10, 0, 10), spacing=0.05)
+        coordinates = bd.grid_coordinates((0, 10, 0, 10), spacing=0.05)
         wavelength_high = 0.5
         wavelength_low = 10
         component_low = (
@@ -707,7 +708,7 @@ def test_gaussian_highpass(self):
         Test gaussian_highpass against a synthetic.
         """
         # Make a synthetic with only 2 known wavelengths
-        coordinates = vd.grid_coordinates((0, 10, 0, 10), spacing=0.05)
+        coordinates = bd.grid_coordinates((0, 10, 0, 10), spacing=0.05)
         wavelength_high = 0.5
         wavelength_low = 10
         component_low = np.sin(2 * np.pi / wavelength_low * coordinates[0]) * np.cos(