Creating a 3D perspective image

Create 3-D perspective image or surface mesh from a grid using pygmt.Figure.grdview.

import pygmt

# Load sample earth relief data
grid = pygmt.datasets.load_earth_relief(resolution="05m", region=[-108, -103, 35, 40])

Out:

grdcut [DEBUG]: gmt_get_filename: In: -108/-103/35/40 Out: -108/-103/35/40
grdcut [DEBUG]: Look for file -108/-103/35/40 in /vercel/.gmt
grdcut [DEBUG]: Look for file -108/-103/35/40 in /vercel/.gmt/cache
grdcut [DEBUG]: Look for file -108/-103/35/40 in /vercel/.gmt/server
grdcut [DEBUG]: Got regular w/e/s/n for region (-108/-103/35/40)
grdcut [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdcut [DEBUG]: Given full path to file /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdcut [DEBUG]: Replace file /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000 with /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdcut [INFORMATION]: Processing input grid
grdcut [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdcut [DEBUG]: Given full path to file /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdcut [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdcut [DEBUG]: Given full path to file /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdcut [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdcut [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdcut [DEBUG]: Found readable file /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdcut [DEBUG]: Object ID 0 : Registered Grid File /vercel/.gmt/sessions/gmt_session.4708/=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 555a653d50a0 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 = 4708]
grdcut [DEBUG]: Revised options: @earth_relief_05m_p/ -R-108/-103/35/40 -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-108/-103/35/40
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: -108/-103/35/40 Out: -108/-103/35/40
grdblend [DEBUG]: Look for file -108/-103/35/40 in /vercel/.gmt
grdblend [DEBUG]: Look for file -108/-103/35/40 in /vercel/.gmt/cache
grdblend [DEBUG]: Look for file -108/-103/35/40 in /vercel/.gmt/server
grdblend [DEBUG]: Got regular w/e/s/n for region (-108/-103/35/40)
grdblend [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdblend [DEBUG]: Given full path to file /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdblend [DEBUG]: Replace file /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000 with /vercel/.gmt/sessions/gmt_session.4708/=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.4708/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdblend [DEBUG]: Given full path to file /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdblend [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdblend [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdblend [DEBUG]: Found readable file /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdblend [DEBUG]: Object ID 2 : Registered Data Table File /vercel/.gmt/sessions/gmt_session.4708/=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.4708/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdblend [DEBUG]: Given full path to file /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdblend [DEBUG]: Found readable file /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdblend [DEBUG]: Found readable file /vercel/.gmt/sessions/gmt_session.4708/=tiled_84_GX.000000
grdblend [INFORMATION]: Reading Data Table from file /vercel/.gmt/sessions/gmt_session.4708/=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 555a652ea730 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 [INFORMATION]: Downloading earth_relief_05m_p/ tile 1 of 1 [S90W180]
grdblend [DEBUG]: Get remote file https://oceania.generic-mapping-tools.org/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.jp2 and write to /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.jp2
grdblend [DEBUG]: Download https://oceania.generic-mapping-tools.org/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.jp2 to /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.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): S90W180
grdblend [INFORMATION]: Downloading file https://oceania.generic-mapping-tools.org/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.jp2 ...
grdblend [DEBUG]: Delete /tmp/S90W180.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/S90W180.earth_relief_05m_p.nc]
grdblend [DEBUG]: Running: grdconvert /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.jp2 -G/vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.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 = 4708]
grdblend [DEBUG]: Revised options: /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.jp2 -G/vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.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-fr1kea5m.nc -R-108/-103/35/40 -Vd /tmp/pygmt-fr1kea5m.nc 16192

The pygmt.Figure.grdview method takes the grid input. The perspective argument changes the azimuth and elevation of the viewpoint; the default is [180, 90], which is looking directly down on the figure and north is “up”. The zsize argument sets how tall the three-dimensional portion appears.

The default grid surface type is mesh plot.

fig = pygmt.Figure()
fig.grdview(
    grid=grid,
    # Sets the view azimuth as 130 degrees, and the view elevation as 30 degrees
    perspective=[130, 30],
    # Sets the x- and y-axis labels, and annotates the west, south, and east axes
    frame=["xa", "ya", "WSnE"],
    # Sets a Mercator projection on a 15-centimeter figure
    projection="M15c",
    # Sets the height of the three-dimensional relief at 1.5 centimeters
    zsize="1.5c",
)
fig.show()

Out:

<IPython.core.display.Image object>

The grid surface type can be set with the surftype parameter.

fig = pygmt.Figure()
fig.grdview(
    grid=grid,
    perspective=[130, 30],
    frame=["xa", "ya", "WSnE"],
    projection="M15c",
    zsize="1.5c",
    # Set the surftype to "surface"
    surftype="s",
)
fig.show()

Out:

<IPython.core.display.Image object>

The default CPT is turbo and can be customized with the cmap parameter.

fig = pygmt.Figure()
fig.grdview(
    grid=grid,
    perspective=[130, 30],
    frame=["xa", "yaf", "WSnE"],
    projection="M15c",
    zsize="1.5c",
    surftype="s",
    # Set the CPT to "geo"
    cmap="geo",
)
fig.show()

Out:

<IPython.core.display.Image object>

The plane argument sets the elevation and color of a plane that provides a fill below the surface relief.

fig = pygmt.Figure()
fig.grdview(
    grid=grid,
    perspective=[130, 30],
    frame=["xa", "yaf", "WSnE"],
    projection="M15c",
    zsize="1.5c",
    surftype="s",
    cmap="geo",
    # Set the plane elevation to 1,000 meters and make the fill "gray"
    plane="1000+ggray",
)
fig.show()

Out:

<IPython.core.display.Image object>

The perspective azimuth can be changed to set the direction that is “up” in the figure.

fig = pygmt.Figure()
fig.grdview(
    grid=grid,
    # Set the azimuth to -130 (230) degrees and the elevation to 30 degrees
    perspective=[-130, 30],
    frame=["xa", "yaf", "WSnE"],
    projection="M15c",
    zsize="1.5c",
    surftype="s",
    cmap="geo",
    plane="1000+ggrey",
)
fig.show()

Out:

<IPython.core.display.Image object>

The contourpen parameter sets the pen used to draw contour lines on the surface.

fig = pygmt.Figure()
fig.grdview(
    grid=grid,
    perspective=[-130, 30],
    frame=["xaf", "yaf", "WSnE"],
    projection="M15c",
    zsize="1.5c",
    surftype="s",
    cmap="geo",
    plane="1000+ggrey",
    # Set the contour pen thickness to "0.5p"
    contourpen="0.5p",
)
fig.show()

Out:

<IPython.core.display.Image object>

pygmt.Figure.colorbar can be used to add a color bar to the figure. The cmap argument does not need to be passed again. To keep the color bar’s alignment similar to the figure, use True as the perspective argument.

fig = pygmt.Figure()
fig.grdview(
    grid=grid,
    perspective=[-130, 30],
    frame=["xaf", "yaf", "WSnE"],
    projection="M15c",
    zsize="1.5c",
    surftype="s",
    cmap="geo",
    plane="1000+ggrey",
    contourpen="0.1p",
)
fig.colorbar(perspective=True, frame=["a500", "x+lElevation", "y+lm"])
fig.show()

Out:

<IPython.core.display.Image object>