# pyGMT
Generic Mapping Tools (GMT) 誕生於1988年，由在 Lamont-Doherty Earth Observatory 的 Pål Wessel 和 Walter H. F. Smith 開始製作，並在 1991 年公開發表第 1 版。在 2019 年 11 月已經更新到第 6 版。過去三十多年間，GMT 已然成為全世界通用的地理製圖工具，在地質學、地球物理學與海洋科學等更可以說是最常接觸到的軟體之一。
隨著時代的演進，Python 的地位逐漸爬升，成為了當代科學分析的龍頭語言。長遠來看，將各種不同的軟體置換成 Python 的語法儼然成為一大趨勢。GMT 也被重新包裝成 Python 的形式。與 GMT 同樣不變的是，pyGMT 也有圖層的概念，先執行的指令會是最底下的圖層，後執行的指令所繪製的圖層會覆蓋在原有的圖層上，一層疊一層的圖層結構可以被 Ai（adobe illustrator）等繪圖軟體識別並拆開圖層作修改。而在編輯 pyGMT 時也別忘了利用圖層的概念去設計自己想展示的圖片喔！
接下來，本筆記將參考 [pyGMT 官方網站](https://www.pygmt.org/latest/)與[鄭懷傑所撰寫的中文手冊](https://gmt-tutorials.org/main/intro_install.html)，將一些常用的指令與功能說明撰寫於本筆記中。主要也會著重在地震學相關的地理繪圖上，若需要更深入的了解 pyGMT 或學習其他學門的製圖，可以查閱本筆記的參考資料。
> @KGJ0717 （柯俊杰）撰寫於 2024 年 7 月 30 日 - 8 月 13 日。
# 下載與安裝 pyGMT
先建立一個新環境（名為：pygmt），並安裝 pygmt 及其他相應的工具
``` linux=
conda create --name pygmt --channel conda-forge pygmt
```
接下來只要啟動該環境，就可以執行 python 的腳本囉！
``` linux=
conda activate pygmt
```
# 繪製海岸線
海岸線的區隔是地理製圖的基本功。
以下將透過一個簡單的範例說明一些基本的製圖函數使用方法。
> 利用 `pygmt.Figure()` 來建立一張空白的底圖並命名為 fig。
>
> 接著利用 `pygmt.Figure().coast()` 來繪製海岸線。其中需要給予以下參數：
> > `region=[119.2, 122.2, 21.8, 25.8]` 依序為起始經度、終止經度、起始緯度、終止緯度，範例為台灣的經緯度範圍。
> >
> > `frame=['a2f0.5g1', 'WSne+t"Taiwan"']` 用於設定圖的數標、格線、刻度及標題，`a` 表示經緯度數值每2度顯示一次；`f` 表示黑白標線每0.5度繪製成一格；`g` 表示中央畫布的格線每1度繪製一條格線。`WSne` 代表圖的左（W）、下（S）顯示經緯度數值標示（大寫表示要標示），而右（e）、上（n）則不顯示（小寫表示不標示）。`+t"[title]"` 用於幫圖加上標題。
> >
> > `projection="M6i"` 為投影方式，除了麥卡托圓柱投影（M[lon/[lat/]]width）之外，常用的投影方式還有繪製遠震分布時會用到的方位等距投影（Elon/lat[/horizon]/width）。6i則是圖的大小（width）為6英吋。
> >
> > `shorelines=['0.1p', 'black']` 為海岸線的繪製，粗細為0.1p，顏色為黑色。
> >
> > `water='cornflowerblue'` 用於設定海洋為矢車菊藍色。
> >
> > `land='lightgreen'` 用於設定陸地顏色為綠色。
> >
> > `rivers='r/0.7p,cornflowerblue'` 用於設定河流資訊，r代表只畫重要大河，0.7p為繪製粗度，cornflowerblue表示用矢車菊藍色來繪製河流。
> >
> > `area_thresh=0.1` 省略面積小於0.1平方公里的小島。
> >
> `pygmt.Figure().savefig('Taiwan.png')` 可以用來儲存圖片。
>
> `pygmt.Figure().show()` 則會將透過外部視窗開啟圖片檢視。
>
完整範例如下：
``` python=
import pygmt
fig = pygmt.Figure()
fig.coast(region=[119.2, 122.2, 21.8, 25.8],
frame=['a2f0.5g1', 'WSne+t"Taiwan"'],
projection="M6i",
shorelines=['0.1p', 'black'],
water='cornflowerblue',
land='lightgreen',
rivers='r/0.7p,cornflowerblue',
area_thresh=0.1)
fig.savefig('Taiwan.png')
fig.show()
```
結果如下圖所示:
# 繪製地形圖
完成了海岸線的區隔之後，陸地與海底地形也必須被加入到圖中。
除了自己本身就有地形資料之外，以下也會簡單介紹如何向 GMT6 的遠端伺服器索取地形資料。
**註：向遠端伺服器下載地形檔時請務必確認網路連線正常。**
> `pygmt.datasets.load_earth_relief` 用來從 GMT6 伺服器下載地形的網格檔（.grd）。
> > 利用 `region` 來指定切網格的範圍
> > `resolution` 用來決定解析度（網格大小），有 `"01d"`, `"30m"`, `"20m"`, `"15m"`, `"10m"`, `"06m"`, `"05m"`, `"04m"`, `"03m"`, `"02m"`, `"01m"`, `"30s"`, `"15s"`, `"03s"` 或 `"01s"` 等 15 種解析度可以選擇，其中d表示角度、m表示角分、s表示角秒。
> >
> `pygmt.makecpt` 用來製作個人化色標，當然也可以用預設即可。
> > `series=[-6000, 6000, 10]` 表示色軸最小值為-6000，最大值為6000，每間隔10就分隔出一種顏色。
> >
> > `cmap='gmt/globe'` 表示利用gmt預設的色軸（colormap）globe來配置個人化色軸。
> > * 除此之外常用的還有 `'gmt/abyss'`, `'gmt/bathy'`, `'gmt/ocean'`, `'SCM/oslo'`, `'SCM/devon'` 可以用來繪製純海洋
> > `'gmt/dem1'`, `'gmt/dem2'`, `'gmt/dem3'` 用來繪製純陸地
> > * 同時有海陸的地形可以用 `'gmt/earth'`, `'gmt/globe'`, `'gmt/geo'`, `'gmt/etopo1'`, `'gmt/relief'`, `'gmt/terra'`, `'gmt/world'`, `'SCM/bukavu'`, `'SCM/oleron'`, `'SCM/fes'`, `'cmocean/topo'` 等。
> > * `'gmt/gray'` 與 `'SCM/grayC'` 則是用灰階處理所有地形起伏。
> > *[色軸表](https://docs.generic-mapping-tools.org/6.5/reference/cpts.html#built-in-color-palette-tables-cpt)中其餘用途的色軸選項之後再介紹。*
> >
> > `background="i"` 表示低於新色軸最小與最大值的數值用cmap的最小與最大值顏色來代替，預設為"o"表示用預設的顏色（最大與最小為黑色；Nan為白色）來繪製。
> >
> > `reverse=False` 表示將色軸對應到的顏色顛倒過來。
> >
> > `transparency=0` 表示透明度，數字可以允許從0到100，100為完全透明。
> >
> > `output='TW_topo.cpt'` 則是將製作好的光柵圖形編輯檔（.cpt）儲存起來。
> >
>
> `pygmt.Figure().grdimage` 用於將網格檔繪製成色彩圖。
> > `grid` 輸入網格檔。
> >
> > `region` 決定繪圖範圍，網格檔可以比繪圖範圍更寬廣，但不會繪製出來。若網格檔所涵蓋的範圍比繪圖範圍窄，那麼沒有涵蓋的部分則會沿用原來底圖的設置。
> >
> > `projection` 表示地圖的投影方式，與上一章節相同。
> >
> > `cmap` 表示繪製色彩時所使用的色軸，可以用自定義的色軸，也可以用內建的色軸。
> > `transparency` 用來調整繪圖透明度。輸入0-100之間的數，數字越高越透明。
> >
> `pygmt.Figure().grdcontour`
> > `grid` 輸入網格檔。
> >
> > `annotation=1000` 表示標示數值的等值線間距為1000。
> >
> > `levels=500` 表示等值線間距為500，也就是每條等高線都間隔500公尺的高差，也可以給一個列表詳細列出每一條等值線所代表的數值。此功能用於繪製等值線但未必會顯示數值。
> > *不過在測試中發現此功能常常會無法正常執行*
> >
> >
> `pygmt.Figure().colorbar`
> > `cmap='TW_topo.cpt'` 表示繪製色彩時所使用的色軸，可以用自定義的色軸，也可以用內建的色軸。
> >
> > `position='n0.35/0.95+w6c/0.2c+h'` 表示色條（colorbar）放置的位置。其中有多種模式可以選擇，如利用地圖座標來定位的g模式、利用參考點定位的j模式、畫布座標來定位的x模式（用公分或英吋輸入），以及歸一化畫布座標來定位的n模式。其中，筆者認為n模式是最好理解跟使用的。將畫布（就是繪製地圖的中間圖框）最左下角設定為0/0，右上角為1/1，數字可以小於0或大於1。`'n0.35/0.95'` 代表將色調繪製在從畫布左下角起算，0.35個畫布寬度、0.95個畫布高度的地方。`'+w6c/0.2c'` 則是設定色條的形狀為6公分長乘上0.2公分寬。`'+h'` 表示色條以水平方向擺放；`'+v'` 表示色條以鉛直方向擺放。
> >
> > `frame=['a2000f1000', 'x+lBathymetry&Topography', 'y+lm']` 用於設定色條的數標、刻度、格線及標籤，`a` 表示數值每2000度顯示一次；`f` 表示刻度每1000度繪製成一刻度；`g` 表示色條中的格線，不加則不畫。若要新增標籤（`+l`），則用 `'x+lBathymetry&Topography'` 與 `'y+lm'` 分別在色條的下方與右方加上標籤（label）。
> >
> > `box='+gwhite@50'` 用來幫色條增加底框。`'+g[color]@[transparency]'` 用來調整底框的顏色跟透明度。
> >
>
完整範例如下：
``` python=
import pygmt
# 取得地形資料
region = [119.2, 122.2, 21.8, 25.8]
grid = pygmt.datasets.load_earth_relief(resolution='01m',
region=region)
pygmt.makecpt(series=[-6000, 6000, 10],
cmap='gmt/globe',
background="i",
reverse=False,
transparency=0,
output='TW_topo.cpt')
# 繪圖
fig = pygmt.Figure()
fig.grdimage(grid, region=region, projection='M6i', cmap='TW_topo.cpt', transparency=0)
fig.coast(frame=['af', '+t"Taiwan"'],
rivers='r/0.7p,cornflowerblue',
area_thresh=0.1)
fig.grdcontour(grid=grid, annotation=1000, levels=500)
fig.colorbar(cmap='TW_topo.cpt',
position='n0.35/0.95+w6c/0.2c+h',
frame=['a2000f1000', 'x+lBathymetry&Topography', 'y+lm'],
box='+gwhite@50')
fig.savefig('Taiwan_topo.png')
fig.show()
```
結果如下圖所示:
如果覺得地圖的網格感太重，表示解析度不夠高，可以從01m改為15s或更高的03s。
# 標示地震測站與構造線
測站位置與構造線屬於基本的繪製點與線的範例。
在繪製好地形圖之後，老闆會叫你把"**我們的測站**"畫上去，也會告訴你"**順便畫個構造線吧**"，這時候你可以偷偷把地形圖的透明度調高到30左右，讓你的圖可以更加強調出測站與構造線的位置，觀眾會更容易看出這張圖的重點，而且不會太多繽紛的色彩而頭昏眼花。
> `fig.plot()` 用於在畫布上繪製線條、圖示，或加上更多裝飾。
>
> > `x=x1` 給定繪圖位置的 X 座標。
> >
> > `y=y1` 給定繪圖位置的 Y 座標。
> >
> > `data='Manila_Trench.txt'` 可以直接用檔案給予繪圖位置的座標。
> >
> > `fill='#FF2222'` 給定標誌內部需塗滿的顏色，可以用 16 進位的色碼（如：#FF0000），也可以直接用顏色的英文（如：red）。
> >
> > `style='t0.7c'` 決定繪製的圖示。`c` 為圓形、`t` 為正三角形、`i` 為倒立正三角形、`s` 為正方形、`d` 為菱形（詳細可以參考[標誌圖例](https://www.pygmt.org/latest/gallery/symbols/basic_symbols.html#sphx-glr-gallery-symbols-basic-symbols-py)）。後方加上的 `0.7c` 指的是圖示的大小為 0.7 cm。
> >
> > `style='f1.4c/0.25c+l+t'` 如果想畫更酷的線，請參考[線段圖示](https://www.pygmt.org/latest/gallery/lines/linefronts.html#sphx-glr-gallery-lines-linefronts-py)。其中的 `f1.4c` 表示圖示的間隔為 1.4 cm，`0.25c` 表示圖示的大小為 0.25 cm，`+l` 表示將圖示繪製於畫線方向的左手邊，`+t` 表示圖示使用正三角形。這是一個逆斷層常見的線段圖示畫法。
> >
> > `pen='0.1c,black,dash'` 用來控制繪製的線段，`0.1c` 表示線段寬度為 0.1 cm，`black` 表示繪製線段為黑色，`dash` 表示用虛線繪製線段。也可以參考[線型種類](https://www.pygmt.org/latest/gallery/lines/linestyles.html#sphx-glr-gallery-lines-linestyles-py)來繪製其他的線段喔！
> >
完整範例如下：
``` python=
import pygmt
import numpy as np
import pandas as pd
# 取得資料與設定繪圖參數
region = [119.2, 122.2, 21.8, 25.8]
grid = pygmt.datasets.load_earth_relief(resolution='01m',
region=region)
pygmt.makecpt(series=[-6000, 6000, 10],
cmap='gmt/globe',
reverse=False,
output='TW_topo.cpt')
SALUTE_realtime = pd.read_csv('stations.csv')
SALUTE_realtime_style = {'style': 't0.7c',
'fill': '#FF2222',
'pen': '0.05c,black'}
# 繪圖
fig = pygmt.Figure()
fig.grdimage(grid, region=region, projection='M6i',
cmap='TW_topo.cpt', transparency=45)
fig.coast(frame=['af', '+t"Taiwan"'],
rivers='r/0.7p,cornflowerblue',
borders='1/0.7p,,--',
area_thresh=0.1)
fig.grdcontour(grid=grid, annotation=1000, levels=500)
fig.colorbar(cmap='TW_topo.cpt', position='n0.35/0.95+w6c/0.2c+h',
frame=['a2000f1000', 'x+lBathymetry&Topography', 'y+lm'],
box='+gwhite@50')
fig.plot(data=np.array([SALUTE_realtime['longitude'], SALUTE_realtime['latitude']]).T,
**SALUTE_realtime_style)
fig.plot(data='Manila_Trench.txt',
style='f1.4c/0.25c+l+t',
fill="black",
pen='0.1c,black,dash')
fig.savefig('Taiwan_station.png')
```
結果如下圖所示:
# 標示地震分布
作為一個地球物理學家，乃至於地質學家，地震分布除了可以告訴我們板塊邊界的輪廓之外，地震活動度與地震規模大小也都是值得被關注的點。
但是本範例僅用於簡單介紹各函數的使用方法與參數給定，並不會深入討論地震分布的意義與構造解釋，因此用於繪製的地震數量極少，選取地震的區域也很小。
> `fig.plot()` 同樣使用該函數，但這時候需要根據資料數值來客製化資料點的繪圖方式，因此會需要用到更多的參數設定。
>
> > `x=x1` 給定繪圖位置的 X 座標。
> >
> > `y=y1` 給定繪圖位置的 Y 座標。
> >
> > `zvalue=z1` 給定繪圖顏色的數值。
> >
> > `cmap='gmt/seis'` 給定繪圖顏色數值所對應的色軸。在此處使用的是 gmt 的 seis 色軸。
> > 常用於表達震源深度的色軸有除了 `gmt/seis` 之外，`gmt/haxby` 、 `gmt/mag` 、 `gmt/rainbow` 、 `google/turbo` 、 `matlab/jet` 也會被用來繪製彩色的震源深度，但是需要將色軸顛倒過來。筆者認為 `matlab/copper` 、 `matplotlib/viridis` 等也可以用於表示震源深度，同樣需要將色軸顛倒過來。
> > *[色軸表](https://docs.generic-mapping-tools.org/6.5/reference/cpts.html#built-in-color-palette-tables-cpt)中其餘用途的色軸選項之後再介紹。*
> >
> > `size=s1` 給定繪製圖示的大小。
> >
> > `data=[[x1, y1, z1, s1], [x2, y2, z2, s2], ...]` 當然也可以利用矩陣的方式給予，對於每一個繪圖點，依序給予 X, Y 座標、顏色數值、大小，可以一次繪製許多點。
> >
> > `data='xyzs.txt'` 若已經整理成一個文字檔，則檔案內的順序一樣採用 X, Y 座標、顏色數值、大小的順序。
> >
完整範例如下：
``` python=
import pygmt
import numpy as np
import pandas as pd
# 取得資料與設定繪圖參數
region = [119.2, 122.2, 21.8, 25.8]
grid = pygmt.datasets.load_earth_relief(resolution='01m',
region=region)
pygmt.makecpt(series=[-6000, 6000, 10],
cmap='gmt/globe',
reverse=False,
output='TW_topo.cpt')
SALUTE_realtime = pd.read_csv('stations.csv')
SALUTE_realtime_style = {'style': 't0.7c',
'fill': '#FF2222',
'pen': '0.05c,black'}
seismicity = pd.read_csv('seismicity.csv')
pygmt.makecpt(series=[0, 60, 0.1],
cmap='gmt/seis',
reverse=False,
output='TW_seis.cpt')
# 繪圖
fig = pygmt.Figure()
fig.grdimage(grid, region=region, projection='M6i',
cmap='TW_topo.cpt', transparency=45)
fig.coast(frame=['af', '+t"Taiwan"'],
rivers='r/0.7p,cornflowerblue',
borders='1/0.7p,,--',
area_thresh=0.1)
fig.grdcontour(grid=grid, annotation=1000, levels=500)
fig.colorbar(cmap='TW_topo.cpt',
position='n0.35/0.97+w6c/0.2c+h',
frame=['a2000f1000', 'x+lBathymetry&Topography', 'y+lm'],
box='+gwhite@50')
fig.plot(data=np.array([seismicity['lon'],
seismicity['lat'],
seismicity['depth'],
(seismicity['mag']/8)**4]).T,
cmap='TW_seis.cpt',
style='c',
pen='0.01c,black',
transparency=50)
fig.colorbar(cmap='TW_seis.cpt',
position='n0.35/0.91+w6c/0.2c+h',
frame=['a10f5', 'x+lDepth', 'y+lkm'],
box='+gwhite@50')
fig.plot(data=np.array([SALUTE_realtime['longitude'], SALUTE_realtime['latitude']]).T,
**SALUTE_realtime_style)
fig.plot(data='Manila_Trench.txt',
style='f1.4c/0.25c+l+t',
fill="black",
pen='0.1c,black,dash')
fig.savefig('Taiwan_seismicity.png')
```
結果如下圖所示:
# 標示震源機制解
# 繪製剖面
# 繪製立體圖
# 參考資料 Reference
* 中文教學手冊（鄭懷傑） https://gmt-tutorials.org/main/intro_install.html
* 英文官方網站 https://www.pygmt.org/latest/
* GMT 色軸表 https://docs.generic-mapping-tools.org/6.5/reference/cpts.html#built-in-color-palette-tables-cpt
* 標誌圖例 https://www.pygmt.org/latest/gallery/symbols/basic_symbols.html#sphx-glr-gallery-symbols-basic-symbols-py
* 線段圖示 https://www.pygmt.org/latest/gallery/lines/linefronts.html#sphx-glr-gallery-lines-linefronts-py
* 線型種類 https://www.pygmt.org/latest/gallery/lines/linestyles.html#sphx-glr-gallery-lines-linestyles-py
Sign in with Wallet
Connect another wallet