Plotting Earth relief

Plotting a map of Earth relief can use the data accessed by the pygmt.datasets.load_earth_relief method. The data can then be plotted using the pygmt.Figure.grdimage method.

import pygmt

Load sample Earth relief data for the entire globe at a resolution of 30 minutes. The other available resolutions are show at https://docs.generic-mapping-tools.org/latest/datasets/remote-data.html#global-earth-relief-grids.

grid = pygmt.datasets.load_earth_relief(resolution="30m")

Out:

gmtwhich [NOTICE]: Remote data courtesy of GMT data server OCEANIA [https://oceania.generic-mapping-tools.org]
gmtwhich [NOTICE]: Earth Relief at 30x30 arc minutes from Gaussian Cartesian filtering (55 km fullwidth) of SRTM15+V2.1 [Tozer et al., 2019].
gmtwhich [NOTICE]:   -> Download grid file [395K]: earth_relief_30m_p.grd

Create a plot

The pygmt.Figure.grdimage method takes the grid input to create a figure. It creates and applies a color palette to the figure based upon the z-values of the data. By default, it plots the map with the turbo CPT, an equidistant cylindrical projection, and with no frame.

fig = pygmt.Figure()
fig.grdimage(grid=grid)
fig.show()

Out:

<IPython.core.display.Image object>

pygmt.Figure.grdimage can take the optional argument projection for the map. In the example below, the projection is set as "R12c" for 12 centimeter figure with a Winkel Tripel projection. For a list of available projections, see https://docs.generic-mapping-tools.org/latest/cookbook/map-projections.html.

fig = pygmt.Figure()
fig.grdimage(grid=grid, projection="R12c")
fig.show()

Out:

<IPython.core.display.Image object>

Set a color map

pygmt.Figure.grdimage takes the cmap argument to set the CPT of the figure. Examples of common CPTs for Earth relief are shown below. A full list of CPTs can be found at https://docs.generic-mapping-tools.org/latest/cookbook/cpts.html.

Using the geo CPT:

fig = pygmt.Figure()
fig.grdimage(grid=grid, projection="R12c", cmap="geo")
fig.show()

Out:

<IPython.core.display.Image object>

Using the relief CPT:

fig = pygmt.Figure()
fig.grdimage(grid=grid, projection="R12c", cmap="relief")
fig.show()

Out:

<IPython.core.display.Image object>

Add a color bar

The pygmt.Figure.colorbar method displays the CPT and the associated Z-values of the figure, and by default uses the same CPT set by the cmap argument for pygmt.Figure.grdimage. The frame argument for pygmt.Figure.colorbar can be used to set the axis intervals and labels. A list is used to pass multiple arguments to frame. In the example below, "a2500" sets the axis interval to 2,500, "x+lElevation" sets the x-axis label, and "y+lm" sets the y-axis label.

fig = pygmt.Figure()
fig.grdimage(grid=grid, projection="R12c", cmap="geo")
fig.colorbar(frame=["a2500", "x+lElevation", "y+lm"])
fig.show()

Out:

<IPython.core.display.Image object>

Create a region map

In addition to providing global data, the region argument for pygmt.datasets.load_earth_relief can be used to provide data for a specific area. The region argument is required for resolutions at 5 minutes or higher, and accepts a list (as in the example below) or a string. The geographic ranges are passed as x-min/x-max/y-min/y-max.

The example below uses data with a 5 minute resolution, and plots it on a 15 centimeter figure with a Mercator projection and a CPT set to geo. frame="a" is used to add a frame to the figure.

grid = pygmt.datasets.load_earth_relief(resolution="05m", region=[-14, 30, 35, 60])
fig = pygmt.Figure()
fig.grdimage(grid=grid, projection="M15c", frame="a", cmap="geo")
fig.colorbar(frame=["a1000", "x+lElevation", "y+lm"])
fig.show()

Out:

grdcut [DEBUG]: gmt_get_filename: In: -14/30/35/60 Out: -14/30/35/60
grdcut [DEBUG]: Look for file -14/30/35/60 in /vercel/.gmt
grdcut [DEBUG]: Look for file -14/30/35/60 in /vercel/.gmt/cache
grdcut [DEBUG]: Look for file -14/30/35/60 in /vercel/.gmt/server
grdcut [DEBUG]: Got regular w/e/s/n for region (-14/30/35/60)
grdcut [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdcut [DEBUG]: Given full path to file /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdcut [DEBUG]: Replace file /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000 with /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdcut [INFORMATION]: Processing input grid
grdcut [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdcut [DEBUG]: Given full path to file /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdcut [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdcut [DEBUG]: Given full path to file /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdcut [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdcut [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdcut [DEBUG]: Found readable file /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdcut [DEBUG]: Object ID 0 : Registered Grid File /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000 as an Input resource with geometry Surface [n_objects = 1]
grdcut [DEBUG]: gmtapi_begin_io: Input resource access is now enabled [container]
grdcut [DEBUG]: gmtapi_import_grid: Passed ID = 0 and mode = 1
grdcut [DEBUG]: Object ID 1 : Registered Grid Memory Copy 55806faaaf30 as an Output resource with geometry Surface [n_objects = 2]
grdcut [DEBUG]: Successfully created a new Grid container
grdcut [DEBUG]: VirtualFile name created: @GMTAPI@-S-O-G-G-G-Y-000001
grdcut [DEBUG]: GMT now running in modern mode [Session ID = 4587]
grdcut [DEBUG]: Revised options: @earth_relief_05m_p/ -R-14/30/35/60 -I05m -rp -G@GMTAPI@-S-O-G-G-G-Y-000001 -fg -Co+n
grdcut (gmtlib_free_tmp_arrays): tried to free unallocated memory
grdblend [DEBUG]: History: Process -R-14/30/35/60
grdblend [DEBUG]: Map distance calculation will be using great circle approximation with authalic auxiliary latitudes and authalic (R_2) radius = 6371007.1809 m, in meter.
grdblend [DEBUG]: gmt_get_filename: In: -14/30/35/60 Out: -14/30/35/60
grdblend [DEBUG]: Look for file -14/30/35/60 in /vercel/.gmt
grdblend [DEBUG]: Look for file -14/30/35/60 in /vercel/.gmt/cache
grdblend [DEBUG]: Look for file -14/30/35/60 in /vercel/.gmt/server
grdblend [DEBUG]: Got regular w/e/s/n for region (-14/30/35/60)
grdblend [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdblend [DEBUG]: Given full path to file /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdblend [DEBUG]: Replace file /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000 with /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdblend [DEBUG]: Map distance calculation will be using great circle approximation with authalic auxiliary latitudes and authalic (R_2) radius = 6371007.1809 m, in meter.
grdblend [DEBUG]: gmtapi_init_import: Passed family = Data Table and geometry = Text
grdblend [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdblend [DEBUG]: Given full path to file /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdblend [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdblend [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdblend [DEBUG]: Found readable file /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdblend [DEBUG]: Object ID 2 : Registered Data Table File /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000 as an Input resource with geometry Text [n_objects = 3]
grdblend (gmtapi_init_import): tried to free unallocated memory
grdblend [DEBUG]: gmtapi_init_import: Added 1 new sources
grdblend [DEBUG]: GMT_Init_IO: Returned first Input object ID = 2
grdblend [DEBUG]: GMT_Begin_IO: Mode value 1 not considered (ignored)
grdblend [DEBUG]: GMT_Begin_IO: Initialize record-by-record access for Input
grdblend [DEBUG]: gmtapi_next_io_source: Selected object 2
grdblend [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdblend [DEBUG]: Given full path to file /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdblend [DEBUG]: Found readable file /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdblend [DEBUG]: Found readable file /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdblend [INFORMATION]: Reading Data Table from file /vercel/.gmt/sessions/gmt_session.4587/=tiled_84_GX.000000
grdblend [DEBUG]: GMT_Begin_IO: Input resource access is now enabled [record-by-record]
grdblend [DEBUG]: Geographic input grid, longitudes span less than 360
grdblend [DEBUG]: Chosen boundary condition for all edges: geographic
grdblend [DEBUG]: Geographic input grid, longitudes span less than 360
grdblend [DEBUG]: Object ID 3 : Registered Grid Memory Reference 5580701beb90 as an Input resource with geometry Surface [n_objects = 4]
grdblend [DEBUG]: Successfully created a new Grid container
grdblend [DEBUG]: gmt_get_filename: In: S90W180.earth_relief_05m_p.nc Out: S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/cache
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/server
grdblend [DEBUG]: gmt_get_filename: In: S90E000.earth_relief_05m_p.nc Out: S90E000.earth_relief_05m_p.nc
grdblend [DEBUG]: Look for file S90E000.earth_relief_05m_p.nc in /vercel/.gmt
grdblend [DEBUG]: Look for file S90E000.earth_relief_05m_p.nc in /vercel/.gmt/cache
grdblend [DEBUG]: Look for file S90E000.earth_relief_05m_p.nc in /vercel/.gmt/server
grdblend [DEBUG]: Remote file @S90W180.earth_relief_05m_p.nc exists locally as /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: gmt_get_filename: In: S90W180.earth_relief_05m_p.nc Out: S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/cache
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/server
grdblend [DEBUG]: Found readable file /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: Replace file S90W180.earth_relief_05m_p.nc with path /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: gmt_get_filename: In: S90W180.earth_relief_05m_p.nc Out: S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/cache
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/server
grdblend [DEBUG]: gmt_get_filename: In: S90W180.earth_relief_05m_p.nc Out: S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/cache
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/server
grdblend [DEBUG]: Found readable file /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: Object ID 4 : Registered Grid File S90W180.earth_relief_05m_p.nc as an Input resource with geometry Surface [n_objects = 5]
grdblend [DEBUG]: gmtapi_import_grid: Passed ID = 4 and mode = 33
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/cache
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/server
grdblend [DEBUG]: Found readable file /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: gmt_get_filename: In: /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc Out: /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: Call gmtgrdio_doctor_geo_increments on a geographic grid
grdblend [DEBUG]: Geographic input grid, longitudes span less than 360
grdblend [INFORMATION]: Downloading earth_relief_05m_p/ tile 1 of 1 [S90E000]
grdblend [DEBUG]: Get remote file https://oceania.generic-mapping-tools.org/server/earth/earth_relief/earth_relief_05m_p/S90E000.earth_relief_05m_p.jp2 and write to /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90E000.earth_relief_05m_p.jp2
grdblend [DEBUG]: Download https://oceania.generic-mapping-tools.org/server/earth/earth_relief/earth_relief_05m_p/S90E000.earth_relief_05m_p.jp2 to /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90E000.earth_relief_05m_p.jp2
grdblend [NOTICE]: Remote data courtesy of GMT data server OCEANIA [https://oceania.generic-mapping-tools.org]
grdblend [NOTICE]: Earth Relief at 5x5 arc minutes from Gaussian Cartesian filtering (9 km fullwidth) of SRTM15+V2.1 [Tozer et al., 2019].
grdblend [NOTICE]:   -> Download 180x180 degree grid tile (earth_relief_05m_p): S90E000
grdblend [INFORMATION]: Downloading file https://oceania.generic-mapping-tools.org/server/earth/earth_relief/earth_relief_05m_p/S90E000.earth_relief_05m_p.jp2 ...
grdblend [DEBUG]: Delete /tmp/S90E000.earth_relief_05m_p.nc.download
grdblend [INFORMATION]: Convert SRTM tile from JPEG2000 to netCDF grid [/vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90E000.earth_relief_05m_p.nc]
grdblend [DEBUG]: Running: grdconvert /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90E000.earth_relief_05m_p.jp2 -G/vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90E000.earth_relief_05m_p.nc=ns+s0.5+o0 -Z+s0.5+o-0 -fg -Vq --IO_NC4_DEFLATION_LEVEL=9 --GMT_HISTORY=falsen
grdblend [DEBUG]: GMT now running in modern mode [Session ID = 4587]
grdblend [DEBUG]: Revised options: /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90E000.earth_relief_05m_p.jp2 -G/vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90E000.earth_relief_05m_p.nc=ns+s0.5+o0 -Z+s0.5+o-0 -fg -Vq --IO_NC4_DEFLATION_LEVEL=9 --GMT_HISTORY=false
grdblend (gmtlib_free_tmp_arrays): tried to free unallocated memory
@earth_relief_05m -G/tmp/pygmt-55ykmfkf.nc -R-14/30/35/60 -Vd
/tmp/pygmt-55ykmfkf.nc
True
274652

<IPython.core.display.Image object>