gmt-plotting

# 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: ```bash 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: ```bash 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): ```bash # --- 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. ```bash #!/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: ```python 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 ```bash -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: ```bash -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 ```bash 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 ```bash 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: ```bash # 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. ```bash 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 ```bash # 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: ```python 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 """ ```