From 9c237c64156bf8adc3b66e8413ce0d50d124e606 Mon Sep 17 00:00:00 2001 From: mdtanker Date: Wed, 19 Nov 2025 09:52:45 +0100 Subject: [PATCH 01/19] use balance for all diverging colormaps, update cpt limits and colorbar settings --- .../block_averaged_sources.py | 18 ++--- .../equivalent_sources/cartesian.py | 17 +++-- .../equivalent_sources/gradient_boosted.py | 16 ++-- .../equivalent_sources/spherical.py | 18 +++-- doc/gallery_src/forward/point_gravity.py | 6 +- doc/gallery_src/forward/prism_layer.py | 2 +- .../forward/prisms_topo_gravity.py | 10 ++- doc/gallery_src/forward/tesseroid_layer.py | 10 ++- doc/gallery_src/gravity_disturbance.py | 9 ++- .../gravity_disturbance_topofree.py | 9 ++- doc/gallery_src/isostatic_moho_airy.py | 7 +- .../transformations/reduction_to_pole.py | 10 +-- doc/gallery_src/transformations/tga.py | 12 +-- doc/gallery_src/transformations/tilt.py | 16 ++-- .../transformations/upward_continuation.py | 9 ++- .../transformations/upward_derivative.py | 8 +- .../equivalent_sources/block-averaged-eqs.rst | 4 +- .../eq-sources-spherical.rst | 2 +- .../eqs-parameters-estimation.rst | 2 +- .../gradient-boosted-eqs.rst | 2 +- doc/user_guide/equivalent_sources/index.rst | 4 +- doc/user_guide/forward_modelling/point.rst | 6 +- .../forward_modelling/tesseroid.rst | 1 + doc/user_guide/gravity_disturbance.rst | 2 +- doc/user_guide/igrf.rst | 76 ++++++++++--------- doc/user_guide/topographic_correction.rst | 6 +- 26 files changed, 161 insertions(+), 121 deletions(-) diff --git a/doc/gallery_src/equivalent_sources/block_averaged_sources.py b/doc/gallery_src/equivalent_sources/block_averaged_sources.py index 0a50e41aa..150630a39 100644 --- a/doc/gallery_src/equivalent_sources/block_averaged_sources.py +++ b/doc/gallery_src/equivalent_sources/block_averaged_sources.py @@ -41,6 +41,7 @@ """ import ensaio +import numpy as np import pandas as pd import pygmt import pyproj @@ -106,14 +107,14 @@ title = "Observed magnetic anomaly data" +# Get the 99.9th percentile of the absolute value of the point data to use as color +# scale limits +cpt_lim = np.quantile(np.abs(data.total_field_anomaly_nt), 0.999) + # 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), + cmap="balance+h0", + series=[-cpt_lim, cpt_lim], background=True, ) @@ -129,8 +130,7 @@ cmap=True, ) -fig.colorbar(cmap=True, frame=["a400f100", "x+lnT"]) - +fig.colorbar(cmap=True, frame=["x+lnT"], position="+e") fig.shift_origin(xshift=fig_width + 1) title = "Gridded and upward-continued" @@ -142,6 +142,6 @@ cmap=True, ) -fig.colorbar(cmap=True, frame=["a400f100", "x+lnT"]) +fig.colorbar(cmap=True, frame=["x+lnT"], position="+e") fig.show() diff --git a/doc/gallery_src/equivalent_sources/cartesian.py b/doc/gallery_src/equivalent_sources/cartesian.py index ad67084df..c789ffcc1 100644 --- a/doc/gallery_src/equivalent_sources/cartesian.py +++ b/doc/gallery_src/equivalent_sources/cartesian.py @@ -31,6 +31,7 @@ """ import ensaio +import numpy as np import pandas as pd import pygmt import pyproj @@ -95,14 +96,14 @@ title = "Observed magnetic anomaly data" +# Get the 99.9th percentile of the absolute value of the point data to use as color +# scale limits +cpt_lim = np.quantile(np.abs(data.total_field_anomaly_nt), 0.999) + # 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), + cmap="balance+h0", + series=[-cpt_lim, cpt_lim], background=True, ) @@ -118,7 +119,7 @@ cmap=True, ) -fig.colorbar(cmap=True, frame=["a400f100", "x+lnT"]) +fig.colorbar(cmap=True, frame=["x+lnT"], position="+e") fig.shift_origin(xshift=fig_width + 1) @@ -131,6 +132,6 @@ cmap=True, ) -fig.colorbar(cmap=True, frame=["a400f100", "x+lnT"]) +fig.colorbar(cmap=True, frame=["x+lnT"], position="+e") fig.show() diff --git a/doc/gallery_src/equivalent_sources/gradient_boosted.py b/doc/gallery_src/equivalent_sources/gradient_boosted.py index 72b5881f6..1a1260d46 100644 --- a/doc/gallery_src/equivalent_sources/gradient_boosted.py +++ b/doc/gallery_src/equivalent_sources/gradient_boosted.py @@ -31,6 +31,7 @@ import boule as bl import ensaio +import numpy as np import pandas as pd import pygmt import pyproj @@ -115,13 +116,14 @@ title = "Observed gravity disturbance data" +# Get the 99.9th percentile of the absolute value of the point data to use as color +# scale limits +cpt_lim = np.quantile(np.abs(data.gravity_disturbance), 0.999) + # Make colormap of data pygmt.makecpt( - cmap="vik", - series=( - -data.gravity_disturbance.quantile(0.99), - data.gravity_disturbance.quantile(0.99), - ), + cmap="balance+h0", + series=[-cpt_lim, cpt_lim], background=True, ) @@ -137,7 +139,7 @@ cmap=True, ) -fig.colorbar(cmap=True, frame=["a50f25", "x+lmGal"]) +fig.colorbar(cmap=True, frame=["x+lmGal"], position="+e") fig.shift_origin(xshift=fig_width + 1) @@ -150,6 +152,6 @@ cmap=True, ) -fig.colorbar(cmap=True, frame=["a50f25", "x+lmGal"]) +fig.colorbar(cmap=True, frame=["x+lmGal"], position="+e") fig.show() diff --git a/doc/gallery_src/equivalent_sources/spherical.py b/doc/gallery_src/equivalent_sources/spherical.py index d41b43906..255329625 100644 --- a/doc/gallery_src/equivalent_sources/spherical.py +++ b/doc/gallery_src/equivalent_sources/spherical.py @@ -80,17 +80,19 @@ # Plot observed and gridded gravity disturbance fig = pygmt.Figure() +# Get the 99.9th percentile of the absolute value of the point data to use as color +# scale limits +cpt_lim = np.quantile(np.abs(gravity_disturbance), 0.999) + # 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), + cmap="balance+h0", + series=[-cpt_lim, cpt_lim], background=True, ) +title = "Observed gravity disturbance data" + fig.plot( projection="M10c", region=region, @@ -102,7 +104,7 @@ cmap=True, ) -fig.colorbar(cmap=True, frame=["a100f50", "x+lmGal"]) +fig.colorbar(cmap=True, frame=["x+lmGal"], position="+e") fig.shift_origin(xshift="w+3c") @@ -112,6 +114,6 @@ cmap=True, ) -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/forward/point_gravity.py b/doc/gallery_src/forward/point_gravity.py index 38c4a2214..b822188b5 100644 --- a/doc/gallery_src/forward/point_gravity.py +++ b/doc/gallery_src/forward/point_gravity.py @@ -56,7 +56,7 @@ maxabs = vd.maxabs(gravity) * 0.80 -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( @@ -70,6 +70,8 @@ fig.plot(x=easting, y=northing, style="c0.2c", fill="grey") -fig.colorbar(cmap=True, position="JMR", frame=["a.6f.2", "x+lmGal"]) +fig.colorbar(cmap=True, position="JMR+e", frame=["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..4cf887b73 100644 --- a/doc/gallery_src/forward/prism_layer.py +++ b/doc/gallery_src/forward/prism_layer.py @@ -62,6 +62,6 @@ 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_layer.py b/doc/gallery_src/forward/tesseroid_layer.py index 4d26564df..5e42984c6 100644 --- a/doc/gallery_src/forward/tesseroid_layer.py +++ b/doc/gallery_src/forward/tesseroid_layer.py @@ -55,13 +55,19 @@ # Plot gravity field fig = pygmt.Figure() -maxabs = vd.maxabs(gravity.g_z) -pygmt.makecpt(cmap="polar", series=(-maxabs, maxabs)) + +# Get the max absolute value to use as color scale limits +cpt_lims = vd.maxabs(gravity.g_z) + +# Make colormap of data +pygmt.makecpt(cmap="balance+h0", series=[-cpt_lims, cpt_lims]) fig.grdimage( gravity.g_z, projection="M15c", nan_transparent=True, + cmap=True, ) + fig.basemap(frame=True) fig.colorbar(frame='af+l"Gravity [mGal]"', position="JCR") fig.coast(shorelines="0.5p,black", borders=["1/0.5p,black"]) diff --git a/doc/gallery_src/gravity_disturbance.py b/doc/gallery_src/gravity_disturbance.py index 398771d00..df48156ae 100644 --- a/doc/gallery_src/gravity_disturbance.py +++ b/doc/gallery_src/gravity_disturbance.py @@ -18,6 +18,7 @@ import boule as bl import ensaio +import numpy as np import pygmt import xarray as xr @@ -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 +cpt_lims = np.quantile(np.abs(disturbance), 0.999) + +# Make colormap of data +pygmt.makecpt(cmap="balance+h0", series=[-cpt_lims, cpt_lims], 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 07601d7c3..94e353090 100644 --- a/doc/gallery_src/gravity_disturbance_topofree.py +++ b/doc/gallery_src/gravity_disturbance_topofree.py @@ -28,6 +28,7 @@ import boule as bl import ensaio +import numpy as np import pygmt import xarray as xr @@ -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 +cpt_lims = np.quantile(np.abs(disturbance_topofree), 0.99) + +# Make colormap of data +pygmt.makecpt(cmap="balance+h0", series=[-cpt_lims, cpt_lims], 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 cc1e0180a..92e6590b1 100644 --- a/doc/gallery_src/isostatic_moho_airy.py +++ b/doc/gallery_src/isostatic_moho_airy.py @@ -56,10 +56,11 @@ 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) +# Make colormap of data +pygmt.makecpt(cmap="viridis", series=[moho.to_numpy().min(), moho.to_numpy().max()]) title = "Airy isostatic Moho depth of Africa" @@ -73,6 +74,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 8a785a36b..d0e5944d2 100644 --- a/doc/gallery_src/transformations/reduction_to_pole.py +++ b/doc/gallery_src/transformations/reduction_to_pole.py @@ -52,9 +52,9 @@ # 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 = 0.5 * vd.maxabs(magnetic_grid, rtp_grid) + pygmt.makecpt(cmap="balance+h0", series=[-cpt_lim, cpt_lim], background=True) with fig.set_panel(panel=0): # Plot magnetic anomaly grid fig.grdimage( @@ -65,7 +65,7 @@ # Add colorbar fig.colorbar( frame='af+l"Magnetic anomaly [nT]"', - position="JBC+h+o0/1c+e", + position="JBC+h+o0/1c+ef", ) with fig.set_panel(panel=1): # Plot upward reduced to the pole grid @@ -73,6 +73,6 @@ # Add colorbar fig.colorbar( frame='af+l"Reduced to the pole [nT]"', - position="JBC+h+o0/1c+e", + position="JBC+h+o0/1c+ef", ) fig.show() diff --git a/doc/gallery_src/transformations/tga.py b/doc/gallery_src/transformations/tga.py index 21c710e5a..d5fe8396e 100644 --- a/doc/gallery_src/transformations/tga.py +++ b/doc/gallery_src/transformations/tga.py @@ -45,8 +45,8 @@ 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) + cpt_lim = vd.maxabs(magnetic_grid) + pygmt.makecpt(cmap="balance+h0", series=[-cpt_lim, cpt_lim], background=True) # Plot magnetic anomaly grid fig.grdimage( grid=magnetic_grid, @@ -56,17 +56,17 @@ # Add colorbar fig.colorbar( frame='af+l"Magnetic anomaly [nT]"', - position="JBC+h+o0/1c+e", + position="JBC+h+o0/1c", ) 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) + cpt_lim = 0.6 * vd.maxabs(tga) + pygmt.makecpt(cmap="balance+h0", series=[0, cpt_lim], background=True) # Plot total gradient amplitude fig.grdimage(grid=tga, projection="X?", cmap=True) # Add colorbar fig.colorbar( frame='af+l"Total Gradient Amplitude [nT/m]"', - position="JBC+h+o0/1c+e", + 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 9f9c1b579..e62234dd7 100644 --- a/doc/gallery_src/transformations/tilt.py +++ b/doc/gallery_src/transformations/tilt.py @@ -79,10 +79,10 @@ sharey="l", margins=["1c", "1c"], ): - scale = 0.5 * vd.maxabs(magnetic_grid, rtp_grid) + cpt_lim = 0.5 * vd.maxabs(magnetic_grid, rtp_grid) 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=[-cpt_lim, cpt_lim], background=True) # Plot magnetic anomaly grid fig.grdimage( grid=magnetic_grid, @@ -92,7 +92,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=[-cpt_lim, cpt_lim], background=True) # Plot reduced to the pole magnetic anomaly grid fig.grdimage( grid=rtp_grid, @@ -104,13 +104,13 @@ label = "nT" fig.colorbar( frame=f"af+l{label}", - position="JMR+o1/-0.25c+e", + position="JMR+o1/-0.25c+ef", ) - scale = 0.6 * vd.maxabs(tilt_grid, tilt_rtp_grid) + cpt_lim = vd.maxabs(tilt_grid, tilt_rtp_grid) 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=[-cpt_lim, cpt_lim], background=True) # Plot tilt fig.grdimage( grid=tilt_grid, @@ -120,7 +120,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=[-cpt_lim, cpt_lim], background=True) # Plot tilt fig.grdimage( grid=tilt_rtp_grid, @@ -132,6 +132,6 @@ label = "rad" fig.colorbar( frame=f"af+l{label}", - position="JMR+o1/-0.25c+e", + position="JMR+o1/-0.25c", ) fig.show() diff --git a/doc/gallery_src/transformations/upward_continuation.py b/doc/gallery_src/transformations/upward_continuation.py index db6ce9829..edd76893d 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 xrft @@ -44,9 +45,9 @@ # 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) * 0.6 + pygmt.makecpt(cmap="balance+h0", series=[-cpt_lim, cpt_lim], background=True) with fig.set_panel(panel=0): # Plot magnetic anomaly grid fig.grdimage( @@ -65,6 +66,6 @@ # Add colorbar fig.colorbar( frame='af+l"Upward continued to 1000m [nT]"', - position="JBC+h+o0/1c+e", + position="JBC+h+o0/1c", ) fig.show() diff --git a/doc/gallery_src/transformations/upward_derivative.py b/doc/gallery_src/transformations/upward_derivative.py index 2d48b89db..30c74618a 100644 --- a/doc/gallery_src/transformations/upward_derivative.py +++ b/doc/gallery_src/transformations/upward_derivative.py @@ -46,8 +46,8 @@ 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) + cpt_lim = vd.maxabs(magnetic_grid) * 0.6 + pygmt.makecpt(cmap="balance+h0", series=[-cpt_lim, cpt_lim], background=True) # Plot magnetic anomaly grid fig.grdimage( grid=magnetic_grid, @@ -61,8 +61,8 @@ ) 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) + cpt_lim = vd.maxabs(deriv_upward) * 0.6 + pygmt.makecpt(cmap="balance+h0", series=[-cpt_lim, cpt_lim], background=True) # Plot upward derivative fig.grdimage(grid=deriv_upward, projection="X?", cmap=True) # Add colorbar diff --git a/doc/user_guide/equivalent_sources/block-averaged-eqs.rst b/doc/user_guide/equivalent_sources/block-averaged-eqs.rst index c27e08dbb..77c5f28f4 100644 --- a/doc/user_guide/equivalent_sources/block-averaged-eqs.rst +++ b/doc/user_guide/equivalent_sources/block-averaged-eqs.rst @@ -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, ) @@ -182,7 +182,7 @@ 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) diff --git a/doc/user_guide/equivalent_sources/eq-sources-spherical.rst b/doc/user_guide/equivalent_sources/eq-sources-spherical.rst index 9d98993a3..32b91a654 100644 --- a/doc/user_guide/equivalent_sources/eq-sources-spherical.rst +++ b/doc/user_guide/equivalent_sources/eq-sources-spherical.rst @@ -167,7 +167,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,)) title = "Block-median reduced gravity disturbance" fig.plot( diff --git a/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst b/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst index 22b5c2c57..615da657f 100644 --- a/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst +++ b/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst @@ -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,)) title = "Gravity disturbance with first guess" diff --git a/doc/user_guide/equivalent_sources/gradient-boosted-eqs.rst b/doc/user_guide/equivalent_sources/gradient-boosted-eqs.rst index 5cb53f00b..2aa6a2784 100644 --- a/doc/user_guide/equivalent_sources/gradient-boosted-eqs.rst +++ b/doc/user_guide/equivalent_sources/gradient-boosted-eqs.rst @@ -173,7 +173,7 @@ And 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,)) title = "Observed gravity disturbance data" with pygmt.config(FONT_TITLE="14p"): diff --git a/doc/user_guide/equivalent_sources/index.rst b/doc/user_guide/equivalent_sources/index.rst index 266a9e3e7..bbdec13d6 100644 --- a/doc/user_guide/equivalent_sources/index.rst +++ b/doc/user_guide/equivalent_sources/index.rst @@ -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]) title="Predicted gravity disturbance" with pygmt.config(FONT_TITLE="14p"): fig.plot( @@ -204,7 +204,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, diff --git a/doc/user_guide/forward_modelling/point.rst b/doc/user_guide/forward_modelling/point.rst index 004058c0d..8f6999474 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/tesseroid.rst b/doc/user_guide/forward_modelling/tesseroid.rst index 21c4a5c47..4db314fcc 100644 --- a/doc/user_guide/forward_modelling/tesseroid.rst +++ b/doc/user_guide/forward_modelling/tesseroid.rst @@ -111,6 +111,7 @@ And finally plot the computed gravitational field .. jupyter-execute:: import pygmt + grid = vd.make_xarray_grid( coordinates, gravity, data_names="gravity", extra_coords_names="extra") diff --git a/doc/user_guide/gravity_disturbance.rst b/doc/user_guide/gravity_disturbance.rst index a27da1c5e..6aef78354 100644 --- a/doc/user_guide/gravity_disturbance.rst +++ b/doc/user_guide/gravity_disturbance.rst @@ -143,7 +143,7 @@ And plot it: maxabs = vd.maxabs(gravity_disturbance) fig = pygmt.Figure() - pygmt.makecpt(series=[-maxabs, maxabs], cmap="polar+h") + pygmt.makecpt(series=[-maxabs, maxabs], cmap="balance+h0") fig.grdimage( gravity_disturbance, projection="W20c", diff --git a/doc/user_guide/igrf.rst b/doc/user_guide/igrf.rst index 9eff60d82..eef99ee2f 100644 --- a/doc/user_guide/igrf.rst +++ b/doc/user_guide/igrf.rst @@ -15,6 +15,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 +148,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 +209,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..0e1fb0080 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, @@ -125,7 +125,7 @@ We can now compute the Bouguer disturbance and plot it: maxabs = vd.maxabs(bouguer_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, @@ -245,7 +245,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, From 15c844a9d13fbc4267f75bb460863c8af4c6f9fc Mon Sep 17 00:00:00 2001 From: mdtanker Date: Wed, 19 Nov 2025 09:54:37 +0100 Subject: [PATCH 02/19] use pygmt instead of matplotlib for documentation --- doc/user_guide/forward_modelling/dipole.rst | 70 ++++- doc/user_guide/forward_modelling/prism.rst | 289 +++++++++++++------ doc/user_guide/transformations.rst | 300 +++++++++++++------- 3 files changed, 471 insertions(+), 188 deletions(-) diff --git a/doc/user_guide/forward_modelling/dipole.rst b/doc/user_guide/forward_modelling/dipole.rst index e7a5e9a6d..3d42046ff 100644 --- a/doc/user_guide/forward_modelling/dipole.rst +++ b/doc/user_guide/forward_modelling/dipole.rst @@ -92,20 +92,70 @@ 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:: + :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") - import matplotlib.pyplot as plt - fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8)) +.. jupyter-execute:: + + import pygmt + + fig = pygmt.Figure() + + maxabs = vd.maxabs(grid.b_e) + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + fig.grdimage( + grid=grid.b_e, + projection="X10c", + cmap=True, + frame=["WSne+tEasting component", "xa","ya"], + ) + fig.colorbar(frame='+lnT') + + fig.shift_origin(xshift="11c") + + maxabs = vd.maxabs(grid.b_n) + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + fig.grdimage( + grid=grid.b_n, + projection="X10c", + cmap=True, + frame=["wSne+tNorthing component", "xa","ya"], + ) + fig.colorbar(frame='+lnT') + + fig.shift_origin(xshift="11c") + + maxabs = vd.maxabs(grid.b_u) + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + fig.grdimage( + grid=grid.b_u, + projection="X10c", + cmap=True, + frame=["wSnE+tUpward component", "xa","ya"], + ) + fig.colorbar(frame='+lnT') - 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]) - 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) - plt.show() + fig.show() ---- diff --git a/doc/user_guide/forward_modelling/prism.rst b/doc/user_guide/forward_modelling/prism.rst index 1794defb0..08a75112f 100644 --- a/doc/user_guide/forward_modelling/prism.rst +++ b/doc/user_guide/forward_modelling/prism.rst @@ -90,54 +90,167 @@ 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:: + :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 matplotlib.pyplot as plt + import pygmt - plt.pcolormesh(coordinates[0], coordinates[1], results["potential"]) - plt.gca().set_aspect("equal") - plt.gca().ticklabel_format(style="sci", scilimits=(0, 0)) - plt.colorbar(label="J/kg") - plt.show() + fig = pygmt.Figure() + + fig.grdimage( + projection="X10c", + grid=grid.potential, + frame=["a", "x+leasting (m)", "y+lnorthing (m)"], + cmap="viridis", + ) + + fig.colorbar(cmap=True, position="JMR", frame=["x+lJ kg@+-1@+"]) + fig.show() .. jupyter-execute:: - fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8)) + fig = pygmt.Figure() + + maxabs = vd.maxabs(grid.g_e) + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + fig.grdimage( + grid=grid.g_e, + projection="X10c", + cmap=True, + frame=["WSne+tEasting component", "xa","ya"], + ) + fig.colorbar(frame='+lnT') - for field, ax in zip(("g_e", "g_n", "g_z"), axes): - tmp = ax.pcolormesh(coordinates[0], coordinates[1], results[field]) - 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) - plt.show() + fig.shift_origin(xshift="11c") + + maxabs = vd.maxabs(grid.g_n) + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + fig.grdimage( + grid=grid.g_n, + projection="X10c", + cmap=True, + frame=["wSne+tNorthing component", "xa","ya"], + ) + fig.colorbar(frame='+lnT') + + fig.shift_origin(xshift="11c") + + maxabs = vd.maxabs(grid.g_z) + pygmt.makecpt(cmap="balance+h0", series=[0, maxabs], background=True) + fig.grdimage( + grid=grid.g_z, + projection="X10c", + cmap=True, + frame=["wSnE+tDownward component", "xa","ya"], + ) + fig.colorbar(frame='+lnT') + + fig.show() .. jupyter-execute:: - fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8)) + fig = pygmt.Figure() - for field, ax in zip(("g_ee", "g_nn", "g_zz"), axes): - tmp = ax.pcolormesh(coordinates[0], coordinates[1], results[field]) - 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) - plt.show() + maxabs = vd.maxabs(grid.g_ee) + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + fig.grdimage( + grid=grid.g_ee, + projection="X10c", + cmap=True, + frame=["WSne+tEasting-easting tensor", "xa","ya"], + ) + fig.colorbar(frame='+lEotvos') + + fig.shift_origin(xshift="11c") + + maxabs = vd.maxabs(grid.g_nn) + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + fig.grdimage( + grid=grid.g_nn, + projection="X10c", + cmap=True, + frame=["wSne+tNorthing-northing tensor", "xa","ya"], + ) + fig.colorbar(frame='+lEotvos') + + fig.shift_origin(xshift="11c") + + maxabs = vd.maxabs(grid.g_zz) + pygmt.makecpt(cmap="balance+h0", series=[0, maxabs], background=True) + fig.grdimage( + grid=grid.g_zz, + projection="X10c", + cmap=True, + frame=["wSnE+tDownward-downward tensor", "xa","ya"], + ) + fig.colorbar(frame='+lEotvos') + + fig.show() .. jupyter-execute:: - fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8)) + fig = pygmt.Figure() + + maxabs = vd.maxabs(grid.g_en) + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + fig.grdimage( + grid=grid.g_en, + projection="X10c", + cmap=True, + frame=["WSne+tEasting-northing tensor", "xa","ya"], + ) + fig.colorbar(frame='+lEotvos') + + fig.shift_origin(xshift="11c") - for field, ax in zip(("g_en", "g_ez", "g_nz"), axes): - tmp = ax.pcolormesh(coordinates[0], coordinates[1], results[field]) - 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) - plt.show() + maxabs = vd.maxabs(grid.g_ez) + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + fig.grdimage( + grid=grid.g_ez, + projection="X10c", + cmap=True, + frame=["wSne+tEasting-downward tensor", "xa","ya"], + ) + fig.colorbar(frame='+lEotvos') + + fig.shift_origin(xshift="11c") + + maxabs = vd.maxabs(grid.g_nz) + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + fig.grdimage( + grid=grid.g_nz, + projection="X10c", + cmap=True, + frame=["wSnE+tNorthing-downward tensor", "xa","ya"], + ) + fig.colorbar(frame='+lEotvos') + + fig.show() Passing multiple prisms @@ -186,20 +299,7 @@ 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() @@ -259,40 +359,57 @@ 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, - ) - - plt.show() + 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 = pygmt.Figure() + + maxabs = vd.maxabs(grid.b_e) + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + fig.grdimage( + grid=grid.b_e, + projection="X10c", + cmap=True, + frame=["WSne+tEasting component", "xa","ya"], + ) + fig.colorbar(frame='+lnT') + + fig.shift_origin(xshift="11c") + + maxabs = vd.maxabs(grid.b_n) + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + fig.grdimage( + grid=grid.b_n, + projection="X10c", + cmap=True, + frame=["wSne+tNorthing component", "xa","ya"], + ) + fig.colorbar(frame='+lnT') + + fig.shift_origin(xshift="11c") + + maxabs = vd.maxabs(grid.b_u) + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + fig.grdimage( + grid=grid.b_u, + projection="X10c", + cmap=True, + frame=["wSnE+tUpward component", "xa","ya"], + ) + fig.colorbar(frame='+lnT') + + fig.show() Alternatively, we can compute just a single component by choosing ``field`` to @@ -318,16 +435,22 @@ generated by these two prisms: .. jupyter-execute:: - maxabs = vd.maxabs(b_u) + fig = pygmt.Figure() + + grid = vd.make_xarray_grid( + coordinates, b_u, data_names="b_u", extra_coords_names="extra") - tmp = plt.pcolormesh( - coordinates[0], coordinates[1], b_u, vmin=-maxabs, vmax=maxabs, cmap="RdBu_r" + maxabs = vd.maxabs(grid.b_u) + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + fig.grdimage( + grid=grid.b_u, + projection="X10c", + cmap=True, + frame=["WSne+tUpward component", "xa","ya"], ) - 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.show() + fig.colorbar(frame='+lnT') + + fig.show() .. _prism_layer: diff --git a/doc/user_guide/transformations.rst b/doc/user_guide/transformations.rst index ce9ac9734..769b717d5 100644 --- a/doc/user_guide/transformations.rst +++ b/doc/user_guide/transformations.rst @@ -25,15 +25,39 @@ We can load the data file using :mod:`xarray`: And plot it: .. jupyter-execute:: + :hide-code: - import matplotlib.pyplot as plt + import pygmt - tmp = magnetic_grid.plot(cmap="seismic", center=0, add_colorbar=False) - 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() + # 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 + import verde as vd + + fig = pygmt.Figure() + + maxabs = vd.maxabs(magnetic_grid) * .6 + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + + fig.grdimage( + magnetic_grid, + projection="X15c", + cmap=True, + frame=["af", "WeSn+tMagnetic anomaly"] + ) + fig.colorbar( + position="JCB+e", + frame=["af", 'x+lnT'], + ) + + fig.show() .. seealso:: @@ -68,12 +92,23 @@ needed by the :func:`xrft.pad` function): .. jupyter-execute:: - tmp = magnetic_grid_padded.plot(cmap="seismic", center=0, add_colorbar=False) - plt.gca().set_aspect("equal") - plt.title("Padded magnetic anomaly grid") - plt.gca().ticklabel_format(style="sci", scilimits=(0, 0)) - plt.colorbar(tmp, label="nT") - plt.show() + fig = pygmt.Figure() + + maxabs = vd.maxabs(magnetic_grid_padded) * .6 + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + + fig.grdimage( + magnetic_grid_padded, + projection="X15c", + cmap=True, + frame=["af", "WeSn+tPadded magnetic anomaly"] + ) + fig.colorbar( + position="JCB+e", + frame=["af", 'x+lnT'], + ) + + fig.show() Now that we have the padded grid, we can apply any grid transformation. @@ -103,12 +138,23 @@ And plot it: .. jupyter-execute:: - tmp = deriv_upward.plot(cmap="seismic", center=0, add_colorbar=False) - plt.gca().set_aspect("equal") - plt.title("Upward derivative of the magnetic anomaly") - plt.gca().ticklabel_format(style="sci", scilimits=(0, 0)) - plt.colorbar(tmp, label="nT/m") - plt.show() + fig = pygmt.Figure() + + maxabs = vd.maxabs(deriv_upward) * .5 + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + + fig.grdimage( + deriv_upward, + projection="X15c", + cmap=True, + frame=["af", "WeSn+tUpward derivative of the magnetic anomaly"] + ) + fig.colorbar( + position="JCB+e", + frame=["af", 'x+lnT/m'], + ) + + fig.show() Horizontal derivatives @@ -132,24 +178,35 @@ And plot them: .. jupyter-execute:: - fig, (ax1, ax2) = plt.subplots( - nrows=1, ncols=2, sharey=True, figsize=(12, 8) + fig = pygmt.Figure() + + maxabs = vd.maxabs(deriv_easting, deriv_northing) * .4 + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + + fig.grdimage( + deriv_easting, + projection="X15c", + cmap=True, + frame=["af", "WeSn+tEasting derivative of the magnetic anomaly"] ) - cbar_kwargs=dict( - label="nT/m", orientation="horizontal", shrink=0.8, pad=0.08, aspect=42 + fig.shift_origin(xshift="16c") + + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + + fig.grdimage( + deriv_northing, + projection="X15c", + cmap=True, + frame=["af", "wESn+tNorthing derivative of the magnetic anomaly"] ) - kwargs = dict(center=0, cmap="seismic", cbar_kwargs=cbar_kwargs) - tmp = deriv_easting.plot(ax=ax1, **kwargs) - tmp = deriv_northing.plot(ax=ax2, **kwargs) + fig.colorbar( + position="JBC+h+e+o-8c/1c+w15c/.8c", + frame=["af", 'x+lnT/m'], + ) - 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() + fig.show() By default, these two functions compute the horizontal derivatives using central finite differences methods. We can choose to use either the finite @@ -172,25 +229,35 @@ frequency domain: .. jupyter-execute:: - fig, (ax1, ax2) = plt.subplots( - nrows=1, ncols=2, sharey=True, figsize=(12, 8) - ) + fig = pygmt.Figure() + + maxabs = vd.maxabs(deriv_easting, deriv_northing) * .4 + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - cbar_kwargs=dict( - label="nT/m", orientation="horizontal", shrink=0.8, pad=0.08, aspect=42 + fig.grdimage( + deriv_easting, + projection="X15c", + cmap=True, + frame=["af", "WeSn+tEasting derivative of the magnetic anomaly"] ) - kwargs = dict(center=0, cmap="seismic", cbar_kwargs=cbar_kwargs) - tmp = deriv_easting.plot(ax=ax1, **kwargs) - tmp = deriv_northing.plot(ax=ax2, **kwargs) + fig.shift_origin(xshift="16c") + + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - 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() + fig.grdimage( + deriv_northing, + projection="X15c", + cmap=True, + frame=["af", "wESn+tNorthing derivative of the magnetic anomaly"] + ) + + fig.colorbar( + position="JBC+h+e+o-8c/1c+w15c/.8c", + frame=["af", 'x+lnT/m'], + ) + fig.show() .. important:: @@ -228,12 +295,23 @@ And plot it: .. jupyter-execute:: - tmp = upward_continued.plot(cmap="seismic", center=0, add_colorbar=False) - 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() + fig = pygmt.Figure() + + maxabs = vd.maxabs(upward_continued) + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + + fig.grdimage( + upward_continued, + projection="X15c", + cmap=True, + frame=["af", "WeSn+tUpward continued magnetic anomaly at 1000m"] + ) + fig.colorbar( + position="JCB", + frame=["af", 'x+lnT'], + ) + + fig.show() Reduction to the pole @@ -280,12 +358,23 @@ And plot it: .. jupyter-execute:: - tmp = rtp_grid.plot(cmap="seismic", center=0, add_colorbar=False) - plt.gca().set_aspect("equal") - plt.title("Magnetic anomaly reduced to the pole") - plt.gca().ticklabel_format(style="sci", scilimits=(0, 0)) - plt.colorbar(tmp, label="nT") - plt.show() + fig = pygmt.Figure() + + maxabs = vd.maxabs(rtp_grid) * .8 + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + + fig.grdimage( + rtp_grid, + projection="X15c", + cmap=True, + frame=["af", "WeSn+tReduced to the pole magnetic anomaly"] + ) + fig.colorbar( + position="JCB+e", + frame=["af", 'x+lnT'], + ) + + fig.show() If on the other hand we have any knowledge about the orientation of the magnetization vector of the sources, we can specify the @@ -309,12 +398,23 @@ magnetization vector of the sources, we can specify the .. jupyter-execute:: - tmp = rtp_grid.plot(cmap="seismic", center=0, add_colorbar=False) - 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() + fig = pygmt.Figure() + + maxabs = vd.maxabs(rtp_grid) * .8 + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + + fig.grdimage( + rtp_grid, + projection="X15c", + cmap=True, + frame=["af", "WeSn+tReduced to the pole with remanence"] + ) + fig.colorbar( + position="JCB+e", + frame=["af", 'x+lnT'], + ) + + fig.show() Gaussian filters @@ -367,34 +467,35 @@ Let's plot the results side by side: .. jupyter-execute:: - import verde as vd + fig = pygmt.Figure() - fig, (ax1, ax2) = plt.subplots( - nrows=1, ncols=2, sharey=True, figsize=(12, 8) + maxabs = vd.maxabs(magnetic_low_freqs, magnetic_high_freqs) * .6 + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + + fig.grdimage( + magnetic_low_freqs, + projection="X15c", + cmap=True, + frame=["af", "WeSn+tLow-pass filtered magnetic anomaly"] ) - maxabs = vd.maxabs(magnetic_low_freqs, magnetic_high_freqs) - kwargs = dict(cmap="seismic", vmin=-maxabs, vmax=maxabs, add_colorbar=False) - - tmp = magnetic_low_freqs.plot(ax=ax1, **kwargs) - tmp = magnetic_high_freqs.plot(ax=ax2, **kwargs) - - ax1.set_title("Magnetic anomaly after low-pass filter") - ax2.set_title("Magnetic anomaly after high-pass filter") - for ax in (ax1, ax2): - ax.set_aspect("equal") - ax.ticklabel_format(style="sci", scilimits=(0, 0)) - - plt.colorbar( - tmp, - ax=[ax1, ax2], - label="nT", - orientation="horizontal", - aspect=42, - shrink=0.8, - pad=0.08, + fig.shift_origin(xshift="16c") + + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + + fig.grdimage( + magnetic_high_freqs, + projection="X15c", + cmap=True, + frame=["af", "wESn+tHigh-pass filtered magnetic anomaly"] ) - plt.show() + + fig.colorbar( + position="JBC+h+e+o-8c/1c+w15c/.8c", + frame=["af", 'x+lnT'], + ) + + fig.show() Total gradient amplitude @@ -433,14 +534,23 @@ And plot it: .. jupyter-execute:: - import verde as vd + fig = pygmt.Figure() + + maxabs = vd.maxabs(tga_grid) + pygmt.makecpt(cmap="viridis", series=[0, maxabs], background=True) + + fig.grdimage( + tga_grid, + projection="X15c", + cmap=True, + frame=["af", "WeSn+tTotal gradient amplitude of the magnetic anomaly"] + ) + fig.colorbar( + position="JCB", + frame=["af", 'x+lnT/m'], + ) - tmp = tga_grid.plot(cmap="viridis", add_colorbar=False) - 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() + fig.show() ---- From bf00da8eb1bebbbe90f12f61d4fc577c444497d6 Mon Sep 17 00:00:00 2001 From: mdtanker Date: Wed, 19 Nov 2025 09:55:31 +0100 Subject: [PATCH 03/19] wrap returns in print in .rst file --- .../equivalent_sources/eqs-parameters-estimation.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst b/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst index 615da657f..e7f9fcb3d 100644 --- a/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst +++ b/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst @@ -39,7 +39,7 @@ Lets start by loading some sample gravity data: fname = ensaio.fetch_bushveld_gravity(version=1) data = pd.read_csv(fname) - data + print(data) And project their coordinates using a Mercator projection: @@ -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: From 38afc87a8f0b068c8f45cad7b8bb348cf6094913 Mon Sep 17 00:00:00 2001 From: mdtanker Date: Wed, 19 Nov 2025 09:56:00 +0100 Subject: [PATCH 04/19] add hidden pygmt import to igrf.rst --- doc/user_guide/igrf.rst | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/doc/user_guide/igrf.rst b/doc/user_guide/igrf.rst index eef99ee2f..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 From e6d05c8020693617f1d5c43b231249625ff36c03 Mon Sep 17 00:00:00 2001 From: mdtanker Date: Wed, 19 Nov 2025 10:00:07 +0100 Subject: [PATCH 05/19] remove quotes from and update and standardize figure titles and labels --- doc/gallery_src/equivalent_sources/spherical.py | 8 +++++--- doc/gallery_src/forward/prism_layer.py | 4 ++-- doc/gallery_src/forward/tesseroid.py | 2 +- doc/gallery_src/forward/tesseroid_layer.py | 6 +++++- .../forward/tesseroid_variable_density.py | 2 +- doc/gallery_src/transformations/reduction_to_pole.py | 12 +++++++++--- doc/gallery_src/transformations/tga.py | 12 +++++++++--- doc/gallery_src/transformations/tilt.py | 6 ++---- .../transformations/upward_continuation.py | 12 +++++++++--- doc/gallery_src/transformations/upward_derivative.py | 12 +++++++++--- doc/user_guide/forward_modelling/point.rst | 3 +-- doc/user_guide/gravity_disturbance.rst | 6 +++--- doc/user_guide/topographic_correction.rst | 2 +- 13 files changed, 57 insertions(+), 30 deletions(-) diff --git a/doc/gallery_src/equivalent_sources/spherical.py b/doc/gallery_src/equivalent_sources/spherical.py index 255329625..e940f271b 100644 --- a/doc/gallery_src/equivalent_sources/spherical.py +++ b/doc/gallery_src/equivalent_sources/spherical.py @@ -94,9 +94,9 @@ 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, @@ -108,8 +108,10 @@ 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, ) diff --git a/doc/gallery_src/forward/prism_layer.py b/doc/gallery_src/forward/prism_layer.py index 4cf887b73..22a39ba9c 100644 --- a/doc/gallery_src/forward/prism_layer.py +++ b/doc/gallery_src/forward/prism_layer.py @@ -51,14 +51,14 @@ # 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", ) diff --git a/doc/gallery_src/forward/tesseroid.py b/doc/gallery_src/forward/tesseroid.py index 9e805598b..65f240740 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( diff --git a/doc/gallery_src/forward/tesseroid_layer.py b/doc/gallery_src/forward/tesseroid_layer.py index 5e42984c6..a760cbcb4 100644 --- a/doc/gallery_src/forward/tesseroid_layer.py +++ b/doc/gallery_src/forward/tesseroid_layer.py @@ -61,14 +61,18 @@ # Make colormap of data pygmt.makecpt(cmap="balance+h0", series=[-cpt_lims, cpt_lims]) + +title = "Gravitational acceleration of topography with tesseroids" + fig.grdimage( gravity.g_z, + frame=f"+t{title}", projection="M15c", nan_transparent=True, cmap=True, ) fig.basemap(frame=True) -fig.colorbar(frame='af+l"Gravity [mGal]"', position="JCR") +fig.colorbar(frame="af+lmGal", position="JCR") 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..7225c1e64 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( diff --git a/doc/gallery_src/transformations/reduction_to_pole.py b/doc/gallery_src/transformations/reduction_to_pole.py index d0e5944d2..8666b984f 100644 --- a/doc/gallery_src/transformations/reduction_to_pole.py +++ b/doc/gallery_src/transformations/reduction_to_pole.py @@ -61,18 +61,24 @@ 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+ef", ) with fig.set_panel(panel=1): # Plot upward reduced to the pole grid - fig.grdimage(grid=rtp_grid, projection="X?", cmap=True) + fig.grdimage( + grid=rtp_grid, + projection="X?", + cmap=True, + frame="+tReduced to the pole", + ) # Add colorbar fig.colorbar( - frame='af+l"Reduced to the pole [nT]"', + frame="af+lnT", position="JBC+h+o0/1c+ef", ) fig.show() diff --git a/doc/gallery_src/transformations/tga.py b/doc/gallery_src/transformations/tga.py index d5fe8396e..049eae1df 100644 --- a/doc/gallery_src/transformations/tga.py +++ b/doc/gallery_src/transformations/tga.py @@ -52,10 +52,11 @@ 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", ) with fig.set_panel(panel=1): @@ -63,10 +64,15 @@ cpt_lim = 0.6 * vd.maxabs(tga) pygmt.makecpt(cmap="balance+h0", series=[0, cpt_lim], 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]"', + 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 e62234dd7..96dad14df 100644 --- a/doc/gallery_src/transformations/tilt.py +++ b/doc/gallery_src/transformations/tilt.py @@ -101,9 +101,8 @@ 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+ef", ) @@ -129,9 +128,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", ) fig.show() diff --git a/doc/gallery_src/transformations/upward_continuation.py b/doc/gallery_src/transformations/upward_continuation.py index edd76893d..84ca04f2c 100644 --- a/doc/gallery_src/transformations/upward_continuation.py +++ b/doc/gallery_src/transformations/upward_continuation.py @@ -54,18 +54,24 @@ grid=magnetic_grid, projection="X?", cmap=True, + frame="+tMagnetic Anomaly at 500m", ) # Add colorbar fig.colorbar( - frame='af+l"Magnetic anomaly at 500m [nT]"', + frame="af+lnT", position="JBC+h+o0/1c+e", ) with fig.set_panel(panel=1): # Plot upward continued grid - fig.grdimage(grid=upward_continued, projection="X?", cmap=True) + fig.grdimage( + grid=upward_continued, + projection="X?", + frame="+tUpward continued to 1000m", + cmap=True, + ) # Add colorbar fig.colorbar( - frame='af+l"Upward continued to 1000m [nT]"', + frame="af+lnT", position="JBC+h+o0/1c", ) fig.show() diff --git a/doc/gallery_src/transformations/upward_derivative.py b/doc/gallery_src/transformations/upward_derivative.py index 30c74618a..a77856dfb 100644 --- a/doc/gallery_src/transformations/upward_derivative.py +++ b/doc/gallery_src/transformations/upward_derivative.py @@ -53,10 +53,11 @@ 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): @@ -64,10 +65,15 @@ cpt_lim = vd.maxabs(deriv_upward) * 0.6 pygmt.makecpt(cmap="balance+h0", series=[-cpt_lim, cpt_lim], 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/user_guide/forward_modelling/point.rst b/doc/user_guide/forward_modelling/point.rst index 8f6999474..cf03b5178 100644 --- a/doc/user_guide/forward_modelling/point.rst +++ b/doc/user_guide/forward_modelling/point.rst @@ -238,7 +238,6 @@ Lets plot these results using :mod:`pygmt`: coordinates_spherical, g_z, data_names="g_z", extra_coords_names="extra") fig = pygmt.Figure() - title = "Gravitational acceleration (downward)" maxabs = vd.maxabs(g_z) pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) @@ -246,7 +245,7 @@ Lets plot these results using :mod:`pygmt`: region=(-72, -68, -46, -42), projection="M10c", grid=grid.g_z, - frame=[f"WSne+t{title}", "x", "y"], + frame=["a", "x", "y"], cmap=True,) fig.colorbar(cmap=True, position="JMR", frame=["a0.000000005", "x+lmGal"]) diff --git a/doc/user_guide/gravity_disturbance.rst b/doc/user_guide/gravity_disturbance.rst index 6aef78354..a8130a56d 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() @@ -153,7 +153,7 @@ And plot it: fig.basemap(frame=["af", "WEsn"]) fig.colorbar( position="JCB+w10c", - frame=["af", 'y+l"mGal"', 'x+l"gravity disturbance"'], + frame=["af", 'y+lmGal', 'x+lGravity disturbance'], ) fig.coast(shorelines=True, resolution="c", area_thresh=1e4) fig.show() diff --git a/doc/user_guide/topographic_correction.rst b/doc/user_guide/topographic_correction.rst index 0e1fb0080..7e980e9db 100644 --- a/doc/user_guide/topographic_correction.rst +++ b/doc/user_guide/topographic_correction.rst @@ -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() From 725169b57c3693b52164f5ace669c849d126a465 Mon Sep 17 00:00:00 2001 From: mdtanker Date: Wed, 19 Nov 2025 10:00:47 +0100 Subject: [PATCH 06/19] standardize forward calculation example colorbars --- doc/gallery_src/forward/tesseroid.py | 2 +- doc/gallery_src/forward/tesseroid_variable_density.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/gallery_src/forward/tesseroid.py b/doc/gallery_src/forward/tesseroid.py index 65f240740..d18a50fc7 100644 --- a/doc/gallery_src/forward/tesseroid.py +++ b/doc/gallery_src/forward/tesseroid.py @@ -61,7 +61,7 @@ cmap="viridis", ) -fig.colorbar(cmap=True, frame=["a200f50", "x+lmGal"]) +fig.colorbar(cmap=True, position="JMR", frame=["a200f50", "x+lmGal"]) fig.coast(shorelines="1p,black") diff --git a/doc/gallery_src/forward/tesseroid_variable_density.py b/doc/gallery_src/forward/tesseroid_variable_density.py index 7225c1e64..3142e3f0f 100644 --- a/doc/gallery_src/forward/tesseroid_variable_density.py +++ b/doc/gallery_src/forward/tesseroid_variable_density.py @@ -89,7 +89,7 @@ def density(radius): cmap="viridis", ) -fig.colorbar(cmap=True, frame=["a200f50", "x+lmGal"]) +fig.colorbar(cmap=True, position="JMR", frame=["a200f50", "x+lmGal"]) fig.coast(shorelines="1p,black") From 1694600215117b946d341566586c7d4b4349706d Mon Sep 17 00:00:00 2001 From: mdtanker Date: Wed, 19 Nov 2025 10:01:03 +0100 Subject: [PATCH 07/19] add labels to point masses in plot --- doc/gallery_src/forward/point_gravity.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/doc/gallery_src/forward/point_gravity.py b/doc/gallery_src/forward/point_gravity.py index b822188b5..9c5be7555 100644 --- a/doc/gallery_src/forward/point_gravity.py +++ b/doc/gallery_src/forward/point_gravity.py @@ -68,7 +68,14 @@ cmap=True, ) -fig.plot(x=easting, y=northing, style="c0.2c", fill="grey") +fig.plot( + x=easting, + y=northing, + style="c0.15c", + fill="white", + pen="1p,black", + label="Point masses", +) fig.colorbar(cmap=True, position="JMR+e", frame=["x+lmGal"]) From 0bfeb31912efa183e23de99270f01607b0ba3f0f Mon Sep 17 00:00:00 2001 From: mdtanker Date: Wed, 19 Nov 2025 10:01:18 +0100 Subject: [PATCH 08/19] make nans transparent in plot --- doc/gallery_src/equivalent_sources/spherical.py | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/gallery_src/equivalent_sources/spherical.py b/doc/gallery_src/equivalent_sources/spherical.py index e940f271b..81551a8fe 100644 --- a/doc/gallery_src/equivalent_sources/spherical.py +++ b/doc/gallery_src/equivalent_sources/spherical.py @@ -114,6 +114,7 @@ frame=[f"ESnw+t{title}", "xa5", "ya4"], grid=grid.gravity_disturbance, cmap=True, + nan_transparent=True, ) fig.colorbar(cmap=True, frame=["x+lmGal"], position="+e") From 1cf65809e43bdce0352eef4c5651f9fd0060a614 Mon Sep 17 00:00:00 2001 From: mdtanker Date: Wed, 19 Nov 2025 14:17:10 +0100 Subject: [PATCH 09/19] format fix --- doc/gallery_src/transformations/reduction_to_pole.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/doc/gallery_src/transformations/reduction_to_pole.py b/doc/gallery_src/transformations/reduction_to_pole.py index 8666b984f..0741eb28c 100644 --- a/doc/gallery_src/transformations/reduction_to_pole.py +++ b/doc/gallery_src/transformations/reduction_to_pole.py @@ -58,10 +58,7 @@ with fig.set_panel(panel=0): # Plot magnetic anomaly grid fig.grdimage( - grid=magnetic_grid, - projection="X?", - cmap=True, - frame="+tMagnetic anomaly" + grid=magnetic_grid, projection="X?", cmap=True, frame="+tMagnetic anomaly" ) # Add colorbar fig.colorbar( From b839ef1fa87f88bd05f02cd319e4e45a138e5505 Mon Sep 17 00:00:00 2001 From: mdtanker Date: Fri, 28 Nov 2025 12:32:05 +0100 Subject: [PATCH 10/19] remove print statement to keep html rendering of dataframe --- doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst b/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst index e7f9fcb3d..2053c2dc5 100644 --- a/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst +++ b/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst @@ -39,7 +39,7 @@ Lets start by loading some sample gravity data: fname = ensaio.fetch_bushveld_gravity(version=1) data = pd.read_csv(fname) - print(data) + data And project their coordinates using a Mercator projection: From db2ad9e407e45bc693242605dd4f686a7ea24c3e Mon Sep 17 00:00:00 2001 From: mdtanker Date: Thu, 30 Apr 2026 07:39:25 +0200 Subject: [PATCH 11/19] use single colorbar in eqs gallery --- .../block_averaged_sources.py | 21 +++++++++--------- .../equivalent_sources/cartesian.py | 22 +++++++++---------- .../equivalent_sources/gradient_boosted.py | 22 +++++++++---------- .../equivalent_sources/spherical.py | 22 +++++++++---------- 4 files changed, 43 insertions(+), 44 deletions(-) diff --git a/doc/gallery_src/equivalent_sources/block_averaged_sources.py b/doc/gallery_src/equivalent_sources/block_averaged_sources.py index 150630a39..e0c8531aa 100644 --- a/doc/gallery_src/equivalent_sources/block_averaged_sources.py +++ b/doc/gallery_src/equivalent_sources/block_averaged_sources.py @@ -107,16 +107,14 @@ title = "Observed magnetic anomaly data" -# Get the 99.9th percentile of the absolute value of the point data to use as color -# scale limits -cpt_lim = np.quantile(np.abs(data.total_field_anomaly_nt), 0.999) +# 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=[-cpt_lim, cpt_lim], - background=True, -) +pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) with pygmt.config(FONT_TITLE="12p"): fig.plot( @@ -130,7 +128,6 @@ cmap=True, ) -fig.colorbar(cmap=True, frame=["x+lnT"], position="+e") fig.shift_origin(xshift=fig_width + 1) title = "Gridded and upward-continued" @@ -142,6 +139,10 @@ cmap=True, ) -fig.colorbar(cmap=True, frame=["x+lnT"], position="+e") +fig.colorbar( + cmap=True, + frame=["a400f100","x+lnT"], + position=f"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 c789ffcc1..37611c795 100644 --- a/doc/gallery_src/equivalent_sources/cartesian.py +++ b/doc/gallery_src/equivalent_sources/cartesian.py @@ -96,16 +96,14 @@ title = "Observed magnetic anomaly data" -# Get the 99.9th percentile of the absolute value of the point data to use as color -# scale limits -cpt_lim = np.quantile(np.abs(data.total_field_anomaly_nt), 0.999) +# 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=[-cpt_lim, cpt_lim], - background=True, -) +pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) with pygmt.config(FONT_TITLE="12p"): fig.plot( @@ -119,8 +117,6 @@ cmap=True, ) -fig.colorbar(cmap=True, frame=["x+lnT"], position="+e") - fig.shift_origin(xshift=fig_width + 1) title = "Gridded and upward-continued" @@ -132,6 +128,10 @@ cmap=True, ) -fig.colorbar(cmap=True, frame=["x+lnT"], position="+e") +fig.colorbar( + cmap=True, + frame=["a400f100","x+lnT"], + position=f"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 6f723e1fe..cfe022b90 100644 --- a/doc/gallery_src/equivalent_sources/gradient_boosted.py +++ b/doc/gallery_src/equivalent_sources/gradient_boosted.py @@ -116,16 +116,14 @@ title = "Observed gravity disturbance data" -# Get the 99.9th percentile of the absolute value of the point data to use as color -# scale limits -cpt_lim = np.quantile(np.abs(data.gravity_disturbance), 0.999) +# 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=[-cpt_lim, cpt_lim], - background=True, -) +pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) with pygmt.config(FONT_TITLE="14p"): fig.plot( @@ -139,8 +137,6 @@ cmap=True, ) -fig.colorbar(cmap=True, frame=["x+lmGal"], position="+e") - fig.shift_origin(xshift=fig_width + 1) title = "Gridded with gradient-boosted equivalent sources" @@ -152,6 +148,10 @@ cmap=True, ) -fig.colorbar(cmap=True, frame=["x+lmGal"], position="+e") +fig.colorbar( + cmap=True, + frame=["x+lmGal"], + position=f"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 4fadeea87..84aa1f7ef 100644 --- a/doc/gallery_src/equivalent_sources/spherical.py +++ b/doc/gallery_src/equivalent_sources/spherical.py @@ -80,16 +80,12 @@ # Plot observed and gridded gravity disturbance fig = pygmt.Figure() -# Get the 99.9th percentile of the absolute value of the point data to use as color -# scale limits -cpt_lim = np.quantile(np.abs(gravity_disturbance), 0.999) +# 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 -pygmt.makecpt( - cmap="balance+h0", - series=[-cpt_lim, cpt_lim], - background=True, -) +pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) + title = "Observed gravity disturbance data" @@ -104,9 +100,7 @@ cmap=True, ) -fig.colorbar(cmap=True, frame=["x+lmGal"], position="+e") - -fig.shift_origin(xshift="w+3c") +fig.shift_origin(xshift="w+1c") title = "Gridded and upward-continued" @@ -117,6 +111,10 @@ nan_transparent=True, ) -fig.colorbar(cmap=True, frame=["x+lmGal"], position="+e") +fig.colorbar( + cmap=True, + frame=["x+lmGal"], + position=f"n0/0+jTC+w10c/0.5c+h+o-0.5c/0.9c", +) fig.show() From 78f89766c941afce962e10281e1be3a77916184c Mon Sep 17 00:00:00 2001 From: mdtanker Date: Thu, 30 Apr 2026 07:39:43 +0200 Subject: [PATCH 12/19] remove dublicate boule from environment.yml --- environment.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/environment.yml b/environment.yml index 956cf3305..2b20636e5 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.* From 4811cc10cf8fc99945f8a64094dfbaa5536ff947 Mon Sep 17 00:00:00 2001 From: mdtanker Date: Thu, 30 Apr 2026 07:40:02 +0200 Subject: [PATCH 13/19] update listed dependencies in install docs --- doc/install.rst | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) 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 `__ * `pandas `__ -* `numba `__ * `scipy `__ -* `xarray `__ * `scikit-learn `__ -* `pooch `__ +* `numba `__ +* `xarray `__ * `verde `__ * `xrft `__ +* `choclo `__ +* `boule `__ Optional: @@ -92,7 +93,7 @@ Optional: The examples in the :ref:`gallery` also use: -* `boule `__ * `ensaio `__ for downloading sample datasets -* `pygmt `__ for plotting maps +* `pygmt `__ and `matplotlib `__ for plotting maps and figures * `pyproj `__ for cartographic projections +* `bordado `__ for working with geographic coordinates \ No newline at end of file From bbf586dcd29730e7cce5c8270f3e716ab20959df Mon Sep 17 00:00:00 2001 From: mdtanker Date: Thu, 30 Apr 2026 07:51:30 +0200 Subject: [PATCH 14/19] plot tesseroid boundaries in forward modeling gallery --- doc/gallery_src/forward/tesseroid.py | 10 ++++++++++ .../forward/tesseroid_variable_density.py | 14 ++++++++++++++ 2 files changed, 24 insertions(+) diff --git a/doc/gallery_src/forward/tesseroid.py b/doc/gallery_src/forward/tesseroid.py index d18a50fc7..0f1c8a01e 100644 --- a/doc/gallery_src/forward/tesseroid.py +++ b/doc/gallery_src/forward/tesseroid.py @@ -61,6 +61,16 @@ cmap="viridis", ) +# 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_variable_density.py b/doc/gallery_src/forward/tesseroid_variable_density.py index 3142e3f0f..20e12e563 100644 --- a/doc/gallery_src/forward/tesseroid_variable_density.py +++ b/doc/gallery_src/forward/tesseroid_variable_density.py @@ -89,6 +89,20 @@ def density(radius): cmap="viridis", ) +# Plot edges of tesseroids +for i, tesseroid in enumerate(tesseroids): + if i == 0: + label="Tesseroid boundaries" + else: + label=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") From d13339b798086abde9aa0dc1026bd60010c95453 Mon Sep 17 00:00:00 2001 From: mdtanker Date: Thu, 30 Apr 2026 08:12:14 +0200 Subject: [PATCH 15/19] use verde percentiles in other gallery docs --- doc/gallery_src/forward/point_gravity.py | 7 +++---- doc/gallery_src/gravity_disturbance.py | 5 +++-- doc/gallery_src/gravity_disturbance_topofree.py | 5 +++-- doc/gallery_src/isostatic_moho_airy.py | 6 +++++- .../transformations/reduction_to_pole.py | 2 +- doc/gallery_src/transformations/tga.py | 10 +++++----- doc/gallery_src/transformations/tilt.py | 16 ++++++++-------- .../transformations/upward_continuation.py | 2 +- .../transformations/upward_derivative.py | 10 +++++----- 9 files changed, 34 insertions(+), 29 deletions(-) diff --git a/doc/gallery_src/forward/point_gravity.py b/doc/gallery_src/forward/point_gravity.py index 9c5be7555..c6b8ff385 100644 --- a/doc/gallery_src/forward/point_gravity.py +++ b/doc/gallery_src/forward/point_gravity.py @@ -54,7 +54,7 @@ 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], background=True) @@ -71,9 +71,8 @@ fig.plot( x=easting, y=northing, - style="c0.15c", - fill="white", - pen="1p,black", + style="t0.15c", + fill="magenta", label="Point masses", ) diff --git a/doc/gallery_src/gravity_disturbance.py b/doc/gallery_src/gravity_disturbance.py index 45b28bdec..787a0fedb 100644 --- a/doc/gallery_src/gravity_disturbance.py +++ b/doc/gallery_src/gravity_disturbance.py @@ -21,6 +21,7 @@ import numpy as np import pygmt import xarray as xr +import verde as vd # Load the global gravity grid fname = ensaio.fetch_earth_gravity(version=1) @@ -38,10 +39,10 @@ fig = pygmt.Figure() # Get the 99.9th percentile of the absolute value to use as color scale limits -cpt_lims = np.quantile(np.abs(disturbance), 0.999) +maxabs = vd.maxabs(disturbance, percentile=99.9) # Make colormap of data -pygmt.makecpt(cmap="balance+h0", series=[-cpt_lims, cpt_lims], background=True) +pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) title = "Gravity disturbance of the Earth" diff --git a/doc/gallery_src/gravity_disturbance_topofree.py b/doc/gallery_src/gravity_disturbance_topofree.py index d6d9add3f..d9bbae76b 100644 --- a/doc/gallery_src/gravity_disturbance_topofree.py +++ b/doc/gallery_src/gravity_disturbance_topofree.py @@ -31,6 +31,7 @@ import numpy as np import pygmt import xarray as xr +import verde as vd import harmonica as hm @@ -64,10 +65,10 @@ fig = pygmt.Figure() # Get the 99th percentile of the absolute value to use as color scale limits -cpt_lims = np.quantile(np.abs(disturbance_topofree), 0.99) +maxabs = vd.maxabs(disturbance_topofree, percentile=99) # Make colormap of data -pygmt.makecpt(cmap="balance+h0", series=[-cpt_lims, cpt_lims], background=True) +pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) title = "Topography-free (Bouguer) gravity disturbance of the Earth" diff --git a/doc/gallery_src/isostatic_moho_airy.py b/doc/gallery_src/isostatic_moho_airy.py index 7fd0b346e..4a1a9c6fe 100644 --- a/doc/gallery_src/isostatic_moho_airy.py +++ b/doc/gallery_src/isostatic_moho_airy.py @@ -29,6 +29,7 @@ import numpy as np import pygmt import xarray as xr +import verde as vd import harmonica as hm @@ -59,8 +60,11 @@ # Make a plot of data using PyGMT fig = pygmt.Figure() +# 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=[moho.to_numpy().min(), moho.to_numpy().max()]) +pygmt.makecpt(cmap="viridis", series=cpt_lims) title = "Airy isostatic Moho depth of Africa" diff --git a/doc/gallery_src/transformations/reduction_to_pole.py b/doc/gallery_src/transformations/reduction_to_pole.py index 176b8c279..babe1bb0e 100644 --- a/doc/gallery_src/transformations/reduction_to_pole.py +++ b/doc/gallery_src/transformations/reduction_to_pole.py @@ -42,7 +42,7 @@ fig = pygmt.Figure() with fig.subplot(nrows=1, ncols=2, figsize=("28c", "15c"), sharey="l"): # Make colormap for both plots - cpt_lim = 0.5 * vd.maxabs(magnetic_grid, rtp_grid) + 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 diff --git a/doc/gallery_src/transformations/tga.py b/doc/gallery_src/transformations/tga.py index e83ac6937..3619b2cae 100644 --- a/doc/gallery_src/transformations/tga.py +++ b/doc/gallery_src/transformations/tga.py @@ -31,8 +31,8 @@ with fig.subplot(nrows=1, ncols=2, figsize=("28c", "15c"), sharey="l"): with fig.set_panel(panel=0): # Make colormap of data - cpt_lim = vd.maxabs(magnetic_grid) - pygmt.makecpt(cmap="balance+h0", series=[-cpt_lim, cpt_lim], 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, @@ -43,12 +43,12 @@ # Add colorbar fig.colorbar( frame="af+lnT", - position="JBC+h+o0/1c", + position="JBC+h+o0/1c+e", ) with fig.set_panel(panel=1): # Make colormap for total gradient amplitude (saturate it a little bit) - cpt_lim = 0.6 * vd.maxabs(tga) - pygmt.makecpt(cmap="balance+h0", series=[0, cpt_lim], 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, diff --git a/doc/gallery_src/transformations/tilt.py b/doc/gallery_src/transformations/tilt.py index 9be9c723e..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"], ): - cpt_lim = 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="balance+h0", series=[-cpt_lim, cpt_lim], 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="balance+h0", series=[-cpt_lim, cpt_lim], 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, @@ -86,13 +86,13 @@ # Add colorbar fig.colorbar( frame="af+lnT", - position="JMR+o1/-0.25c+ef", + position="JMR+o1/-0.25c+e", ) - cpt_lim = 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="balance+h0", series=[-cpt_lim, cpt_lim], background=True) + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) # Plot tilt fig.grdimage( grid=tilt_grid, @@ -102,7 +102,7 @@ ) with fig.set_panel(panel=3): # Make colormap for tilt rtp (saturate it a little bit) - pygmt.makecpt(cmap="balance+h0", series=[-cpt_lim, cpt_lim], background=True) + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) # Plot tilt fig.grdimage( grid=tilt_rtp_grid, @@ -113,6 +113,6 @@ # Add colorbar fig.colorbar( frame="af+lradians", - position="JMR+o1/-0.25c", + 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 5b8c4c161..dd6daa10e 100644 --- a/doc/gallery_src/transformations/upward_continuation.py +++ b/doc/gallery_src/transformations/upward_continuation.py @@ -32,7 +32,7 @@ fig = pygmt.Figure() with fig.subplot(nrows=1, ncols=2, figsize=("28c", "15c"), sharey="l"): # Make colormap based on original data - cpt_lim = vd.maxabs(magnetic_grid) * 0.6 + 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 diff --git a/doc/gallery_src/transformations/upward_derivative.py b/doc/gallery_src/transformations/upward_derivative.py index 341c410cd..38e90f21b 100644 --- a/doc/gallery_src/transformations/upward_derivative.py +++ b/doc/gallery_src/transformations/upward_derivative.py @@ -32,8 +32,8 @@ with fig.subplot(nrows=1, ncols=2, figsize=("28c", "15c"), sharey="l"): with fig.set_panel(panel=0): # Make colormap of data - cpt_lim = vd.maxabs(magnetic_grid) * 0.6 - pygmt.makecpt(cmap="balance+h0", series=[-cpt_lim, cpt_lim], 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, @@ -47,9 +47,9 @@ position="JBC+h+o0/1c+e", ) with fig.set_panel(panel=1): - # Make colormap for upward derivative (saturate it a little bit) - cpt_lim = vd.maxabs(deriv_upward) * 0.6 - pygmt.makecpt(cmap="balance+h0", series=[-cpt_lim, cpt_lim], 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, From 6f4567a557b5e36711a8561b0e543bd5d0dbfd1a Mon Sep 17 00:00:00 2001 From: mdtanker Date: Thu, 30 Apr 2026 08:12:39 +0200 Subject: [PATCH 16/19] combine colorbars for compare RTP and up continued grids --- .../transformations/reduction_to_pole.py | 17 +++++++---------- .../transformations/upward_continuation.py | 17 +++++++---------- 2 files changed, 14 insertions(+), 20 deletions(-) diff --git a/doc/gallery_src/transformations/reduction_to_pole.py b/doc/gallery_src/transformations/reduction_to_pole.py index babe1bb0e..5ce685e72 100644 --- a/doc/gallery_src/transformations/reduction_to_pole.py +++ b/doc/gallery_src/transformations/reduction_to_pole.py @@ -49,11 +49,6 @@ fig.grdimage( grid=magnetic_grid, projection="X?", cmap=True, frame="+tMagnetic anomaly" ) - # Add colorbar - fig.colorbar( - frame="af+lnT", - position="JBC+h+o0/1c+ef", - ) with fig.set_panel(panel=1): # Plot upward reduced to the pole grid fig.grdimage( @@ -62,9 +57,11 @@ cmap=True, frame="+tReduced to the pole", ) - # Add colorbar - fig.colorbar( - frame="af+lnT", - position="JBC+h+o0/1c+ef", - ) + # Add colorbar + fig.colorbar( + cmap=True, + frame=["a1000f500", "x+lnT"], + position=f"n0/0+jTC+w12c/0.5c+h+o-0.5c/0.9c+e", + ) + fig.show() diff --git a/doc/gallery_src/transformations/upward_continuation.py b/doc/gallery_src/transformations/upward_continuation.py index dd6daa10e..95b147f63 100644 --- a/doc/gallery_src/transformations/upward_continuation.py +++ b/doc/gallery_src/transformations/upward_continuation.py @@ -42,11 +42,6 @@ cmap=True, frame="+tMagnetic Anomaly at 500m", ) - # Add colorbar - fig.colorbar( - frame="af+lnT", - position="JBC+h+o0/1c+e", - ) with fig.set_panel(panel=1): # Plot upward continued grid fig.grdimage( @@ -55,9 +50,11 @@ frame="+tUpward continued to 1000m", cmap=True, ) - # Add colorbar - fig.colorbar( - frame="af+lnT", - position="JBC+h+o0/1c", - ) + # Add colorbar + fig.colorbar( + cmap=True, + frame=["a1000f500", "x+lnT"], + position=f"n0/0+jTC+w12c/0.5c+h+o-0.5c/0.9c+e", + ) + fig.show() From 1cbf976df09e539d3a9cde1c475a9e8b4e0c33c5 Mon Sep 17 00:00:00 2001 From: mdtanker Date: Thu, 30 Apr 2026 14:01:43 +0200 Subject: [PATCH 17/19] add min version for verde for docs --- environment.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/environment.yml b/environment.yml index 2b20636e5..5a010b430 100644 --- a/environment.yml +++ b/environment.yml @@ -39,6 +39,7 @@ dependencies: - matplotlib - bordado>=0.4.0 - ensaio>=0.5.0 + - verde>=1.9.0 - netcdf4 - pygmt - gmt From bada65a4cb0ce1cadc850fd7842113c408999a65 Mon Sep 17 00:00:00 2001 From: mdtanker Date: Thu, 30 Apr 2026 14:02:36 +0200 Subject: [PATCH 18/19] formatting --- doc/gallery_src/equivalent_sources/block_averaged_sources.py | 5 ++--- doc/gallery_src/equivalent_sources/cartesian.py | 5 ++--- doc/gallery_src/equivalent_sources/gradient_boosted.py | 3 +-- doc/gallery_src/equivalent_sources/spherical.py | 2 +- doc/gallery_src/forward/tesseroid_variable_density.py | 5 +---- doc/gallery_src/gravity_disturbance.py | 3 +-- doc/gallery_src/gravity_disturbance_topofree.py | 3 +-- doc/gallery_src/isostatic_moho_airy.py | 2 +- doc/gallery_src/transformations/reduction_to_pole.py | 2 +- doc/gallery_src/transformations/upward_continuation.py | 2 +- 10 files changed, 12 insertions(+), 20 deletions(-) diff --git a/doc/gallery_src/equivalent_sources/block_averaged_sources.py b/doc/gallery_src/equivalent_sources/block_averaged_sources.py index e0c8531aa..5e4d9a6d1 100644 --- a/doc/gallery_src/equivalent_sources/block_averaged_sources.py +++ b/doc/gallery_src/equivalent_sources/block_averaged_sources.py @@ -41,7 +41,6 @@ """ import ensaio -import numpy as np import pandas as pd import pygmt import pyproj @@ -141,8 +140,8 @@ fig.colorbar( cmap=True, - frame=["a400f100","x+lnT"], - position=f"n0/0+jTC+w8c/0.25c+h+o-0.5c/0.8c", + 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 37611c795..8d71806ca 100644 --- a/doc/gallery_src/equivalent_sources/cartesian.py +++ b/doc/gallery_src/equivalent_sources/cartesian.py @@ -31,7 +31,6 @@ """ import ensaio -import numpy as np import pandas as pd import pygmt import pyproj @@ -130,8 +129,8 @@ fig.colorbar( cmap=True, - frame=["a400f100","x+lnT"], - position=f"n0/0+jTC+w8c/0.25c+h+o-0.5c/0.8c", + 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 cfe022b90..5e3d58b7a 100644 --- a/doc/gallery_src/equivalent_sources/gradient_boosted.py +++ b/doc/gallery_src/equivalent_sources/gradient_boosted.py @@ -31,7 +31,6 @@ import boule as bl import ensaio -import numpy as np import pandas as pd import pygmt import pyproj @@ -151,7 +150,7 @@ fig.colorbar( cmap=True, frame=["x+lmGal"], - position=f"n0/0+jTC+w10c/0.5c+h+o-0.5c/0.9c", + 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 84aa1f7ef..3c3c9cb4b 100644 --- a/doc/gallery_src/equivalent_sources/spherical.py +++ b/doc/gallery_src/equivalent_sources/spherical.py @@ -114,7 +114,7 @@ fig.colorbar( cmap=True, frame=["x+lmGal"], - position=f"n0/0+jTC+w10c/0.5c+h+o-0.5c/0.9c", + position="n0/0+jTC+w10c/0.5c+h+o-0.5c/0.9c", ) fig.show() diff --git a/doc/gallery_src/forward/tesseroid_variable_density.py b/doc/gallery_src/forward/tesseroid_variable_density.py index 20e12e563..8ce08634b 100644 --- a/doc/gallery_src/forward/tesseroid_variable_density.py +++ b/doc/gallery_src/forward/tesseroid_variable_density.py @@ -91,10 +91,7 @@ def density(radius): # Plot edges of tesseroids for i, tesseroid in enumerate(tesseroids): - if i == 0: - label="Tesseroid boundaries" - else: - label=None + 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]], diff --git a/doc/gallery_src/gravity_disturbance.py b/doc/gallery_src/gravity_disturbance.py index 787a0fedb..393deeecb 100644 --- a/doc/gallery_src/gravity_disturbance.py +++ b/doc/gallery_src/gravity_disturbance.py @@ -18,10 +18,9 @@ import boule as bl import ensaio -import numpy as np import pygmt -import xarray as xr import verde as vd +import xarray as xr # Load the global gravity grid fname = ensaio.fetch_earth_gravity(version=1) diff --git a/doc/gallery_src/gravity_disturbance_topofree.py b/doc/gallery_src/gravity_disturbance_topofree.py index d9bbae76b..fe989aa44 100644 --- a/doc/gallery_src/gravity_disturbance_topofree.py +++ b/doc/gallery_src/gravity_disturbance_topofree.py @@ -28,10 +28,9 @@ import boule as bl import ensaio -import numpy as np import pygmt -import xarray as xr import verde as vd +import xarray as xr import harmonica as hm diff --git a/doc/gallery_src/isostatic_moho_airy.py b/doc/gallery_src/isostatic_moho_airy.py index 4a1a9c6fe..c6fb1b29f 100644 --- a/doc/gallery_src/isostatic_moho_airy.py +++ b/doc/gallery_src/isostatic_moho_airy.py @@ -28,8 +28,8 @@ import ensaio import numpy as np import pygmt -import xarray as xr import verde as vd +import xarray as xr import harmonica as hm diff --git a/doc/gallery_src/transformations/reduction_to_pole.py b/doc/gallery_src/transformations/reduction_to_pole.py index 5ce685e72..2088a8cfa 100644 --- a/doc/gallery_src/transformations/reduction_to_pole.py +++ b/doc/gallery_src/transformations/reduction_to_pole.py @@ -61,7 +61,7 @@ fig.colorbar( cmap=True, frame=["a1000f500", "x+lnT"], - position=f"n0/0+jTC+w12c/0.5c+h+o-0.5c/0.9c+e", + position="n0/0+jTC+w12c/0.5c+h+o-0.5c/0.9c+e", ) fig.show() diff --git a/doc/gallery_src/transformations/upward_continuation.py b/doc/gallery_src/transformations/upward_continuation.py index 95b147f63..4058d08bf 100644 --- a/doc/gallery_src/transformations/upward_continuation.py +++ b/doc/gallery_src/transformations/upward_continuation.py @@ -54,7 +54,7 @@ fig.colorbar( cmap=True, frame=["a1000f500", "x+lnT"], - position=f"n0/0+jTC+w12c/0.5c+h+o-0.5c/0.9c+e", + position="n0/0+jTC+w12c/0.5c+h+o-0.5c/0.9c+e", ) fig.show() From 7fb99244aa8ba04bf47fad36d64df2ba89a631eb Mon Sep 17 00:00:00 2001 From: mdtanker Date: Thu, 30 Apr 2026 14:04:10 +0200 Subject: [PATCH 19/19] update colorbars and plots and revert some back to xarray --- .../equivalent_sources/block-averaged-eqs.rst | 18 +- .../eq-sources-spherical.rst | 14 +- .../eqs-parameters-estimation.rst | 9 +- .../gradient-boosted-eqs.rst | 15 +- doc/user_guide/equivalent_sources/index.rst | 9 +- doc/user_guide/forward_modelling/dipole.rst | 68 ++--- .../forward_modelling/ellipsoid.rst | 220 ++++++++------ doc/user_guide/forward_modelling/point.rst | 3 +- doc/user_guide/forward_modelling/prism.rst | 280 +++++------------- .../forward_modelling/tesseroid.rst | 3 +- doc/user_guide/gravity_disturbance.rst | 8 +- doc/user_guide/topographic_correction.rst | 37 ++- doc/user_guide/transformations.rst | 251 ++++++++-------- 13 files changed, 428 insertions(+), 507 deletions(-) diff --git a/doc/user_guide/equivalent_sources/block-averaged-eqs.rst b/doc/user_guide/equivalent_sources/block-averaged-eqs.rst index 77c5f28f4..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 @@ -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. @@ -184,7 +184,8 @@ we are efectivelly upward continuing the data. pygmt.makecpt( 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 43d05685c..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="balance+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 2053c2dc5..bc361a302 100644 --- a/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst +++ b/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst @@ -216,7 +216,7 @@ Lets plot it: fig = pygmt.Figure() # Make colormap of data - pygmt.makecpt(cmap="balance+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 94dd03789..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="balance+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 bbdec13d6..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="balance+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() @@ -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 3d42046ff..5555d8539 100644 --- a/doc/user_guide/forward_modelling/dipole.rst +++ b/doc/user_guide/forward_modelling/dipole.rst @@ -104,58 +104,24 @@ We can reshape the results into variables of an dataset: ) .. 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 - - fig = pygmt.Figure() - - maxabs = vd.maxabs(grid.b_e) - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - fig.grdimage( - grid=grid.b_e, - projection="X10c", - cmap=True, - frame=["WSne+tEasting component", "xa","ya"], - ) - fig.colorbar(frame='+lnT') - - fig.shift_origin(xshift="11c") - - maxabs = vd.maxabs(grid.b_n) - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - fig.grdimage( - grid=grid.b_n, - projection="X10c", - cmap=True, - frame=["wSne+tNorthing component", "xa","ya"], - ) - fig.colorbar(frame='+lnT') - - fig.shift_origin(xshift="11c") - - maxabs = vd.maxabs(grid.b_u) - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - fig.grdimage( - grid=grid.b_u, - projection="X10c", - cmap=True, - frame=["wSnE+tUpward component", "xa","ya"], - ) - fig.colorbar(frame='+lnT') - fig.show() + import matplotlib.pyplot as plt + + 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', + ) + 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() ---- 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 108f3e7cf..66bcd3637 100644 --- a/doc/user_guide/forward_modelling/point.rst +++ b/doc/user_guide/forward_modelling/point.rst @@ -238,6 +238,7 @@ Lets plot these results using :mod:`pygmt`: coordinates_spherical, g_z, data_names="g_z", extra_coords_names="extra") fig = pygmt.Figure() + title = "Gravitational acceleration (downward)" maxabs = vd.maxabs(g_z) pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) @@ -245,7 +246,7 @@ Lets plot these results using :mod:`pygmt`: region=(-72, -68, -46, -42), projection="M10c", grid=grid.g_z, - frame=["a", "x", "y"], + frame=[f"WSne+t{title}", "x", "y"], cmap=True,) fig.colorbar(cmap=True, position="JMR", frame=["a0.000000005", "x+lmGal"]) diff --git a/doc/user_guide/forward_modelling/prism.rst b/doc/user_guide/forward_modelling/prism.rst index 8ae048179..8a7d8e8eb 100644 --- a/doc/user_guide/forward_modelling/prism.rst +++ b/doc/user_guide/forward_modelling/prism.rst @@ -105,152 +105,71 @@ We can reshape the results into variables of an dataset: Plot the results: .. 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 + import matplotlib.pyplot as plt - fig = pygmt.Figure() - - fig.grdimage( - projection="X10c", - grid=grid.potential, - frame=["a", "x+leasting (m)", "y+lnorthing (m)"], - cmap="viridis", - ) - - fig.colorbar(cmap=True, position="JMR", frame=["x+lJ kg@+-1@+"]) - fig.show() + grid.potential.plot(cbar_kwargs={"label":"J/kg"}) + plt.gca().set_aspect("equal") + plt.show() .. jupyter-execute:: - fig = pygmt.Figure() - - maxabs = vd.maxabs(grid.g_e) - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - fig.grdimage( - grid=grid.g_e, - projection="X10c", - cmap=True, - frame=["WSne+tEasting component", "xa","ya"], - ) - fig.colorbar(frame='+lnT') - - fig.shift_origin(xshift="11c") + fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8)) - maxabs = vd.maxabs(grid.g_n) - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - fig.grdimage( - grid=grid.g_n, - projection="X10c", - cmap=True, - frame=["wSne+tNorthing component", "xa","ya"], - ) - fig.colorbar(frame='+lnT') - - fig.shift_origin(xshift="11c") - - maxabs = vd.maxabs(grid.g_z) - pygmt.makecpt(cmap="balance+h0", series=[0, maxabs], background=True) - fig.grdimage( - grid=grid.g_z, - projection="X10c", - cmap=True, - frame=["wSnE+tDownward component", "xa","ya"], - ) - fig.colorbar(frame='+lnT') - - fig.show() + 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) + fig.colorbar(tmp, ax=axes.ravel().tolist(), orientation="horizontal", label="mGal", fraction=.03, pad=0.08) + plt.show() .. jupyter-execute:: - fig = pygmt.Figure() - - maxabs = vd.maxabs(grid.g_ee) - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - fig.grdimage( - grid=grid.g_ee, - projection="X10c", - cmap=True, - frame=["WSne+tEasting-easting tensor", "xa","ya"], - ) - fig.colorbar(frame='+lEotvos') - - fig.shift_origin(xshift="11c") - - maxabs = vd.maxabs(grid.g_nn) - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - fig.grdimage( - grid=grid.g_nn, - projection="X10c", - cmap=True, - frame=["wSne+tNorthing-northing tensor", "xa","ya"], - ) - fig.colorbar(frame='+lEotvos') - - fig.shift_origin(xshift="11c") + fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8)) - maxabs = vd.maxabs(grid.g_zz) - pygmt.makecpt(cmap="balance+h0", series=[0, maxabs], background=True) - fig.grdimage( - grid=grid.g_zz, - projection="X10c", - cmap=True, - frame=["wSnE+tDownward-downward tensor", "xa","ya"], - ) - fig.colorbar(frame='+lEotvos') + 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) + fig.colorbar(tmp, ax=axes.ravel().tolist(), orientation="horizontal", label="Eotvos", fraction=.03, pad=0.08) + plt.show() - fig.show() .. jupyter-execute:: - fig = pygmt.Figure() - - maxabs = vd.maxabs(grid.g_en) - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - fig.grdimage( - grid=grid.g_en, - projection="X10c", - cmap=True, - frame=["WSne+tEasting-northing tensor", "xa","ya"], - ) - fig.colorbar(frame='+lEotvos') - - fig.shift_origin(xshift="11c") + fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8)) - maxabs = vd.maxabs(grid.g_ez) - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - fig.grdimage( - grid=grid.g_ez, - projection="X10c", - cmap=True, - frame=["wSne+tEasting-downward tensor", "xa","ya"], - ) - fig.colorbar(frame='+lEotvos') - - fig.shift_origin(xshift="11c") - - maxabs = vd.maxabs(grid.g_nz) - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - fig.grdimage( - grid=grid.g_nz, - projection="X10c", - cmap=True, - frame=["wSnE+tNorthing-downward tensor", "xa","ya"], - ) - fig.colorbar(frame='+lEotvos') - - fig.show() + 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) + fig.colorbar(tmp, ax=axes.ravel().tolist(), orientation="horizontal", label="Eotvos", fraction=.03, pad=0.08) + plt.show() Passing multiple prisms @@ -277,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 ) @@ -302,15 +219,7 @@ Lets plot this gravitational field: 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 @@ -373,43 +282,26 @@ We can reshape the results into variables of an dataset: .. jupyter-execute:: - fig = pygmt.Figure() + fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8)) - maxabs = vd.maxabs(grid.b_e) - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - fig.grdimage( - grid=grid.b_e, - projection="X10c", - cmap=True, - frame=["WSne+tEasting component", "xa","ya"], - ) - fig.colorbar(frame='+lnT') - - fig.shift_origin(xshift="11c") - - maxabs = vd.maxabs(grid.b_n) - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - fig.grdimage( - grid=grid.b_n, - projection="X10c", - cmap=True, - frame=["wSne+tNorthing component", "xa","ya"], - ) - fig.colorbar(frame='+lnT') - - fig.shift_origin(xshift="11c") - - maxabs = vd.maxabs(grid.b_u) - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - fig.grdimage( - grid=grid.b_u, - projection="X10c", - cmap=True, - frame=["wSnE+tUpward component", "xa","ya"], - ) - fig.colorbar(frame='+lnT') - - fig.show() + 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() Alternatively, we can compute just a single component by choosing ``field`` to @@ -435,22 +327,16 @@ generated by these two prisms: .. jupyter-execute:: - fig = pygmt.Figure() + maxabs = vd.maxabs(b_u) - grid = vd.make_xarray_grid( - coordinates, b_u, data_names="b_u", extra_coords_names="extra") - - maxabs = vd.maxabs(grid.b_u) - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - fig.grdimage( - grid=grid.b_u, - projection="X10c", - cmap=True, - frame=["WSne+tUpward component", "xa","ya"], + tmp = plt.pcolormesh( + coordinates[0], coordinates[1], b_u, vmin=-maxabs, vmax=maxabs, cmap="RdBu_r" ) - fig.colorbar(frame='+lnT') - - fig.show() + 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, shrink=0.8) + plt.show() .. _prism_layer: @@ -538,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 49f781fac..2250b95ee 100644 --- a/doc/user_guide/forward_modelling/tesseroid.rst +++ b/doc/user_guide/forward_modelling/tesseroid.rst @@ -116,7 +116,6 @@ And finally plot the computed gravitational field .. jupyter-execute:: import pygmt - grid = vd.make_xarray_grid( coordinates, gravity, data_names="gravity", extra_coords_names="extra") @@ -133,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 4d9b41203..93ecdcf02 100644 --- a/doc/user_guide/gravity_disturbance.rst +++ b/doc/user_guide/gravity_disturbance.rst @@ -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="balance+h0") + 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+lmGal', 'x+lGravity 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/topographic_correction.rst b/doc/user_guide/topographic_correction.rst index 7e980e9db..f85abbb1c 100644 --- a/doc/user_guide/topographic_correction.rst +++ b/doc/user_guide/topographic_correction.rst @@ -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="balance+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() @@ -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 142833589..4515e540b 100644 --- a/doc/user_guide/transformations.rst +++ b/doc/user_guide/transformations.rst @@ -29,19 +29,19 @@ We can load the data file using :mod:`xarray`: And plot it: .. jupyter-execute:: - :hide-code: - -.. 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:: In case we have a regular grid defined in geographic coordinates (longitude, @@ -65,23 +65,17 @@ And plot it: .. jupyter-execute:: - fig = pygmt.Figure() - - maxabs = vd.maxabs(deriv_upward) * .5 - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - - fig.grdimage( - deriv_upward, - projection="X15c", - cmap=True, - frame=["af", "WeSn+tUpward derivative of the magnetic anomaly"] - ) - fig.colorbar( - position="JCB+e", - frame=["af", 'x+lnT/m'], + maxabs = vd.maxabs(deriv_upward, percentile=99.9) + tmp = deriv_upward.plot( + cmap="RdBu_r", + vmin=-maxabs, + vmax=maxabs, + cbar_kwargs={"label":"nT/m"}, ) - - fig.show() + plt.gca().set_aspect("equal") + plt.title("Magnetic anomaly grid") + plt.gca().ticklabel_format(style="sci", scilimits=(0, 0)) + plt.show() Horizontal derivatives @@ -105,35 +99,37 @@ And plot them: .. jupyter-execute:: - fig = pygmt.Figure() - - maxabs = vd.maxabs(deriv_easting, deriv_northing) * .4 - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - - fig.grdimage( - deriv_easting, - projection="X15c", - cmap=True, - frame=["af", "WeSn+tEasting derivative of the magnetic anomaly"] + fig, (ax1, ax2) = plt.subplots( + nrows=1, ncols=2, sharey=True, figsize=(12, 8) ) + maxabs = vd.maxabs(deriv_easting, deriv_northing, percentile=99) - fig.shift_origin(xshift="16c") - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - - fig.grdimage( - deriv_northing, - projection="X15c", - cmap=True, - frame=["af", "wESn+tNorthing derivative of the magnetic anomaly"] + tmp = deriv_easting.plot( + ax=ax1, + add_colorbar=False, + vmin=-maxabs, + vmax=maxabs, + cmap='RdBu_r', ) - fig.colorbar( - position="JBC+h+e+o-8c/1c+w15c/.8c", - frame=["af", 'x+lnT/m'], + tmp = deriv_northing.plot( + ax=ax2, + add_colorbar=False, + vmin=-maxabs, + vmax=maxabs, + cmap='RdBu_r', ) - fig.show() + 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") + 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 central finite differences methods. We can choose to use either the finite @@ -154,35 +150,37 @@ frequency domain: .. jupyter-execute:: - fig = pygmt.Figure() - - maxabs = vd.maxabs(deriv_easting, deriv_northing) * .4 - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - - fig.grdimage( - deriv_easting, - projection="X15c", - cmap=True, - frame=["af", "WeSn+tEasting derivative of the magnetic anomaly"] + fig, (ax1, ax2) = plt.subplots( + nrows=1, ncols=2, sharey=True, figsize=(12, 8) ) + maxabs = vd.maxabs(deriv_easting, deriv_northing, percentile=99) - fig.shift_origin(xshift="16c") - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - - fig.grdimage( - deriv_northing, - projection="X15c", - cmap=True, - frame=["af", "wESn+tNorthing derivative of the magnetic anomaly"] + tmp = deriv_easting.plot( + ax=ax1, + add_colorbar=False, + vmin=-maxabs, + vmax=maxabs, + cmap='RdBu_r', ) - fig.colorbar( - position="JBC+h+e+o-8c/1c+w15c/.8c", - frame=["af", 'x+lnT/m'], + tmp = deriv_northing.plot( + ax=ax2, + add_colorbar=False, + vmin=-maxabs, + vmax=maxabs, + cmap='RdBu_r', ) - fig.show() + 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") + 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:: @@ -226,23 +224,11 @@ Now we can plot it: .. jupyter-execute:: - fig = pygmt.Figure() - - maxabs = vd.maxabs(upward_continued) - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - - fig.grdimage( - upward_continued, - projection="X15c", - cmap=True, - frame=["af", "WeSn+tUpward continued magnetic anomaly at 1000m"] - ) - fig.colorbar( - position="JCB", - frame=["af", 'x+lnT'], - ) - - fig.show() + 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.show() Reduction to the pole @@ -305,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 @@ -345,23 +349,14 @@ example of what happens to the reduction using arbitrary values: .. jupyter-execute:: - fig = pygmt.Figure() - - maxabs = vd.maxabs(rtp_grid) * .8 - pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs], background=True) - - fig.grdimage( - rtp_grid, - projection="X15c", - cmap=True, - frame=["af", "WeSn+tReduced to the pole with remanence"] + maxabs = vd.maxabs(rtp_grid, percentile=99.9) + tmp = rtp_grid.plot( + cmap="RdBu_r", vmin=-maxabs, vmax=maxabs, cbar_kwargs={'label':'nT'}, ) - fig.colorbar( - position="JCB+e", - frame=["af", 'x+lnT'], - ) - - fig.show() + plt.gca().set_aspect("equal") + plt.title("Reduced to the pole with remanence") + plt.gca().ticklabel_format(style="sci", scilimits=(0, 0)) + plt.show() Gaussian filters @@ -404,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() @@ -467,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() ----