GMT5 Tutorial

GMT5 教程 | 編者: Po-Chin Tseng

等高線圖及剖面

目錄

  1. 總覽
  2. GMT介紹及安裝
  3. 網路資源及配套軟體
  4. 第零章: 基本概念及默認值
  5. 第一章: 製作地圖(地理投影法)
  6. 第二章: XY散佈圖(其他投影法)
  7. 第三章: 等高線圖及剖面
  8. 第四章: 地形圖與色階
  9. 第五章: 地震活動性與機制解
  10. 第六章: 向量與速度場
  11. 第七章: 台灣地理資訊
  12. 第八章: 直方、圓餅、三元圖
  13. 第九章: 三維空間視圖
  14. 第十章: 地質圖

7. 等高線圖及剖面

本章將會先說明GMT所使用的網格檔是什麼、如何製作，接著如何將數值高程以等高線的方式呈現， 最後，當要看一條特定線上的高程變化時，如何利用GMT來繪製該線段的高度變化剖面。

7.1 目的

本章將學習如何繪製

  1. 簡介網格檔(Grid data)
  2. 等高線圖(Contour)
  3. 地形剖面(Profile)

7.2 學習的指令與概念

7.3 簡介網格檔

在開始學習繪製等高線圖之前，要先介紹GMT在繪製網格圖時所使用的檔案格式.grd， 它是屬於netCDF檔案格式，全名被稱做網路通用數據格式(Net Common Data Form)， 由於一般的10進位檔案格式(ascii)，相當地占用硬碟空間，因此常會將資料儲存成2進位的格式(binary)， 這樣可以有效地壓縮檔案所占用的空間。

用實際例子來說明檔案格式的差異，本章將會使用到政府資料開放平台中的20公尺網格數值地形模組資料， 其中的不分福_全台及澎湖的資料，下載下來的.tif檔，透過一般的編譯器打開，如下圖:

顯示出一推看不懂的亂碼，但如何想簡單的得知這個檔案理面的資料資訊，可以使用

gmt grdinfo 檔名

讀取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，透過指令，

gdal_translate -of XYZ [input.tif] [output.xyz]

將.tif檔轉換成.xyz檔。第二，由於內政部提供的檔案其座標系是TWD97， 而GMT繪製地圖慣用在WGS84座標系，因此需要做座標系轉換，透過Python，將檔案轉成WGS84座標系， 完成後再將.xyz檔轉換成.grd檔。

gmt xyz2grd [input.xyz] -R119.0/119.9/21.8/25.4 -I0.6s/0.6s -G[output.grd]

其中

同樣的，也可以透過

gmt grd2xyz input.grd > output.xyz

將網格檔轉回.xyz檔。

透過上述方式，成功地將內政部提供的數值高程檔轉換成GMT繪圖使用的格式，在這個範例中， 一般.xyz的檔案大小約為420MB，轉換成.grd的格式後，變成4.61MB，大大地節省硬碟空間。 使用一樣的做法，把台灣本島的數值高程也轉換成功後，不希望之後畫圖都要分開讀取台灣本島以及澎湖，所以使用:

gmt grdpaste input1.grd input2.grd -Goutput.grd -fg

把兩個擁有共同邊界的網格檔給合併，-f是指定輸入及輸出資料的格式，

當然一次讀取這麼大的數值資料，會耗費很長的時間，因此可以利用GMT提供的網格切割指令:

gmt grdcut input.grd -Goutput.grd -R範圍

將網格檔切割成較小的區域來做使用，接下來，繪製等高線圖的部份，將會使用以花東縱谷做為範圍切割的網格檔。 資料壓縮的方式非常多種，無法在此一一介紹，編者以使用過的經驗作為分享，請見諒。

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

本節學習的新指令:

希望透過兩張圖的比較，能更加知道各指令變化後，繪製出來的差異。從左右兩張的差異來看， 最明顯的莫過於海的部份，由於在數值地形在海的部份，並無資料，而當初在轉檔時，皆將此視作為零， 導致了大量的等高線繪製於此，為了避免此現象發生，可以再畫一次海的顏色，也可以透過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 當在製作小區域地圖時，往往需要利用全區域的小張地圖來框繪出小區域的範圍 ，有兩種方式。

6 製作AA’剖面:

7 同6，換成BB`剖面。

首先AA’剖面切過中央山脈、花東縱谷、海岸山脈，從高層剖面可以清楚地看出三個不同的區域，你也可以嘗試在圖上， 利用pstext加上文字做說明，在約8km處，有一個凹陷處，從等高線來判識，應該是一河谷地形(木瓜溪)。 另一條BB`剖面線是做為對比，展示一般山區的地形。

在一開始設立一個輸入區域，可以方便你做更改，以這範例來說，如果我想換切別條剖面線， 只要在一開始更改經緯度數值，就可以展示出不同的高程剖面，可以自己玩看看。

sed的使用是為了讓剖面高程的多邊形關閉，以方便著色。關於sed以及linux其他指令在windows上如何使用， 可參考5-5，編者另有安裝Cygwin。

7.6 習題

玉山為東亞第一高峰，其海拔3952公尺，其景色優美、氣勢磅礡，吸引大量登山客前往。利用維基百科中玉山主峰以及 其他四個至高點(北峰、東峰、南峰、西峰)，來繪製玉山區域的等高線圖，以及主峰到其餘四峰的高度剖面。

使用的資料檔:

完成圖如下:

7.7 參考批次檔

列出本章節使用的批次檔，供讀者參考使用，檔案路經可能會有些許不同，再自行修改。

上一章 – 下一章