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. 出力結果

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







GMT 5.3.1のインストール (Cygwin編)

2012年のGMT 4.5.7のインストール (Cygwin編) に多くのアクセスがありますが、現在のGMTの最新バージョンは5系列です。こちらは過去の記事の内容を最新バージョン用にアップデートした記事となります。

GMT (Generic Mapping Tools)とはコマンドラインベースの地図作成ツールで、地球科学分野でよく用いられるツールで、綺麗な空間分布図を描画することができます。クロスプラットフォームなソフトなのでWindows、Mac、Linuxで使うことができます。


今回はWindows 10でのGMTのインストール(&初期設定)について説明します。


1. はじめに


WindowsでGMTを使う場合は、いくつかの選択肢があります。

  1. Bash on Windows でGMTをコンパイル・インストール
  2. 昔ながらのCygwinでGMTを使用

今回は2番目のCygwinを用いた方法を紹介します。
以下のツールをインストールする必要があります。


必須
GSView
GMTで出力したps/epsファイルの表示や変換処理に必要 (※実際に変換しているのはGhostscriptやpstoedit、GSViewはインターフェース)。
SumatraPDF
GMTで出力したps/eps/pdfファイルの表示に便利。Adobe readerで表示すると排他アクセス処理が働くので、こちらのほうが良い。
推奨
Cygwin
シェルスクリプトでGMTコマンドを使う場合必要。WindowsバッチでもGMTは扱えるがCygwinの方が多機能で便利。
Ghostscript
GMTで出力したps/epsファイルをビットマップ画像に変換するのに必要。
PStoedit
GMTで出力したps/epsファイルをベクター形式のままSVG形式に変換できる。Inkscapeなどのイラストソフトを使っている人には便利。
改行コードやエンコーディングが変更できるエディタ
Notopad++EUC-JP対応版がおすすめ、自分はemacs modifiedを使ってるけど…)
異なるOSや他人とのファイルの移動を行う場合は、環境の違いによって日本語コメントなどが文字化けする場合がある。
受け渡し先に合わせたエンコーディングや改行コードに変換したり、人からもらったスクリプトファイルのエンコーディング&改行コードが自分の環境に合わない場合などにあると便利(でも本当はコメントは英語で書いておくのが無難)

2. インストール

GMTサイトのwindows用ダウンロードのページから、最新バージョンのインストーラをダウンロードして、インストールを行います。

ちなみに2016/12/22現在の最新バージョンは5.3.1なので
  • Windows (64 bit)の人は、 gmt-5.3.1_install64.exe
  • Windows (32 bit)の人は、 gmt-5.3.1_install32.exe


GMTの実行ファイルは、 C:\programs\gmt5 フォルダにあります。 昔のGMTではインストール後(というかインストーラが無かったので、圧縮された実行ファイルをダウンロード後、解凍して任意の場所に置くだけ)、パスの設定を行う必要がありましたが、最近のGMTでは自動でパスの設定をしてくれます。

なのでコマンドプロンプトやCygwinターミナルで何かGMTのコマンド、例えば、

gmt psxy

※GMT4まではpsxyなどのコマンドだけで良かったのですが、GMT5よりgmtを最初に付ける必要があります。

を打ち込むと、大量のメッセージが出てきたらインストール成功です。
もし「'gmt psxy'は、内部コマンド~(中略)~として認識されていません」とか「gmt psxy:コマンドが見つかりません」などと出てきたらインストール失敗です。そのようなことは通常無いと思うのですが(少なくとも自分は見たことないです)、万が一そのようなことが起こった場合は自分で環境変数の設定をしてパスを通す必要があります。 必要なユーザー環境変数は、以下のとおりです。

変数
GMT_SHAREDIRC:\programs\gmt5\share
pathC:\programs\gmt5\bin

※設定後は再起動が必要です(たぶん)。


上記のGMT&周辺ツールのインストールが終わったら、Cygwinターミナルで試しに

gmt pscoast -Jm1:30000000 -R120/150/20/50 -Ba10f5g5 -Gtan > sample.ps

と打ってみて下さい。
※表示の都合で2行にまたがっているかもしれませんが、1行で入力して下さい

カレントディレクトリにsample.psというPSファイルができているはずです。このPSファイルをWindowsで見るために、SumatraPDFやGsviewが必要となります。

ちなみにPSファイルをPDFファイルに変換したい場合は
gmt psconvert -A -Tf world.ps

ちなみにPSファイルをPNGファイルに変換したい場合は
gmt psconvert -A -Tg world.ps

を実行するとそれぞれ変換できます。


From Research_public
pscoastコマンドで作成した図

2016年12月19日月曜日

QGIS processing でのRアルゴリズムの自作例 (気象庁の気象観測データのGISへのインポート)

この記事はFOSS4G Advent Calendar 2016 の18日目の記事として書きました。

色んな場所のQGISハンズオンにて、QGIS processing機能による他のソフトウェアとの連携処理 (R言語とかGRASS GISとか)のお話を最近ちょこちょこしています。

Slideshareに資料をあげているものだと


※他にもちょこちょこ


Rの初期設定やRSXファイルとは何かについてはSliedshareを参照していただくとして、ここでは気象庁で公開されている気象観測データをGISに読み込むRアルゴリズムを作ってみたので、RSXの自作の仕方の例として順を追って紹介していきます。


使用データ(気象庁の気象観測データ)について

過去の気象データ検索から日本全国の気象台・アメダス、日時を指定することで該当の気象観測データのページが表示されます。(例. 銚子の2016/12/18の1時間ごとのデータ

似たような過去の気象データ・ダウンロードとは異なり、URLに規則性があったり、10分毎の観測値があったりとプログラミング上、使い勝手が良いので今回は過去の気象データ検索を使っています。


まずは普通にRスクリプトの作成


ざっくりとした処理の流れは、
  1. Webページのテーブル情報の読み取り
  2. 予め調べておいた観測所ごとの緯度経度の紐付け
  3. SpatialPointsDataFrameクラスに変換
  4. GISデータとして出力
となります。Step2の観測所ごとの緯度経度を調べたり、Webページのテーブルの列数が気象台とアメダス(アメダス内でも観測項目の数はいろいろ)で異なり、場合分け処理の記述がめんどいので、今回は関東地方の気象台のみの観測データに対応させてあります。





RSXファイル独自の記述

(詳細はQGIS Training Manual 17.31, 17.32 を参照)


Rスクリプトの冒頭に以下の記述をしています。
スクリプトでは指定した年月日時の観測値を読み取りGISdataに変換しているので、年月日時を入力するように設定しています。

year の2016は初期入力値で変更可能です。
month,day,hour は以下のような記述でドロップボックスから選択するようにすることができます。

##Data downloading=group
##year=number 2016
##month=selection;1;2;3;4;5;6;7;8;9;10;11;12
##day=selection;1;2;3;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18;19;20;21;22;23;24;25;26;27;28;29;30;31
##hour=selection;1;2;3;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18;19;20;21;22;23;24
##Output=output vector


また、スクリプトの最後ではWriteOGR関数によるGISデータの出力はせずに、以下のように冒頭で指定したOutputという名前の変数(SpatialPointsDataFrameクラス)を作成して終わっています。

Output <- SpatialPointsDataFrame(sp,data=sp.dfm)


出力結果

出力されるポイントデータは属性情報として、気象庁の観測データの各項目が入っています(気温、降水量、気圧や風速、風向など)。
※風向は北を0とした時計回りの角度に変換してます。

以下の図はOpenStreetMapの背景図に気温分布のラスターデータと観測地点のシンボルを矢印に変えて風向の値に応じて回転させたものが載せてあります。



とまあRでやっている処理をQGISでGUIで操作できるようにできるというのがprocessingのメリットです。

ぶっちゃけRのプログラミングができる人はそちらで完結したほうが楽なのですが…。同じ処理をプログラミングのできない人とも共有できるのがpurocessingのメリットかなと思います。