跳到主要內容

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 操作不求人系列 - 貳章之貳 - Shell Script 程式設計(II) - BASH 與 TCSH / CSH

     在上章,我們介紹很多bash shell指令的應用方式,並讓它們變成 script,在這章此節,我們要承襲上節,繼續討論 bash shell script 的其它程式設計概念,與讓程式可重複使用的方法,就是利用函式(function)。      首先,我們來創作一個判斷是否為閏年的函數,script的名稱就取為 check_year.sh ,請自行將其設為可執行。程式碼如圖2-6,為了解說方便,筆者利用指令 nl check_year.sh 將程式碼包含行數印出(圖2-7),其他除顏色外,都與圖2-6同。以下 export 宣告環境變數,自訂變數無需加宣告保留字。(2021.10.07更正export說明,以下自訂變數宣告請直接去除export開頭) (圖2-6) (圖2-7)        圖2-7第2-6行,與之前的範例相似,皆有防堵參數個數輸入錯誤的判斷。第7-11行為接著判斷輸入的年分,是否為真的正整數,也就是大於零的數字。其中第7行可解釋為,利用正規表示式搜尋 $2 字串值得頭至尾部的字元,皆由 0-9 組成,若有,則為真(True)會進入 if 內的陳述執行, 但我們想要的,應該是僅要字串其中一字元為非正整數,便進入if 內的警告並跳出 。故,筆者在判斷式前多加一個 ! ,代表著當字元完全是正整數時,就不要執行  if  內陳述,而直接往第12行執行,但若其中有一個字元為非正整數,則會進入 if 內印出錯誤訊息並跳出 Script 。在此例使用者輸入非正整數等字串(如:12ab、cde、1a1b),便會出現錯誤訊息"Error Value",並跳出 Script 。而第8行的判斷式,效果跟第7行相同,但只能在BASH 3.0才能支援,故筆者故意保留,讓讀者可以學到另一種表示方式。         第12-28行,為宣告一個函式 leapyr () ,在 BASH Script內若要使用含式,必須在使用之前先建立函式的功能,如函式建立在第12-28行,則若要呼叫使用(Call)函式,則必須在第29行之後才能呼叫,並且可重複呼叫。第13行為定義函式呼叫時,一定會有一個外部參數, 此外部參數非該 Script 的外部參數 ,而是由程式撰寫者給入,故在此無設定防呆判斷。第14-

Linux操作不求人 - 伍章之陸 - Intel PXE 與利用 tftp + anaconda kickstart來自動部署系統

     現代的個人電腦與筆電爲了節能省碳,往往皆省略了光碟機的設置。雖然說 USB隨身碟亦可應用於安裝作業系統,但若遇到機房的大量機器需安裝之問題,亦不適合此種用光碟片或隨身碟逐臺安裝的方式,或是需準備多片光碟或隨身碟。故爲了大量安裝與部屬 Linux作業系統的機器,我們便可透過具有PXE功能之網路卡(https://en.wikipedia.org/wiki/Preboot_Execution_Environment),其所具備的網路開機的能力,來作透過網路方式來安裝與大量部署 CentOS Linux作業系統。其原理 wikipedia  的PXE說明,都解釋得很清楚,筆者不需要再 " 掉書包 ",在開發者大神面前班門弄斧,我們就直接來做做看。      首先,我們一樣透過 yum來安裝所需的 tftp伺服器軟體(http://www.jounin.net/tftpd32.html),利用PXE支援 tftp 透過 udp 埠號 69連線,來提供檔案的能力,給利用 PXE開機後安裝作業系統的機器,圖5-60   (圖5-60) 接著修改 tftp 設定檔 /etc/xinetd.d/tftp,如圖5-61,請記得 disable  要改為 no (圖5-61)      若要修改 server_args 參數為自訂的 tftp root 路徑,記得要變更 SELinux 的權限設定,如要改為/tftpboot,則使用指令 chcon  來調整,如以下: $ mkdir /tftpboot $ chcon  --reference /var/lib/tftpboot  /tftpboot 接著將 xinetd 與 tftp 加入開機啟動,並重啟動 xinetd,如下: $ chkconfig  --level  235  xinetd on $ chkconfig  --level  235  tftp  on $ service xinetd restart 開放防火牆通行 $ iptables -A INPUT -p udp --dport 69 -j ACCEPT $ service iptables re

Linux操作不求人- 肆章之壹 - 伺服器架設(I) - SSH(SFTP、SCP)、FTP伺服器與遠端連線

@ ssh, sftp      通常安裝好 CentOS6_x64 作業系統後,sshd, ssh daemon 的服務功能是預設開啟的,如筆者的前面篇章所述,預設的 iptables 防火牆設定,亦是開啟讓 ssh 的連線是可通過的,不僅可以連出,也可以被連入。若要確認是否有安裝 sshd 套件以及在啟動時的 runlevel 2 3 5 是否有被載入,可以使用以下指令搜尋: $   rpm  -qa  |   grep  openssh     #  ssh 與 sshd 連線服務皆由 openssh 應用軟體提供。 或使用以下 $   rpm  -qa  |   grep  ssh         #  比用關鍵字 openssh 搜尋更模糊,故符合的條件更多。      如圖4-1可以查詢到有關於 openssh-client 與 openssh-server 的套件, openssh-server 便是提供連入服務的軟體,openssh-client 為提供可以連出的工具。若無以上套件,則使用 yum install openssh ,則可下載安裝。 (圖4-1) 再接著輸入以下指令查詢到 sshd 這個服務,是否有再開機程序內載入,如圖4-2 $ chkconfig | grep ssh   # 查詢 sshd 是否有於 開機 runlevel 啟動 再利用以下指令,來查詢是否防火牆有允許連線 $ iptables  -L  |  grep ssh   # 出現如圖4-2 允許通過之條件 (圖4-2)      接下來,我們要先來調整 sshd 的設定檔,利用 vim  /etc/sshd/sshd_conf ,如圖4-3-1與4-3-2。因為設定檔參數很多,筆者為方便說明,將 /etc/sshd/sshd_conf 檔案內容分成兩張圖。 (圖4-3-1)