diff --git a/doc/gallery_src/equivalent_sources/block_averaged_sources.py b/doc/gallery_src/equivalent_sources/block_averaged_sources.py
index 0a50e41aa..5e4d9a6d1 100644
--- a/doc/gallery_src/equivalent_sources/block_averaged_sources.py
+++ b/doc/gallery_src/equivalent_sources/block_averaged_sources.py
@@ -106,17 +106,15 @@
 
 title = "Observed magnetic anomaly data"
 
-# Make colormap of data
-# Get the 95 percentile of the maximum absolute value between the original and
-# gridded data so we can use the same color scale for both plots and have 0
-# centered at the white color.
-maxabs = vd.maxabs(data.total_field_anomaly_nt, grid.magnetic_anomaly.values) * 0.95
-pygmt.makecpt(
-    cmap="vik",
-    series=(-maxabs, maxabs),
-    background=True,
+# Get the 99th percentile of the absolute value to use as color scale limits
+maxabs = vd.maxabs(
+    data.total_field_anomaly_nt,
+    percentile=99,
 )
 
+# Make colormap of data
+pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True)
+
 with pygmt.config(FONT_TITLE="12p"):
     fig.plot(
         projection=fig_proj,
@@ -129,8 +127,6 @@
         cmap=True,
     )
 
-fig.colorbar(cmap=True, frame=["a400f100", "x+lnT"])
-
 fig.shift_origin(xshift=fig_width + 1)
 
 title = "Gridded and upward-continued"
@@ -142,6 +138,10 @@
         cmap=True,
     )
 
-fig.colorbar(cmap=True, frame=["a400f100", "x+lnT"])
+fig.colorbar(
+    cmap=True,
+    frame=["a400f100", "x+lnT"],
+    position="n0/0+jTC+w8c/0.25c+h+o-0.5c/0.8c",
+)
 
 fig.show()
diff --git a/doc/gallery_src/equivalent_sources/cartesian.py b/doc/gallery_src/equivalent_sources/cartesian.py
index ad67084df..8d71806ca 100644
--- a/doc/gallery_src/equivalent_sources/cartesian.py
+++ b/doc/gallery_src/equivalent_sources/cartesian.py
@@ -95,17 +95,15 @@
 
 title = "Observed magnetic anomaly data"
 
-# Make colormap of data
-# Get the 95 percentile of the maximum absolute value between the original and
-# gridded data so we can use the same color scale for both plots and have 0
-# centered at the white color.
-maxabs = vd.maxabs(data.total_field_anomaly_nt, grid.magnetic_anomaly.values) * 0.95
-pygmt.makecpt(
-    cmap="vik",
-    series=(-maxabs, maxabs),
-    background=True,
+# Get the 99th percentile of the absolute value to use as color scale limits
+maxabs = vd.maxabs(
+    data.total_field_anomaly_nt,
+    percentile=99,
 )
 
+# Make colormap of data
+pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True)
+
 with pygmt.config(FONT_TITLE="12p"):
     fig.plot(
         projection=fig_proj,
@@ -118,8 +116,6 @@
         cmap=True,
     )
 
-fig.colorbar(cmap=True, frame=["a400f100", "x+lnT"])
-
 fig.shift_origin(xshift=fig_width + 1)
 
 title = "Gridded and upward-continued"
@@ -131,6 +127,10 @@
         cmap=True,
     )
 
-fig.colorbar(cmap=True, frame=["a400f100", "x+lnT"])
+fig.colorbar(
+    cmap=True,
+    frame=["a400f100", "x+lnT"],
+    position="n0/0+jTC+w8c/0.25c+h+o-0.5c/0.8c",
+)
 
 fig.show()
diff --git a/doc/gallery_src/equivalent_sources/gradient_boosted.py b/doc/gallery_src/equivalent_sources/gradient_boosted.py
index 9de8143e6..5e3d58b7a 100644
--- a/doc/gallery_src/equivalent_sources/gradient_boosted.py
+++ b/doc/gallery_src/equivalent_sources/gradient_boosted.py
@@ -115,16 +115,15 @@
 
 title = "Observed gravity disturbance data"
 
-# Make colormap of data
-pygmt.makecpt(
-    cmap="vik",
-    series=(
-        -data.gravity_disturbance.quantile(0.99),
-        data.gravity_disturbance.quantile(0.99),
-    ),
-    background=True,
+# Get the 99th percentile of the absolute value to use as color scale limits
+maxabs = vd.maxabs(
+    data.gravity_disturbance,
+    percentile=99,
 )
 
+# Make colormap of data
+pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True)
+
 with pygmt.config(FONT_TITLE="14p"):
     fig.plot(
         projection=fig_proj,
@@ -137,8 +136,6 @@
         cmap=True,
     )
 
-fig.colorbar(cmap=True, frame=["a50f25", "x+lmGal"])
-
 fig.shift_origin(xshift=fig_width + 1)
 
 title = "Gridded with gradient-boosted equivalent sources"
@@ -150,6 +147,10 @@
         cmap=True,
     )
 
-fig.colorbar(cmap=True, frame=["a50f25", "x+lmGal"])
+fig.colorbar(
+    cmap=True,
+    frame=["x+lmGal"],
+    position="n0/0+jTC+w10c/0.5c+h+o-0.5c/0.9c",
+)
 
 fig.show()
diff --git a/doc/gallery_src/equivalent_sources/spherical.py b/doc/gallery_src/equivalent_sources/spherical.py
index d972de748..3c3c9cb4b 100644
--- a/doc/gallery_src/equivalent_sources/spherical.py
+++ b/doc/gallery_src/equivalent_sources/spherical.py
@@ -80,21 +80,19 @@
 # Plot observed and gridded gravity disturbance
 fig = pygmt.Figure()
 
+# Get the 99th percentile of the absolute value to use as color scale limits
+maxabs = vd.maxabs(gravity_disturbance, percentile=99)
+
 # Make colormap of data
-# Get the 90% of the maximum absolute value between the original and gridded
-# data so we can use the same color scale for both plots and have 0 centered
-# at the white color.
-maxabs = vd.maxabs(gravity_disturbance, grid.gravity_disturbance.values) * 0.90
-pygmt.makecpt(
-    cmap="vik",
-    series=(-maxabs, maxabs),
-    background=True,
-)
+pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True)
+
+
+title = "Observed gravity disturbance data"
 
 fig.plot(
-    projection="M10c",
+    projection="M12c",
     region=region,
-    frame=["WSne", "xa5", "ya4"],
+    frame=[f"WSne+t{title}", "xa5", "ya4"],
     x=longitude,
     y=latitude,
     fill=gravity_disturbance,
@@ -102,16 +100,21 @@
     cmap=True,
 )
 
-fig.colorbar(cmap=True, frame=["a100f50", "x+lmGal"])
+fig.shift_origin(xshift="w+1c")
 
-fig.shift_origin(xshift="w+3c")
+title = "Gridded and upward-continued"
 
 fig.grdimage(
-    frame=["ESnw", "xa5", "ya4"],
+    frame=[f"ESnw+t{title}", "xa5", "ya4"],
     grid=grid.gravity_disturbance,
     cmap=True,
+    nan_transparent=True,
 )
 
-fig.colorbar(cmap=True, frame=["a100f50", "x+lmGal"])
+fig.colorbar(
+    cmap=True,
+    frame=["x+lmGal"],
+    position="n0/0+jTC+w10c/0.5c+h+o-0.5c/0.9c",
+)
 
 fig.show()
diff --git a/doc/gallery_src/forward/point_gravity.py b/doc/gallery_src/forward/point_gravity.py
index 38c4a2214..c6b8ff385 100644
--- a/doc/gallery_src/forward/point_gravity.py
+++ b/doc/gallery_src/forward/point_gravity.py
@@ -54,9 +54,9 @@
 
 title = "Gravitational acceleration (downward)"
 
-maxabs = vd.maxabs(gravity) * 0.80
+maxabs = vd.maxabs(gravity, percentile=100)
 
-pygmt.makecpt(cmap="vik", series=(-maxabs, maxabs, 0.3))
+pygmt.makecpt(cmap="vik", series=[-maxabs, maxabs, 0.3], background=True)
 
 with pygmt.config(FONT_TITLE="16p"):
     fig.grdimage(
@@ -68,8 +68,16 @@
         cmap=True,
     )
 
-fig.plot(x=easting, y=northing, style="c0.2c", fill="grey")
+fig.plot(
+    x=easting,
+    y=northing,
+    style="t0.15c",
+    fill="magenta",
+    label="Point masses",
+)
+
+fig.colorbar(cmap=True, position="JMR+e", frame=["x+lmGal"])
 
-fig.colorbar(cmap=True, position="JMR", frame=["a.6f.2", "x+lmGal"])
+fig.legend()
 
 fig.show()
diff --git a/doc/gallery_src/forward/prism_layer.py b/doc/gallery_src/forward/prism_layer.py
index bec45c826..22a39ba9c 100644
--- a/doc/gallery_src/forward/prism_layer.py
+++ b/doc/gallery_src/forward/prism_layer.py
@@ -51,17 +51,17 @@
 # Plot gravity field
 fig = pygmt.Figure()
 
-title = "Gravitational acceleration of a layer of prisms"
+title = "Gravitational acceleration of topography with prisms"
 
 with pygmt.config(FONT_TITLE="14p"):
     fig.grdimage(
         region=region,
         projection="X10c/10c",
         grid=grid.gravity,
-        frame=["a", f"+t{title}", 'x+l"easting (m)"', 'y+l"northing (m)"'],
+        frame=["a", f"+t{title}", "x+leasting (m)", "y+lnorthing (m)"],
         cmap="viridis",
     )
 
-fig.colorbar(cmap=True, position="JMR", frame=["a2f1", "x+lmGal"])
+fig.colorbar(cmap=True, position="JMR", frame=["x+lmGal"])
 
 fig.show()
diff --git a/doc/gallery_src/forward/prisms_topo_gravity.py b/doc/gallery_src/forward/prisms_topo_gravity.py
index d71d7f651..de4046264 100644
--- a/doc/gallery_src/forward/prisms_topo_gravity.py
+++ b/doc/gallery_src/forward/prisms_topo_gravity.py
@@ -87,15 +87,21 @@
 
 title = "Gravitational acceleration of the topography"
 
+# Get the max absolute value to use as color scale limits
+cpt_lims = vd.maxabs(grid.gravity)
+
+# Make colormap of data
+pygmt.makecpt(cmap="balance+h0", series=[-cpt_lims, cpt_lims])
+
 with pygmt.config(FONT_TITLE="14p"):
     fig.grdimage(
         region=xy_region,
         projection=fig_proj,
         grid=grid.gravity,
         frame=["ag", f"+t{title}"],
-        cmap="vik",
+        cmap=True,
     )
 
-fig.colorbar(cmap=True, frame=["a100f50", "x+lmGal"])
+fig.colorbar(cmap=True, position="JMR", frame=["a100f50", "x+lmGal"])
 
 fig.show()
diff --git a/doc/gallery_src/forward/tesseroid.py b/doc/gallery_src/forward/tesseroid.py
index 9e805598b..0f1c8a01e 100644
--- a/doc/gallery_src/forward/tesseroid.py
+++ b/doc/gallery_src/forward/tesseroid.py
@@ -50,7 +50,7 @@
 # Plot the gravitational field
 fig = pygmt.Figure()
 
-title = "Downward component of gravitational acceleration"
+title = "Gravitational acceleration of a tesseroid"
 
 with pygmt.config(FONT_TITLE="16p"):
     fig.grdimage(
@@ -61,7 +61,17 @@
         cmap="viridis",
     )
 
-fig.colorbar(cmap=True, frame=["a200f50", "x+lmGal"])
+# Plot edges of tesseroid
+fig.plot(
+    x=[tesseroid[0], tesseroid[1], tesseroid[1], tesseroid[0], tesseroid[0]],
+    y=[tesseroid[2], tesseroid[2], tesseroid[3], tesseroid[3], tesseroid[2]],
+    pen="1p,red",
+    label="Tesseroid boundaries",
+)
+
+fig.legend()
+
+fig.colorbar(cmap=True, position="JMR", frame=["a200f50", "x+lmGal"])
 
 fig.coast(shorelines="1p,black")
 
diff --git a/doc/gallery_src/forward/tesseroid_layer.py b/doc/gallery_src/forward/tesseroid_layer.py
index 6697a7ba7..065b5e159 100644
--- a/doc/gallery_src/forward/tesseroid_layer.py
+++ b/doc/gallery_src/forward/tesseroid_layer.py
@@ -86,9 +86,11 @@
     gravity.g_z,
     projection="M15c",
     nan_transparent=True,
+    cmap=True,
     frame="+tForward gravity",
 )
+
 fig.basemap(frame=True)
-fig.colorbar(frame="af+lGravity (mGal)")
+fig.colorbar(frame="af+lgravity (mGal)")
 fig.coast(shorelines="0.5p,black", borders=["1/0.5p,black"])
 fig.show()
diff --git a/doc/gallery_src/forward/tesseroid_variable_density.py b/doc/gallery_src/forward/tesseroid_variable_density.py
index 16be3382d..8ce08634b 100644
--- a/doc/gallery_src/forward/tesseroid_variable_density.py
+++ b/doc/gallery_src/forward/tesseroid_variable_density.py
@@ -78,7 +78,7 @@ def density(radius):
 # Plot the gravitational field
 fig = pygmt.Figure()
 
-title = "Downward component of gravitational acceleration"
+title = "Gravitational acceleration of variable density tesseroids"
 
 with pygmt.config(FONT_TITLE="16p"):
     fig.grdimage(
@@ -89,7 +89,18 @@ def density(radius):
         cmap="viridis",
     )
 
-fig.colorbar(cmap=True, frame=["a200f50", "x+lmGal"])
+# Plot edges of tesseroids
+for i, tesseroid in enumerate(tesseroids):
+    label = "Tesseroid boundaries" if i == 0 else None
+    fig.plot(
+        x=[tesseroid[0], tesseroid[1], tesseroid[1], tesseroid[0], tesseroid[0]],
+        y=[tesseroid[2], tesseroid[2], tesseroid[3], tesseroid[3], tesseroid[2]],
+        pen="1p,red",
+        label=label,
+    )
+fig.legend()
+
+fig.colorbar(cmap=True, position="JMR", frame=["a200f50", "x+lmGal"])
 
 fig.coast(shorelines="1p,black")
 
diff --git a/doc/gallery_src/gravity_disturbance.py b/doc/gallery_src/gravity_disturbance.py
index 2f038793c..393deeecb 100644
--- a/doc/gallery_src/gravity_disturbance.py
+++ b/doc/gallery_src/gravity_disturbance.py
@@ -19,6 +19,7 @@
 import boule as bl
 import ensaio
 import pygmt
+import verde as vd
 import xarray as xr
 
 # Load the global gravity grid
@@ -36,7 +37,11 @@
 # Make a plot of data using PyGMT
 fig = pygmt.Figure()
 
-pygmt.grd2cpt(grid=disturbance, cmap="polar", continuous=True)
+# Get the 99.9th percentile of the absolute value to use as color scale limits
+maxabs = vd.maxabs(disturbance, percentile=99.9)
+
+# Make colormap of data
+pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True)
 
 title = "Gravity disturbance of the Earth"
 
@@ -50,6 +55,6 @@
 
 fig.coast(shorelines="0.5p,black", resolution="crude")
 
-fig.colorbar(cmap=True, frame=["a100f50", "x+lmGal"])
+fig.colorbar(cmap=True, frame=["x+lmGal"], position="+e")
 
 fig.show()
diff --git a/doc/gallery_src/gravity_disturbance_topofree.py b/doc/gallery_src/gravity_disturbance_topofree.py
index 092a5ff58..fe989aa44 100644
--- a/doc/gallery_src/gravity_disturbance_topofree.py
+++ b/doc/gallery_src/gravity_disturbance_topofree.py
@@ -29,6 +29,7 @@
 import boule as bl
 import ensaio
 import pygmt
+import verde as vd
 import xarray as xr
 
 import harmonica as hm
@@ -62,7 +63,11 @@
 # Make a plot of data using PyGMT
 fig = pygmt.Figure()
 
-pygmt.grd2cpt(grid=disturbance_topofree, cmap="vik+h0", continuous=True)
+# Get the 99th percentile of the absolute value to use as color scale limits
+maxabs = vd.maxabs(disturbance_topofree, percentile=99)
+
+# Make colormap of data
+pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True)
 
 title = "Topography-free (Bouguer) gravity disturbance of the Earth"
 
@@ -77,6 +82,6 @@
 
 fig.coast(shorelines="0.5p,black", resolution="crude")
 
-fig.colorbar(cmap=True, frame=["a200f50", "x+lmGal"])
+fig.colorbar(cmap=True, frame=["x+lmGal"], position="+e")
 
 fig.show()
diff --git a/doc/gallery_src/isostatic_moho_airy.py b/doc/gallery_src/isostatic_moho_airy.py
index c5aab75e2..c6fb1b29f 100644
--- a/doc/gallery_src/isostatic_moho_airy.py
+++ b/doc/gallery_src/isostatic_moho_airy.py
@@ -28,6 +28,7 @@
 import ensaio
 import numpy as np
 import pygmt
+import verde as vd
 import xarray as xr
 
 import harmonica as hm
@@ -56,10 +57,14 @@
 print("\nMoho depth grid:")
 print(moho)
 
-# Draw the maps
+# Make a plot of data using PyGMT
 fig = pygmt.Figure()
 
-pygmt.grd2cpt(grid=moho, cmap="viridis", reverse=True, continuous=True)
+# Get the min and max values to use as color scale limits
+cpt_lims = vd.minmax(moho)
+
+# Make colormap of data
+pygmt.makecpt(cmap="viridis", series=cpt_lims)
 
 title = "Airy isostatic Moho depth of Africa"
 
@@ -73,6 +78,6 @@
 
 fig.coast(shorelines="0.5p,black", resolution="crude")
 
-fig.colorbar(cmap=True, frame=["a10000f2500", "x+lmeters"])
+fig.colorbar(cmap=True, frame=["x+lmeters"])
 
 fig.show()
diff --git a/doc/gallery_src/transformations/reduction_to_pole.py b/doc/gallery_src/transformations/reduction_to_pole.py
index 64d7aa3a7..2088a8cfa 100644
--- a/doc/gallery_src/transformations/reduction_to_pole.py
+++ b/doc/gallery_src/transformations/reduction_to_pole.py
@@ -41,27 +41,27 @@
 # Plot original magnetic anomaly and the reduced to the pole
 fig = pygmt.Figure()
 with fig.subplot(nrows=1, ncols=2, figsize=("28c", "15c"), sharey="l"):
-    # Make colormap for both plots (saturate it a little bit)
-    scale = 0.5 * vd.maxabs(magnetic_grid, rtp_grid)
-    pygmt.makecpt(cmap="polar+h", series=[-scale, scale], background=True)
+    # Make colormap for both plots
+    cpt_lim = vd.maxabs(magnetic_grid, rtp_grid, percentile=99.9)
+    pygmt.makecpt(cmap="balance+h0", series=[-cpt_lim, cpt_lim], background=True)
     with fig.set_panel(panel=0):
         # Plot magnetic anomaly grid
         fig.grdimage(
-            grid=magnetic_grid,
-            projection="X?",
-            cmap=True,
-        )
-        # Add colorbar
-        fig.colorbar(
-            frame='af+l"Magnetic anomaly [nT]"',
-            position="JBC+h+o0/1c+e",
+            grid=magnetic_grid, projection="X?", cmap=True, frame="+tMagnetic anomaly"
         )
     with fig.set_panel(panel=1):
         # Plot upward reduced to the pole grid
-        fig.grdimage(grid=rtp_grid, projection="X?", cmap=True)
-        # Add colorbar
-        fig.colorbar(
-            frame='af+l"Reduced to the pole [nT]"',
-            position="JBC+h+o0/1c+e",
+        fig.grdimage(
+            grid=rtp_grid,
+            projection="X?",
+            cmap=True,
+            frame="+tReduced to the pole",
         )
+    # Add colorbar
+    fig.colorbar(
+        cmap=True,
+        frame=["a1000f500", "x+lnT"],
+        position="n0/0+jTC+w12c/0.5c+h+o-0.5c/0.9c+e",
+    )
+
 fig.show()
diff --git a/doc/gallery_src/transformations/tga.py b/doc/gallery_src/transformations/tga.py
index 24d3b1bbe..3619b2cae 100644
--- a/doc/gallery_src/transformations/tga.py
+++ b/doc/gallery_src/transformations/tga.py
@@ -31,28 +31,34 @@
 with fig.subplot(nrows=1, ncols=2, figsize=("28c", "15c"), sharey="l"):
     with fig.set_panel(panel=0):
         # Make colormap of data
-        scale = 2500
-        pygmt.makecpt(cmap="polar+h", series=[-scale, scale], background=True)
+        maxabs = vd.maxabs(magnetic_grid, percentile=99.9)
+        pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True)
         # Plot magnetic anomaly grid
         fig.grdimage(
             grid=magnetic_grid,
             projection="X?",
             cmap=True,
+            frame="+tMagnetic Anomaly",
         )
         # Add colorbar
         fig.colorbar(
-            frame='af+l"Magnetic anomaly [nT]"',
+            frame="af+lnT",
             position="JBC+h+o0/1c+e",
         )
     with fig.set_panel(panel=1):
         # Make colormap for total gradient amplitude (saturate it a little bit)
-        scale = 0.6 * vd.maxabs(tga)
-        pygmt.makecpt(cmap="polar+h", series=[0, scale], background=True)
+        maxabs = vd.maxabs(tga, percentile=99.9)
+        pygmt.makecpt(cmap="balance+h0", series=[0, maxabs], background=True)
         # Plot total gradient amplitude
-        fig.grdimage(grid=tga, projection="X?", cmap=True)
+        fig.grdimage(
+            grid=tga,
+            projection="X?",
+            cmap=True,
+            frame="+tTotal Gradient Amplitude",
+        )
         # Add colorbar
         fig.colorbar(
-            frame='af+l"Total Gradient Amplitude [nT/m]"',
-            position="JBC+h+o0/1c+e",
+            frame="af+lnT/m",
+            position="JBC+h+o0/1c+ef",
         )
 fig.show()
diff --git a/doc/gallery_src/transformations/tilt.py b/doc/gallery_src/transformations/tilt.py
index e7efdbb76..b3081a545 100644
--- a/doc/gallery_src/transformations/tilt.py
+++ b/doc/gallery_src/transformations/tilt.py
@@ -62,10 +62,10 @@
     sharey="l",
     margins=["1c", "1c"],
 ):
-    scale = 0.5 * vd.maxabs(magnetic_grid, rtp_grid)
+    maxabs = vd.maxabs(magnetic_grid, rtp_grid, percentile=99)
     with fig.set_panel(panel=0):
         # Make colormap of data
-        pygmt.makecpt(cmap="polar+h", series=[-scale, scale], background=True)
+        pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True)
         # Plot magnetic anomaly grid
         fig.grdimage(
             grid=magnetic_grid,
@@ -75,7 +75,7 @@
         )
     with fig.set_panel(panel=1):
         # Make colormap of data
-        pygmt.makecpt(cmap="polar+h", series=[-scale, scale], background=True)
+        pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True)
         # Plot reduced to the pole magnetic anomaly grid
         fig.grdimage(
             grid=rtp_grid,
@@ -84,16 +84,15 @@
             frame=["a", "+tReduced to the pole (RTP)"],
         )
         # Add colorbar
-        label = "nT"
         fig.colorbar(
-            frame=f"af+l{label}",
+            frame="af+lnT",
             position="JMR+o1/-0.25c+e",
         )
 
-    scale = 0.6 * vd.maxabs(tilt_grid, tilt_rtp_grid)
+    maxabs = vd.maxabs(tilt_grid, tilt_rtp_grid, percentile=99)
     with fig.set_panel(panel=2):
         # Make colormap for tilt (saturate it a little bit)
-        pygmt.makecpt(cmap="polar+h", series=[-scale, scale], background=True)
+        pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True)
         # Plot tilt
         fig.grdimage(
             grid=tilt_grid,
@@ -103,7 +102,7 @@
         )
     with fig.set_panel(panel=3):
         # Make colormap for tilt rtp (saturate it a little bit)
-        pygmt.makecpt(cmap="polar+h", series=[-scale, scale], background=True)
+        pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True)
         # Plot tilt
         fig.grdimage(
             grid=tilt_rtp_grid,
@@ -112,9 +111,8 @@
             frame=["a", "+tTilt of RTP grid"],
         )
         # Add colorbar
-        label = "rad"
         fig.colorbar(
-            frame=f"af+l{label}",
+            frame="af+lradians",
             position="JMR+o1/-0.25c+e",
         )
 fig.show()
diff --git a/doc/gallery_src/transformations/upward_continuation.py b/doc/gallery_src/transformations/upward_continuation.py
index 182b29261..4058d08bf 100644
--- a/doc/gallery_src/transformations/upward_continuation.py
+++ b/doc/gallery_src/transformations/upward_continuation.py
@@ -11,6 +11,7 @@
 
 import ensaio
 import pygmt
+import verde as vd
 import xarray as xr
 
 import harmonica as hm
@@ -30,27 +31,30 @@
 # Plot original magnetic anomaly and the upward continued grid
 fig = pygmt.Figure()
 with fig.subplot(nrows=1, ncols=2, figsize=("28c", "15c"), sharey="l"):
-    # Make colormap for both plots data
-    scale = 2500
-    pygmt.makecpt(cmap="polar+h", series=[-scale, scale], background=True)
+    # Make colormap based on original data
+    cpt_lim = vd.maxabs(magnetic_grid, percentile=99.9)
+    pygmt.makecpt(cmap="balance+h0", series=[-cpt_lim, cpt_lim], background=True)
     with fig.set_panel(panel=0):
         # Plot magnetic anomaly grid
         fig.grdimage(
             grid=magnetic_grid,
             projection="X?",
             cmap=True,
-        )
-        # Add colorbar
-        fig.colorbar(
-            frame='af+l"Magnetic anomaly at 500m [nT]"',
-            position="JBC+h+o0/1c+e",
+            frame="+tMagnetic Anomaly at 500m",
         )
     with fig.set_panel(panel=1):
         # Plot upward continued grid
-        fig.grdimage(grid=upward_continued, projection="X?", cmap=True)
-        # Add colorbar
-        fig.colorbar(
-            frame='af+l"Upward continued to 1000m [nT]"',
-            position="JBC+h+o0/1c+e",
+        fig.grdimage(
+            grid=upward_continued,
+            projection="X?",
+            frame="+tUpward continued to 1000m",
+            cmap=True,
         )
+    # Add colorbar
+    fig.colorbar(
+        cmap=True,
+        frame=["a1000f500", "x+lnT"],
+        position="n0/0+jTC+w12c/0.5c+h+o-0.5c/0.9c+e",
+    )
+
 fig.show()
diff --git a/doc/gallery_src/transformations/upward_derivative.py b/doc/gallery_src/transformations/upward_derivative.py
index b432a0b57..38e90f21b 100644
--- a/doc/gallery_src/transformations/upward_derivative.py
+++ b/doc/gallery_src/transformations/upward_derivative.py
@@ -32,28 +32,34 @@
 with fig.subplot(nrows=1, ncols=2, figsize=("28c", "15c"), sharey="l"):
     with fig.set_panel(panel=0):
         # Make colormap of data
-        scale = 2500
-        pygmt.makecpt(cmap="polar+h", series=[-scale, scale], background=True)
+        maxabs = vd.maxabs(magnetic_grid, percentile=99)
+        pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True)
         # Plot magnetic anomaly grid
         fig.grdimage(
             grid=magnetic_grid,
             projection="X?",
             cmap=True,
+            frame="+tMagnetic Anomaly",
         )
         # Add colorbar
         fig.colorbar(
-            frame='af+l"Magnetic anomaly [nT]"',
+            frame="af+lnT",
             position="JBC+h+o0/1c+e",
         )
     with fig.set_panel(panel=1):
-        # Make colormap for upward derivative (saturate it a little bit)
-        scale = 0.6 * vd.maxabs(deriv_upward)
-        pygmt.makecpt(cmap="polar+h", series=[-scale, scale], background=True)
+        # Make colormap for upward derivative
+        maxabs = vd.maxabs(deriv_upward, percentile=99)
+        pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True)
         # Plot upward derivative
-        fig.grdimage(grid=deriv_upward, projection="X?", cmap=True)
+        fig.grdimage(
+            grid=deriv_upward,
+            projection="X?",
+            cmap=True,
+            frame="+tUpward Derivative",
+        )
         # Add colorbar
         fig.colorbar(
-            frame='af+l"Upward derivative [nT/m]"',
+            frame="af+lnT/m",
             position="JBC+h+o0/1c+e",
         )
 fig.show()
diff --git a/doc/install.rst b/doc/install.rst
index c91c9ff8f..ed379989b 100644
--- a/doc/install.rst
+++ b/doc/install.rst
@@ -68,13 +68,14 @@ Required:
 
 * `numpy <http://www.numpy.org/>`__
 * `pandas <http://pandas.pydata.org/>`__
-* `numba <https://numba.pydata.org/>`__
 * `scipy <https://www.scipy.org/>`__
-* `xarray <https://xarray.pydata.org/>`__
 * `scikit-learn <https://scikit-learn.org>`__
-* `pooch <http://www.fatiando.org/pooch/>`__
+* `numba <https://numba.pydata.org/>`__
+* `xarray <https://xarray.pydata.org/>`__
 * `verde <http://www.fatiando.org/verde/>`__
 * `xrft <https://xrft.readthedocs.io/>`__
+* `choclo <http://www.fatiando.org/choclo/>`__
+* `boule <http://www.fatiando.org/boule/>`__
 
 Optional:
 
@@ -92,7 +93,7 @@ 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
+* `pygmt <https://www.pygmt.org/>`__ and `matplotlib <https://matplotlib.org/>`__ for plotting maps and figures
 * `pyproj <https://jswhit.github.io/pyproj/>`__ for cartographic projections
+* `bordado <http://www.fatiando.org/bordado/>`__ for working with geographic coordinates
\ No newline at end of file
diff --git a/doc/user_guide/equivalent_sources/block-averaged-eqs.rst b/doc/user_guide/equivalent_sources/block-averaged-eqs.rst
index c27e08dbb..acf9c21fd 100644
--- a/doc/user_guide/equivalent_sources/block-averaged-eqs.rst
+++ b/doc/user_guide/equivalent_sources/block-averaged-eqs.rst
@@ -79,7 +79,7 @@ And project the geographic coordinates to plain Cartesian ones:
 
     import pygmt
 
-    maxabs = vd.maxabs(data.total_field_anomaly_nt)*.8
+    maxabs = vd.maxabs(data.total_field_anomaly_nt, percentile=99)
 
     # Set figure properties
     w, e, s, n = xy_region
@@ -94,7 +94,7 @@ And project the geographic coordinates to plain Cartesian ones:
     title = "Observed total-field magnetic anomaly"
 
     pygmt.makecpt(
-        cmap="polar+h0",
+        cmap="balance+h0",
         series=(-maxabs, maxabs),
         background=True,
     )
@@ -110,14 +110,14 @@ And project the geographic coordinates to plain Cartesian ones:
             style="c0.1c",
             cmap=True,
         )
-    fig.colorbar(cmap=True, position="JMR", frame=["a200f100", "x+lnT"])
+    fig.colorbar(cmap=True, position="JMR+e", frame=["a200f100", "y+lnT"])
     fig.show()
 
 
-Most airborne surveys like this one present an anysotropic distribution of the
+Most airborne surveys like this one present an anisotropic distribution of the
 data: there are more observation points along the flight lines that goes west
 to east than the ones going south to north.
-Placing a single source beneath each observation point generates an anysotropic
+Placing a single source beneath each observation point generates an anisotropic
 distribution of the equivalent sources, which might lead to aliases on the
 generated outputs.
 
@@ -182,9 +182,10 @@ we are efectivelly upward continuing the data.
 
     title = "Observed magnetic anomaly data"
     pygmt.makecpt(
-        cmap="polar+h0",
+        cmap="balance+h0",
         series=(-maxabs, maxabs),
-        background=True)
+        background=True,
+    )
 
     with pygmt.config(FONT_TITLE="14p"):
         fig.plot(
@@ -197,7 +198,6 @@ we are efectivelly upward continuing the data.
             style="c0.1c",
             cmap=True,
         )
-    fig.colorbar(cmap=True, frame=["a200f100", "x+lnT"])
 
     fig.shift_origin(xshift=fig_width + 1)
 
@@ -209,7 +209,11 @@ we are efectivelly upward continuing the data.
             grid=grid.magnetic_anomaly,
             cmap=True,
         )
-    fig.colorbar(cmap=True, frame=["a200f100", "x+lnT"])
+    fig.colorbar(
+        cmap=True,
+        frame=["a200f50", "x+lnT"],
+        position=f"n0/0+jTC+w{fig_width*.75}c/0.5c+h+o-0.5c/1c+e",
+    )
 
     fig.show()
 
diff --git a/doc/user_guide/equivalent_sources/eq-sources-spherical.rst b/doc/user_guide/equivalent_sources/eq-sources-spherical.rst
index a27d7f529..fe8461691 100644
--- a/doc/user_guide/equivalent_sources/eq-sources-spherical.rst
+++ b/doc/user_guide/equivalent_sources/eq-sources-spherical.rst
@@ -162,12 +162,12 @@ Lets plot it:
 
     import pygmt
 
-    maxabs = vd.maxabs(gravity_disturbance, grid.gravity_disturbance.values)
+    maxabs = vd.maxabs(gravity_disturbance, grid.gravity_disturbance.values, percentile=99)
 
     fig = pygmt.Figure()
 
     # Make colormap of data
-    pygmt.makecpt(cmap="polar+h0",series=(-maxabs, maxabs,))
+    pygmt.makecpt(cmap="balance+h0", series=(-maxabs, maxabs), background=True)
 
     title = "Block-median reduced gravity disturbance"
     fig.plot(
@@ -181,9 +181,8 @@ Lets plot it:
         cmap=True,
     )
     fig.coast(shorelines="0.5p,black", area_thresh=1e4)
-    fig.colorbar(cmap=True, frame=["a50f25", "x+lmGal"])
 
-    fig.shift_origin(xshift='w+3c')
+    fig.shift_origin(xshift='w+1c')
 
     title = "Gridded gravity disturbance"
     fig.grdimage(
@@ -193,7 +192,12 @@ Lets plot it:
         nan_transparent=True,
     )
     fig.coast(shorelines="0.5p,black", area_thresh=1e4)
-    fig.colorbar(cmap=True, frame=["a50f25", "x+lmGal"])
+
+    fig.colorbar(
+      cmap=True,
+      frame=["a50f25", "x+lmGal"],
+      position=f"n0/0+jTC+w10c/.5c+h+o-0.5c/1c+e",
+    )
 
     fig.show()
 
diff --git a/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst b/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst
index 22b5c2c57..bc361a302 100644
--- a/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst
+++ b/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst
@@ -80,7 +80,7 @@ value of the score obtained after each cross validation.
            data.gravity_disturbance_mgal,
        )
    )
-   score_first_guess
+   print(score_first_guess)
 
 The resulting score corresponds to the R^2. It represents how well the
 equivalent sources can reproduce the variation of our data. As closer it gets
@@ -145,7 +145,7 @@ And now we can actually ran one cross validation for each pair of parameters:
            )
        )
        scores.append(score)
-   scores
+   print(scores)
 
 Once every score has been computed, we can obtain the best score and the
 corresponding parameters that generate it:
@@ -216,7 +216,7 @@ Lets plot it:
    fig = pygmt.Figure()
 
    # Make colormap of data
-   pygmt.makecpt(cmap="polar+h0",series=(-maxabs, maxabs,))
+   pygmt.makecpt(cmap="balance+h0",series=(-maxabs, maxabs), background=True)
 
    title = "Gravity disturbance with first guess"
 
@@ -227,7 +227,6 @@ Lets plot it:
       grid=grid_first_guess.scalars,
       cmap=True,
    )
-   fig.colorbar(cmap=True, frame=["a50f25", "x+lmGal"])
 
    fig.shift_origin(xshift=fig_width + 1)
 
@@ -238,7 +237,11 @@ Lets plot it:
       grid=grid.scalars,
       cmap=True,
    )
-   fig.colorbar(cmap=True, frame=["a50f25", "x+lmGal"])
+   fig.colorbar(
+      cmap=True,
+      frame=["a50f25", "x+lmGal"],
+      position=f"n0/0+jTC+w{fig_width*.75}c/0.5c+h+o-0.5c/1c+e",
+   )
 
    fig.show()
 
diff --git a/doc/user_guide/equivalent_sources/gradient-boosted-eqs.rst b/doc/user_guide/equivalent_sources/gradient-boosted-eqs.rst
index e44a3d218..4affdcef0 100644
--- a/doc/user_guide/equivalent_sources/gradient-boosted-eqs.rst
+++ b/doc/user_guide/equivalent_sources/gradient-boosted-eqs.rst
@@ -168,12 +168,16 @@ And plot it:
     fig_ratio = (n - s) / (fig_height / 100)
     fig_proj = f"x1:{fig_ratio}"
 
-    maxabs = vd.maxabs(disturbance, grid_masked.gravity_disturbance)
+    maxabs = vd.maxabs(disturbance, grid_masked.gravity_disturbance, percentile=99)
 
     fig = pygmt.Figure()
 
     # Make colormap of data
-    pygmt.makecpt(cmap="polar+h0",series=(-maxabs, maxabs,))
+    pygmt.makecpt(
+        cmap="balance+h0",
+        series=(-maxabs, maxabs),
+        background=True,
+    )
 
     title = "Observed gravity disturbance data"
     with pygmt.config(FONT_TITLE="14p"):
@@ -187,7 +191,6 @@ And plot it:
             style="c0.1c",
             cmap=True,
         )
-    fig.colorbar(cmap=True, frame=["a50f25", "x+lmGal"])
 
     fig.shift_origin(xshift=fig_width + 1)
 
@@ -199,7 +202,11 @@ And plot it:
             cmap=True,
         )
 
-    fig.colorbar(cmap=True, frame=["a50f25", "x+lmGal"])
+    fig.colorbar(
+      cmap=True,
+      frame=["a50f25", "x+lmGal"],
+      position=f"n0/0+jTC+w{fig_width*.75}c/0.5c+h+o-0.5c/1c+e",
+    )
 
     fig.show()
 
diff --git a/doc/user_guide/equivalent_sources/index.rst b/doc/user_guide/equivalent_sources/index.rst
index 266a9e3e7..457815d0c 100644
--- a/doc/user_guide/equivalent_sources/index.rst
+++ b/doc/user_guide/equivalent_sources/index.rst
@@ -133,7 +133,7 @@ And plot it:
    import pygmt
 
    # Get max absolute value for the observed gravity disturbance
-   maxabs = vd.maxabs(data.gravity_disturbance_mgal)
+   maxabs = vd.maxabs(disturbance, data.gravity_disturbance_mgal)
 
    # Set figure properties
    w, e, s, n = region
@@ -143,7 +143,7 @@ And plot it:
    fig_proj = f"x1:{fig_ratio}"
 
    fig = pygmt.Figure()
-   pygmt.makecpt(cmap="polar+h0", series=[-maxabs, maxabs])
+   pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True)
    title="Predicted gravity disturbance"
    with pygmt.config(FONT_TITLE="14p"):
       fig.plot(
@@ -156,7 +156,6 @@ And plot it:
          region=region,
          frame=['ag', f"+t{title}"],
       )
-   fig.colorbar(cmap=True, position="JMR", frame=["a50f25", "y+lmGal"])
 
    fig.shift_origin(yshift=fig_height + 2)
 
@@ -170,7 +169,7 @@ And plot it:
          style="c3p",
          frame=['ag', f"+t{title}"],
       )
-   fig.colorbar(cmap=True, position="JMR", frame=["a50f25", "y+lmGal"])
+   fig.colorbar(cmap=True, position=f"JMR+o1c/{-(fig_height+2)/2}+w12c/.5c+e", frame=["a50f25", "y+lmGal"])
 
    fig.show()
 
@@ -204,7 +203,7 @@ And plot it
    maxabs = vd.maxabs(grid.gravity_disturbance)
 
    fig = pygmt.Figure()
-   pygmt.makecpt(cmap="polar+h0", series=[-maxabs, maxabs])
+   pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True)
    fig.grdimage(
       frame=['af', 'WSen'],
       grid=grid.gravity_disturbance,
@@ -212,7 +211,7 @@ And plot it
       projection=fig_proj,
       cmap=True,
    )
-   fig.colorbar(cmap=True, frame=["a50f25", "x+lgravity disturbance", "y+lmGal"])
+   fig.colorbar(cmap=True, frame=["a50f25", "x+lGravity disturbance", "y+lmGal"])
 
    fig.show()
 
diff --git a/doc/user_guide/forward_modelling/dipole.rst b/doc/user_guide/forward_modelling/dipole.rst
index e7a5e9a6d..5555d8539 100644
--- a/doc/user_guide/forward_modelling/dipole.rst
+++ b/doc/user_guide/forward_modelling/dipole.rst
@@ -92,19 +92,35 @@ every observation point:
 
    b_e, b_n, b_u = hm.dipole_magnetic(coordinates, dipoles, magnetic_moments, field="b")
 
+We can reshape the results into variables of an dataset:
+
+.. jupyter-execute::
+
+   grid = vd.make_xarray_grid(
+      coordinates,
+      data=(b_e, b_n, b_u),
+      data_names=["b_e", "b_n", "b_u"],
+      extra_coords_names="extra"
+   )
+
 .. jupyter-execute::
 
    import matplotlib.pyplot as plt
 
    fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8))
 
-   fields = {"b_e": b_e, "b_n": b_n, "b_u": b_u}
-   for field, ax in zip(fields, axes):
-      tmp = ax.pcolormesh(coordinates[0], coordinates[1], fields[field])
+   maxabs = vd.maxabs(*[grid[v] for v in grid.variables], percentile=99)
+   for field, ax in zip(grid.variables, axes):
+      tmp = grid[field].plot(
+         ax=ax,
+         add_colorbar=False,
+         vmin=-maxabs,
+         vmax=maxabs,
+         cmap='RdBu_r',
+      )
       ax.set_aspect("equal")
       ax.set_title(field)
-      ax.ticklabel_format(style="sci", scilimits=(0, 0), axis="both")
-      plt.colorbar(tmp, ax=ax, orientation="horizontal", label="nT", pad=0.008)
+   fig.colorbar(tmp, ax=axes.ravel().tolist(), orientation="horizontal", label="nT", fraction=.03, pad=0.08)
    plt.show()
 
 ----
diff --git a/doc/user_guide/forward_modelling/ellipsoid.rst b/doc/user_guide/forward_modelling/ellipsoid.rst
index 9b3e3f3f2..a25e6290d 100644
--- a/doc/user_guide/forward_modelling/ellipsoid.rst
+++ b/doc/user_guide/forward_modelling/ellipsoid.rst
@@ -95,30 +95,32 @@ We can then forward model the ellipsoid:
    import verde as vd
    import matplotlib.pyplot as plt
 
+   maxabs = vd.maxabs(ge, gn, gz, percentile=99)
+
    _, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(10, 6))
    for ax, g_component, label in zip(
       axes, (ge, gn, gz), ("ge", "gn", "gz"), strict=True
    ):
-       maxabs = vd.maxabs(g_component)
-       tmp = ax.pcolormesh(
-           coordinates[0],
-           coordinates[1],
-           g_component,
-           vmin=-maxabs,
-           vmax=maxabs,
-           cmap="RdBu_r",
-       )
-       plt.colorbar(
-           tmp,
-           ax=ax,
-           label=f"{label} [mGal]",
-           orientation="horizontal",
-           pad=0.12,
-           format='%.0e',
-       )
-       ax.set_xlabel("Easting [m]")
-       ax.set_ylabel("Northing [m]")
-       ax.set_aspect("equal")
+      tmp = ax.pcolormesh(
+         coordinates[0],
+         coordinates[1],
+         g_component,
+         vmin=-maxabs,
+         vmax=maxabs,
+         cmap="RdBu_r",
+      )
+      ax.set_title(label)
+      ax.set_xlabel("Easting [m]")
+      ax.set_ylabel("Northing [m]")
+      ax.set_aspect("equal")
+   plt.colorbar(
+      tmp,
+      ax=axes,
+      label="mGal",
+      orientation="horizontal",
+      pad=0.12,
+      fraction=.03,
+   )
    plt.show()
 
 
@@ -138,30 +140,32 @@ We can also forward model multiple ellipsoid by creating a list of them:
 
 .. jupyter-execute::
 
+   maxabs = vd.maxabs(ge, gn, gz, percentile=95)
+
    _, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(10, 6))
    for ax, g_component, label in zip(
       axes, (ge, gn, gz), ("ge", "gn", "gz"), strict=True
    ):
-       maxabs = vd.maxabs(g_component)
-       tmp = ax.pcolormesh(
-           coordinates[0],
-           coordinates[1],
-           g_component,
-           vmin=-maxabs,
-           vmax=maxabs,
-           cmap="RdBu_r",
-       )
-       plt.colorbar(
-           tmp,
-           ax=ax,
-           label=f"{label} [mGal]",
-           orientation="horizontal",
-           pad=0.12,
-           format='%.0e',
-       )
-       ax.set_xlabel("Easting [m]")
-       ax.set_ylabel("Northing [m]")
-       ax.set_aspect("equal")
+      tmp = ax.pcolormesh(
+         coordinates[0],
+         coordinates[1],
+         g_component,
+         vmin=-maxabs,
+         vmax=maxabs,
+         cmap="RdBu_r",
+      )
+      ax.set_title(label)
+      ax.set_xlabel("Easting [m]")
+      ax.set_ylabel("Northing [m]")
+      ax.set_aspect("equal")
+   plt.colorbar(
+      tmp,
+      ax=axes,
+      label="mGal",
+      orientation="horizontal",
+      pad=0.12,
+      fraction=.03,
+   )
    plt.show()
 
 
@@ -213,22 +217,32 @@ And forward the magnetic field in the grid of observation points:
 
 .. jupyter-execute::
 
+   maxabs = vd.maxabs(be, bn, bu, percentile=99)
+
    _, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(10, 6))
    for ax, b_component, label in zip(
       axes, (be, bn, bu), ("Be", "Bn", "Bu"), strict=True
    ):
-       tmp = ax.pcolormesh(coordinates[0], coordinates[1], b_component)
-       plt.colorbar(
-           tmp,
-           ax=ax,
-           label=f"{label} [nT]",
-           orientation="horizontal",
-           pad=0.12,
-           format='%.0e',
-       )
-       ax.set_xlabel("Easting [m]")
-       ax.set_ylabel("Northing [m]")
-       ax.set_aspect("equal")
+      tmp = ax.pcolormesh(
+         coordinates[0],
+         coordinates[1],
+         b_component,
+         vmin=-maxabs,
+         vmax=maxabs,
+         cmap="RdBu_r",
+      )
+      ax.set_title(label)
+      ax.set_xlabel("Easting [m]")
+      ax.set_ylabel("Northing [m]")
+      ax.set_aspect("equal")
+   plt.colorbar(
+      tmp,
+      ax=axes,
+      label="nT",
+      orientation="horizontal",
+      pad=0.12,
+      fraction=.03,
+   )
    plt.show()
 
 
@@ -257,22 +271,32 @@ through the ``remanent_mag`` physical property:
 
 .. jupyter-execute::
 
+   maxabs = vd.maxabs(be, bn, bu, percentile=99)
+
    _, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(10, 6))
    for ax, b_component, label in zip(
       axes, (be, bn, bu), ("Be", "Bn", "Bu"), strict=True
    ):
-       tmp = ax.pcolormesh(coordinates[0], coordinates[1], b_component)
-       plt.colorbar(
-           tmp,
-           ax=ax,
-           label=f"{label} [nT]",
-           orientation="horizontal",
-           pad=0.12,
-           format='%.0e',
-       )
-       ax.set_xlabel("Easting [m]")
-       ax.set_ylabel("Northing [m]")
-       ax.set_aspect("equal")
+      tmp = ax.pcolormesh(
+         coordinates[0],
+         coordinates[1],
+         b_component,
+         vmin=-maxabs,
+         vmax=maxabs,
+         cmap="RdBu_r",
+      )
+      ax.set_title(label)
+      ax.set_xlabel("Easting [m]")
+      ax.set_ylabel("Northing [m]")
+      ax.set_aspect("equal")
+   plt.colorbar(
+      tmp,
+      ax=axes,
+      label="nT",
+      orientation="horizontal",
+      pad=0.12,
+      fraction=.03,
+   )
    plt.show()
 
 
@@ -302,22 +326,32 @@ ellipsoid:
 
 .. jupyter-execute::
 
+   maxabs = vd.maxabs(be, bn, bu, percentile=99)
+
    _, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(10, 6))
    for ax, b_component, label in zip(
       axes, (be, bn, bu), ("Be", "Bn", "Bu"), strict=True
    ):
-       tmp = ax.pcolormesh(coordinates[0], coordinates[1], b_component)
-       plt.colorbar(
-           tmp,
-           ax=ax,
-           label=f"{label} [nT]",
-           orientation="horizontal",
-           pad=0.12,
-           format='%.0e',
-       )
-       ax.set_xlabel("Easting [m]")
-       ax.set_ylabel("Northing [m]")
-       ax.set_aspect("equal")
+      tmp = ax.pcolormesh(
+         coordinates[0],
+         coordinates[1],
+         b_component,
+         vmin=-maxabs,
+         vmax=maxabs,
+         cmap="RdBu_r",
+      )
+      ax.set_title(label)
+      ax.set_xlabel("Easting [m]")
+      ax.set_ylabel("Northing [m]")
+      ax.set_aspect("equal")
+   plt.colorbar(
+      tmp,
+      ax=axes,
+      label="nT",
+      orientation="horizontal",
+      pad=0.12,
+      fraction=.03,
+   )
    plt.show()
 
 
@@ -356,20 +390,30 @@ We can also forward a collection of ellipsoids with mixed physical properties:
 
 .. jupyter-execute::
 
+   maxabs = vd.maxabs(be, bn, bu, percentile=99)
+
    _, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(10, 6))
    for ax, b_component, label in zip(
       axes, (be, bn, bu), ("Be", "Bn", "Bu"), strict=True
    ):
-       tmp = ax.pcolormesh(coordinates[0], coordinates[1], b_component)
-       plt.colorbar(
-           tmp,
-           ax=ax,
-           label=f"{label} [nT]",
-           orientation="horizontal",
-           pad=0.12,
-           format='%.0e',
-       )
-       ax.set_xlabel("Easting [m]")
-       ax.set_ylabel("Northing [m]")
-       ax.set_aspect("equal")
+      tmp = ax.pcolormesh(
+         coordinates[0],
+         coordinates[1],
+         b_component,
+         vmin=-maxabs,
+         vmax=maxabs,
+         cmap="RdBu_r",
+      )
+      ax.set_title(label)
+      ax.set_xlabel("Easting [m]")
+      ax.set_ylabel("Northing [m]")
+      ax.set_aspect("equal")
+   plt.colorbar(
+      tmp,
+      ax=axes,
+      label="nT",
+      orientation="horizontal",
+      pad=0.12,
+      fraction=.03,
+   )
    plt.show()
diff --git a/doc/user_guide/forward_modelling/point.rst b/doc/user_guide/forward_modelling/point.rst
index 1b4942c9e..66bcd3637 100644
--- a/doc/user_guide/forward_modelling/point.rst
+++ b/doc/user_guide/forward_modelling/point.rst
@@ -113,7 +113,7 @@ Lets plot this gravitational field:
    fig = pygmt.Figure()
 
    maxabs = vd.maxabs(g_z)
-   pygmt.makecpt(cmap="polar", series=(-maxabs, maxabs), background=True)
+   pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True)
 
    fig.grdimage(
       region=(-250, 1250, -250, 1250),
@@ -239,8 +239,8 @@ Lets plot these results using :mod:`pygmt`:
 
    fig = pygmt.Figure()
    title = "Gravitational acceleration (downward)"
-   maxabs = vd.maxabs(g_z)*.95
-   pygmt.makecpt(cmap="polar", series=(-maxabs, maxabs), background=True)
+   maxabs = vd.maxabs(g_z)
+   pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True)
 
    fig.grdimage(
       region=(-72, -68, -46, -42),
diff --git a/doc/user_guide/forward_modelling/prism.rst b/doc/user_guide/forward_modelling/prism.rst
index de520f41b..8a7d8e8eb 100644
--- a/doc/user_guide/forward_modelling/prism.rst
+++ b/doc/user_guide/forward_modelling/prism.rst
@@ -90,16 +90,26 @@ point:
    for field in fields:
       results[field] = hm.prism_gravity(coordinates, prism, density, field=field)
 
+We can reshape the results into variables of an dataset:
+
+.. jupyter-execute::
+
+   grid = vd.make_xarray_grid(
+      coordinates,
+      tuple(results.values()),
+      data_names=results.keys(),
+      extra_coords_names="extra",
+   )
+   print(grid)
+
 Plot the results:
 
 .. jupyter-execute::
 
    import matplotlib.pyplot as plt
 
-   plt.pcolormesh(coordinates[0], coordinates[1], results["potential"])
+   grid.potential.plot(cbar_kwargs={"label":"J/kg"})
    plt.gca().set_aspect("equal")
-   plt.gca().ticklabel_format(style="sci", scilimits=(0, 0))
-   plt.colorbar(label="J/kg")
    plt.show()
 
 
@@ -107,36 +117,58 @@ Plot the results:
 
    fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8))
 
-   for field, ax in zip(("g_e", "g_n", "g_z"), axes):
-      tmp = ax.pcolormesh(coordinates[0], coordinates[1], results[field])
+   fields = ["g_e", "g_n", "g_z"]
+   maxabs = vd.maxabs(*[grid[f] for f in fields], percentile=99)
+   for field, ax in zip(fields, axes):
+      tmp = grid[field].plot(
+         ax=ax,
+         add_colorbar=False,
+         vmin=-maxabs,
+         vmax=maxabs,
+         cmap='RdBu_r',
+      )
       ax.set_aspect("equal")
       ax.set_title(field)
-      ax.ticklabel_format(style="sci", scilimits=(0, 0))
-      plt.colorbar(tmp, ax=ax, label="mGal", orientation="horizontal", pad=0.08)
+   fig.colorbar(tmp, ax=axes.ravel().tolist(), orientation="horizontal", label="mGal", fraction=.03, pad=0.08)
    plt.show()
 
 .. jupyter-execute::
 
    fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8))
 
-   for field, ax in zip(("g_ee", "g_nn", "g_zz"), axes):
-      tmp = ax.pcolormesh(coordinates[0], coordinates[1], results[field])
+   fields = ["g_ee", "g_nn", "g_zz"]
+   maxabs = vd.maxabs(*[grid[f] for f in fields], percentile=99)
+   for field, ax in zip(fields, axes):
+      tmp = grid[field].plot(
+         ax=ax,
+         add_colorbar=False,
+         vmin=-maxabs,
+         vmax=maxabs,
+         cmap='RdBu_r',
+      )
       ax.set_aspect("equal")
       ax.set_title(field)
-      ax.ticklabel_format(style="sci", scilimits=(0, 0))
-      plt.colorbar(tmp, ax=ax, label="Eotvos", orientation="horizontal", pad=0.08)
+   fig.colorbar(tmp, ax=axes.ravel().tolist(), orientation="horizontal", label="Eotvos", fraction=.03, pad=0.08)
    plt.show()
 
+
 .. jupyter-execute::
 
    fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8))
 
-   for field, ax in zip(("g_en", "g_ez", "g_nz"), axes):
-      tmp = ax.pcolormesh(coordinates[0], coordinates[1], results[field])
+   fields = ["g_en", "g_ez", "g_nz"]
+   maxabs = vd.maxabs(*[grid[f] for f in fields], percentile=99)
+   for field, ax in zip(fields, axes):
+      tmp = grid[field].plot(
+         ax=ax,
+         add_colorbar=False,
+         vmin=-maxabs,
+         vmax=maxabs,
+         cmap='RdBu_r',
+      )
       ax.set_aspect("equal")
       ax.set_title(field)
-      ax.ticklabel_format(style="sci", scilimits=(0, 0))
-      plt.colorbar(tmp, ax=ax, label="Eotvos", orientation="horizontal", pad=0.08)
+   fig.colorbar(tmp, ax=axes.ravel().tolist(), orientation="horizontal", label="Eotvos", fraction=.03, pad=0.08)
    plt.show()
 
 
@@ -164,8 +196,6 @@ height:
 
 .. jupyter-execute::
 
-   import verde as vd
-
    coordinates = vd.grid_coordinates(
        region=(0, 10e3, 0, 10e3), shape=(40, 40), extra_coords=0
    )
@@ -185,32 +215,11 @@ generated by the whole set of prisms on every computation point:
 
 Lets plot this gravitational field:
 
-.. jupyter-execute::
-   :hide-code:
-
-    import pygmt
-
-    # Needed so that displaying works on jupyter-sphinx and sphinx-gallery at
-    # the same time. Using PYGMT_USE_EXTERNAL_DISPLAY="false" in the Makefile
-    # for sphinx-gallery to work means that fig.show won't display anything here
-    # either.
-    pygmt.set_display(method="notebook")
-
-
 .. jupyter-execute::
 
-   import pygmt
    grid = vd.make_xarray_grid(
       coordinates, g_z, data_names="g_z", extra_coords_names="extra")
-   fig = pygmt.Figure()
-   fig.grdimage(
-      region=(0, 10e3, 0, 10e3),
-      projection="X10c",
-      grid=grid.g_z,
-      frame=["WSne", "x+leasting (m)", "y+lnorthing (m)"],
-      cmap='viridis',)
-   fig.colorbar(cmap=True, position="JMR", frame=["a2", "x+lmGal"])
-   fig.show()
+   grid.g_z.plot(cbar_kwargs={"label":"mGal"})
 
 
 Magnetic fields
@@ -259,39 +268,39 @@ points by choosing ``field="b"``:
 
    b_e, b_n, b_u = hm.prism_magnetic(coordinates, prisms, magnetization, field="b")
 
+
+We can reshape the results into variables of an dataset:
+
 .. jupyter-execute::
 
-   fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 6))
-
-   for ax, mag_component, title in zip(axes, (b_e, b_n, b_u), ("Be", "Bn", "Bu")):
-       maxabs = vd.maxabs(mag_component)
-       tmp = ax.pcolormesh(
-           coordinates[0],
-           coordinates[1],
-           mag_component,
-           vmin=-maxabs,
-           vmax=maxabs,
-           cmap="RdBu_r",
-       )
-       ax.contour(
-           coordinates[0],
-           coordinates[1],
-           mag_component,
-           colors="k",
-           linewidths=0.5,
-       )
-       ax.set_title(title)
-       ax.set_aspect("equal")
-       plt.colorbar(
-           tmp,
-           ax=ax,
-           orientation="horizontal",
-           label="nT",
-           pad=0.08,
-           aspect=42,
-           shrink=0.8,
-       )
+   grid = vd.make_xarray_grid(
+      coordinates,
+      data=(b_e, b_n, b_u),
+      data_names=["b_e", "b_n", "b_u"],
+      extra_coords_names="extra"
+   )
 
+.. jupyter-execute::
+
+   fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8))
+
+   maxabs = vd.maxabs(*[grid[v] for v in grid.variables], percentile=99)
+   for field, ax in zip(grid.variables, axes):
+      tmp = grid[field].plot(
+         ax=ax,
+         add_colorbar=False,
+         vmin=-maxabs,
+         vmax=maxabs,
+         cmap='RdBu_r',
+      )
+      grid[field].plot.contour(
+         ax=ax,
+         colors="k",
+         linewidths=0.5,
+      )
+      ax.set_aspect("equal")
+      ax.set_title(field)
+   fig.colorbar(tmp, ax=axes.ravel().tolist(), orientation="horizontal", label="nT", fraction=.03, pad=0.08)
    plt.show()
 
 
@@ -321,12 +330,12 @@ generated by these two prisms:
    maxabs = vd.maxabs(b_u)
 
    tmp = plt.pcolormesh(
-       coordinates[0], coordinates[1], b_u, vmin=-maxabs, vmax=maxabs, cmap="RdBu_r"
+      coordinates[0], coordinates[1], b_u, vmin=-maxabs, vmax=maxabs, cmap="RdBu_r"
    )
    plt.contour(coordinates[0], coordinates[1], b_u, colors="k", linewidths=0.5)
    plt.title("Bu")
    plt.gca().set_aspect("equal")
-   plt.colorbar(tmp, label="nT", pad=0.03, aspect=42, shrink=0.8)
+   plt.colorbar(tmp, label="nT", pad=0.03, shrink=0.8)
    plt.show()
 
 
@@ -415,16 +424,8 @@ Finally, let's plot the gravitational field:
    grid = vd.make_xarray_grid(
       coordinates, gravity, data_names="gravity", extra_coords_names="extra")
 
-   fig = pygmt.Figure()
-   title = "Gravitational acceleration of a layer of prisms"
-   fig.grdimage(
-      region=region_pad,
-      projection="X10c",
-      grid=grid.gravity,
-      frame=[f"WSne+t{title}", "x+leasting (m)", "y+lnorthing (m)"],
-      cmap='viridis',)
-   fig.colorbar(cmap=True, position="JMR", frame=["a.02", "x+lmGal"])
-   fig.show()
+   grid.gravity.plot(cbar_kwargs={"label":"mGal"})
+
 
 ----
 
diff --git a/doc/user_guide/forward_modelling/tesseroid.rst b/doc/user_guide/forward_modelling/tesseroid.rst
index 292ab2683..2250b95ee 100644
--- a/doc/user_guide/forward_modelling/tesseroid.rst
+++ b/doc/user_guide/forward_modelling/tesseroid.rst
@@ -132,7 +132,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]],
diff --git a/doc/user_guide/gravity_disturbance.rst b/doc/user_guide/gravity_disturbance.rst
index 2b2775822..93ecdcf02 100644
--- a/doc/user_guide/gravity_disturbance.rst
+++ b/doc/user_guide/gravity_disturbance.rst
@@ -89,7 +89,7 @@ Lets plot it:
    fig.basemap(frame=["af", "WEsn"])
    fig.colorbar(
        position="JCB+w10c",
-       frame=["af", 'y+l"mGal"', 'x+l"observed gravity"'],
+       frame=["af", 'y+lmGal', 'x+lObserved gravity'],
    )
    fig.coast(shorelines=True, resolution="c", area_thresh=1e4)
    fig.show()
@@ -122,7 +122,7 @@ And plot it:
    fig.basemap(frame=["af", "WEsn"])
    fig.colorbar(
        position="JCB+w10c",
-       frame=["af", 'y+l"mGal"', 'x+l"normal gravity"'],
+       frame=["af", 'y+lmGal', 'x+lNormal gravity'],
    )
    fig.coast(shorelines=True, resolution="c", area_thresh=1e4)
    fig.show()
@@ -140,10 +140,10 @@ And plot it:
 
    import verde as vd
 
-   maxabs = vd.maxabs(gravity_disturbance)
+   maxabs = vd.maxabs(gravity_disturbance, percentile=98)
 
    fig = pygmt.Figure()
-   pygmt.makecpt(series=[-maxabs, maxabs], cmap="polar+h")
+   pygmt.makecpt(series=[-maxabs, maxabs], cmap="balance+h0", background=True)
    fig.grdimage(
        gravity_disturbance,
        projection="W20c",
@@ -152,8 +152,8 @@ And plot it:
    )
    fig.basemap(frame=["af", "WEsn"])
    fig.colorbar(
-       position="JCB+w10c",
-       frame=["af", 'y+l"mGal"', 'x+l"gravity disturbance"'],
+       position="JCB+w10c+e",
+       frame=["a20f10", 'y+lmGal', 'x+lGravity disturbance'],
    )
    fig.coast(shorelines=True, resolution="c", area_thresh=1e4)
    fig.show()
diff --git a/doc/user_guide/igrf.rst b/doc/user_guide/igrf.rst
index 9eff60d82..ce0623ee0 100644
--- a/doc/user_guide/igrf.rst
+++ b/doc/user_guide/igrf.rst
@@ -8,6 +8,17 @@ magnetic field from 1900 until 5 years after the latest release (based on
 predictions of the secular variation). Harmonica allows calculating the 14th
 generation IGRF field with :class:`harmonica.IGRF14`. Here's how it works.
 
+.. jupyter-execute::
+   :hide-code:
+
+    import pygmt
+
+    # Needed so that displaying works on jupyter-sphinx and sphinx-gallery at
+    # the same time. Using PYGMT_USE_EXTERNAL_DISPLAY="false" in the Makefile
+    # for sphinx-gallery to work means that fig.show won't display anything here
+    # either.
+    pygmt.set_display(method="notebook")
+    
 .. jupyter-execute::
 
    import datetime
@@ -15,6 +26,7 @@ generation IGRF field with :class:`harmonica.IGRF14`. Here's how it works.
    import pygmt
    import boule as bl
    import harmonica as hm
+   import verde as vd
 
 All of the functionality is wrapped in the :class:`~harmonica.IGRF14` class.
 When creating an instance of it, we need to provide the date on which we want
@@ -147,24 +159,31 @@ We can plot the three components using :mod:`pygmt` on a nice map:
 .. jupyter-execute::
 
    fig = pygmt.Figure()
+
    for c in ["b_east", "b_north", "b_up"]:
-       fig.grdimage(
-           grid[c], cmap="polar+h", projection="W15c", region="g",
-       )
-       fig.coast(shorelines=True)
-       fig.colorbar(
-           position="JMR+ml+o0.5c",
-           frame=[
+      max_abs = vd.maxabs(grid[c])
+      pygmt.makecpt(
+         cmap="balance+h0",
+         series=[-max_abs, max_abs],
+         background=True,
+      )
+      fig.grdimage(
+         grid[c], cmap=True, projection="W15c", region="g",
+      )
+      fig.coast(shorelines=True)
+      fig.colorbar(
+         position="JMR+ml+o0.5c",
+         frame=[
                "a10000",
                f"x+l{grid[c].attrs['long_name']}",
                f"y+l{grid[c].attrs['units']}",
-           ]
-       )
-       if c == "b_east":
-           fig.basemap(frame=["a", f"+t{grid.attrs['title']}"])
-       else:
-           fig.basemap(frame="a")
-       fig.shift_origin(yshift="-h-0.5c")
+         ]
+      )
+      if c == "b_east":
+         fig.basemap(frame=["a", f"+t{grid.attrs['title']}"])
+      else:
+         fig.basemap(frame="a")
+      fig.shift_origin(yshift="-h-0.5c")
    fig.show()
 
 The grid spacing was calculated automatically to match the maximum degree of
@@ -201,33 +220,33 @@ with PyGMT the same way we did the vector components:
 
    fig = pygmt.Figure()
    cmaps = {
-       "intensity": "viridis",
-       "inclination": "polar+h",
-       "declination": "polar+h",
+      "intensity": "viridis",
+      "inclination": "balance+h0",
+      "declination": "balance+h0",
    }
    cb_annot = {
-       "intensity": "a10000",
-       "inclination": "a20+u\\260",
-       "declination": "a40+u\\260",
+      "intensity": "a10000",
+      "inclination": "a20+u\\260",
+      "declination": "a40+u\\260",
    }
    for c in ["intensity", "inclination", "declination"]:
-       fig.grdimage(
-           grid[c], cmap=cmaps[c], projection="W0/15c",
-       )
-       fig.coast(shorelines=True)
-       fig.colorbar(
-           position="JMR+ml+o0.5c",
-           frame=[
+      fig.grdimage(
+         grid[c], cmap=cmaps[c], projection="W0/15c",
+      )
+      fig.coast(shorelines=True)
+      fig.colorbar(
+         position="JMR+ml+o0.5c",
+         frame=[
                cb_annot[c],
                f"x+l{grid[c].attrs['long_name']}",
                f"y+l{grid[c].attrs['units']}",
-           ]
-       )
-       if c == "intensity":
-           fig.basemap(frame=["a", f"+t{grid.attrs['title']}"])
-       else:
-           fig.basemap(frame="a")
-       fig.shift_origin(yshift="-h-0.5c")
+         ]
+      )
+      if c == "intensity":
+         fig.basemap(frame=["a", f"+t{grid.attrs['title']}"])
+      else:
+         fig.basemap(frame="a")
+      fig.shift_origin(yshift="-h-0.5c")
    fig.show()
 
 We can clearly see the `South Atlantic Magnetic Anomaly
diff --git a/doc/user_guide/topographic_correction.rst b/doc/user_guide/topographic_correction.rst
index ec7b60c07..f85abbb1c 100644
--- a/doc/user_guide/topographic_correction.rst
+++ b/doc/user_guide/topographic_correction.rst
@@ -72,7 +72,7 @@ And plot it:
    maxabs = vd.maxabs(data.gravity_disturbance_mgal)
 
    fig = pygmt.Figure()
-   pygmt.makecpt(cmap="polar+h0", series=[-maxabs, maxabs])
+   pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs])
    fig.plot(
       x=data.longitude,
       y=data.latitude,
@@ -82,7 +82,7 @@ And plot it:
       projection="M15c",
       frame=['ag', 'WSen'],
    )
-   fig.colorbar(cmap=True, frame=["a50f25", "x+lgravity disturbance", "y+lmGal"])
+   fig.colorbar(cmap=True, frame=["a50f25", "x+lGravity disturbance", "y+lmGal"])
    fig.show()
 
 
@@ -122,10 +122,10 @@ We can now compute the Bouguer disturbance and plot it:
 
 .. jupyter-execute::
 
-   maxabs = vd.maxabs(bouguer_disturbance)
+   maxabs = vd.maxabs(bouguer_disturbance, percentile=98)
 
    fig = pygmt.Figure()
-   pygmt.makecpt(cmap="polar+h0", series=[-maxabs, maxabs])
+   pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True)
    fig.plot(
       x=data.longitude,
       y=data.latitude,
@@ -135,7 +135,11 @@ We can now compute the Bouguer disturbance and plot it:
       projection="M15c",
       frame=['ag', 'WSen'],
    )
-   fig.colorbar(cmap=True, frame=["a50f25", "x+lBouguer disturbance (with simple Bouguer correction)", "y+lmGal"])
+   fig.colorbar(
+      cmap=True,
+      frame=["a50f25", "x+lBouguer disturbance (with simple Bouguer correction)", "y+lmGal"],
+      position="+e",
+   )
    fig.show()
 
 
@@ -245,7 +249,7 @@ And plot it:
    maxabs = vd.maxabs(topo_free_disturbance)
 
    fig = pygmt.Figure()
-   pygmt.makecpt(cmap="polar+h0", series=[-maxabs, maxabs])
+   pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs])
    fig.plot(
       x=data.longitude,
       y=data.latitude,
@@ -258,6 +262,33 @@ And plot it:
    fig.colorbar(cmap=True, frame=["a50f25", "x+lTopography-free gravity disturbance", "y+lmGal"])
    fig.show()
 
+Compare the Bouguer and Topography-free disturbances
+----------------------------------------------------
+
+.. jupyter-execute::
+
+   difference = topo_free_disturbance-bouguer_disturbance
+
+   maxabs = vd.maxabs(difference, percentile=99)
+
+   fig = pygmt.Figure()
+   pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True)
+   fig.plot(
+      x=data.longitude,
+      y=data.latitude,
+      fill=difference,
+      cmap=True,
+      style="c3p",
+      projection="M15c",
+      frame=['ag', 'WSen'],
+   )
+   fig.colorbar(
+      cmap=True,
+      frame=["af", "x+lDifference between Bouguer and Topography-free disturbances", "y+lmGal"],
+      position="+e",
+   )
+   fig.show()
+
 ----
 
 .. grid:: 2
diff --git a/doc/user_guide/transformations.rst b/doc/user_guide/transformations.rst
index 15bef7def..4515e540b 100644
--- a/doc/user_guide/transformations.rst
+++ b/doc/user_guide/transformations.rst
@@ -30,11 +30,16 @@ And plot it:
 
 .. jupyter-execute::
 
-    tmp = magnetic_grid.plot(cmap="seismic", center=0, add_colorbar=False)
+    maxabs = vd.maxabs(magnetic_grid, percentile=99.9)
+    tmp = magnetic_grid.plot(
+        cmap="RdBu_r",
+        vmin=-maxabs,
+        vmax=maxabs,
+        cbar_kwargs={"label":"nT"},
+    )
     plt.gca().set_aspect("equal")
     plt.title("Magnetic anomaly grid")
     plt.gca().ticklabel_format(style="sci", scilimits=(0, 0))
-    plt.colorbar(tmp, label="nT")
     plt.show()
 
 .. seealso::
@@ -60,11 +65,16 @@ And plot it:
 
 .. jupyter-execute::
 
-    tmp = deriv_upward.plot(cmap="seismic", center=0, add_colorbar=False)
+    maxabs = vd.maxabs(deriv_upward, percentile=99.9)
+    tmp = deriv_upward.plot(
+        cmap="RdBu_r",
+        vmin=-maxabs,
+        vmax=maxabs,
+        cbar_kwargs={"label":"nT/m"},
+    )
     plt.gca().set_aspect("equal")
-    plt.title("Upward derivative of the magnetic anomaly")
+    plt.title("Magnetic anomaly grid")
     plt.gca().ticklabel_format(style="sci", scilimits=(0, 0))
-    plt.colorbar(tmp, label="nT/m")
     plt.show()
 
 
@@ -92,20 +102,33 @@ And plot them:
     fig, (ax1, ax2) = plt.subplots(
         nrows=1, ncols=2, sharey=True, figsize=(12, 8)
     )
+    maxabs = vd.maxabs(deriv_easting, deriv_northing, percentile=99)
 
-    cbar_kwargs=dict(
-        label="nT/m", orientation="horizontal", shrink=0.8, pad=0.08, aspect=42
+
+    tmp = deriv_easting.plot(
+        ax=ax1,
+        add_colorbar=False,
+        vmin=-maxabs,
+        vmax=maxabs,
+        cmap='RdBu_r',
     )
-    kwargs = dict(center=0, cmap="seismic", cbar_kwargs=cbar_kwargs)
 
-    tmp = deriv_easting.plot(ax=ax1, **kwargs)
-    tmp = deriv_northing.plot(ax=ax2, **kwargs)
+    tmp = deriv_northing.plot(
+        ax=ax2,
+        add_colorbar=False,
+        vmin=-maxabs,
+        vmax=maxabs,
+        cmap='RdBu_r',
+    )
 
+    ax1.set_aspect("equal")
+    ax2.set_aspect("equal")
     ax1.set_title("Easting derivative of the magnetic anomaly")
     ax2.set_title("Northing derivative of the magnetic anomaly")
-    for ax in (ax1, ax2):
-        ax.set_aspect("equal")
-        ax.ticklabel_format(style="sci", scilimits=(0, 0))
+    ax1.ticklabel_format(style="sci", scilimits=(0, 0))
+    ax2.ticklabel_format(style="sci", scilimits=(0, 0))
+
+    fig.colorbar(tmp, ax=(ax1, ax2), orientation="horizontal", label="nT/m", fraction=.03, pad=0.08)
     plt.show()
 
 By default, these two functions compute the horizontal derivatives using
@@ -130,22 +153,34 @@ frequency domain:
     fig, (ax1, ax2) = plt.subplots(
         nrows=1, ncols=2, sharey=True, figsize=(12, 8)
     )
+    maxabs = vd.maxabs(deriv_easting, deriv_northing, percentile=99)
+
 
-    cbar_kwargs=dict(
-        label="nT/m", orientation="horizontal", shrink=0.8, pad=0.08, aspect=42
+    tmp = deriv_easting.plot(
+        ax=ax1,
+        add_colorbar=False,
+        vmin=-maxabs,
+        vmax=maxabs,
+        cmap='RdBu_r',
     )
-    kwargs = dict(center=0, cmap="seismic", cbar_kwargs=cbar_kwargs)
 
-    tmp = deriv_easting.plot(ax=ax1, **kwargs)
-    tmp = deriv_northing.plot(ax=ax2, **kwargs)
+    tmp = deriv_northing.plot(
+        ax=ax2,
+        add_colorbar=False,
+        vmin=-maxabs,
+        vmax=maxabs,
+        cmap='RdBu_r',
+    )
 
+    ax1.set_aspect("equal")
+    ax2.set_aspect("equal")
     ax1.set_title("Easting derivative of the magnetic anomaly")
     ax2.set_title("Northing derivative of the magnetic anomaly")
-    for ax in (ax1, ax2):
-        ax.set_aspect("equal")
-        ax.ticklabel_format(style="sci", scilimits=(0, 0))
-    plt.show()
+    ax1.ticklabel_format(style="sci", scilimits=(0, 0))
+    ax2.ticklabel_format(style="sci", scilimits=(0, 0))
 
+    fig.colorbar(tmp, ax=(ax1, ax2), orientation="horizontal", label="nT/m", fraction=.03, pad=0.08)
+    plt.show()
 
 .. important::
 
@@ -189,11 +224,10 @@ Now we can plot it:
 
 .. jupyter-execute::
 
-    tmp = upward_continued.plot(cmap="seismic", center=0, add_colorbar=False)
+    tmp = upward_continued.plot(cmap="RdBu_r", center=0, cbar_kwargs={'label':'nT'})
     plt.gca().set_aspect("equal")
     plt.title("Upward continued magnetic anomaly to 1000m")
     plt.gca().ticklabel_format(style="sci", scilimits=(0, 0))
-    plt.colorbar(tmp, label="nT")
     plt.show()
 
 
@@ -257,18 +291,36 @@ And plot it:
 
 .. jupyter-execute::
 
-    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 8))
-    cbar_kwargs=dict(
-        label="nT", orientation="horizontal", shrink=0.8, pad=0.08, aspect=42
+    fig, (ax1, ax2) = plt.subplots(
+        nrows=1, ncols=2, sharey=True, figsize=(12, 8)
+    )
+    maxabs = vd.maxabs(magnetic_grid, rtp_grid, percentile=99)
+
+
+    tmp = magnetic_grid.plot(
+        ax=ax1,
+        add_colorbar=False,
+        vmin=-maxabs,
+        vmax=maxabs,
+        cmap='RdBu_r',
+    )
+
+    tmp = rtp_grid.plot(
+        ax=ax2,
+        add_colorbar=False,
+        vmin=-maxabs,
+        vmax=maxabs,
+        cmap='RdBu_r',
     )
-    magnetic_grid.plot(ax=ax1, cmap="seismic", center=0, cbar_kwargs=cbar_kwargs)
-    rtp_grid.plot(ax=ax2, cmap="seismic", center=0, cbar_kwargs=cbar_kwargs)
+
     ax1.set_aspect("equal")
     ax2.set_aspect("equal")
     ax1.set_title("Magnetic anomaly")
     ax2.set_title("Reduced to the pole")
     ax1.ticklabel_format(style="sci", scilimits=(0, 0))
     ax2.ticklabel_format(style="sci", scilimits=(0, 0))
+
+    fig.colorbar(tmp, ax=(ax1, ax2), orientation="horizontal", label="nT", fraction=.03, pad=0.08)
     plt.show()
 
 The reduction concentrated the Lightning Creek anomaly and the negative
@@ -297,11 +349,13 @@ example of what happens to the reduction using arbitrary values:
 
 .. jupyter-execute::
 
-    tmp = rtp_grid.plot(cmap="seismic", center=0, add_colorbar=False)
+    maxabs = vd.maxabs(rtp_grid, percentile=99.9)
+    tmp = rtp_grid.plot(
+        cmap="RdBu_r", vmin=-maxabs, vmax=maxabs, cbar_kwargs={'label':'nT'},
+    )
     plt.gca().set_aspect("equal")
     plt.title("Reduced to the pole with remanence")
     plt.gca().ticklabel_format(style="sci", scilimits=(0, 0))
-    plt.colorbar(tmp, label="nT")
     plt.show()
 
 
@@ -345,33 +399,23 @@ Let's plot the results side by side:
 
 .. jupyter-execute::
 
-    fig, (ax1, ax2, ax3) = plt.subplots(
-        nrows=1, ncols=3, sharey=True, figsize=(14, 8)
-    )
-
-    maxabs = vd.maxabs(magnetic_grid, magnetic_low, magnetic_high)
-    kwargs = dict(cmap="seismic", vmin=-maxabs, vmax=maxabs, add_colorbar=False)
-
-    tmp = magnetic_grid.plot(ax=ax1, **kwargs)
-    tmp = magnetic_low.plot(ax=ax2, **kwargs)
-    tmp = magnetic_high.plot(ax=ax3, **kwargs)
-
-    ax1.set_title("Original")
-    ax2.set_title("After low-pass filter")
-    ax3.set_title("After high-pass filter")
-    for ax in (ax1, ax2, ax3):
+    fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8))
+
+    grids = [magnetic_grid, magnetic_low, magnetic_high]
+    maxabs = vd.maxabs(*grids, percentile=99)
+    for grid, ax in zip(grids, axes):
+        tmp = grid.plot(
+            ax=ax,
+            add_colorbar=False,
+            vmin=-maxabs,
+            vmax=maxabs,
+            cmap='RdBu_r',
+        )
         ax.set_aspect("equal")
-        ax.ticklabel_format(style="sci", scilimits=(0, 0))
-
-    plt.colorbar(
-        tmp,
-        ax=[ax1, ax2, ax3],
-        label="nT",
-        orientation="horizontal",
-        aspect=42,
-        shrink=0.8,
-        pad=0.08,
-    )
+    axes[0].set_title("Original")
+    axes[1].set_title("After low-pass filter")
+    axes[2].set_title("After high-pass filter")
+    fig.colorbar(tmp, ax=axes.ravel().tolist(), orientation="horizontal", label="nT", fraction=.03, pad=0.08)
     plt.show()
 
 
@@ -408,11 +452,11 @@ And plot it:
 
 .. jupyter-execute::
 
-    tmp = tga_grid.plot(cmap="viridis", add_colorbar=False)
+    cpt_lims = vd.minmax(tga_grid, min_percentile=1, max_percentile=99)
+    tga_grid.plot(vmin=cpt_lims[0], vmax=cpt_lims[1],cbar_kwargs={'label':'nT/m'})
     plt.gca().set_aspect("equal")
     plt.title("Total gradient amplitude of the magnetic anomaly")
     plt.gca().ticklabel_format(style="sci", scilimits=(0, 0))
-    plt.colorbar(tmp, label="nT/m")
     plt.show()
 
 ----
diff --git a/environment.yml b/environment.yml
index 956cf3305..5a010b430 100644
--- a/environment.yml
+++ b/environment.yml
@@ -26,7 +26,6 @@ dependencies:
   - pytest
   - pytest-cov
   - coverage
-  - boule
   # Documentation requirements
   - sphinx==7.2.*
   - numpydoc==1.7.*
@@ -40,6 +39,7 @@ dependencies:
   - matplotlib
   - bordado>=0.4.0
   - ensaio>=0.5.0
+  - verde>=1.9.0
   - netcdf4
   - pygmt
   - gmt