gmt-plotting

name	gmt_plotting
category	visualization
description	Use GMT / Generic Mapping Tools for publication-style geoscience maps, topography, coastlines, focal mechanisms, cross-sections, and shell-based GMT plotting. Choose this skill whenever the user explicitly says GMT.
keywords	GMT, map, seismicity, epicenter distribution, station distribution, fault, topography, terrain, terrain background, 地形底图, focal mechanism, coastline, contour, cross-section, travel time, earthquake catalog, basemap, coast, pscoast, meca, coupe, grdimage, gmt begin, gmt end, run_gmt, 地图, 震中分布, 台站分布, 地形, 地形底图, 震源机制, 海岸线
prefer_when	GMT, gmt, Generic Mapping Tools, 使用GMT, 用GMT, GMT绘制, gmt begin, gmt coast, run_gmt
related_skills	_gen_gmt_docs_zh, _gen_gmt_docs_6_5
workflow	gmt_terrain_map

GMT Map Plotting

Description

Use GMT (Generic Mapping Tools) to create professional seismology maps: epicenter distribution maps, station location maps, topographic maps, focal mechanism beachballs, cross-sections, etc.

Use this skill whenever the user explicitly asks for GMT, gmt, Generic Mapping Tools, gmt begin, gmt coast, run_gmt, or shell-based GMT plotting.

If the user explicitly requests 地形底图 or terrain background, always include a terrain grid (gmt grdimage) before coastlines and plot layers.

---

⚠️ Critical Rules

• Output a bash code block — the executor runs it directly, no Python wrapper needed
• Script must cd "${SAGE_OUTDIR}" at the top so output files land in the right place
• Script must use gmt begin <name> PNG ... gmt end (GMT6 modern mode)
• For tasks that need Python to prepare data first, call the pre-injected run_gmt(script_str, outname) from within Python code
• Requires GMT >= 6.0 installed on the system

Point Symbol Input Rule

For single-point symbols such as a red five-point star, epicenter, station marker, or label anchor, do not generate echo "lon lat" | gmt plot .... That syntax is valid GMT/Bash, but users often read it as printing text instead of plotting data. Prefer an explicit coordinate file or heredoc so the script clearly looks like drawing.

Preferred coordinate-file form:

cat > star_point.txt << 'EOF'
104 36
EOF
gmt plot star_point.txt -R${R} -J${J} -Sa1.5c -Gred -W0.5p,red

Preferred heredoc form:

gmt plot -R${R} -J${J} -Sa1.5c -Gred -W0.5p,red << 'EOF'
104 36
EOF

-Sa is GMT's star symbol. The input row is still plotting data: longitude latitude.

⚠️ TOPOGRAPHY RULES — READ FIRST

ALWAYS add terrain background when plotting geographic maps. The user almost always wants topography. Do NOT skip it.

Layer order (mandatory):

1. gmt grdimage — terrain background (FIRST, before everything)
2. gmt coast — coastlines / borders on top of terrain (NO -G fill when using grdimage)
3. gmt plot / gmt meca etc. — data on top
4. gmt colorbar — color scale

Robust topography snippet (copy-paste this every time):

# --- Download relief grid (skip 01m — too slow; chain 02m→05m) ---
# Note: always check the file was actually created (grdcut may exit 0 without output)
if gmt grdcut @earth_relief_02m -R${R} -Gtopo.grd 2>/dev/null && [ -f topo.grd ]; then
  echo "02m OK"
else
  gmt grdcut @earth_relief_05m -R${R} -Gtopo.grd
fi

# --- Verify grid exists before proceeding ---
if [ ! -f topo.grd ]; then
  echo "ERROR: topo.grd not created. Check region bounds: R=${R}" >&2
  exit 1
fi

# --- ALWAYS build the CPT explicitly — never rely on -Cetopo1 alone ---
Z_MIN=$(gmt grdinfo topo.grd -C 2>/dev/null | awk '{print $6}')
Z_MAX=$(gmt grdinfo topo.grd -C 2>/dev/null | awk '{print $7}')
if [ -z "${Z_MIN}" ] || [ -z "${Z_MAX}" ] || ! printf '%s\n' "${Z_MIN} ${Z_MAX}" | grep -Eq '^-?[0-9]+(\.[0-9]+)?\s+-?[0-9]+(\.[0-9]+)?$'; then
  echo "ERROR: invalid elevation range (${Z_MIN}, ${Z_MAX})" >&2
  exit 1
fi
gmt makecpt -Cgeo -T${Z_MIN}/${Z_MAX} -Z > topo.cpt
if [ ! -s topo.cpt ]; then
  echo "ERROR: topo.cpt empty" >&2
  exit 1
fi

# --- Render terrain with hillshading ---
gmt grdimage topo.grd -J${J} -R${R} -Ctopo.cpt -I+d

# --- Coast on top of terrain (NO -G fill — would hide terrain) ---
gmt coast -R${R} -J${J} -W0.6p,gray30 -N1/0.8p,gray50 -A500

# --- Elevation colorbar ---
gmt colorbar -DJBC+w7c/0.35c+o0/0.5c -Ctopo.cpt -Baf+l"Elevation"

Why explicit makecpt is required:

• ❌ -Cetopo1 — GMT looks for built-in CPT; if not found, falls back to grayscale → land appears solid black
• ✅ gmt makecpt -Cgeo -T${Z_MIN}/${Z_MAX} -Z > topo.cpt then -Ctopo.cpt — always works, correct colors
• -Cgeo is built into every GMT6 installation (land=tan/green/brown, ocean=blue/cyan)
• Use gmt grdinfo -C to get real min/max so the color range exactly fits the data

Other common mistakes:

• ❌ gmt coast -Gtan / -Gwhite / -G<ANY color> — SOLID fill completely hides terrain. NEVER use -G with grdimage.
• ❌ gmt basemap before gmt grdimage — basemap frame gets buried under terrain
• ❌ Missing -I+d on grdimage — flat, washed-out colors without hillshading
• ❌ Hard-coding -T-6000/6000 when region is all land — white/grey scale, looks wrong
• ❌ gmt coast -Slightblue alone without grdimage — ocean will be colored but land has no terrain

---

Core Pattern

Pure bash output (preferred for all GMT-only tasks)

Output a ```bash code block. The engine's execute_bash() runs it with SAGE_OUTDIR set to the temp work directory.

#!/bin/bash
cd "${SAGE_OUTDIR}"      # all output files go here automatically

R="73/136/18/54"        # west/east/south/north
J="M15c"                # Mercator, 15 cm wide

# ... your GMT commands ...

gmt begin mymap PNG
  gmt grdimage topo.grd -J${J} -R${R} -Ctopo.cpt -I+d
  gmt coast -R${R} -J${J} -W0.6p,gray30 -N1/0.8p,gray50 -A500 -Baf
gmt end

Normal bash syntax throughout — ${VAR}, $(cmd), awk '{print $6}' — no escaping needed.

Mixed Python + GMT (when CSV data prep is required)

Call the pre-injected run_gmt(script_str, outname) from Python code:

import numpy as np, os

pts_file = os.path.join(os.environ.get('SAGE_OUTDIR', '/tmp'), 'points.txt')
np.savetxt(pts_file, np.column_stack([lon, lat]), fmt='%.6f')

R = f"{lon.min()-1:.2f}/{lon.max()+1:.2f}/{lat.min()-1:.2f}/{lat.max()+1:.2f}"
J = "M15c"

# Inside f-string: Python vars → {R}, bash vars → ${{Z_MIN}}, awk → {{print $6}}
script = f"""
cd "${{SAGE_OUTDIR}}"
gmt begin mymap PNG
  gmt plot {pts_file} -R{R} -J{J} -Sc0.2c -Gred -W0.3p,black -Baf
gmt end
"""
run_gmt(script, outname="mymap")

---

Legend / 图例

Use gmt legend to draw a boxed legend after all data layers. Never use gmt basemap for legend placement.

Quick rule of thumb

Want to plot	Legend spec line
Circle point	S 0.3c c <fill> <pen> Label
Triangle (station)	S 0.3c t <fill> <pen> Label
Square	S 0.3c s <fill> <pen> Label
Line	S 0.5c - - <pen> Label
Dashed line	S 0.5c - - <pen>,- Label

Legend spec file format — all directives

H <fontsize> [<font>] <heading>       ← bold title inside legend box
D [gap] <pen>                          ← horizontal divider line (gap optional e.g. 0.1c)
N <ncols>                              ← number of columns for following S lines
S [dx1] <symbol> <size> <fill> <pen> [dx2] <label>
G <gap>                                ← extra vertical space (e.g. 0.1c)
T <text>                               ← plain text line
L <size> <justification> <text>        ← left/center/right-justified text (L/C/R)

S line field widths:

S  dx1  symbol  size  fill     pen        [dx2]  label
S  0.3c   c     0.2c   red   0.5p,black         Earthquake

• dx1: horizontal offset of symbol from left edge of box (default 0.2c is fine)
• symbol letters: c circle, t triangle, s square, d diamond, i inverted-triangle, - line, f fault, v vector
• fill and pen must match the gmt plot command exactly; use - for "none"
• dx2 (optional): gap between symbol and label; omit to use default

Placement -D options

-DjBR+w5c+o0.2c/0.2c    # bottom-right, box width 5c, 0.2c offset from edge
-DjBL+w5c+o0.2c/0.2c    # bottom-left
-DjTR+w5c+o0.2c/0.2c    # top-right
-DjBC+w8c+o0/0.5c        # bottom-center (e.g. above the colorbar)
-DjTL+w5c+o0.2c/0.2c    # top-left

Box background -F options

Always use -F on terrain maps so the legend is readable:

-F+p0.8p,black+gwhite          # white background, black border
-F+p0.8p,gray30+gwhite@20      # semi-transparent white
-F+p0.5p,gray50+glightgray     # gray background

---

Example 1 — Station-only legend

gmt plot stations.txt -St0.3c -Gred -W0.5p,black

cat > legend.txt << 'EOF'
H 11 Legend
D 0.1c 0.5p
S 0.2c t 0.3c red 0.5p,black 0.3c Seismic station
EOF

gmt legend legend.txt -DjBR+w4.5c+o0.2c/0.2c -F+p0.8p,black+gwhite

---

Example 2 — Multi-layer legend: stations + earthquakes + faults

gmt plot faults.txt    -W1.2p,firebrick
gmt plot earthquakes.txt -Sc0.2c -Gblue   -W0.3p,black
gmt plot stations.txt    -St0.3c -Gred    -W0.5p,black

cat > legend.txt << 'EOF'
H 11 Legend
D 0.1c 0.5p
S 0.5c - - - 1.2p,firebrick 0.3c Fault
S 0.2c c 0.2c blue 0.3p,black 0.3c Earthquake
S 0.2c t 0.3c red 0.5p,black  0.3c Station
EOF

gmt legend legend.txt -DjBR+w5.5c+o0.2c/0.2c -F+p0.8p,black+gwhite

---

Example 3 — Magnitude-scaled earthquake legend

When earthquake circles are scaled by magnitude (-Sc with variable size), use several representative rows to communicate the scale:

# Variable-size circles: size = mag * 0.08c
awk '{print $1,$2,$3*0.08"c"}' catalog.txt | gmt plot -Sc -Gblue@40 -W0.3p,gray30

cat > legend.txt << 'EOF'
H 11 Magnitude
D 0.1c 0.5p
S 0.25c c 0.24c blue@40 0.3p,gray30 0.4c M 3
S 0.25c c 0.40c blue@40 0.3p,gray30 0.4c M 5
S 0.25c c 0.56c blue@40 0.3p,gray30 0.4c M 7
EOF

gmt legend legend.txt -DjBR+w3.8c+o0.2c/0.2c -F+p0.8p,black+gwhite

---

Example 4 — Focal mechanism legend (gmt meca)

gmt meca symbols use type letter matching the -S flag: a = Aki-Richards, d = double-couple.

gmt meca focal.txt -Sa0.5c -Gred -W0.5p,black

cat > legend.txt << 'EOF'
H 11 Focal mechanism
D 0.1c 0.5p
S 0.3c a 0.5c red 0.5p,black 0.3c Focal mechanism (Mw≥4)
EOF

gmt legend legend.txt -DjTR+w5c+o0.2c/0.2c -F+p0.8p,black+gwhite

---

Example 5 — Depth-colored earthquakes with separate legend + colorbar

# Color by depth using a CPT
gmt makecpt -Cjet -T0/100/10 > depth.cpt
awk '{print $1,$2,$3}' catalog.txt | gmt plot -Sc0.15c -Cdepth.cpt -W0.2p,gray30

# Colorbar for depth
gmt colorbar -DJBC+w7c/0.35c+o0/0.5c -Cdepth.cpt -Baf+l"Depth (km)"

# Separate legend for symbol type
cat > legend.txt << 'EOF'
H 11 Data
D 0.1c 0.5p
S 0.2c c 0.15c gray50 0.2p,gray30 0.3c Earthquake (color = depth)
EOF

gmt legend legend.txt -DjTR+w5c+o0.2c/0.2c -F+p0.8p,black+gwhite

---

f-string safe legend patterns (Python-generated scripts)

When the GMT script is built inside a Python f-string, use printf + pipe or a heredoc with 'EOF' (single-quoted = no variable expansion). Both are safe:

gmt_script = f"""
gmt begin mymap PNG

  # ... terrain, coast, plot layers ...

  gmt plot {pts_file} -R{R} -J{J} -Sc0.2c -Gblue -W0.3p,black

  printf 'H 11 Legend\\nD 0.1c 0.5p\\nS 0.2c c 0.2c blue 0.3p,black 0.3c Data point\\n' \\
    | gmt legend -DjBR+w4.5c+o0.2c/0.2c -F+p0.8p,black+gwhite

gmt end
"""