| 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):
gmt grdimage — terrain background (FIRST, before everything)gmt coast — coastlines / borders on top of terrain (NO -G fill when using grdimage)gmt plot / gmt meca etc. — data on topgmt colorbar — color scale
Robust topography snippet (copy-paste this every time):
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
if [ ! -f topo.grd ]; then
echo "ERROR: topo.grd not created. Check region bounds: R=${R}" >&2
exit 1
fi
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
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
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}"
R="73/136/18/54"
J="M15c"
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"
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
-DjBL+w5c+o0.2c/0.2c
-DjTR+w5c+o0.2c/0.2c
-DjBC+w8c+o0/0.5c
-DjTL+w5c+o0.2c/0.2c
Box background -F options
Always use -F on terrain maps so the legend is readable:
-F+p0.8p,black+gwhite
-F+p0.8p,gray30+gwhite@20
-F+p0.5p,gray50+glightgray
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:
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
gmt makecpt -Cjet -T0/100/10 > depth.cpt
awk '{print $1,$2,$3}' catalog.txt | gmt plot -Sc0.15c -Cdepth.cpt -W0.2p,gray30
gmt colorbar -DJBC+w7c/0.35c+o0/0.5c -Cdepth.cpt -Baf+l"Depth (km)"
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
"""
Or, using a single-quoted heredoc (no ${{}} escaping needed):
gmt_script = f"""
gmt begin mymap PNG
gmt plot {pts_file} -R{R} -J{J} -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
gmt end
"""
Rule: inside an f-string, use 'EOF' (single-quoted) for the legend heredoc — this prevents the shell from expanding any $ inside the legend spec, and you don't need to escape anything.
Complete working example: terrain map with stations + earthquakes + legend
#!/bin/bash
cd "${SAGE_OUTDIR}"
STATIONS="/path/to/stations.txt"
CATALOG="/path/to/catalog.txt"
R="95/115/25/42"
J="M15c"
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
Z_MIN=$(gmt grdinfo topo.grd -C | awk '{print $6}')
Z_MAX=$(gmt grdinfo topo.grd -C | awk '{print $7}')
gmt makecpt -Cgeo -T${Z_MIN}/${Z_MAX} -Z > topo.cpt
gmt begin seismic_map 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 \
-Bxaf+l"Longitude (E)" -Byaf+l"Latitude (N)" -BWSne+t"Seismicity Map"
gmt plot ${CATALOG} -R${R} -J${J} -Sc0.15c -Gblue@40 -W0.3p,gray20
gmt plot ${STATIONS} -R${R} -J${J} -St0.35c -Gred -W0.6p,black
gmt colorbar -DJBC+w7c/0.35c+o0/0.5c -Ctopo.cpt -Baf+l"Elevation (m)"
cat > legend.txt << 'EOF'
H 12 Legend
D 0.1c 0.5p
S 0.2c c 0.15c blue@40 0.3p,gray20 0.3c Earthquake
S 0.2c t 0.35c red 0.6p,black 0.3c Station
EOF
gmt legend legend.txt -DjBR+w5c+o0.3c/0.3c -F+p0.8p,black+gwhite
gmt end
Common legend mistakes
- ❌ Placing
gmt legend before gmt grdimage or before data layers — legend gets buried
- ❌ Missing
-F — legend text invisible on top of terrain
- ❌ Symbol letter in spec doesn't match
gmt plot -S flag (e.g. plotting -St but writing c in spec)
- ❌ Fill/pen in legend spec doesn't match the plot command — use copy-paste to avoid mismatch
- ❌ Using
<< EOF (unquoted) inside f-string — shell expands $ in spec lines, causing errors
- ❌ Legend box too narrow (
+w too small) — label text gets clipped; use +w5c or wider as default
GMT Basemap Skill
Objective
Add and configure a basemap in GMT scripts to control map region, projection, axes, ticks, annotations, title, and overall layout.
When to Use
Use this skill when the user needs to:
- Define map boundaries (
-R)
- Set projection (
-J)
- Add axes, ticks, gridlines (
-B)
- Add map title or labels
- Standardize GMT plotting structure
- Prepare a base layer before plotting data (e.g., earthquakes, stations, faults)
Core Command
gmt basemap -R<west>/<east>/<south>/<north> -J<projection> -B<axis_settings>
Typical Example
gmt basemap -R100/110/20/30 -JM12c -Bxa2f1 -Bya2f1 -BWSen
Meaning
-R100/110/20/30 Region: lon 100–110, lat 20–30
-JM12c Mercator projection, 12 cm width
-Bxa2f1 X-axis: major ticks every 2°, minor every 1°
-Bya2f1 Y-axis: major ticks every 2°, minor every 1°
-BWSen Draw west & south with annotations; east & north frame only
Standard Script Template
#!/bin/bash
gmt begin map png
gmt basemap \
-R100/110/20/30 \
-JM12c \
-Bxa2f1+l"Longitude" \
-Bya2f1+l"Latitude" \
-BWSen+t"Study Area"
gmt coast \
-W0.5p,black \
-Df \
-N1/0.5p,gray \
-A1000
gmt end show
Recommended Workflow
Always plot basemap first, then overlay data layers:
gmt begin seismic_map pdf
gmt basemap -R100/110/20/30 -JM15c -Bxa2f1 -Bya2f1 -BWSen+t"Seismicity Map"
gmt coast -W0.5p -Df -N1/0.5p,gray -A1000
gmt plot faults.txt -W1p,red
gmt plot stations.txt -St0.25c -Gblue -W0.3p,black
gmt plot earthquakes.txt -Sc0.08c -Gred -W0.1p,black
gmt end show
Key Parameters
1. Region (-R)
-Rwest/east/south/north
Example:
-R95/110/20/35
Cartesian:
-R0/100/0/50
2. Projection (-J)
Common options:
-JM12c
-JX12c/8c
-JL...
Typical usage:
- Maps →
-JM
- Profiles / sections →
-JX
3. Axes & Frame (-B)
Example:
-Bxa1f0.5 -Bya1f0.5 -BWSen
Meaning:
a1 major ticks = 1
f0.5 minor ticks = 0.5
W S draw + annotate west & south
e n draw east & north frame only
Add labels:
-Bxa2f1+l"Longitude"
-Bya2f1+l"Latitude"
Add title:
-BWSen+t"Earthquake Distribution"
Geographic Map Template
gmt basemap \
-R${lon_min}/${lon_max}/${lat_min}/${lat_max} \
-JM${width}c \
-Bxa${x_major}f${x_minor}+l"Longitude" \
-Bya${y_major}f${y_minor}+l"Latitude" \
-BWSen+t"${title}"
Cross-section (Depth Profile) Template
gmt basemap \
-R0/100/0/40 \
-JX12c/-6c \
-Bxa20f10+l"Distance (km)" \
-Bya10f5+l"Depth (km)" \
-BWSen+t"Depth Section"
Note:
-JX12c/-6c
Negative height flips the Y-axis (depth increasing downward).
Complete Working Example: Topographic Basemap
Use this template directly. It handles all common pitfalls: terrain download, CPT creation, hillshading, coastlines, and colorbar.
#!/usr/bin/env bash
region=95/110/20/35
projection=M15c
if gmt grdcut @earth_relief_01m -R${region} -Gtopo.grd 2>/dev/null && [ -f topo.grd ]; then
echo "✓ Terrain resolution 01m downloaded"
elif gmt grdcut @earth_relief_02m -R${region} -Gtopo.grd 2>/dev/null && [ -f topo.grd ]; then
echo "✓ Terrain resolution 02m downloaded"
else
echo "⚠ Using fallback resolution 05m"
gmt grdcut @earth_relief_05m -R${region} -Gtopo.grd
fi
if [ ! -f topo.grd ]; then
echo "ERROR: topo.grd not created. Check region: R=${region}" >&2
exit 1
fi
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
echo "Elevation range: ${Z_MIN} — ${Z_MAX} m"
gmt makecpt -Cgeo -T${Z_MIN}/${Z_MAX}/100 > topo.cpt
if [ ! -s topo.cpt ]; then
echo "ERROR: topo.cpt empty" >&2
exit 1
fi
gmt begin topo_map png
gmt grdimage topo.grd -R${region} -J${projection} \
-Ctopo.cpt -I+d
gmt coast -R${region} -J${projection} \
-W0.6p,gray30 \
-N1/0.8p,gray50 \
-A500
gmt basemap -R${region} -J${projection} \
-Bxa2f1+l"Longitude (°E)" \
-Bya2f1+l"Latitude (°N)" \
-BWSen+t"Topographic Basemap"
gmt colorbar -DJBC+w7c/0.35c+o0/0.5c \
-Ctopo.cpt \
-Baf+l"Elevation (m)"
gmt end show
Key Features
| Feature | How It Works |
|---|
| Fallback resolution | Tries 01m → 02m → 05m to ensure download succeeds |
| Validation | Checks topo.grd exists and elevation range is valid before continuing |
| Explicit CPT | Uses gmt makecpt -Cgeo -T${Z_MIN}/${Z_MAX} instead of relying on built-in colormap names |
| Hillshading | -I+d adds relief shading for 3D effect |
| No terrain override | gmt coast has NO -G fill; coastlines draw cleanly over terrain |
| Colorbar | Specifies -D position and -C topo.cpt explicitly |
| Frame | gmt basemap adds axes, ticks, labels, and title after terrain is rendered |
Common Customizations
Different region:
region=100/115/25/40
region=70/100/10/30
Different projection:
projection=M20c
projection=L100/25/20/30/12c
Add data points (stations, earthquakes):
gmt plot stations.txt -Sc0.2c -Gblue -W0.3p,black
gmt plot earthquakes.txt -Sc0.08c -Gred -W0.1p,darkred
Color and legend example:
gmt plot stations.txt -Sc0.25c -Gblue -W0.3p,black
cat > legend.txt << EOF
S 0.25c c blue 0.3p,black Seismic station
EOF
gmt legend legend.txt -DjTR+w4c/1.5c+o0.2c/0.2c -F+p0.8p,black+gwhite
Different color scheme:
gmt makecpt -Cturbo -T${Z_MIN}/${Z_MAX} -Z > topo.cpt
gmt makecpt -Cbilbao -T${Z_MIN}/${Z_MAX} -Z > topo.cpt
gmt makecpt -Chot -T${Z_MIN}/${Z_MAX} -Z > topo.cpt
Minimal Working Example
gmt begin test_basemap png
gmt basemap \
-R100/110/20/30 \
-JM12c \
-Bxa2f1+l"Longitude" \
-Bya2f1+l"Latitude" \
-BWSen+t"GMT Basemap Example"
gmt coast -W0.5p -Df -N1/0.5p,gray -A1000
gmt end show
Common Mistakes
1. Missing axis definition
❌
-B
✔
-Bxa2f1 -Bya2f1 -BWSen
2. Region mismatch
If data fall outside -R, nothing is plotted.
Check:
gmt info data.txt
3. Wrong depth direction
Fix with:
-JX12c/-6c
4. Wrong projection choice
- Maps →
-JM
- Profiles →
-JX
Output Requirement (for agents)
When the user asks to “add basemap”:
Note
If you want, I can extend this into a full GMT seismic plotting skill set (basemap + coast + fault + beachball + colorbar + legend), which is closer to your actual workflow.
Example Collection
0. Plot Points from CSV with Non-Standard Column Names (e.g. lon1/lat1)
When your CSV has coordinate columns named lon1, lat1 (not lon, lat), extract
them first in Python, write a temp GMT-format file, then call run_gmt().
import pandas as pd, numpy as np, os
df = pd.read_csv("/path/to/data.csv")
print("columns:", df.columns.tolist()); print(df.head(3))
lon_col = 'lon1'
lat_col = 'lat1'
lon = df[lon_col].values
lat = df[lat_col].values
assert -180 <= lon.min() and lon.max() <= 180, \
f"lon out of range {lon.min():.2f}~{lon.max():.2f} — column '{lon_col}' is not longitude"
assert -90 <= lat.min() and lat.max() <= 90, \
f"lat out of range {lat.min():.2f}~{lat.max():.2f} — column '{lat_col}' is not latitude!"
print(f"lon: {lon.min():.4f}~{lon.max():.4f} lat: {lat.min():.4f}~{lat.max():.4f}")
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"
print(f"Region: {R}")
script = f"""
gmt begin location_map PNG
if gmt grdcut @earth_relief_02m -R{R} -Gtopo.grd 2>/dev/null && [ -f topo.grd ]; then
echo "02m OK"
elif gmt grdcut @earth_relief_05m -R{R} -Gtopo.grd 2>/dev/null && [ -f topo.grd ]; then
echo "05m OK"
else
echo "Terrain download failed" >&2; exit 1
fi
Z_MIN=$(gmt grdinfo topo.grd -C | awk '{{print $6}}')
Z_MAX=$(gmt grdinfo topo.grd -C | awk '{{print $7}}')
gmt makecpt -Cgeo -T${{Z_MIN}}/${{Z_MAX}} -Z > topo.cpt
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 -BWSne+t"Location Map"
gmt plot {pts_file} -R{R} -J{J} -Sc0.2c -Gred -W0.3p,black
gmt colorbar -DJBC+w7c/0.35c -Ctopo.cpt -Baf+l"Elevation (m)"
gmt end
"""
run_gmt(script, outname="location_map", title="Location map from CSV")
1. Epicenter Distribution Map (Most Common) — with topography
#!/bin/bash
cd "${SAGE_OUTDIR}"
R="70/140/15/55"
J="M15c"
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
Z_MIN=$(gmt grdinfo topo.grd -C | awk '{print $6}')
Z_MAX=$(gmt grdinfo topo.grd -C | awk '{print $7}')
gmt makecpt -Cgeo -T${Z_MIN}/${Z_MAX} -Z > topo.cpt
gmt begin seismicity PNG
gmt grdimage topo.grd -J${J} -R${R} -Ctopo.cpt -I+d
gmt coast -R${R} -J${J} -W0.5p,gray40 -N1/0.8p,gray60 -A500 \
-Baf -BWSne+t"Seismicity Map"
gmt end
2. Station Location Map
#!/bin/bash
cd "${SAGE_OUTDIR}"
gmt begin station_map PNG
gmt basemap -R100/115/25/42 -JM14c -Baf -BWSne+t"Station Map"
gmt coast -Gwheat -Slightblue -W0.5p -N1/0.8p,gray
gmt plot stations.txt -i0,1 -St0.5c -Gred -W1p,black
gmt text stations.txt -i0,1,2 -F+f7p,Helvetica,black+jBL -D0.1c/0.1c
gmt end
3. Focal Mechanism Beachballs
#!/bin/bash
cd "${SAGE_OUTDIR}"
gmt begin focal_map PNG
gmt basemap -R95/105/28/38 -JM12c -Baf -BWSne+t"Focal Mechanisms"
gmt coast -Gtan -Slightblue -W0.5p -N1/0.5p
cat << 'EOF' | gmt meca -Sa1c -Gred -W0.5p,black
100.5 33.2 15 120 60 -90 5.8 1 1
102.1 30.5 20 45 75 30 4.9 1 1
103.4 31.8 10 90 45 180 5.2 1 1
EOF
gmt end
4. Topographic Map (ETOPO / SRTM) — Robust Version
#!/bin/bash
cd "${SAGE_OUTDIR}"
R="100/115/25/40"
J="M14c"
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
Z_MIN=$(gmt grdinfo topo.grd -C | awk '{print $6}')
Z_MAX=$(gmt grdinfo topo.grd -C | awk '{print $7}')
gmt makecpt -Cgeo -T${Z_MIN}/${Z_MAX} -Z > topo.cpt
gmt begin topo_map 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 \
-Bxaf+l"Longitude" -Byaf+l"Latitude" -BWSne+t"Topographic Map"
gmt colorbar -DJBC+w8c/0.4c+e -Ctopo.cpt -Baf+l"Elevation (m)"
gmt end
5. Cross-Section
#!/bin/bash
cd "${SAGE_OUTDIR}"
gmt begin cross_section PNG
gmt basemap -R0/1200/0/100 -JX15c/-8c \
-Bxaf+l"Distance (km)" -Byaf+l"Depth (km)" -BWSne+t"Cross Section A-B"
printf '0 35\n1200 35\n' | gmt plot -W1p,red,-
gmt end
Complete Chained Example: CSV Catalog to Epicenter Map
import os
catalog_file = "/data/catalog/eq_catalog.csv"
gmt_script = """
gmt begin seismicity_final PNG
gmt basemap -R70/140/15/55 -JM16c -Baf+g245/245/240 -BWSne+t"2020-2024 Seismicity"
gmt coast -Gtan -Slightblue -W0.5p,gray50 -N1/0.8p,gray70 -A500
# Small earthquakes (M<4): small gray dots
awk -F',' '$4<4 {print $1,$2}' """ + catalog_file + """ | gmt plot -Sc0.08c -Ggray60 -t50
# Moderate earthquakes (4≤M<5): orange dots
awk -F',' '$4>=4 && $4<5 {print $1,$2}' """ + catalog_file + """ | gmt plot -Sc0.18c -Gorange -W0.3p,black
# Large earthquakes (M≥5): red dots
awk -F',' '$4>=5 {print $1,$2,($4*0.15)"c"}' """ + catalog_file + """ | gmt plot -Sc -Gred -W0.5p,black
gmt colorbar -DJBC+w6c/0.35c -Baf+l"Depth (km)"
gmt end
"""
run_gmt(gmt_script, outname="seismicity_final", title="Epicenter distribution map")
print("Earthquake count: " + str(sum(1 for _ in open(catalog_file)) - 1))
Complete Chained Example: Read CSV Data and Plot Locations with GMT + Topography
This is the canonical template for any "read CSV → plot locations on map" request.
Always follow this pattern when the user provides a data file with lon/lat columns.
import pandas as pd
import numpy as np
import os, tempfile
data_file = "/path/to/your/data.csv"
df = pd.read_csv(data_file)
print("Columns:", list(df.columns))
print("Shape:", df.shape)
lon_col = next((c for c in df.columns if 'lon' in c.lower()), df.columns[0])
lat_col = next((c for c in df.columns if 'lat' in c.lower()), df.columns[1])
print(f"Using lon={lon_col}, lat={lat_col}")
lons = pd.concat([df[lon_col]] + ([df['lon2']] if 'lon2' in df.columns else [])).dropna()
lats = pd.concat([df[lat_col]] + ([df['lat2']] if 'lat2' in df.columns else [])).dropna()
pad = 0.05
lon_min = lons.min(); lon_max = lons.max()
lat_min = lats.min(); lat_max = lats.max()
dlon = max(lon_max - lon_min, 1.0); dlat = max(lat_max - lat_min, 1.0)
R = f"{lon_min - pad*dlon:.2f}/{lon_max + pad*dlon:.2f}/{lat_min - pad*dlat:.2f}/{lat_max + pad*dlat:.2f}"
J = "M15c"
print(f"Region: {R}")
tmp = tempfile.mktemp(suffix='.txt')
points = pd.DataFrame({'lon': df[lon_col], 'lat': df[lat_col]}).drop_duplicates()
points.to_csv(tmp, sep=' ', index=False, header=False)
print(f"Plotting {len(points)} unique locations")
gmt_script = f"""
gmt begin locations_map PNG
# ── Topography (try 01m → 02m → 05m) ──────────────────────────────────
if gmt grdcut @earth_relief_01m -R{R} -Gtopo.grd 2>/dev/null; then
echo "terrain 01m loaded"
elif gmt grdcut @earth_relief_02m -R{R} -Gtopo.grd 2>/dev/null; then
echo "terrain 02m loaded"
else
gmt grdcut @earth_relief_05m -R{R} -Gtopo.grd
echo "terrain 05m loaded"
fi
# Render terrain with hillshading
Z_MIN=$(gmt grdinfo topo.grd -C | awk '{{print $6}}')
Z_MAX=$(gmt grdinfo topo.grd -C | awk '{{print $7}}')
gmt makecpt -Cgeo -T${{Z_MIN}}/${{Z_MAX}} -Z > topo.cpt
gmt grdimage topo.grd -J{J} -R{R} -Ctopo.cpt -I+d
# ── Coastlines + borders on top (no -G fill — preserves terrain) ──────
gmt coast -R{R} -J{J} \\
-W0.6p,gray30 -N1/0.8p,gray50 -A500 \\
-Bxaf+l"Longitude (°E)" -Byaf+l"Latitude (°N)" \\
-BWSne+t"Seismic Station Locations"
# ── Plot data points ───────────────────────────────────────────────────
gmt plot {tmp} -R{R} -J{J} -St0.35c -Gred -W0.6p,black
# ── Elevation color bar ────────────────────────────────────────────────
gmt colorbar -DJBC+w7c/0.35c+o0/0.5c -Baf+l"Elevation (m)" -F+g255/255/255@40
gmt end
"""
run_gmt(gmt_script, outname="seismic_locations_topo", title="Seismic locations on topographic map")
os.remove(tmp)
print("Map with topography saved.")
Complete Chained Example: Ray Path Visualization with Topography
import pandas as pd
import os, tempfile
data_file = "/path/to/your/data.csv"
df = pd.read_csv(data_file)
all_lons = list(df['lon1']) + list(df['lon2'])
all_lats = list(df['lat1']) + list(df['lat2'])
pad = 0.5
R = f"{min(all_lons)-pad:.2f}/{max(all_lons)+pad:.2f}/{min(all_lats)-pad:.2f}/{max(all_lats)+pad:.2f}"
J = "M15c"
rays_file = tempfile.mktemp(suffix='_rays.txt')
pts_file = tempfile.mktemp(suffix='_pts.txt')
ray_ids = df['ray_id'].unique()[:500]
with open(rays_file, 'w') as fr, open(pts_file, 'w') as fp:
seen = set()
for rid in ray_ids:
sub = df[df['ray_id'] == rid].iloc[0]
fr.write(f">\n{sub.lon1} {sub.lat1}\n{sub.lon2} {sub.lat2}\n")
for pt in [(sub.lon1, sub.lat1), (sub.lon2, sub.lat2)]:
if pt not in seen:
fp.write(f"{pt[0]} {pt[1]}\n")
seen.add(pt)
gmt_script = f"""
gmt begin ray_paths PNG
# ── Topography ──────────────────────────────────────────────────────────
if gmt grdcut @earth_relief_01m -R{R} -Gtopo.grd 2>/dev/null; then
echo "01m OK"
elif gmt grdcut @earth_relief_02m -R{R} -Gtopo.grd 2>/dev/null; then
echo "02m OK"
else
gmt grdcut @earth_relief_05m -R{R} -Gtopo.grd
fi
Z_MIN=$(gmt grdinfo topo.grd -C | awk '{{print $6}}')
Z_MAX=$(gmt grdinfo topo.grd -C | awk '{{print $7}}')
gmt makecpt -Cgeo -T${{Z_MIN}}/${{Z_MAX}} -Z > topo.cpt
gmt grdimage topo.grd -J{J} -R{R} -Ctopo.cpt -I+d
# ── Coast + frame ───────────────────────────────────────────────────────
gmt coast -R{R} -J{J} -W0.5p,gray30 -N1/0.6p,gray50 -A500 \\
-Bxaf+l"Longitude" -Byaf+l"Latitude" -BWSne+t"Ray Paths"
# ── Ray paths (semi-transparent blue) ──────────────────────────────────
gmt plot {rays_file} -R{R} -J{J} -W0.4p,dodgerblue@60
# ── Station / source points ─────────────────────────────────────────────
gmt plot {pts_file} -R{R} -J{J} -Sc0.12c -Gred -W0.4p,black
gmt colorbar -DJBC+w7c/0.35c+o0/0.5c -Baf+l"Elevation (m)"
gmt end
"""
run_gmt(gmt_script, outname="ray_paths_topo", title="Ray path visualization with topography")
os.remove(rays_file); os.remove(pts_file)
print(f"Plotted {len(ray_ids)} ray paths, {len(seen)} unique endpoints")
Projection Rule: Rectangular Map vs Orthographic Projection
If the user asks for a normal rectangular geographic map, DO NOT use orthographic or oblique projections.
Use:
-JM15c
or, for regional maps:
-JL${lon0}/${lat0}/${lat1}/${lat2}/15c
Avoid these projections unless the user explicitly requests a globe-style or oblique map:
-JG
-JO
-JA
Important:
- A normal longitude-latitude regional map should look rectangular.
- If the plotted map appears as a tilted parallelogram or curved globe view, the projection is probably wrong.
- For most seismic station maps, use
-JM by default.
Stable station map template:
```python
gmt_script = f"""
gmt begin station_map PNG
gmt set MAP_FRAME_TYPE plain
gmt set FORMAT_GEO_MAP ddd.x
REGION="-R{lon_min}/{lon_max}/{lat_min}/{lat_max}"
PROJ="-JM15c"
# Topography with fallback + explicit CPT (avoids black-terrain bug)
if gmt grdcut @earth_relief_01m $REGION -Gtopo.grd 2>/dev/null; then
echo "01m OK"
elif gmt grdcut @earth_relief_02m $REGION -Gtopo.grd 2>/dev/null; then
echo "02m OK"
else
gmt grdcut @earth_relief_05m $REGION -Gtopo.grd
fi
Z_MIN=$(gmt grdinfo topo.grd -C | awk '{{print $6}}')
Z_MAX=$(gmt grdinfo topo.grd -C | awk '{{print $7}}')
gmt makecpt -Cgeo -T${{Z_MIN}}/${{Z_MAX}} -Z > topo.cpt
gmt grdimage topo.grd $REGION $PROJ -Ctopo.cpt -I+d
# Coastline and frame (no -G fill)
gmt coast $REGION $PROJ -W0.5p,gray40 -N1/0.5p,gray60 -A500 \\
-Bxa2f1+l"Longitude" \\
-Bya2f1+l"Latitude" \\
-BWSne+t"Station Map"
# Plot stations
gmt plot stations.txt $REGION $PROJ -St0.25c -Gred -W0.4p,black
gmt colorbar -DJBC+w6c/0.3c+o0/0.4c -Ctopo.cpt -Baf+l"Elevation (m)"
gmt end
"""
Common GMT6 Command Reference
| Command | Purpose |
|---|
gmt basemap | Draw base map frame and axes |
gmt coast | Coastlines, national borders, terrain fill |
gmt plot | Plot points, lines, polygons |
gmt text | Text annotation |
gmt grdimage | Raster data rendering (topography, velocity models) |
gmt colorbar | Color bar |
gmt meca | Focal mechanism beachball |
gmt coupe | Focal mechanism cross-section |
gmt project | Coordinate projection (for cross-sections) |
gmt surface | Data interpolation to grid |
gmt contour | Contour lines |
Notes
- GMT6 modern mode:
gmt begin <name> PNG ... gmt end (recommended)
- Script execution directory is
SAGE_OUTDIR (temporary directory); data files must use absolute paths
run_gmt automatically replaces the file name after gmt begin with outname- Output
.sh script can be re-run directly in terminal, fully reproducible
