等高線圖及剖面
7. 等高線圖及剖面
本章將會先說明GMT所使用的網格檔是什麼、如何製作，接著如何將數值高程以等高線的方式呈現， 最後，當要看一條特定線上的高程變化時，如何利用GMT來繪製該線段的高度變化剖面。
7.1 目的
本章將學習如何繪製
- 簡介網格檔(Grid data)
- 等高線圖(Contour)
- 地形剖面(Profile)
7.2 學習的指令與概念
grdcontour: 繪製等高線圖grdcut: 切割網格檔grdinfo: 從網格檔中獲取資訊grdpaste: 合併兩個有共同邊界的網格檔grdtrack: 對網格檔取樣grd2xyz: 轉換網格檔至ascii檔案格式gmtinfo: 從表格資料中讀取資訊project: 將資料投影至線、大圓或轉換座標系psbasemap: 繪製圖框、刻度、標籤等等pscoast: 繪製海岸線xyz2grd: 轉換ascii檔案格式至網格檔gdal_translate轉檔程式
7.3 簡介網格檔
在開始學習繪製等高線圖之前，要先介紹GMT在繪製網格圖時所使用的檔案格式.grd， 它是屬於netCDF檔案格式，全名被稱做網路通用數據格式(Net Common Data Form)， 由於一般的10進位檔案格式(ascii)，相當地占用硬碟空間，因此常會將資料儲存成2進位的格式(binary)， 這樣可以有效地壓縮檔案所占用的空間。
用實際例子來說明檔案格式的差異，本章將會使用到政府資料開放平台中的20公尺網格數值地形模組資料， 其中的不分幅_全台及澎湖的資料，下載下來的.tif檔，透過一般的編譯器打開，如下圖:
顯示出一推看不懂的亂碼，但如何想簡單的得知這個檔案理面的資料資訊，可以使用
讀取Penghu_20m.tif，可得到資訊如下:
Penghu_20m.tif: Title: Grid imported via GDAL
Penghu_20m.tif: Command:
Penghu_20m.tif: Remark:
Penghu_20m.tif: Pixel node registration used [Cartesian grid]
Penghu_20m.tif: Grid file format: gd = Import/export through GDAL
Penghu_20m.tif: x_min: 76661.4764748 x_max: 121061.476475 x_inc: 20 name: x nx: 2220
Penghu_20m.tif: y_min: 2563952.79393 y_max: 2633852.79393 y_inc: 20 name: y ny: 3495
Penghu_20m.tif: z_min: -999 z_max: 70.5 name: z
Penghu_20m.tif: scale_factor: 1 add_offset: 0
+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=GRS80 +units=m +no_defs
從這些資訊可以得知標題、X、Y、Z軸的範圍和間距，放大倍率(scale_factor)等等的訊息。 但這檔案要提供給GMT使用，會遇到兩個問題，第一，檔案格式為Geotiff， 利用GDAL(地理空間資料抽象化文庫, Geospatial Data Abstraction Library) 提供的轉檔程式gdal_translate，透過指令，
將.tif檔轉換成.xyz檔。第二，由於內政部提供的檔案其座標系是TWD97， 而GMT繪製地圖慣用在WGS84座標系，因此需要做座標系轉換，透過Python，將檔案轉成WGS84座標系， 完成後再將.xyz檔轉換成.grd檔。
其中
-R是給定檔案的X、Y軸範圍。-I給與X、Y軸的間隔，後面的英文字代表單位，可參考4-4距離的單位。- m對應角分。
- s對應角秒。
- e對應公尺。
- k對應公里。
-G給予輸出檔名稱。- 另外常用到的就是
-A處理多數據重疊在某點時的方法: - -Af選第一個落在重複座標上的數值做為該座標值。
- -As選最後一個落在重複座標上的數值做為該座標值。
- -Al選最小...
- -Au選最大...
- -Ad選最大最小的差值...
- -Am選平均...
- -Ar選RMS(方均根)...
- -As選標準差...
- -An選出現次數...
- -Az選總合...
- 以及
-D寫檔頭資訊: - +xx軸名稱，格式為名稱 [單位]，例如distance [km]。
- +y、+z同上。
- +s讀取網格數據後的放大倍率，預設值為1。
- +o讀取網格數據後要加上的常數，預設值為0。
- +n指定特定數值為無效值，預設為NaN。
- +t網格檔的標題。
- +r網格檔的註解。
同樣的，也可以透過
將網格檔轉回.xyz檔。
透過上述方式，成功地將內政部提供的數值高程檔轉換成GMT繪圖使用的格式，在這個範例中， 一般.xyz的檔案大小約為420MB，轉換成.grd的格式後，變成4.61MB，大大地節省硬碟空間。 使用一樣的做法，把台灣本島的數值高程也轉換成功後，不希望之後畫圖都要分開讀取台灣本島以及澎湖，所以使用:
把兩個擁有共同邊界的網格檔給合併，-f是指定輸入及輸出資料的格式，
* -ff給予浮點數格式
* -fg地理座標系統
* -fc卡式座標系統
當然一次讀取這麼大的數值資料，會耗費很長的時間，因此可以利用GMT提供的網格切割指令:
將網格檔切割成較小的區域來做使用，接下來，繪製等高線圖的部份，將會使用以花東縱谷做為範圍切割的網格檔。 資料壓縮的方式非常多種，無法在此一一介紹，編者以使用過的經驗作為分享，請見諒。
7.4 等高線圖
等高線或稱等值線圖，是將資料中等值的資料點連接起來，借此了解地勢的高低起伏，其中每條線的疏密程度，
則可用來判斷坡度的緩急，辨識河谷、山稜線、峭壁等等的地形。設定範圍在121.33/121.68/23.55/24.1，
使用grdcut將網格檔切割成較小的區域。
花東縱谷由菲律賓海板塊與歐亞板塊擠壓而形成，原本的谷地受到河川的侵蝕，形成廣闊的沖積扇平原， 成為了花東一帶，最富庶的地區。在平原中，寬廣的視野，觀看東西兩側(中央山脈、海岸山脈)的高聳入雲， 不由得讚嘆大自然的神奇，就讓我們來看看其高程的變化吧！
使用的資料檔: - 花東縱谷數值地形檔
成果圖
批次檔
set ps=7_4_east_rift_valley.ps
# left figure without clip
gmt pscoast -R121.33/121.68/23.55/24.1 -JM10 -BWeSn -Bxa.2 -Bya.2 ^
-Df -W1 -G194/250/216 -S175/243/255 -K > %ps%
gmt grdcontour east_rift_valley.grd -R -JM -C250 -A1000+f12p -K -O >> %ps%
# right figure with clip
gmt pscoast -R121.33/121.68/23.55/24.1 -JM10 -BWeSn -Bxa.2 -Bya.2 ^
-Df -W1 -G194/250/216 -S175/243/255 -X13 -K -O >> %ps%
gmt pscoast -R -JM -Df -Gc -K -O >> %ps%
gmt grdcontour east_rift_valley.grd -R -JM -C250 -A1000+an+f12p+g255/153/199 ^
-Gd15c -Wc.5,255/110/110 -Wa1,180/13/13 -Q180 -L0/2000 -K -O >> %ps%
gmt pscoast -R -JM -Q -K -O >> %ps%
gmt psxy -R -JM -T -O >> %ps%
gmt psconvert %ps% -Tg -A -P
本節學習的新指令:
- 第6,12行:
grdcontour透過網格檔繪製等高線，語法是grdcontour 網格檔 -C間隔值，其中最常用到的選項: - -A繪製有註解的等高線:
- +a註解的角度，+an垂直等高線；+ap平行等高線(預設值)
- +f註解字體的屬性
- +g字體框的底色
- +p字體框框線的屬性
- -C繪製無註解的等高線，有三種方式:
- 給一個數字，會以此繪製固定間隔的等高線
- 使用.cpt檔，等高線會依照色階檔的設定來繪製，未來將介紹
- 使用一個檔案，第一欄是高度，第二欄是C或A
- -G註解的間隔設定:
- 依照距離來取決間隔，d用圖紙的單位；D用地圖的單位
- f給一個檔案ffile.d，依照檔案中的位置，繪製註解
- 依照線段來取決，l起始/結束，用,來區分每條線，L使用大圓路徑
- n給一條線中等距標籤的數量
- -L低值/高值，限制範圍。
- -W等高線的屬性，-Wa設定有註解的等高線；-Wc設定無註解的等高線。
- -Q省略資料點小於多少的等高線。
- 第11行:
pscoast -Gc開啟海岸線切割模式。 - 第13行:
pscoast -Q關閉海岸線切割模式。
希望透過兩張圖的比較，能更加知道各指令變化後，繪製出來的差異。從左右兩張的差異來看，
最明顯的莫過於海的部份，由於在數值地形在海的部份，並無資料，而當初在轉檔時，皆將此視作為零，
導致了大量的等高線繪製於此，為了避免此現象發生，可以再畫一次海的顏色，也可以透過pscoast -Gc及-Q，
來限制等高線繪製的範圍。
-A、-C、-W，其詳細的設定方式，像是+f、+p都和前兩章提到的方式一樣，透過上述的簡介
配合指令檔，就可以理解左右兩圖的差異。唯-Q是較特別的選項，有時後等高線的資料點數不多，
造成繪製在圖上是一點一點的，有礙觀瞻，透過限制最小資料點數，將可避免這現象的發生。
7.5 地形剖面
在上一節中，介紹了呈現在平面上的等高線圖(X、Y方向)，那如果想看特定某一條線段，其X、Z方向的變化， 本節將介紹如何利用GMT來完成繪製地形剖面。
成果圖
批次檔
# 1. set two profile start point and end point
set ps=7_5_elevation_profile.ps
set Alon1=121.38
set Alat1=24.07
set Alon2=121.5
set Alat2=23.6
set Blon1=121.54
set Blat1=24.06
set Blon2=121.47
set Blat2=24.0
# 2. contour basemap
gmt pscoast -R121.33/121.68/23.55/24.1 -JM10 -BWeSn -Bxa.2 -Bya.2 ^
-Df -W1 -G194/250/216 -S175/243/255 -K --MAP_FRAME_TYPE=plain > %ps%
gmt pscoast -R -JM -Df -Gc -K -O >> %ps%
gmt grdcontour east_rift_valley.grd -R -JM -C250 -A1000+f10p ^
-Wc.5,150 -Wa1,30 -Q180 -K -O >> %ps%
gmt pscoast -R -JM -Q -K -O >> %ps%
# 3. AA' line
echo %Alon1% %Alat1% > tmp
echo %Alon2% %Alat2% >> tmp
gmt psxy tmp -R -JM -W2,238/91/78 -K -O >> %ps%
echo %Alon1% %Alat1% A > tmp
gmt psxy tmp -R -JM -Sc.8 -G238/91/78 -W1 -K -O >> %ps%
gmt pstext tmp -R -JM -F+f16p,0,blue -K -O >> %ps%
echo %Alon2% %Alat2% A' > tmp
gmt psxy tmp -R -JM -Sc.8 -G238/91/78 -W1 -K -O >> %ps%
gmt pstext tmp -R -JM -F+f16p,0,blue -K -O >> %ps%
# 4. BB' line
echo %Blon1% %Blat1% > tmp
echo %Blon2% %Blat2% >> tmp
gmt psxy tmp -R -JM -W2,234/235/128 -K -O >> %ps%
echo %Blon1% %Blat1% B > tmp
gmt psxy tmp -R -JM -Sc.8 -G234/235/128 -W1 -K -O >> %ps%
gmt pstext tmp -R -JM -F+f16p -K -O >> %ps%
echo %Blon2% %Blat2% B' > tmp
gmt psxy tmp -R -JM -Sc.8 -G234/235/128 -W1 -K -O >> %ps%
gmt pstext tmp -R -JM -F+f16p -K -O >> %ps%
# 5. insert map
gmt pscoast -R119.9/122.1/21.8/25.4 -JM3 -Bwesn -Ba -Df -W1 -S255 -G230 ^
-X7 -K -O --MAP_FRAME_TYPE=inside >> %ps%
gmt psbasemap -R -JM -D121.33/121.68/23.55/24.1 -F+p2 -K -O >> %ps%
# 6. AA' profile
gmt project -C%Alon1%/%Alat1% -E%Alon2%/%Alat2% -G0.1 -Q | ^
gmt grdtrack -Geast_rift_valley.grd > tmp
gmtinfo tmp -i2,3 -I1/10 > tmp1
set /p pr=<tmp1
gmtinfo tmp -i2 -C -o1 > tmp1
set /p md=<tmp1
sed -i '1i %Alon1% %Alat1% 0 0' tmp
sed -i '$a %Alon2% %Alat2% %md% 0' tmp
awk "{print $3, $4}" tmp | gmt psxy %pr% -JX12/6 -W2 -G238/91/78 -X5 -Y.5 -K -O >> %ps%
gmt psbasemap -R -JX -BwESn+t"AA' Profile" -Bxa+l"Distance (km)" ^
-Bya+l"Elevation (m)" -K -O --FONT_TITLE=24p,0,blue >> %ps%
# 7. BB' profile
gmt project -C%Blon1%/%Blat1% -E%Blon2%/%Blat2% -G0.1 -Q | ^
gmt grdtrack -Geast_rift_valley.grd > tmp
gmtinfo tmp -i2,3 -I1/10 > tmp1
set /p pr=<tmp1
gmtinfo tmp -i2 -C -o1 > tmp1
set /p md=<tmp1
sed -i '1i %Blon1% %Blat1% 0 0' tmp
sed -i '$a %Blon2% %Blat2% %md% 0' tmp
awk "{print $3, $4}" tmp | gmt psxy %pr% -JX -W2 -G234/235/128 -Y9 -K -O >> %ps%
gmt psbasemap -R -JX -BwESn+t"BB' Profile" -Bxa -Bya+l"Elevation (m)" -K -O >> %ps%
gmt psxy -R -JX -T -O >> %ps%
gmt psconvert %ps% -Tg -A -P
del tmp*
本次批次檔篇幅較長，將利用註解號碼來依序解說，本節學習的新指令:
1 設定兩條剖面線的經緯度，漸漸去習慣把一些常需要更改的值，設定為變數，並在批次檔的開頭， 規劃一個區域去編輯這些變數。
2 透過上一節學到的技巧，製作一張等高線圖底圖。
3 呼叫AA'線段經緯度，畫出線及加上註解。
4 呼叫BB'線段經緯度，畫出線及加上註解。
5 當在製作小區域地圖時，往往需要利用全區域的小張地圖來框繪出小區域的範圍 ，有兩種方式。
psbasemap-Dx軸最小/x軸最大/y軸最小/y軸最大。-D參考點模式，可參考6-6極軸中pslegend -D。- +wx軸寬度/y軸寬度，外框的寬度
- +j指定insert map的參考點
- +ox軸偏移量/y軸偏移量
- +s檔名，輸出左下角點位及長寬
- +t將繪圖原點改至insert map的左下角
-Fd、l、t，對應地圖邊界(-D)、比例尺(-L)、方向標(-T)， 預設值是三者設定。- +c小圖與外框間隔
- +g填滿顏色
- +p外框線屬性
- +r圓角矩形外框
- +s加陰影
6 製作AA'剖面:
project將表格式資料投影到一條線上、大圓上、別的座標系上，此區只介紹投影至線上。-C起始點座標-E結束點座標-F輸出的選量，輸出可選xyzpqrs(預設值)，代表X、Y、Z軸(以度為單位)；如果搭配使用-G， 輸出成rsp(以公里為單位)-G取樣距離-Q地圖的單位，轉換成KMgrdtrack對網格檔做指定位置的取樣-G被取樣的網格檔- 透過
project和grdtrack搭配使用，將資料(x軸 y軸 距離 高度)輸出成暫存檔(tmp)。 gmtinfo從表格資料中抓取訊息。-A當有多個檔案的使用:- -Aa輸出全檔案的範圍
- -Af每個檔案分開輸出範圍
- -As當檔案屬於多重區塊檔案時，輸出每個區塊的範圍
-C輸出最大最小值，配合-o來指定輸出的欄位。-D配合-I修改輸出的範圍-E輸出最大或最小值:- i最小值；I絕對值後的最小值
- h最大值；H絕對值後的最大值
-I類似四捨五入的方式取-R需要的範圍，給定x軸/y軸取捨的標準，1表示個位數，10表示10位數等等。-i選擇輸入的欄位，0表示第一欄-o選擇輸出的欄位，0表示第一欄- 透過
gmtinfo輸出範圍訊息至暫存檔(tmp1)，接著使用set /p pr=<tmp1將範圍訊息設為變數pr使用。 sed資料處理工具，詳細介紹請參考鳥哥的Linux私房菜。- -i: 直接修改讀取的檔案內容
- a: 插入字串在目前的下一行
- i: 插入字串在目前的上一行$符號表示行尾， 以及最後一行中的*符號，用來表示任意一個字元與重覆字元， 都可以從鳥哥對正規表示法(Regular Expression)的介紹中學習。
--FONT_TITLE改變標題的文字屬性
7 同6，換成BB`剖面。
首先AA'剖面切過中央山脈、花東縱谷、海岸山脈，從高層剖面可以清楚地看出三個不同的區域，你也可以嘗試在圖上，
利用pstext加上文字做說明，在約8km處，有一個凹陷處，從等高線來判識，應該是一河谷地形(木瓜溪)。
另一條BB`剖面線是做為對比，展示一般山區的地形。
在一開始設立一個輸入區域，可以方便你做更改，以這範例來說，如果我想換切別條剖面線， 只要在一開始更改經緯度數值，就可以展示出不同的高程剖面，可以自己玩看看。
sed的使用是為了讓剖面高程的多邊形關閉，以方便著色。關於sed以及linux其他指令在windows上如何使用，
可參考5-5，編者另有安裝Cygwin。
7.6 習題
玉山為東亞第一高峰，其海拔3952公尺，其景色優美、氣勢磅礡，吸引大量登山客前往。利用維基百科中玉山主峰以及 其他四個至高點(北峰、東峰、南峰、西峰)，來繪製玉山區域的等高線圖，以及主峰到其餘四峰的高度剖面。
使用的資料檔: - 玉山區域網格檔的製作可參考7-3的方式，或者直接使用前面連結下載。
完成圖如下:
7.7 參考批次檔
列出本章節使用的批次檔，供讀者參考使用，檔案路經可能會有些許不同，再自行修改。