2016年12月22日木曜日

GMTによる地図作成の基礎

毎年この時期になるとGMTやLaTeX関連の記事のアクセス数が増加してきて、卒論・修論・D論が架橋に差し掛かっているのだなぁと季節を感じます。(いや、D論はもう方がついていないとやばいですね)


0. 前準備

まず、GMTをインストールします。Bash on WindowsでもCygwinでもLinuxでも構いません。この記事ではとりあえずWindows10でCygwinという環境ということで説明していきます。

データの前処理としてQGISも使いますのでそちらもインストールしてください。


1. データダウンロード

今回はUSGS(米国地質調査所)のSRTM3データ(3秒メッシュ≒100 mメッシュ)を用います。USGSのダウンロードページから、SRTM3/Eurasia とたどり、以下のファイルをダウンロード&解凍して同じフォルダ(例. srtm といったわかりやすい名前のフォルダ)に入れておいてください。

N35E139.hgt.zip
N36E139.hgt.zip
N35E140.hgt.zip
N36E140.hgt.zip


2. QGISでデータの変換
先ほど解凍して同じフォルダに入れたhgtファイルをQGIS上で結合して、GMTで読み込める形式(NetCDF形式)に変換します。

メニューの「ラスター」→「その他」→「結合」と選択して、以下の結合Windowを立ち上げます。

「ファイルの代わりに入力ディレクトリを選択する」にチェックを入れて、入力ディレクトリ(hgtファイルが入っているフォルダ)を指定、出力ファイル名(GeoTIFF形式で)を指定します。
図1. 結合Window

図2. 結合結果(色は見やすく変えてあります。0の値を透過設定)

メニューの「ラスター」→「変換」→「変換(形式変換)」を選択し、出力形式に「GMT NetCDF Grid Format」を選択して保存します。ここではsrtm2_dem.nc というファイル名で保存をしました。


3. GMTで地形表示


Cygwinを立ち上げて、データの置いてあるディレクトリ(=フォルダ)に移動します。

今回使用する地形データは海域に0の値が入っていますので、まずは以下のコマンドでNo data値が正しく設定されたNetCDFファイルを作成します。

gmt grdmath srtm3_dem.nc 0 NAN = srtm3_dem_withnodata.nc

上のコマンドはファイル変換に使うだけですので、最初1回だけ実行すればOKです。
srtm3_dem_withnodata.nc というファイルが作成されます。


GMT上で下記の3行のコマンドを実行します。

gmt makecpt -Cdem1 -T0/3000/100 -Z > sample.cpt
gmt psbasemap -R139/141/35/37 -Jm1:2000000 -Ba1f1g1 -K > output.ps
gmt grdimage srtm3_dem_withnodata.nc -Csample.cpt -Q -R -J -O >> output.ps

図3. 出力結果


4. GMTで地形に他の情報の重ね合わせ

緯度経度の列を持つ下記のようなデータを用意します。(今回は適当に作りました)

図4. ダミーデータ

このようなデータをCSV形式で保存します。


GMT上で下記の4行のコマンドを実行します。
gmt makecpt -Cdem1 -T0/3000/100 -Z > sample.cpt
gmt psbasemap -R139/141/35/37 -Jm1:2000000 -Ba1f1g1 -K > output.ps
gmt grdimage srtm3_dem_withnodata.nc -Csample.cpt -Q -R -J -O -K >> output.ps
gmt psxy data.csv -Sc0.6 -W1,black -Gred -R -J -O >> output.ps


図3. 出力結果

するとこのような図が作成できます。







0 件のコメント:

コメントを投稿