5. 為地形圖上色#

注意

本教學適用於 GMT 6 的現代模式。如須參閱 GMT 6 (傳統模式) 與 GMT 4-5 繪製相同地圖的教程,請至這裡

把地圖依照不同屬性、在不同區域著色,可以把龐大的資料有效視覺化。最常用來著色的屬性,高度絕對是其中之一。在典型的地形圖中,不同的海拔高度會套上不同的顏色。本章中將會詳細說明如何利用地表的不同高度,製作彩色的地形圖。

5.1. 目標#

製作一張的的喀喀湖 (Lago Titicaca) 的分層設色地形圖,如下圖所示。的的喀喀湖位於秘魯和玻利維亞交界,面積達 8372 平方公里。特殊的是,湖面位置的海拔高達 3812 公尺,是全世界最高的商業通航湖泊1。圖中的藍色線條表示河流,注意右側的河流其實是亞馬遜河的源頭,而且並沒有和的的喀喀湖連接。黑色的虛線是兩國國界。地形圖的顏色從海拔 0 至 5000 公尺,依循以下顏色過渡:綠色 -> 黃色 -> 紅色 -> 灰色 -> 白色,其中 2800 公尺以上即開始以灰色系展示。這麼高的湖泊,位於一片白色系背景當中,更增添了些許雪國的氣氛 (實際上,的的喀喀湖冬天非常乾燥,幾乎不會下雪2)。

../_images/titicaca_gmt6.png

直接觀看指令稿

5.2. 使用的指令與概念#

  • grdinfo - 查閱網格檔的基本資訊

  • grdimage - 繪製著色影像

  • coast - 繪製河流、湖泊與國界

  • 使用 GMT 伺服器的地形資料繪圖

  • cpt-city 尋找適合地圖的色階檔

  • 如何疊加輸出圖層

  • 使用腳本、gmt begingmt end,一次下達所有繪圖指令

5.3. 取得資料#

這次我們必須要先取得繪圖區域的地形資料。GMT 6 本身有提供全球的地形資料,不過是放在遠端伺服器上,當使用者需要時才會下載到本機電腦中。這些地形資料的檔名為:

@earth_relief_rru

其中 rr 是解析度的大小,而 u 是單位 (有度 d、分 m、秒 s 可選)。我們在這邊要使用的是全球解析度 1 弧分的地形高程,因此檔名為 @earth_relief_01m。其他可供選擇的資料,皆在 GMT 官方手冊中有說明。如果你想查閱地形資料的更多資訊,可以透過終端機輸入 grdinfo 指令來達成。執行畫面會像如下所示:

$ gmt grdinfo @earth_relief_01m
earth_relief_01m: Download file from the GMT data server [data set size is 214M].
earth_relief_01m: Earth Relief at 1x1 arc minutes obtained by Gaussian Cartesian filtering (1.9 km fullwidth) of SRTM15+V2 [Tozer et al., 2019].

/home/whyj/.gmt/server/earth_relief_01m.grd: Title: Earth Relief at 01 arc minute
/home/whyj/.gmt/server/earth_relief_01m.grd: Command: grdfilter SRTM15+V2.nc -Fg1.9 -D1 -I01m -rg -Gearth_relief_01m.grd=ns --IO_NC4_DEFLATION_LEVEL=9 --PROJ_ELLIPSOID=Sphere
/home/whyj/.gmt/server/earth_relief_01m.grd: Remark: Obtained by Gaussian Cartesian filtering (1.9 km fullwidth) from SRTM15+V2.nc [Tozer et al., 2019; http://dx.doi.org/10.1029/2019EA000658]
/home/whyj/.gmt/server/earth_relief_01m.grd: Gridline node registration used [Geographic grid]
/home/whyj/.gmt/server/earth_relief_01m.grd: Grid file format: ns = GMT netCDF format (16-bit integer), CF-1.7
/home/whyj/.gmt/server/earth_relief_01m.grd: x_min: -180 x_max: 180 x_inc: 0.0166666666667 (1 min) name: longitude n_columns: 21601
/home/whyj/.gmt/server/earth_relief_01m.grd: y_min: -90 y_max: 90 y_inc: 0.0166666666667 (1 min) name: latitude n_rows: 10801
/home/whyj/.gmt/server/earth_relief_01m.grd: z_min: -10907 z_max: 8170 name: elevation (m)
/home/whyj/.gmt/server/earth_relief_01m.grd: scale_factor: 1 add_offset: 0
/home/whyj/.gmt/server/earth_relief_01m.grd: format: netCDF-4 chunk_size: 129,129 shuffle: on deflation_level: 9

訊息的前兩行在有需要下載時才會出現。下載下來的檔案預設是放在 ~/.gmt/server/ 這個位置,之後 GMT 會先查看這個位置有沒有已經下載的檔案,如果沒有的話才會再次下載。下載檔檔名為 earth_relief_01m.grd,副檔名 .grd 表明了它是某種「網格檔 (grid)」,亦即地球的二維表面被切割成網格狀,每個「格子」內,都有一個相對應的高度數據。(非常類似 GeoTiff 的概念,只不過在 GeoTiff 中,相對應的術語稱為「像素」而非「網格點」。) .grd 檔是 netCDF 格式常見的附檔名,也就是說它可被 GMT 讀取並作圖。另外,從 grdinfo 的輸出中你還可以看到其他的基本資訊,像是

  • 網格檔是從 SRTM15+V2 這個資料庫重新取樣而來。

  • X 從 -180 度到 180 度,Y 從 -90 度到 90 度,意味著此網格檔涵蓋全球的範圍。

  • 網格大小是 0.0166666666667,也就是 1/60 度 (1 弧分)。

  • 總共有 21601 x 10801 = 233312401 個格點 (也就是兩億三千三百萬像素)。

  • Z 值 (高度) 從 -10907 到 8170 公尺。為什麼最高的地方不是 8848 公尺 (聖母峰的高度)? 那是因為網格的解析度還是不夠高,1 弧分其實只相當於 1 至 2 公里才取一個點而已。

5.4. 操作流程#

從本章開始,我們需要使用兩個以上的 GMT 指令製作地圖。在上一章介紹過的 coast 用來畫向量式資料 (河流、湖泊以及國界),而在本章中新介紹的 grdimage 則用來畫地形起伏。因為上一章介紹的 -png 選項只適用於單一指令,所以在本章我們要用另一種更通用的方式指定輸出檔。具體做法是使用 gmt begingmt end 這兩個關鍵字串把其他所有繪圖指令包起來,像這樣:

$ gmt begin [輸出地圖主檔名,不用加附檔名] [格式]      # 輸入後按下 enter
$ [繪圖指令1]    # 輸入後按下 enter
$ [繪圖指令2]    # 輸入後按下 enter
...
$ gmt end

你可以一行一行的在命令列中打好指令再送出執行,不過為了修改方便,我們通常會選擇把所有的指令先寫在腳本內,再執行此腳本以一次送出所有指令給電腦。在本章中,我們就要來介紹這種更爲便捷的方式。

畫地圖的第一步就是決定繪製區域。如同前一章「製作第一張地圖」所述,我們可以使用 Google Map 等工具決定區域。這裡我們考慮到高度的呈現效果,選擇經度從 -70.8 到 -66.56 (負值表示西經),緯度從 -17.17 到 -14.42 (負值表示南緯。)

../_images/coloring_topography_fig1.png

圖 5.1 本次作圖區域在 Google Map 上呈現的樣子。#

使用 grdimage 呈現地形數據,最少需要提供兩個參數,分別是 輸入檔-J (投影方法和地圖大小)。在這裡我們也需要使用 -R 指定繪製區域,如果不加 -R 的話,GMT 會把全網格檔的資料話在地圖上 (因爲我們使用的是全球地形資料,所以輸出會是全球地圖)。就讓我們來畫一張跟前一章一樣寬度的地圖吧。請建立新文字檔,檔名可以為 titicaca.bash 或任何你喜歡、且符合作業系統格式慣例的名字。然後在檔案中輸入:

gmt begin titicaca png
    gmt grdimage @earth_relief_01m -R-70.8/-17.17/-66.56/-14.42r -JM6i
gmt end

你會注意到我們在 gmt grdimage 前面多加了四個空格,這純粹只是製造縮排效果,讓閱讀腳本的人知道這行被包在 gmt begingmt end 之間,增加腳本的可讀性。縮排與否和空格數量不會影響腳本的執行。如需 -R-J 選項的語法說明,請至「製作第一張地圖」中查閱。存檔後,在終端機模式下前往 titicaca.bash 所在的資料夾,然後執行腳本 (以 Linux 系統爲例):

$ bash titicaca.bash

備註

不同的作業系統的腳本會有不同的額外設定,可以讓你的腳本能夠更方便的執行,詳情請參閱「基本概念與名詞」。

輸出的圖檔 titicaca.png 開啟後應該會是如下的圖案:

../_images/coloring_topography_gmt6_fig2.png

有點難以看懂?這是因為 GMT 預設的套色沒那麼適合這個地區的緣故。在未指定地圖的套色下繪製 @earth_relief 地形資料時,GMT 會用一套名為 geo 的色階來呈色。這個色階的顏色如下所示

../_images/coloring_topography_gmt6_fig3.png

圖中的 01 的顏色,也就是深藍色和白色,會對應到地圖範圍中最低和最高的的地方。有沒有辦法換成另一種比較好看的色階?當然可以,GMT 有數種方法可以製造或調整色階,但這邊我們要先來看一種最簡單的方法:直接下載已經做好的色階來用。

cpt-city 網站蒐集了不少常用的色階,其中也包括 GMT 格式的色階檔。GMT 色階檔通常使用副檔名 .cpt,是 color palette table (色階表) 的縮寫。你可以在網站中搜尋你想要的色階,例如這個網頁就列了許多最常被下載的 .cpt 檔。我們這次來用其中的一個叫做 mby 的色階看看。按下 mby 的色階表,你應該就會看到如下頁面

../_images/coloring_topography_fig4.png

注意圖中的 -8000 … 5000,這是指此色階已經被指定在 -8000 到 5000 公尺之內呈現了。按一下網頁上的 cpt 連結就可以把此色階檔下載下來。

為了方便起見,這邊也提供了 mby.cpt 的下載點。

使用 grdimage 中的 -C 選項,就可以指定要用的色階,只要在 -C 後接續色階檔檔名就行了。所以我們更新版的腳本為:

gmt begin titicaca png
    gmt grdimage @earth_relief_01m -R-70.8/-17.17/-66.56/-14.42r -JM6i -Cmby.cpt
gmt end

使用相同的方式再執行一次腳本,這次地圖會變成這樣

../_images/coloring_topography_gmt6_fig5.png

看起來已經有模有樣了。

備註

有關於如何自訂色階表,讓地圖的顏色可以隨心所欲變化的技巧,請參閱「畫筆與調色技巧」與「色階檔的編修與呈現」。

這個地圖目前有個大問題,就是我們雖然知道地圖中大部份的區域都很高 (灰白色),但是卻看不到我們的主角的的喀喀湖。這是因為地形資料只有記載高度,並沒有記載地表的狀態,所以光靠 grdimage 是無從得知哪邊是湖泊、海洋等等地貌的。因此,為了標上的的喀喀湖的位置,我們需要呼叫 coast 這個指令,利用內建的水體資料幫我們畫上湖泊,順便再畫上邊框。

當一張地圖包含兩個以上的作圖指令時,越上方的圖層要越後畫。由於我們希望 coast 畫的湖泊水體疊加在 grdimage 畫的地形資料上,所以 coast 指令要排在 grdimage 之後。再次打開腳本,把內容修改爲:

gmt begin titicaca png
    gmt grdimage @earth_relief_01m -R-70.8/-17.17/-66.56/-14.42r -JM6i -Cmby.cpt
    gmt coast -Scornflowerblue -Baf -B+t"Lago Titicaca"
gmt end

關於 coast 的指令解說,請至「製作第一張地圖」查看。另外,你或許會發現在 coast 中,沒有 -R-J 選項。當不指定這些選項時,GMT 會自動沿用上次的設定,也就是 grdimage-R-J 資訊,省下重複輸入相同內容的麻煩。

小訣竅

  1. coast-S 選項會套用在所有水域。如果你的地圖中有海又有湖,而你想分別為海跟湖指定不同顏色,可以試試看換成 -C 選項。因為本圖沒有海域,-Scornflowerblue-Cl/cornflowerblue 具有相同效果,不過後者其實只會為湖泊填上顏色。

  2. 讓我們看一下 coast 中的兩個 -B 選項。實際上它們也可以加在 grdimage 中,因為邊框屬於哪個圖層在本例中不影響地圖疊加,所以效果會一模一樣。

在終端機中執行腳本,然後開啟輸出檔,應該會看到下圖:

../_images/coloring_topography_gmt6_fig6.png

這次總算有個令人滿意的輸出結果。接下來,我們還可以利用 coast 的其餘功能來美化這張地圖。首先是 -I,它可以幫我們標上河流,語法為

-I編號/畫筆

其中的編號可以為不同的英數字,對應到不同等級的河流;而畫筆則用來指定河流的繪圖樣式。這裡我們要使用 r 這個編號,它對應的是所有的永久河流。畫筆樣式則使用 0.7p,cornflowerblue

另外一個 coast 的選項是 -N,它可以在地圖上標上政治性的邊界。語法同樣為

-N編號/畫筆

這裡的編號一樣可以為不同的英數字,對應到不同等級政治邊界 (國界、州界等)。我們要用的是 1,亦即標上國界。至於畫筆樣式,則使用 0.7p,,--,注意這裡我們空下了「顏色」的欄位,因此 GMT 會使用預設值「黑色」進行繪圖。在「線條樣式」欄位中,我們選擇 --,這樣國界就會以虛線呈現。

備註

有關畫筆的進階樣式設定,例如設定線條樣式、利用 RGB 值自訂顏色等,會在「畫筆與調色技巧」詳加描述。

-N-I 以及他們的設定參數加入我們的 coast 指令中,就完成了這次的地圖設計。

5.5. 指令稿#

本地圖的最終指令稿如下:

gmt begin titicaca png
    gmt grdimage @earth_relief_01m -R-70.8/-17.17/-66.56/-14.42r -JM6i -Cmby.cpt
    gmt coast -Ir/0.7p,cornflowerblue -N1/0.7p,,-- -Scornflowerblue -Baf -B+t"Lago Titicaca"
gmt end

備註

「繪製地形設色圖:以 70.8E - 66.56E,17.17S - 14.42S 為邊界,麥卡托投影,地圖橫向寬 6 吋,套色使用 mby.cpt 內的色階表。使用 cornflowerblue 把湖泊和河流上色,畫上國界,四邊繪製邊框與刻度數值,自動選擇數值間隔與刻度間隔,不畫格線,標題為 Lago Titicaca,存檔為 titicaca.png。」

觀看最終版地圖

5.6. PyGMT 指令稿#

使用 PyGMT,可以如下腳本繪製本地圖:

import pygmt

# 取得地形資料
region = [-70.8, -66.56, -17.17, -14.42]
grid = pygmt.datasets.load_earth_relief(resolution='01m', region=region)

# 繪圖
fig = pygmt.Figure()
fig.grdimage(grid, region=region, projection='M6i', cmap='mby.cpt')
fig.coast(rivers='r/0.7p,cornflowerblue', borders='1/0.7p,,--', water='cornflowerblue', frame=['af', '+t"Lago Titicaca"'])
fig.show()
fig.savefig('titicaca_pygmt.png')

你可以使用以下的 Binder 連結嘗試此程式碼:

https://mybinder.org/badge_logo.svg

5.7. 習題#

  1. 畫出貝加爾湖 (Lake BaikalBaygal nuur) 與週邊區域的地形設色圖。貝加爾湖是全世界最深的湖泊,也是世界第七大湖,位於俄羅斯,南端接近蒙古的邊界。

  2. 從地圖中可看到有幾條河流注入或流出貝加爾湖?

  3. 貝加爾湖週邊的海拔為何?湖底的海拔為何?你的地圖分別用什麼顏色表示這兩個海拔?

1

Data Summary: Lago Titicaca (Lake Titicaca). International Lake Environment Committee Foundation - ILEC. Retrieved 2009-01-03.

2

ClClimatological Information for Juliaca, Peru. Hong Kong Observatory.