Linux操作不求人 - 伍章之壹 - Generic Mapping Tools (GMT) 繪圖於 Linux 環境下

      本篇的GMT(Generic Mapping Tool,http://gmt.soest.hawaii.edu/)是地球科學學術領域相當常見的地圖繪圖工具,並非常看到的格林威治時間(GMT +8)。GMT科學繪圖軟體(以下簡稱gmt),在很多地球科學相關的科學文件內圖片,可常見到他繪圖的身影。當然有很多套裝軟體亦可達到相同的目的,如 MathWorks 的 MATLAB (https://www.mathworks.com/products/matlab.html)或是 HARRIS 的 IDL(http://www.harrisgeospatial.com/ProductsandSolutions/GeospatialProducts/IDL.aspx) 也是相當常見的,或是使用esri ArcGIS等相關繪製地圖的軟體,但這些軟體若要安裝在所有電腦,應當所費不貲。故口袋不深的研究人員,皆偏好 Open Source 的繪圖軟體,可能是 Python (http://matplotlib.org/basemap/)或是 R 環境內使用繪圖Library或免費的GIS軟體,如OpenStreetMap、Quantum GIS(QGIS)來製作。但筆者僅想透過本篇,介紹筆者在研究生時期曾經做過的 gmt 繪圖,以期能對地球科學研究的同好有些許幫助。
    首先,我們先回想第二章的 bash 或 csh 所撰寫的 Shell Scripts 方法,接下來所要展示的 gmt 繪圖指令,皆會透過 Shell Scripts 的撰寫方式帶入,來畫出 PS (PostScript,http://www.adobe.com/products/postscript/,往後簡稱 ps )格式的向量型圖檔,在 Linux 環境內,可用 Ghostscript 或Ghostview (hhttp://pages.cs.wisc.edu/~ghost/,http://www.ghostscript.com/download/gsdnld.html)來觀看圖。地形深度與高度網格數值檔 ETOPO,皆可透過 NOAA(https://www.ngdc.noaa.gov/mgg/global/global.html)下載,在此先感謝這些研究單位無私地提供資料,與這些 Open Source 的創作者們,分享這些免費 "滿滿的大平台"。
     在CentOS_6 下安裝GMT軟體,我們一樣使用指令 yum install gmt,因現在epel repository 裡(用指令 yum install epel-6* 安裝) 已經包含 gmt 套件,執行後如圖5-1























(圖5-1)


安裝完成後,要再下載 NOAA 提供的全球海岸線資料(GSHHS Shoreline Database,https://www.soest.hawaii.edu/pwessel/gshhg/)(https://www.ngdc.noaa.gov/mgg/shorelines/data/gshhg/latest/),解開gshhg-gmt-2.3.5.tar.gz檔案後,再複製到 /usr/share/GMT/coast/ 內覆蓋舊檔。接著可以在Terminal介面,利用 gmtdefaults -L 來查詢預設的 gmt 環境變數與其值,如圖5-2,因變數較多,筆者將其分成圖5-2內的左中右。


























(圖5-2)

接著,我們利用 touch plot_stn.gmt ,並將權限利用 chmod u+x plot_stn.gmt 改成可以執行,來新增一個 gmt 繪圖步驟檔,並利用指令 vim plot_stn.gmt 來編輯它。筆者此篇內之gmt 的繪圖 scripts皆利用 gmt 早前版本 4.X 版,而最新版本為 5.X 版,若有一些落後的用法,再請讀者見諒,合先敘明。將 plot_stn.gmt 輸入圖5-3-1以下指令,筆者方便說明,使用指令 nl plot_stn.gmt 來印出行數,如圖5-3-2,執行 ./plot_stn.gmt 後就會出現如圖5-4的繪圖:



(圖5-3-1)







(圖5-3-2)


























(圖5-4)

     圖5-4所使用的資料為中央研究院地球科學研究所(IES,SINICA)之台灣寬頻地震網(BATS)的測站(http://bats.earth.sinica.edu.tw/Station/)(ref. 1),並與發生的一個地震事件相對地圖位置圖。第1行為宣告以下為bash shell 的指令語法,第2-4 的指令 gmtset 為要改變gmtdefaults 內的預設值,參數說明如以下
ANNOT_FONT_SIZE 代表邊框標記的經緯度字體大小,若是x-y分佈圖,則為x-y軸刻度標記的大小。
HEADER_FONT_SIZE 代表圖表主題的字體大小,本例子為"BATS Station vs Event"。
LABEL_FONT_SIZE 代表垂直與水平軸的說明標記字體大小。
接著第5-7行定義再來繪圖步驟會常用到的變數,range 代表投影的範圍,本例為投影地圖,故設定最小經度/最大經度/最小緯度/最大緯度(經度正號代表東經,經度負號代表西經,緯度正號代表北緯,緯度負號代表南緯)。proj 變數代表投影的方法,如果單純的 x-y 圖,則使用 X 表 Cartesian 座標,此例的 M 代表 Mercator 投影(各種地理投影的方式,可參考Wikipedia 說明https://en.wikipedia.org/wiki/Map_projection,gmt支援的投影方式,可參考http://gmt.soest.hawaii.edu/gmt4/gmt/html/man/psbasemap.html),M後面的數字代表投影地圖繪出的大小。psfile變數為要產生的 ps 檔,亦可透過指令 ps2epsi  stn.ps  stn_reformat.eps,來將 ps 檔案轉換為 eps 檔案格式。
第8行的指令 psbasemap 代表先繪出帶座標的基底圖,可能是x-y軸,可能是經-緯度,端看投影方法參數 -J  來決定。投影範圍參數 -R 為帶入 $range 數值,-B (Border)則為自訂邊框的參數,可分別選擇如下
-Ba1f1:."BATS station vs Event":WNse -P -G100 -K -U "station map"
a1: 為自訂每1度一個標記(annotation)
f1: 為自訂每1度隔一個邊框,如例圖黑色的長度為1度
以上若無分別,則四邊界共同設定值,分別以 : : 來隔開值,底圖的題目為 :. : 內的字串,後利用 W[w]N[n]S[s]E[e]Z[z] 來表示大寫的軸是顯示標註值,小寫則否。
若為 x-y 繪圖,則 -B 的參數值,可給x-y軸不同的說明,格式如下
-Ba1f1:"x axis"::,"time seconds"/a0.1f0.1:"y axis"::,"displacement":."Graph Subject":wnse

     參數 -P 表起始的底圖繪圖為直式(Portrait),參數 -G 為灰階顏色(白-黑-白,0-255),本例為100。參數 -K為允許還可以再執行下一個繪圖指令,參數 -U 為畫上時間戳章,預設位置在圖的左下角,此範例戳章說明為 station map。
第9行的指令 pscoast 為帶入 gshhg 資料來繪海岸線圖一樣給定 -J -R -B 表示不更改繼承 psbasemap 的 -J -R -B 值,參數 -G 一樣為灰階顏色,-A 為海岸線的精細程度,預設全繪出為0/0/4,本例給1000為小於1000 km^2 地區不繪出,愈精細花愈多繪圖時間,但小地區範圍亦較精準。參數 -D 為解析度,可接( full, high, intermediate, low, crude)的第一個字母。參數-K如psbasemap參數-K,僅在繪圖最後一個指令無須加上。參數 -O 為覆蓋圖層,表示一直將圖層堆疊上去作圖。接著要利用接續寫入轉向子 >> 來寫入ps檔。
第10至12行為繪出event 在地圖上的位置,第13至15行則為對event的說明
指令 psxy 可以繪出點,線等標記,而我們可以給定多點座標放在讀入範圍 << EOF 至 EOF(本例只給一點),或是直接給定座標清單檔如 psxy coordinate.xy  ,便可以繪出數點位置。參數
-W 代表繪圖線的寬度,-S 代表符號種類( c 表圓形,d 表菱形,s 表正方形,a 表星形,t 表三角型,其他形狀符號請參閱 man psxy )。其他參數 -G,-K,-O 皆同於 pscoast。
第14行指令 pstext 主要將字串繪成圖型物件,讀入之欄位格式由左而右為座標,接著為字的大小,傾斜角度(0-90,水平到垂直),字型(圖5-5),對齊文字框的([L | C | R](左|中|右) 與 [T | M | B](上|中|下)),最後為要畫的說明字串。
第16-23行如第10-15行的方式,繪入圖形物件,第24-25行僅僅為了要結束gmt 的script 的結尾描述,記得一定要去掉參數 -K ,因為後面沒有繪圖指令了。



















(圖5-5,引自GMT Technical Reference and Cookbook)

     以上皆可利用指令 man psbasemapman pscoastman psxyman pstext 來查詢。想想看,一個完整的繪圖可能需要哪些資訊,才能讓讀圖的人"秒懂",是否要有比例尺?圖的標記代表意義?圖的顏色代表意義?是否需要標示出高度/深度圖?當然,這樣看起來,第一個範例是很陽春的,接著,我們要用 gmt 來表示這些閱讀圖表所需之工具。
為了繪出地形高低,我們接著從美國國家海洋與大氣總署(NOAA) 的https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data(ref. 2),下載1 arc-minute的格點資料(取樣點)(約1.85km,https://en.wikipedia.org/wiki/Minute_and_second_of_arc)ETOPO1的網格數值檔
,我們要利用gmt來繪出地形的樣貌,如圖5-6。圖5-6上圖為修改圖5-4,並給上等高線,與地形並以顏色劃分,加上顏色比例尺,並加入事件(event)發生日期與時間,並給上震源機制解(qCMT)。圖5-6下圖為將上圖轉個視角,來看島周圍的水深深度,筆者爲了要表現清楚的水深,不合等比例來放大Z軸投影。此例地形資料來源為NOAA的ETOPO1(ref. 2),gmt script 內容如圖5-7。








































(圖5-6)







































(圖5-7)

    第2至6行爲設定gmt基本參數,便不再贅述。第7至11行是先設定會用到的變數值,此例多了viewaz與cptfile,分別用來改變視角與塗色色盤。各種 cpt 顏色檔,亦可以從 cpt-city 網站 http://soliton.vm.bytemark.co.uk/pub/cpt-city/views/totp-cpt.html 來取得,或是如18,19行的指令皆可產生。筆者爲了展示,故已將 ETOPO1.cpt 的選項註解了,且色盤是由網格資料來產生後(指令grd2cpt 參數 -C 爲使用下載或是GMT內建的塗色盤,如圖5-8,位置爲/usr/share/GMT/cpt ),再手動編輯自己要表現的顏色,當然也可以把臺灣着色成”藍瘦香菇“。各顏色R/G/B數字對應,可從 w3school (http://www.w3schools.com/colors/colors_rgb.asp)來求得。
    第12行利用指令 grdcut 來截取要繪圖的範圍,本例爲臺灣南方東經119.5至121.5度,北緯22至24度,使原先的ETOPO_bed_g.g98僅保留該範圍,轉存成 myloc.grd 檔,藉此增快繪圖速度。第14指令grdgradient,爲將要用的grd檔求出光影的角度(參數-A)與強度(參數-N),來產生一個myloc.int的光影強度檔,有興趣的讀者可以利用指令man grdgradient來瞭解細節,用於調整繪grdview與grdimage的光線與陰影。第15至16行,則是利用指令grdhisteqgrdmath,來將grd檔資料梯度變化變較平緩。第20行則宣告cptfile變數值,以利接後面的cpt檔案名稱只需引用cptfile變數值。接着第21至54行,便爲主要的繪圖。






















(圖5-8)

     接下來,筆者想要知道這二十多年來,所觀測到的臺灣附近區域的地震活動度,並繪出投影在地表上的地震事件與投影出地表下的震源位態與機制(Harvard CMT,http://www.globalcmt.org/CMTsearch.html )。如圖5-9-1 與 圖5-9-2。首先我們先至美國地質調查所(USGS)的地震目錄搜尋網站,https://earthquake.usgs.gov/earthquakes/search/ ,來搜尋我們想要知道的年代範圍,規模範圍,區域範圍的地震清單,並下載搜尋結果。再將結果整理成我們所需要的欄位格式,如 緯度  經度  震源深度  規模,各欄位間留空格,或如csv格式給分欄符號 ','(comma) ,此過程就不再贅述。因繪圖版面限制,筆者將與剖面圖分開畫,圖5-9-1之 gmt script 爲圖5-10,圖5-9-2之gmt script爲圖5-11。






































(圖5-9-1)





(圖5-9-2)







































(圖5-10)








未完,待續...
To be continued...


GMT功能相當多樣,也有符合各種地圖繪圖需求的參數,或科學符號特殊字元的參數,若想知道更多GMT操作指令,可參考http://gmt.soest.hawaii.edu/projects/gmt/wiki/Documentation與指令速查http://gmt.soest.hawaii.edu/doc/latest/


References
1. Broadband Array in Taiwan for Seismology, Institude of Earth Science, Sinica (http://bats.earth.sinica.edu.tw/)
2. Amante, C. and B.W. Eakins, 2009. ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis. NOAA Technical Memorandum NESDIS NGDC-24. National Geophysical Data Center, NOAA. doi:10.7289/V5C8276M [access date].
3. USGS Earthquake Hazard Program Website



If you have any feedback or question, please go to my forum to discuss.

這個網誌中的熱門文章

Linux操作不求人 - 伍章之伍 - make 巨集式編譯器

Linux 操作不求人系列 - 貳章之壹 - Shell Script 程式設計(I) - BASH

Linux操作不求人 - 肆章之貳 - 伺服器架設(II) - 郵件伺服器 - postfix 與 dovecot