2017年1月19日木曜日

GMTによる散布図作成の基礎

GMTというと地図作成ソフトで知られますが、普通のグラフも作ることができます。
今日は簡単な散布図作成につい紹介します。

GMTのインストールは既に完了しているものとします。


1. ファイルの説明

以下のファイルをダウンロードして同じフォルダにおきます。


sample_xy.sh はGMTのコマンドが並べられたシェルスクリプトで、メモ帳などのエディタで開くことができます。ただし改行コードがUnix(LF)となっているためWindowsのメモ帳では1行に表示されてしまいます。メモ帳ではなくNotepad++などの正しく様々な改行コードを認識できるエディタを使うことを進めます。

data1.csv、data2.csv は適当に作成した2変数のスペース区切りデータで、以下のような内容です。
0 1000
-5 900
-18 800
-14 700
-15 600


2. GMTコマンドの説明

sample_xy.shの中身は以下のとおりです。

順番に説明をすると、
1行目はシェルスクリプト(bash)のおまじない※無くても良い

4行目はグラフの枠とdata1.csvのプロット、それぞれのオプションの意味は
-Jはグラフの横と縦の長さ、縦をマイナスにすると数値の大小が逆向きになる。
-Rは値の範囲、X軸は-30から30、Y軸は500から1000
-BxはX軸の目盛間隔とラベル
-ByはY軸の目盛間隔とラベル
3つめの-Bは目盛をふる側(WとSだけ大文字なので左と下の目盛をふる)とタイトル
-Wは線の太さと色(ここでは黒色)と種類(実線)
-Vはターミナルに情報表示
-Pはグラフの向きを縦に
-KはGMT描画コマンド(psxy,pscoast,psbasemapなど、psconvertは関係ない)が2つ以上の場合、最後のコマンド以外につける必要がある。「まだ続きがあるよ」という意味

6行目はdata2.csvのプロット、それぞれのオプションの意味は
-Jは2回目なので省略可能
-Rも2回目なので省略可能
-Wで線の太さと色(ここでは赤色)と種類(破線)
-Vはターミナルに情報表示
-Pはグラフの向きを縦に
-OはGMT描画コマンド(psxy,pscoast,psbasemapなど、psconvertは関係ない)が2つ以上の場合、最初のコマンド以外につける必要がある。「前のコマンドの続きだよ」という意味
-Kは先程と同じく「まだ続きがあるよ」という意味

8-14行目で凡例の指定

17行目でpsファイルをPDFファイルに変換
18行目でpsファイルをPNGファイルに変換


3. 実行
ファイルのある場所に移動してから、
bash sample_xy.sh

とターミナルにコマンドを入力します。
sample_xy.shをテキストエディタで開いて一行一行ターミナルにコピペしてもできますが、上のようにコマンドを入力したほうが簡単です。


出来上がったグラフはこんな感じ




※2017/1/21 線の種類指定の追加と凡例のプロットを追加
※2017/2/6 ℃記号を修正しました(^_^;)。

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のメリットかなと思います。






2016年9月2日金曜日

GMT 5.2.1のインストール(Bash on Ubuntu on Windows 編)

Bash on Ubuntu on Windows (長い…、以下BoUoW)の登場でWindows上でのUnix系の処理のやり方が大きく変わりつつあるのでGMTのインストール方法の情報を更新。
(以前掲載したGMT 4.5.7のインストールのしかたをベースに作成しています。)


GMT (Generic Mapping Tools)とはコマンドラインベースの地図作成ツールで、地球科学分野でよく用いられるツールで、綺麗な空間分布図を描画することができます。クロスプラットフォームなソフトなのでWindows、Mac、Linuxで使うことができます。
追記:2016年9月2日現在GMT本体の最新バージョンは5.2.1です。海岸線データ(GSHHG)の最新バージョンは2.3.6です。後に出てくるインストールスクリプトのバージョン番号を自分が使用するバージョン番号に置き換えてください。

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


1. はじめに

Bash on Ubuntu on Windows のインストールは下記サイトを参照
Bash on Ubuntu on Windowsをインストールしてみよう!

作成したps/pdfファイルの表示はSumatraPDFがおすすめ。
Adobe Readerとかだとファイルを開くとロックが掛かり更新ができない。

また、改行コードやエンコーディングが変更できるエディタ(Notepad++など)も必要。

また、GMTインストールにあたって以下のパッケージが必要となりますので、apt-get でインストールしてください。

build-essential
make
cmake
netcdf-bin
libnetcdf-dev
gdal-bin
libgdal-dev


2. GMTのインストール

BoUoWにて、下記スクリプトを実行してください。



処理内容の説明

Line 4-5:
バージョン番号は最新のものに置き換えてください。

Line 10-14:
ファイルをダウンロードして展開。

Line 19-22:
海岸線データをインストールするための処理

Line 25-30:
インストール処理(時間かかるかも)



3. GMTの動作確認

上記のGMTのインストールが終わったら、BoUoWで試しに

gmt pscoast -R0/360/-80/80 -Jm1:300000000 -Ba30df15dWeSn -A10000 -Dl  -W0.5 -P > world.ps

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

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

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

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

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



pscoastコマンドで作成した図


※BoUoWのホームディレクトリは以下の場所となります。
(nuimuraをあなたのユーザー名に置き換えてください。)
C:\Users\nuimura\AppData\Local\lxss\home\nuimura


また、別のやり方としてはsources.list を強制的にxenialに書き換えて、16.04にアップグレードしてからソースからではなくapt-getからgmt5をインストールするやり方もあるそうですが、そのやり方は私もまだ試したことがありません。試される場合は自己責任でお願いします。


※GMTコマンドや説明の誤り、16.04にアップグレードの情報などを @Say_no さんにいただき修正・追記いたしました。どうもありがとうございましたー。


2016年2月3日水曜日

SRTM GL1 (1秒メッシュ) の欠損値処理

他大の学生さんから解析手法について問い合わせがあったのでついでにブログ投稿。


地形データを利用する様々な分野の多くの研究者によって使われているSRTM DEMの前処理についてです。

2000年2月のスペースシャトルエンデバーにより計測されたSRTM DEM は3秒メッシュ (約90 m) 解像度 (SRTM3) については当初から全世界のDEMが公開されていましたが、1秒メッシュ(約30 m) 解像度は米国内の領域しか公開されていませんでした。

2014年には米国以外の領域の1秒メッシュ解像度DEMも公開され (SRTM GL1) ダウンロードができるようになっていますが、利用者の利便性を考えてなのか、本来計測できていない部分の標高値も他のDEM (ASTER GDEM、GMTED、NED) を使って補間してしまっています。

手っ取り早く高解像度のDEMを使いたい人にとっては問題ないのですが、地形や氷河体積の時間変化などを明らかにしたい人にとっては、2000年時点 (スナップショット) での標高値ではない標高 (例えばASTER GDEM はグリッドによって2000年~2011年と異なる) が混ざるため困ります。

そこで標高の時間変化計算に用いるためにはSRTM GL1から外来DEMのグリッドを除去する処理が不可欠となります。SRTM GL1には対応するタイルごとにSRTMGL1Nというフラグデータがあり、SRTM GL1のグリッドごとに、外来DEMを使ったかどうか使った場合はどの種類かがわかる情報があります。この情報を元にオリジナル のSRTMを使っているグリッドのみを抽出すればOKです。

QGISでの読み込み

以上の処理はRやMATLABなどのプログラミングでももちろんできますが、ここではQGISを使ってやり方を紹介します。

SRTM GL1Nの仕様のページを見ると、値は8 bit (0-255) で格納されており、201以上の値を使えば大丈夫そうです。SRTM GL1のフォーマット (*.hgt) は単純な raw binary形式ですが、QGISはSRTMのファイル名(ex. NaaEbbb.hgtの場合は北緯aa度東経bbb度を北西端)から自動で認識して正しい位置に表示してくれます。

SRTM GL1Nも単純な raw binary 形式となっていますがこちらはQGISが自動認識はしてくれないためヘッダーファイル (*.hdr) を用意して位置情報などのデータに関する上をQGISに教えてやる必要があります。


hdrファイルとは

テキスト形式でデータについての情報を記述したものです。 raw binaryファイルと拡張子より前の部分の名前を揃える必要があります (ex. データがN27E086.numの場合はN27E086.hdr)。

詳しくはこちらのページのEHdr -- ESRI .hdr labelled を参照
http://www.gdal.org/frmt_various.html


N27E086.num の場合は、N27E086.hdrという名前で以下の内容が記述してあるテキストファイルを用意すればOKです。


BYTEORDER      I
LAYOUT         BIL
NROWS          3601
NCOLS          3601
NBANDS         1
NBITS          8
BANDROWBYTES   7202
TOTALROWBYTES  7202
PIXELTYPE      UNSIGNEDINT
ULXMAP         86
ULYMAP         28
XDIM           0.000277777777777778
YDIM           0.000277777777777778
NODATA         0


一応こちらからダウンロードできます。適当なエディタで開いてみてください。
https://dl.dropboxusercontent.com/u/870568/blogger/N27E086.hdr


準備ができたら N27E086.num をQGISにドラッグ&ドロップすれば表示されるはずです。あとはQGISのラスターメニューの変換コマンド (または右クリックで別名保存) でGeoTIFF形式に変換すると他のソフトウェアでも簡単に読み込むことができます。
SRTM GL1のファイル名と同じだとややこしいので N27E086_num.tifとして保存します。


SRTM GL1からオリジナルのSRTMのみの抽出

QGISにSRTM GL1 (N27E086.hgt) とSRTM GL1N (N27E086_num.tif) の2つを読み込んで置いて下さい。

まずラスタ計算機 (Raster calculater) で以下の内容を図のようにウィンドウ下部の演算欄に入力します。

"N27E086@1" * ("N27E086_num@1" >= 201)



フラグデータで201以上の値が入っている場所のSRTM GL1のみを抽出するという意味です。ウィンドウ右上の出力ファイル名も指定しておいて下さい。
出力ファイル名を N27E086_temp.tif とします。

変換して出力されたファイルは以下の図のようにSRTMじゃない部分 (つまり欠損値) に数値の0が入っています。



標高変化の計算のためには0じゃなくてNo Dataだと認識されて欲しいのでもうひと手間かけます。

N27E086_temp.tif を右クリックで別名保存で、ウィンドウ下部のNo data valueで0を指定します。出力ファイル名を指定してOKをクリックします。
ここでは出力ファイル名をN27E086_void.tif とします。



以下のような感じでできあがりです。

2016年2月2日火曜日

R言語(ggplot2)による気温変動グラフの作成

毎年(といっても去年から)後期授業も一段落した2月ぐらいに、ゼミの学生たち+α対象に勉強会を開いています。

去年はGISの勉強会統計(Excel)の勉強会を行いましたが、今年はそれに加えてプログラミング(R言語)の勉強会も追加しようと考えています。

でいろいろネタ探しをしていたら、気象庁のWebページで以下のような報道発表があり、グラフに使われたデータもありました(CSVなどではなくWebページ上の表ですが)。



報道発表のリンク先にある世界の年平均気温の偏差のページにあるような図をちょっとRのggplot2ライブラリで作成してみました。





以下のような図が出来上がります。
 気象庁のページに似せて作った図。
世界全体での年平均値(黒丸)、回帰直線(青)および5年移動平均(赤)


世界全体(赤)、北半球(緑)および南半球(青)の5年移動平均


Githubリポジトリはこちら