Tutorial 5 - 3-D Topography (Planetary / Antarctic maps) 🏔️

Tutorial 5 - 3-D Topography (Planetary / Antarctic maps) 🏔️#

In this tutorial, let’s use PyGMT to create 3-D perspective plots of Digital Elevation Models (DEM) over Mars (the planet) and Antarctica (the continent)!

Note

This tutorial is part of the AGU24 annual meeting GMT/PyGMT pre-conference workshop (PREWS9) Mastering Geospatial Visualizations with GMT/PyGMT

Website: https://www.generic-mapping-tools.org/agu24workshop
GitHub: GenericMappingTools/agu24workshop
Conference: https://agu.confex.com/agu/agu24/meetingapp.cgi/Session/226736

History

Authors: Wei Ji Leong and André Belém
Created: November-December 2024
Recommended versions: PyGMT v0.13.0 with GMT 6.5.0

References

https://www.generic-mapping-tools.org/remote-datasets/mars-relief.html
andrebelem/PlanetaryMaps
Neumann et al. [2003]

Fee free to play around with these code examples 🚀. In case you found any kind of error, just report it by opening an issue or provide a fix via a pull request. Please use the GMT forum to ask questions.

import pygmt
import rioxarray
import rioxarray.merge

1️⃣ Using `grdview` for 3-D Visualization#

The grdview method in PyGMT is a powerful tool for creating 3-D perspective views of gridded data. By adjusting azimuth and elevation parameters, you can change the viewpoint, helping you to highlight specific terrain features or data patterns. Let’s go through how these parameters affect the visualization.

Setting the Perspective: Azimuth and Elevation

Azimuth (azimuth): Controls the horizontal rotation of the view, ranging from 0° to 360°. Lower values (close to 0°) represent a viewpoint from the north, while 90° gives a view from the east, 180° from the south, and 270° from the west. Experimenting with azimuth can help showcase specific aspects of the data from different angles.
Elevation (elevation): Controls the vertical angle of the view, with 90° representing a top-down perspective and lower values providing more of a side view. Typically, values between 20° and 60° create a balanced 3-D effect.

Example: Using perspective=[150, 45] means azimuth=150 and elevation=45, which gives you a southeast view, tilted at a moderate angle to capture both horizontal and vertical details.

fig = pygmt.Figure()

fig.grdview(
    grid=da_olympus,
    cmap="SCM/managua+h0",
    region=[-151, -117, 12, 25, -5000, 23000],  # xmin, xmax, ymin, ymax, zmin, zmax
    projection="M12c",
    perspective=[150, 45],  # azimuth bearing, and elevation angle
    zsize="4c",  # vertical exaggeration
    surftype="s",  # surface plot
    shading=True,
    frame=["xaf", "yaf", "z5000+lmeters"],
)

fig.colorbar(frame=["a5000", "x+lElevation", "y+lm"], perspective=True, shading=True)
fig.show()

_images/14bcf9767aede971398cbc56fc2bd4e6aba04b49f61497063707a3ada60abe02.png

Note that there are other things we have configured such as:

zsize - vertical exaggeration as z-axis size, we used 4c here for 4 centimeters
surftype - using s will just create a regular surface
shading - set to True to get the default hillshading effect
frame - A proper 3-D map frame that consists of:
- automatic tick marks on x and y axis (e.g., xaf and yaf)
- z-axis tick marks every 5000m, plus a label (z5000+lLabel)

Tip

When choosing azimuth and elevation, always consider how the scene is illuminated. Azimuth angles that align with typical light directions (e.g., from the northwest) often provide the most natural and visually appealing shadows. Elevations between 20° and 45° typically create a good balance, highlighting terrain features without flattening or obscuring them. You can experiment with different combinations to best reveal the data’s structure.

3️⃣ Draping RGB image on 3-D topography#

We have our RGB imagery from Sentinel-2, and a DEM from REMA, and now we can learn how to render the color image on top of the 3-D topography! Once again, we will be using grdview, but pass in some extra arguments:

drapegrid: This is the image that will be overlaid on top of the relief grid.
surftype: We are plotting an RGB image with 3-bands, so we will use i for an image plot. It is also possible to append a number (e.g. 600) as the dots-per-unit resolution for the rasterization
zscale: vertical exaggeration factor, usually given as a fractional number.

Tip

When setting the zscale for vertical exaggeration, choose a value that balances clarity and realism. For subtle topographies, higher exaggeration (e.g., smaller zscale values) can emphasize elevation differences, making features more visible. However, for steep terrains, lower exaggeration helps maintain a natural appearance. You may use shading=True to add hillshading, which enhances the 3-D effect by simulating light and shadows, making terrain features easier to interpret. Experiment with both parameters to find the best combination for your dataset and visualization goals.

This is how the code will look like. We’ll also use pygmt.config to set PS_PAGE_COLOR (the background color) to an off-white color instead of the default black to better match the polar landscape.

fig = pygmt.Figure()
with pygmt.config(PS_PAGE_COLOR="#a9aba5"):
    fig.grdview(
        grid=dem_mosaic,  # DEM layer
        drapegrid=rgb_image,  # Sentinel-2 image layer
        surftype="i600",  # image draping with 600dpi resolution
        perspective=[170, 20],  # view azimuth and angle
        zscale="0.0005",  # vertical exaggeration
        shading=True,  # hillshading
        # frame="af",
    )
fig.show()

grdview [WARNING]: (w - x_min) must equal (NX + eps) * x_inc), where NX is an integer and |eps| <= 0.0001.

grdview [WARNING]: w reset from 250016 to 250000

grdview [WARNING]: (e - x_min) must equal (NX + eps) * x_inc), where NX is an integer and |eps| <= 0.0001.

grdview [WARNING]: e reset from 329984 to 330000

grdview [WARNING]: (s - y_min) must equal (NY + eps) * y_inc), where NY is an integer and |eps| <= 0.0001.

grdview [WARNING]: s reset from -1370000 to -1370016

grdview [WARNING]: (n - y_min) must equal (NY + eps) * y_inc), where NY is an integer and |eps| <= 0.0001.

grdview [WARNING]: n reset from -1300016 to -1300000

_images/0a6fc097b1fcaa99eff8632dab714787757564ed60e07f20825101c6a8e045ae.png

Tutorial 5 - 3-D Topography (Planetary / Antarctic maps) 🏔️

Contents

Tutorial 5 - 3-D Topography (Planetary / Antarctic maps) 🏔️#

0️⃣ Mars relief data#

0.1 2-D map view#

0.2 Zoomed in view#

1️⃣ Using `grdview` for 3-D Visualization#

2️⃣ Antarctic Digital Elevation Model#

2.1 Getting a DEM mosaic#

2.2 Getting RGB imagery#

3️⃣ Draping RGB image on 3-D topography#

Tutorial 5 - 3-D Topography (Planetary / Antarctic maps) 🏔️

Contents

Tutorial 5 - 3-D Topography (Planetary / Antarctic maps) 🏔️#

0️⃣ Mars relief data#

0.1 2-D map view#

0.2 Zoomed in view#

1️⃣ Using grdview for 3-D Visualization#

2️⃣ Antarctic Digital Elevation Model#

2.1 Getting a DEM mosaic#

2.2 Getting RGB imagery#

3️⃣ Draping RGB image on 3-D topography#

1️⃣ Using `grdview` for 3-D Visualization#