Replace Cartopy with PyGMT

data/examples/britain_magnetic.py

@@ -29,9 +29,8 @@
 If the file isn't already in your data directory, it will be downloaded
 automatically.
 """
-import cartopy.crs as ccrs
-import matplotlib.pyplot as plt
 import numpy as np
+import pygmt
 import verde as vd
 
 import harmonica as hm
@@ -40,31 +39,33 @@
 data = hm.datasets.fetch_britain_magnetic()
 print(data)
 
-# Plot the observations in a Mercator map using Cartopy
-fig = plt.figure(figsize=(7.5, 10))
-ax = plt.axes(projection=ccrs.Mercator())
-ax.set_title("Magnetic data from Great Britain", pad=25)
-maxabs = np.percentile(data.total_field_anomaly_nt, 99)
-tmp = ax.scatter(
- data.longitude,
- data.latitude,
- c=data.total_field_anomaly_nt,
- s=0.001,
- cmap="seismic",
- vmin=-maxabs,
- vmax=maxabs,
- transform=ccrs.PlateCarree(),
+# Get the region of the data
+region = vd.get_region((data.longitude, data.latitude))
+
+# Plot the observations in a Mercator map using PyGMT
+fig = pygmt.Figure()
+
+title = "Magnetic data from Great Britain"
+
+# Make colormap scaled to 97 percentile of data.
+maxabs = np.percentile(data.total_field_anomaly_nt, 97)
+pygmt.makecpt(cmap="polar", series=(-maxabs, maxabs), background=True)
+
+fig.plot(
+ region=region,
+ projection="M10c",
+ frame=["ag", f"+t{title}"],
+ x=data.longitude,
+ y=data.latitude,
+ color=data.total_field_anomaly_nt,
+ style="c0.02c",
+ cmap=True,
 )
-plt.colorbar(
- tmp,
- ax=ax,
- label="total field magnetic anomaly [nT]",
- orientation="vertical",
- aspect=50,
- shrink=0.7,
- pad=0.1,
+
+fig.colorbar(
+ cmap=True, position="JMR", frame=["a100f50", "x+ltotal field magnetic anomaly [nT]"]
 )
-ax.set_extent(vd.get_region((data.longitude, data.latitude)))
-ax.gridlines(draw_labels=True)
-ax.coastlines(resolution="50m")
-plt.show()
+
+fig.coast(shorelines="1p,black")
+
+fig.show()

data/examples/earth_geoid.py

@@ -16,22 +16,29 @@
 with 0.5 degree spacing and was generated from the spherical harmonic model
 EIGEN-6C4 [Forste_etal2014]_.
 """
-import cartopy.crs as ccrs
-import matplotlib.pyplot as plt
+import pygmt
 
 import harmonica as hm
 
 # Load the geoid grid
 data = hm.datasets.fetch_geoid_earth()
 print(data)
 
-# Make a plot of data using Cartopy
-plt.figure(figsize=(10, 10))
-ax = plt.axes(projection=ccrs.Orthographic(central_longitude=100))
-pc = data.geoid.plot.pcolormesh(ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False)
-plt.colorbar(
- pc, label="meters", orientation="horizontal", aspect=50, pad=0.01, shrink=0.6
+# Make a plot using PyGMT
+fig = pygmt.Figure()
+
+title = "Geoid heights (EIGEN-6C4)"
+
+fig.grdimage(
+ region="g",
+ projection="G100/0/15c",
+ frame=f"+t{title}",
+ grid=data.geoid,
+ cmap="vik",
 )
-ax.set_title("Geoid heights (EIGEN-6C4)")
-ax.coastlines()
-plt.show()
+
+fig.coast(shorelines="0.5p,black", resolution="crude")
+
+fig.colorbar(cmap=True, frame=["a25f10", "x+lmeters"])
+
+fig.show()

data/examples/earth_gravity.py

@@ -13,24 +13,29 @@
 spacing at 10km ellipsoidal height. It was generated from the spherical
 harmonic model EIGEN-6C4 [Forste_etal2014]_.
 """
-import cartopy.crs as ccrs
-import matplotlib.pyplot as plt
+import pygmt
 
 import harmonica as hm
 
 # Load the gravity grid
 data = hm.datasets.fetch_gravity_earth()
 print(data)
 
-# Make a plot of data using Cartopy
-plt.figure(figsize=(10, 10))
-ax = plt.axes(projection=ccrs.Orthographic(central_longitude=150))
-pc = data.gravity.plot.pcolormesh(
- ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False
-)
-plt.colorbar(
- pc, label="mGal", orientation="horizontal", aspect=50, pad=0.01, shrink=0.6
+# Make a plot using Pygmt
+fig = pygmt.Figure()
+
+title = "Gravity of the Earth (EIGEN-6C4)"
+
+fig.grdimage(
+ region="g",
+ projection="G150/0/15c",
+ frame=f"+t{title}",
+ grid=data.gravity,
+ cmap="viridis",
 )
-ax.set_title("Gravity of the Earth (EIGEN-6C4)")
-ax.coastlines()
-plt.show()
+
+fig.coast(shorelines="0.5p,black", resolution="crude")
+
+fig.colorbar(cmap=True, frame=["a1000f250", "x+lmGal"])
+
+fig.show()

data/examples/earth_topography.py

@@ -13,24 +13,29 @@
 we downsampled to 0.5 degree grid spacing to save space and download times.
 Heights are referenced to sea level.
 """
-import cartopy.crs as ccrs
-import matplotlib.pyplot as plt
+import pygmt
 
 import harmonica as hm
 
 # Load the topography grid
 data = hm.datasets.fetch_topography_earth()
 print(data)
 
-# Make a plot of data using Cartopy
-plt.figure(figsize=(10, 10))
-ax = plt.axes(projection=ccrs.Orthographic(central_longitude=-30))
-pc = data.topography.plot.pcolormesh(
- ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False, cmap="terrain"
-)
-plt.colorbar(
- pc, label="meters", orientation="horizontal", aspect=50, pad=0.01, shrink=0.6
+# Make a plot using PyGMT
+fig = pygmt.Figure()
+
+title = "Topography of the Earth (ETOPO1)"
+
+fig.grdimage(
+ region="g",
+ projection="G-30/0/15c",
+ frame=f"+t{title}",
+ grid=data.topography,
+ cmap="globe",
 )
-ax.set_title("Topography of the Earth (ETOPO1)")
-ax.coastlines()
-plt.show()
+
+fig.coast(shorelines="0.5p,black", resolution="crude")
+
+fig.colorbar(cmap=True, frame=["a2000f500", "x+lmeters"])
+
+fig.show()

data/examples/south_africa_gravity.py

@@ -18,8 +18,7 @@
 documentation for :func:`harmonica.datasets.fetch_south_africa_gravity` for
 more information.
 """
-import cartopy.crs as ccrs
-import matplotlib.pyplot as plt
+import pygmt
 import verde as vd
 
 import harmonica as hm
@@ -28,22 +27,29 @@
 data = hm.datasets.fetch_south_africa_gravity()
 print(data)
 
-# Plot the observations in a Mercator map using Cartopy
-fig = plt.figure(figsize=(6.5, 5))
-ax = plt.axes(projection=ccrs.Mercator())
-ax.set_title("Observed gravity data from South Africa", pad=25)
-tmp = ax.scatter(
- data.longitude,
- data.latitude,
- c=data.gravity,
- s=0.8,
- cmap="viridis",
- transform=ccrs.PlateCarree(),
-)
-plt.colorbar(
- tmp, ax=ax, label="observed gravity [mGal]", aspect=50, pad=0.1, shrink=0.92
+# Get the region of the grid
+region = vd.get_region((data.longitude.values, data.latitude.values))
+
+# Make a plot using PyGMT
+fig = pygmt.Figure()
+
+title = "Observed gravity data from South Africa"
+
+pygmt.makecpt(cmap="viridis", series=(data.gravity.min(), data.gravity.max()))
+
+fig.plot(
+ region=region,
+ projection="M15c",
+ frame=["ag", f"+t{title}"],
+ x=data.longitude,
+ y=data.latitude,
+ color=data.gravity,
+ style="c0.1c",
+ cmap=True,
 )
-ax.set_extent(vd.get_region((data.longitude, data.latitude)))
-ax.gridlines(draw_labels=True)
-ax.coastlines()
-plt.show()
+
+fig.coast(shorelines="1p,black")
+
+fig.colorbar(cmap=True, frame=["a200f50", "x+lobserved gravity [mGal]"], position="JMR")
+
+fig.show()

data/examples/south_africa_topography.py

@@ -13,24 +13,33 @@
 we downsampled to 0.1 degree grid spacing to save space and download times.
 Heights are referenced to sea level.
 """
-import cartopy.crs as ccrs
-import matplotlib.pyplot as plt
+import pygmt
+import verde as vd
 
 import harmonica as hm
 
 # Load the topography grid
 data = hm.datasets.fetch_south_africa_topography()
 print(data)
 
-# Make a plot of data using Cartopy
-plt.figure(figsize=(10, 10))
-ax = plt.axes(projection=ccrs.Mercator())
-pc = data.topography.plot.pcolormesh(
- ax=ax, transform=ccrs.PlateCarree(), add_colorbar=False, cmap="terrain"
-)
-plt.colorbar(
- pc, label="meters", orientation="horizontal", aspect=50, pad=0.01, shrink=0.6
+# Get the region of the grid
+region = vd.get_region((data.longitude.values, data.latitude.values))
+
+# Make a plot using PyGMT
+fig = pygmt.Figure()
+
+title = "Topography of South africa (ETOPO1)"
+
+fig.grdimage(
+ region=region,
+ projection="M15c",
+ grid=data.topography,
+ frame=["ag", f"+t{title}"],
+ cmap="earth",
 )
-ax.set_title("Topography of South africa (ETOPO1)")
-ax.coastlines()
-plt.show()
+
+fig.colorbar(cmap=True, frame=["a2000f500", "x+lmeters"])
+
+fig.coast(shorelines="1p,black")
+
+fig.show()

doc/Makefile

@@ -42,7 +42,8 @@ html: api
  @echo
  @echo "Building HTML files."
  @echo
- $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html
+ # Set PYGMT_USE_EXTERNAL_DISPLAY to "false" to disable external display
+ PYGMT_USE_EXTERNAL_DISPLAY="false" $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html
  @echo
  @echo "Build finished. The HTML pages are in $(BUILDDIR)/html."
 

doc/conf.py

@@ -6,6 +6,7 @@
 #
 import datetime
 
+import pygmt.sphinx_gallery
 import pyvista
 from sphinx_gallery.sorting import FileNameSortKey
 
@@ -49,14 +50,17 @@
  "pandas": ("http://pandas.pydata.org/pandas-docs/stable/", None),
  "xarray": ("http://xarray.pydata.org/en/stable/", None),
  "xrft": ("https://xrft.readthedocs.io/en/stable/", None),
- "cartopy": ("https://scitools.org.uk/cartopy/docs/latest/", None),
  "pooch": ("https://www.fatiando.org/pooch/latest/", None),
  "verde": ("https://www.fatiando.org/verde/latest/", None),
  "boule": ("https://www.fatiando.org/boule/latest/", None),
  "matplotlib": ("https://matplotlib.org/", None),
  "pyproj": ("https://pyproj4.github.io/pyproj/stable/", None),
  "pyvista": ("https://docs.pyvista.org", None),
- "numba_progress": ("https://pypi.org/project/numba-progress/", None),
+ "numba_progress": (
+ "https://pypi.org/project/numba-progress/",
+ None,
+ ),
+ "pygmt": ("https://www.pygmt.org/latest/", None),
 }
 
 # Autosummary pages will be generated by sphinx-autogen instead of sphinx-build
@@ -103,7 +107,11 @@
  # Insert links to documentation of objects in the examples
  "reference_url": {"harmonica": None},
  # Add pyvista to the image scrapers
- "image_scrapers": ("pyvista", "matplotlib"),
+ "image_scrapers": (
+ "pyvista",
+ "matplotlib",
+ pygmt.sphinx_gallery.PyGMTScraper(),
+ ),
 }
 
 

doc/install.rst

@@ -47,8 +47,7 @@ Optional:
 The examples in the :ref:`gallery` also use:
 
 * `boule <http://www.fatiando.org/boule/>`__
-* `matplotlib <https://matplotlib.org/>`__
-* `cartopy <https://scitools.org.uk/cartopy/>`__ for plotting maps
+* `pygmt <https://www.pygmt.org/>`__ for plotting maps
 * `pyproj <https://jswhit.github.io/pyproj/>`__ for cartographic projections