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..81551a8fe 100644 --- a/doc/gallery_src/equivalent_sources/spherical.py +++ b/doc/gallery_src/equivalent_sources/spherical.py @@ -80,21 +80,23 @@ # 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", + projection="M12c", region=region, - frame=["WSne", "xa5", "ya4"], + frame=[f"WSne+t{title}", "xa5", "ya4"], x=longitude, y=latitude, fill=gravity_disturbance, @@ -102,16 +104,19 @@ 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") +title = "Gridded and upward-continued" + fig.grdimage( - frame=["ESnw", "xa5", "ya4"], + frame=[f"ESnw+t{title}", "xa5", "ya4"], grid=grid.gravity_disturbance, cmap=True, + nan_transparent=True, ) -fig.colorbar(cmap=True, frame=["a100f50", "x+lmGal"]) +fig.colorbar(cmap=True, frame=["x+lmGal"], position="+e") fig.show() diff --git a/doc/gallery_src/forward/point_gravity.py b/doc/gallery_src/forward/point_gravity.py index 38c4a2214..9c5be7555 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( @@ -68,8 +68,17 @@ 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"]) -fig.colorbar(cmap=True, position="JMR", frame=["a.6f.2", "x+lmGal"]) +fig.legend() fig.show() diff --git a/doc/gallery_src/forward/prism_layer.py b/doc/gallery_src/forward/prism_layer.py index bec45c826..22a39ba9c 100644 --- a/doc/gallery_src/forward/prism_layer.py +++ b/doc/gallery_src/forward/prism_layer.py @@ -51,17 +51,17 @@ # Plot gravity field fig = pygmt.Figure() -title = "Gravitational acceleration of a layer of prisms" +title = "Gravitational acceleration of topography with prisms" with pygmt.config(FONT_TITLE="14p"): fig.grdimage( region=region, projection="X10c/10c", grid=grid.gravity, - frame=["a", f"+t{title}", 'x+l"easting (m)"', 'y+l"northing (m)"'], + frame=["a", f"+t{title}", "x+leasting (m)", "y+lnorthing (m)"], cmap="viridis", ) -fig.colorbar(cmap=True, position="JMR", frame=["a2f1", "x+lmGal"]) +fig.colorbar(cmap=True, position="JMR", frame=["x+lmGal"]) fig.show() diff --git a/doc/gallery_src/forward/prisms_topo_gravity.py b/doc/gallery_src/forward/prisms_topo_gravity.py index d71d7f651..de4046264 100644 --- a/doc/gallery_src/forward/prisms_topo_gravity.py +++ b/doc/gallery_src/forward/prisms_topo_gravity.py @@ -87,15 +87,21 @@ title = "Gravitational acceleration of the topography" +# Get the max absolute value to use as color scale limits +cpt_lims = vd.maxabs(grid.gravity) + +# Make colormap of data +pygmt.makecpt(cmap="balance+h0", series=[-cpt_lims, cpt_lims]) + with pygmt.config(FONT_TITLE="14p"): fig.grdimage( region=xy_region, projection=fig_proj, grid=grid.gravity, frame=["ag", f"+t{title}"], - cmap="vik", + cmap=True, ) -fig.colorbar(cmap=True, frame=["a100f50", "x+lmGal"]) +fig.colorbar(cmap=True, position="JMR", frame=["a100f50", "x+lmGal"]) fig.show() diff --git a/doc/gallery_src/forward/tesseroid.py b/doc/gallery_src/forward/tesseroid.py index 9e805598b..d18a50fc7 100644 --- a/doc/gallery_src/forward/tesseroid.py +++ b/doc/gallery_src/forward/tesseroid.py @@ -50,7 +50,7 @@ # Plot the gravitational field fig = pygmt.Figure() -title = "Downward component of gravitational acceleration" +title = "Gravitational acceleration of a tesseroid" with pygmt.config(FONT_TITLE="16p"): fig.grdimage( @@ -61,7 +61,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_layer.py b/doc/gallery_src/forward/tesseroid_layer.py index 4d26564df..a760cbcb4 100644 --- a/doc/gallery_src/forward/tesseroid_layer.py +++ b/doc/gallery_src/forward/tesseroid_layer.py @@ -55,14 +55,24 @@ # 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]) + +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..3142e3f0f 100644 --- a/doc/gallery_src/forward/tesseroid_variable_density.py +++ b/doc/gallery_src/forward/tesseroid_variable_density.py @@ -78,7 +78,7 @@ def density(radius): # Plot the gravitational field fig = pygmt.Figure() -title = "Downward component of gravitational acceleration" +title = "Gravitational acceleration of variable density tesseroids" with pygmt.config(FONT_TITLE="16p"): fig.grdimage( @@ -89,7 +89,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") 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..0741eb28c 100644 --- a/doc/gallery_src/transformations/reduction_to_pole.py +++ b/doc/gallery_src/transformations/reduction_to_pole.py @@ -52,27 +52,30 @@ # 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( - grid=magnetic_grid, - projection="X?", - cmap=True, + grid=magnetic_grid, projection="X?", cmap=True, frame="+tMagnetic anomaly" ) # Add colorbar fig.colorbar( - frame='af+l"Magnetic anomaly [nT]"', - position="JBC+h+o0/1c+e", + 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]"', - position="JBC+h+o0/1c+e", + 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 21c710e5a..049eae1df 100644 --- a/doc/gallery_src/transformations/tga.py +++ b/doc/gallery_src/transformations/tga.py @@ -45,28 +45,34 @@ with fig.subplot(nrows=1, ncols=2, figsize=("28c", "15c"), sharey="l"): with fig.set_panel(panel=0): # Make colormap of data - scale = 2500 - pygmt.makecpt(cmap="polar+h", series=[-scale, scale], background=True) + 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, projection="X?", cmap=True, + frame="+tMagnetic Anomaly", ) # Add colorbar fig.colorbar( - frame='af+l"Magnetic anomaly [nT]"', - position="JBC+h+o0/1c+e", + frame="af+lnT", + 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) + fig.grdimage( + grid=tga, + projection="X?", + cmap=True, + frame="+tTotal Gradient Amplitude", + ) # Add colorbar fig.colorbar( - frame='af+l"Total Gradient Amplitude [nT/m]"', - position="JBC+h+o0/1c+e", + frame="af+lnT/m", + position="JBC+h+o0/1c+ef", ) fig.show() diff --git a/doc/gallery_src/transformations/tilt.py b/doc/gallery_src/transformations/tilt.py index 9f9c1b579..96dad14df 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, @@ -101,16 +101,15 @@ frame=["a", "+tReduced to the pole (RTP)"], ) # Add colorbar - label = "nT" fig.colorbar( - frame=f"af+l{label}", - position="JMR+o1/-0.25c+e", + frame="af+lnT", + 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 +119,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, @@ -129,9 +128,8 @@ frame=["a", "+tTilt of RTP grid"], ) # Add colorbar - label = "rad" fig.colorbar( - frame=f"af+l{label}", - position="JMR+o1/-0.25c+e", + 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 db6ce9829..84ca04f2c 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,27 +45,33 @@ # 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( 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]"', - position="JBC+h+o0/1c+e", + 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 2d48b89db..a77856dfb 100644 --- a/doc/gallery_src/transformations/upward_derivative.py +++ b/doc/gallery_src/transformations/upward_derivative.py @@ -46,28 +46,34 @@ with fig.subplot(nrows=1, ncols=2, figsize=("28c", "15c"), sharey="l"): with fig.set_panel(panel=0): # Make colormap of data - scale = 2500 - pygmt.makecpt(cmap="polar+h", series=[-scale, scale], background=True) + 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, projection="X?", cmap=True, + frame="+tMagnetic Anomaly", ) # Add colorbar fig.colorbar( - frame='af+l"Magnetic anomaly [nT]"', + frame="af+lnT", position="JBC+h+o0/1c+e", ) with fig.set_panel(panel=1): # Make colormap for upward derivative (saturate it a little bit) - scale = 0.6 * vd.maxabs(deriv_upward) - pygmt.makecpt(cmap="polar+h", series=[-scale, scale], background=True) + 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/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..2053c2dc5 100644 --- a/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst +++ b/doc/user_guide/equivalent_sources/eqs-parameters-estimation.rst @@ -80,7 +80,7 @@ value of the score obtained after each cross validation. data.gravity_disturbance_mgal, ) ) - score_first_guess + print(score_first_guess) The resulting score corresponds to the R^2. It represents how well the equivalent sources can reproduce the variation of our data. As closer it gets @@ -145,7 +145,7 @@ And now we can actually ran one cross validation for each pair of parameters: ) ) scores.append(score) - scores + print(scores) Once every score has been computed, we can obtain the best score and the corresponding parameters that generate it: @@ -216,7 +216,7 @@ Lets plot it: fig = pygmt.Figure() # Make colormap of data - pygmt.makecpt(cmap="polar+h0",series=(-maxabs, maxabs,)) + pygmt.makecpt(cmap="balance+h0",series=(-maxabs, maxabs,)) 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/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/point.rst b/doc/user_guide/forward_modelling/point.rst index 004058c0d..cf03b5178 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), @@ -238,15 +238,14 @@ 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)*.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), 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/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/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..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() @@ -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", @@ -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/igrf.rst b/doc/user_guide/igrf.rst index 9eff60d82..ce0623ee0 100644 --- a/doc/user_guide/igrf.rst +++ b/doc/user_guide/igrf.rst @@ -8,6 +8,17 @@ magnetic field from 1900 until 5 years after the latest release (based on predictions of the secular variation). Harmonica allows calculating the 14th generation IGRF field with :class:`harmonica.IGRF14`. Here's how it works. +.. jupyter-execute:: + :hide-code: + + import pygmt + + # Needed so that displaying works on jupyter-sphinx and sphinx-gallery at + # the same time. Using PYGMT_USE_EXTERNAL_DISPLAY="false" in the Makefile + # for sphinx-gallery to work means that fig.show won't display anything here + # either. + pygmt.set_display(method="notebook") + .. jupyter-execute:: import datetime @@ -15,6 +26,7 @@ generation IGRF field with :class:`harmonica.IGRF14`. Here's how it works. import pygmt import boule as bl import harmonica as hm + import verde as vd All of the functionality is wrapped in the :class:`~harmonica.IGRF14` class. When creating an instance of it, we need to provide the date on which we want @@ -147,24 +159,31 @@ We can plot the three components using :mod:`pygmt` on a nice map: .. jupyter-execute:: fig = pygmt.Figure() + for c in ["b_east", "b_north", "b_up"]: - fig.grdimage( - grid[c], cmap="polar+h", projection="W15c", region="g", - ) - fig.coast(shorelines=True) - fig.colorbar( - position="JMR+ml+o0.5c", - frame=[ + max_abs = vd.maxabs(grid[c]) + pygmt.makecpt( + cmap="balance+h0", + series=[-max_abs, max_abs], + background=True, + ) + fig.grdimage( + grid[c], cmap=True, projection="W15c", region="g", + ) + fig.coast(shorelines=True) + fig.colorbar( + position="JMR+ml+o0.5c", + frame=[ "a10000", f"x+l{grid[c].attrs['long_name']}", f"y+l{grid[c].attrs['units']}", - ] - ) - if c == "b_east": - fig.basemap(frame=["a", f"+t{grid.attrs['title']}"]) - else: - fig.basemap(frame="a") - fig.shift_origin(yshift="-h-0.5c") + ] + ) + if c == "b_east": + fig.basemap(frame=["a", f"+t{grid.attrs['title']}"]) + else: + fig.basemap(frame="a") + fig.shift_origin(yshift="-h-0.5c") fig.show() The grid spacing was calculated automatically to match the maximum degree of @@ -201,33 +220,33 @@ with PyGMT the same way we did the vector components: fig = pygmt.Figure() cmaps = { - "intensity": "viridis", - "inclination": "polar+h", - "declination": "polar+h", + "intensity": "viridis", + "inclination": "balance+h0", + "declination": "balance+h0", } cb_annot = { - "intensity": "a10000", - "inclination": "a20+u\\260", - "declination": "a40+u\\260", + "intensity": "a10000", + "inclination": "a20+u\\260", + "declination": "a40+u\\260", } for c in ["intensity", "inclination", "declination"]: - fig.grdimage( - grid[c], cmap=cmaps[c], projection="W0/15c", - ) - fig.coast(shorelines=True) - fig.colorbar( - position="JMR+ml+o0.5c", - frame=[ + fig.grdimage( + grid[c], cmap=cmaps[c], projection="W0/15c", + ) + fig.coast(shorelines=True) + fig.colorbar( + position="JMR+ml+o0.5c", + frame=[ cb_annot[c], f"x+l{grid[c].attrs['long_name']}", f"y+l{grid[c].attrs['units']}", - ] - ) - if c == "intensity": - fig.basemap(frame=["a", f"+t{grid.attrs['title']}"]) - else: - fig.basemap(frame="a") - fig.shift_origin(yshift="-h-0.5c") + ] + ) + if c == "intensity": + fig.basemap(frame=["a", f"+t{grid.attrs['title']}"]) + else: + fig.basemap(frame="a") + fig.shift_origin(yshift="-h-0.5c") fig.show() We can clearly see the `South Atlantic Magnetic Anomaly diff --git a/doc/user_guide/topographic_correction.rst b/doc/user_guide/topographic_correction.rst index ec7b60c07..7e980e9db 100644 --- a/doc/user_guide/topographic_correction.rst +++ b/doc/user_guide/topographic_correction.rst @@ -72,7 +72,7 @@ And plot it: maxabs = vd.maxabs(data.gravity_disturbance_mgal) fig = pygmt.Figure() - pygmt.makecpt(cmap="polar+h0", series=[-maxabs, maxabs]) + pygmt.makecpt(cmap="balance+h0", series=[-maxabs, maxabs]) fig.plot( x=data.longitude, y=data.latitude, @@ -82,7 +82,7 @@ And plot it: projection="M15c", frame=['ag', 'WSen'], ) - fig.colorbar(cmap=True, frame=["a50f25", "x+lgravity disturbance", "y+lmGal"]) + fig.colorbar(cmap=True, frame=["a50f25", "x+lGravity disturbance", "y+lmGal"]) fig.show() @@ -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, 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() ----