Tutorial:Plotting features on a 3-D surface

examples/tutorials/advanced/features_on_3d_surface.py

@@ -0,0 +1,152 @@
+"""
+Plotting features on a 3-D surface
+==================================
+
+In addition to draping a dataset (grid or image) on top of a topographic surface,
+you may want to add additional features like coastlines, symbols, and text
+annotations. This tutorial shows how to use :meth:`pygmt.Figure.coast`,
+:meth:`pygmt.Figure.plot3d`, and :meth:`pygmt.Figure.text` to add these features
+on a 3-D surface created by :meth:`pygmt.Figure.grdview`.
+
+This tutorial builds a 3-D map with additional features in four steps:
+
+1. Creating a 3-D surface
+2. Adding coastlines on a 3-D surface
+3. Adding symbols on a 3-D surface
+4. Adding text annotations on a 3-D surface
+"""
+
+# %%
+
+import pandas as pd
+import pygmt
+from pygmt.params import Axis, Frame
+
+# %%
+# 1. Creating a 3-D surface
+# -------------------------
+#
+# In the first step, we create a 3-D topographic map using :meth:`pygmt.Figure.grdview`.
+# We use a region around Taiwan to demonstrate adding features on a 3-D surface.
+
+# Define the study area in degrees East or North
+region_2d = [119, 123, 21, 26]  # [lon_min, lon_max, lat_min, lat_max]
+
+# Download elevation grid for the study region with a resolution of 5 arc-minutes.
+grd_relief = pygmt.datasets.load_earth_relief(resolution="05m", region=region_2d)
+
+# Determine the 3-D region from the minimum and maximum values of the relief grid
+region_3d = [*region_2d, grd_relief.min().to_numpy(), grd_relief.max().to_numpy()]
+
+fig = pygmt.Figure()
+
+# Set up a colormap for topography and bathymetry
+pygmt.makecpt(cmap="gmt/globe", series=[-6000, 3000])
+
+# Create a 3-D surface
+fig.grdview(
+    projection="M12c",  # Mercator projection with a width of 12 cm
+    region=region_3d,
+    grid=grd_relief,
+    cmap=True,
+    surftype="surface",
+    shading="+a0/270+ne0.6",
+    perspective=[157.5, 30],  # Azimuth and elevation for the 3-D plot
+    zsize="1.5c",
+    facade_fill="darkgray",
+    frame=Frame(axes="WSnE", axis=Axis(annot=True, tick=True)),
+)
+
+# Add a colorbar
+fig.colorbar(perspective=True, annot=1000, tick=500, label="Elevation", unit="m")
+
+fig.show()
+
+# %%
+# 2. Adding coastlines on a 3-D surface
+# -------------------------------------
+#
+# Next, we add coastlines using :meth:`pygmt.Figure.coast` with a matching
+# ``perspective`` setting. Here we set the z-level to 0 so coastlines are drawn
+# at sea level.
+
+# Add coastlines on top of the 3-D surface
+# Use an explicit perspective to match grdview (azimuth=157.5, elevation=30)
+# and set the z-level to 0 so the coastlines are drawn at sea level.
+fig.coast(
+    perspective=[157.5, 30, 0],
+    resolution="full",
+    shorelines="1/2p,black",
+)
+
+fig.show()
+
+# %%
+# 3. Adding symbols on a 3-D surface
+# ----------------------------------
+#
+# In the third step, we add star symbols on top of the same 3-D map. To plot
+# symbols on a 3-D surface, use :meth:`pygmt.Figure.plot3d`. The z-coordinate should be
+# set to a value at or above the maximum elevation to ensure the symbols are visible.
+# Note that 3-D rendering in GMT/PyGMT uses a painter's algorithm (depth sorting)
+# rather than true 3-D occlusion. From some viewpoints, symbols that should be
+# hidden behind terrain may still appear visible.
+
+# Sample point data: five coastal cities around Taiwan
+cities = pd.DataFrame(
+    {
+        "longitude": [
+            121.74,
+            121.61,
+            121.14,
+            120.30,
+            120.53,
+        ],  # Keelung, Hualien, Taitung, Kaohsiung, Taichung Port
+        "latitude": [25.13, 23.99, 22.76, 22.63, 24.27],
+        "name": ["Keelung", "Hualien", "Taitung", "Kaohsiung", "Taichung Port"],
+    }
+)
+
+# Use one common z-level so all stars share the same shape and size.
+max_elevation = float(grd_relief.max().to_numpy())
+z_stars = [max_elevation + 1500] * len(cities.index)
+
+# Add five identical star symbols on top of the 3-D surface
+fig.plot3d(
+    x=cities.longitude,
+    y=cities.latitude,
+    z=z_stars,
+    style="a0.65c",
+    fill="gold",
+    pen="0.8p,black",
+    perspective=True,
+    zsize="1.5c",  # Use the same zsize as grdview
+    no_clip=True,
+)
+
+fig.show()
+
+# %%
+# 4. Adding text annotations on a 3-D surface
+# -------------------------------------------
+#
+# In the final step, we add text labels to the same 3-D map. To add text
+# annotations on a 3-D surface, use :meth:`pygmt.Figure.text` with
+# ``perspective=True``. Note that the current implementation of ``text`` does not
+# support a ``z`` parameter for controlling the vertical position of text labels.
+# The text will be placed at the base of the 3-D plot (z=0).
+
+# Add text labels for cities
+# Note: text is placed at z=0 (base level) since z parameter is not yet supported
+fig.text(
+    x=cities.longitude,
+    y=cities.latitude,
+    text=cities.name,
+    perspective=True,
+    font="11p,Times-Bold,red",
+    no_clip=True,  # Prevent text from being clipped at the frame boundaries
+)
+
+fig.show()
+
+# sphinx_gallery_thumbnail_number = 3