2012年12月27日木曜日

GDAL/OGRコマンドをWindowsでもLinuxライクに使う手順

研究室で購入した衛星データを他のメンバーに使いやすくするために(HDF->GeoTIFF変換や各種派生データ(コンポジット、地形指標など)の作成)、その処理の多くをGDAL/OGRコマンドを使ったシェルスクリプトに記述しています。
それらの処理を研究室の後輩への引き継ぐために、一連の手順の中でも汎用的な部分 (Cygwinのインストール&初期設定) をブログに載せます。

  1. Cygwinのインストール
  2. OSGeo4Wのインストール
  3. Cygwin.batの編集

1. Cygwinのインストール
CygwinはWindows上でLinuxの各種コマンドを使えるようにするツールです。後述のOSGeo4WでインストールされるGDAL/OGRコマンドはCygwinをインストールしなくてもWindowsコマンドプロンプトからでも使用出来るのですが、非常に使いづらいので(※個人の感想です)、Cygwinをインストールします。

以下はインストールの手順です。
Cygwinのサイト ( http://www.cygwin.com/ ) からインストーラをダウンロードして実行します。


インターネット経由でダウンロードが選択されているのを確認して次へ(つまりインストール時にはインタネット環境が必要)。


インストールディレクトリを変更することもできます。が特別な理由がなければデフォルトのままでいいでしょう。日本語フォルダをインストール先にするなどもってのほか。
※Windowsのユーザー名に日本語やスペースが入っていると以前はいろいろエラーが出たらしいです。最近はどうだか知りません。



インターネットからダウンロードしたインストーラの一時的な置き場所。インストール完了後はこのファイルは削除してもOK。このファイルを残しておくと次からはインターネット環境がなくてもこのファイルを呼んでCygwinインストールができる。まあよくわかんない人は消してOK。



で次へ次へと進んでいくとダウンロード元のサーバーを選択する画面がでます。適当に好きなとこを選んで下さい。日本国内のサーバーだとなんとなく早いような気がします(※個人の感想です)。選んだサーバーの状態が良くない場合は別の所を選んでと言われるので、その場合は別の所を選択して次へ進んで下さい。



次にインストールしたいソフトウェアを選択する場面になります。カテゴリ別に多くのソフトウェアが含まれているのでインストールしたいソフトウェアがあれば"Skip"となっているところをクリックして"Install"に変更してください。デフォルトの状態で基本的なソフトウェアはすでに選択されているので特に追加インストールしたいソフトウェアが無ければそのまま次へをクリックして下さい。
※大量のソフトウェアがあり探すのが大変なので検索窓にソフトウェア名を打ち込んで探すと楽だと思います。例えば Imagemagick と打ち込むとGraphicカテゴリだけ表示されます。今回はためしに Imagemagick を追加でインストールしてみましょう。


インストールが完了すると下のような画面が出てきます。アイコンをデスクトップに置きたい場合はチェックを付けてください。


2. OSGeo4Wのインストール
OSGeo4Wのサイト ( http://trac.osgeo.org/osgeo4w/ ) からインストーラをダウンロードして実行します。

インストールモードの選択画面が出てくるので、よくわかんない人は「エクスプレスデスクトップインストール」が選択されているのを確認して次へ進みましょう。



インストールするソフトウェアを選択する画面が出てきます。よくわかんない人は全て選択して次へ進んで下さい。
※GDAL/OGRコマンドだけを使いたい場合はGDALだけにチェックを入れればOKなのですが、QuantumGISはあると何かと便利ですし、QuantumGISを入れると依存関係にあるGRASS GISもインストールされます。uDigとOpenEVはあまり使う機会はないかも(※個人の感想です)。



3. Cygwin.batの編集
実は先ほどまでの説明はググればいくらでも説明しているサイトはあります。ここからが重要な部分です。

Cygwinのアイコンのプロパティを開いてリンク先を見るとわかるのですが、Cygwinの起動はC:\cygwinにあるCygwin.batというバッチファイルで行われています。
CygwinのなかでGDAL/OGRコマンドを使う場合はそのバッチファイル中にパスを通すコマンドを書く必要があります。
※環境変数をいじってもいいのですが別のソフトウェアに悪影響を与える可能性もあるので(特にpythonまわり)、バッチファイルに記述しておくのが安全です。


まずCドライブのcygwinフォルダに移動して、Cygwin.batというバッチファイルで右クリックをして"編集"をクリックします。

バッチファイルには以下のように記述されていると思います。

@echo off

C:
chdir C:\cygwin\bin

bash --login -i



ここに以下の3行を@echoのあとに加えます。
set OSGEO4W_ROOT=C:\OSGeo4W
PATH=%OSGEO4W_ROOT%\bin;%PATH%
for %%f in (%OSGEO4W_ROOT%\etc\ini\*.bat) do call %%f

追加後は以下のような感じになると思います。
※追加行は赤で表示。
@echo off

set OSGEO4W_ROOT=C:\OSGeo4W
PATH=%OSGEO4W_ROOT%\bin;%PATH%
for %%f in (%OSGEO4W_ROOT%\etc\ini\*.bat) do call %%f


C:
chdir C:\cygwin\bin

bash --login -i

で保存すれば、次にCygwinを起動したときからGDAL/OGRのコマンドが使えるようになります。

2012年12月13日木曜日

AGU2012 Fall Meetingで行われたOSGeo Workshop参加報告

FOSS4G Advent Calendar 2012 の13日目の記事です。
http://atnd.org/events/34052

今年のFOSS4G Advent Calendarは今のところなかなかの猛者揃いで、個人的には満足です。自分も当初は「クロスプラットフォームなスクリプト作成法(GDAL、GRASS、GMT、R)」というネタを密かに考えていたのですが、AGUのOSGeoのWorkshopが面白かったのでその紹介をします。当初予定していたネタはまた別の機会に。

  1. はじめに 
  2. ハンズオンの概要 
  3. WELD (Web-Enabled Landsat Data)について 
  4. ハンズオンの内容 
  5. おわりに 

1. はじめに

12月3日(月)から7日(金)の間サンフランシスコで行われたAGU (American Geophysical Union: 米国地球惑星科学連合) Fall Meetingという、地球惑星科学系では最大の研究集会(参加者数20000名以上)に参加してきました。自分はそちらでポスター発表を行なってきたのですが、今回は個人的に注目すべきイベントとしてアメリカのOSGeoメンバーによる「オープンソースソフトウェアを用いたLandsatデータ解析」というハンズオンがありました。通常のセッション修了後の6:00 pm~10:00 pmと遅い時間だったのですが、早めの夕食を取り後輩と2人で参加してきました。

ハンズオンの会場は、AGUのメイン会場であるMoscone Centerではなく、すぐ近くのMarriott Marquisという高級そうなホテルの会議室で行われました。講師を務めたのはサウスダコタ大学の人(David RoyGiuseppe Amatulli)で、アシスタントは4名、参加者は意外と多く20数名はいました。WELDの紹介はDavid Roy教授が行なっていましたが、その他は全てDr. Giuseppe Amatulli(←じーだる発音派)が行なっていました。
 

2. ハンズオンの概要

このハンズオンは以下のような構成になっていました。

<プレゼンテーション>
6:00 pm - 6:10 pm ワークショップの概要説明とUSBドライブでのデータ(VirtualBoxイメージとインストーラ(Windows、Mac))配布。
6:10 pm - 6:30 pm WELDプロジェクトの説明 (Roy)
6.30 pm - 7.00 pm VirtualBoxのインストールと機能の説明、Linuxのコマンドの説明 (Amatulli)
7.00 pm - 7.30 pm 使用するオープンソースソフトウェアの紹介:GDAL/OGR、R、OFT、PKTOOLS (Amatulli)
7.30 pm - 7.45 pm 休憩

 <ハンズオン>
7.45 pm - 8.15 pm OpenEVでWELDデータの可視化 (Amatulli)
8.15 pm - 8.30 pm WELDデータのグリッド単位での時系列変化について (Amatulli)
8.30 pm - 9.30 pm WELDデータの森林火災エリアのマッピング:自動分類や変化の検出 (Amatulli)
9.30 pm - 9.45 pm WELDデータの大気補正前と後の比較 (Amatulli)
9.45 pm - 10.00 pm オープンディスカッション (Roy)

会場に入るとまずVirtualBoxインストーラー&イメージの入ったUSBディスクのコピーをするように言われます。VirtualBoxイメージはzip圧縮された状態で約3GB(解凍後は約8GB)で、このイメージを起動すると今回使用するオープンソースソフトウェア、使用データ、マニュアルの含まれたXubuntuが起動します。自分はこのzipファイルをWindows 7 上で解凍する際に何故か「600PB(ペタバイト)以上の空き領域が必要です」とのエラーメッセージが出て解凍できず(無理w)。デュアルブートのUbuntuで解凍して事なきを得ました。その他にも小さい危機(xubuntuのキーボードがUS設定なので記号入力に最初戸惑うなど)がありましたが、無事最後まで脱落せずについていくことができました。

前半(プレゼンテーション)では、WELDプロダクトの説明や、使用するオープンソースソフトウェアの紹介をしていました。使用したソフトウェアについては、GDAL/OGRRについてはよく知っていたのですが、OFTPKTOOLSは今回のハンズオンで初めて知りました。また、後半のハンズオンではLinuxのコマンドを多く使用するため、初歩的なLinuxコマンドの紹介(ls、less、head、awkなど)もしていました。さらりとawkを紹介している辺りに後半は荒れるなという予感がしました。

後半(ハンズオン)の大まかな流れは、WELDデータを使って森林火災エリアを含むシーンの教師付き分類をSVM(サポートベクターマシン)法とRandomforest法で行い、各プロセスが終わるごとにOpenEVあるいはQGISで結果の確認を行い、最後におまけとしてWELD version 1.5(大気補正前)とversion 2.0(大気補正後)の違いを紹介していました。


3. WELD (Web-Enabled Landsat Data)について

WELDはLandsat 7 ETM+ オルソ補正画像をモザイクしてシームレスに作成した画像プロダクトで、土地被覆分類や、様々な地表面解析を目的として、USGSとサウスダコタ大学の共同プロジェクトで作成されたプロダクトです。対象領域はアメリカ本土とアラスカのみで、30m解像度で提供されています。

複数時期のモザイクしたプロダクトのため、各グリッドに用いられた元データの時期がわかる補助データも付随しています。特定の時間での地表面状態を知りたいという研究にはあまり向いていないのですが、地表面被覆の大まかな分布がわかれば良い研究分野にとってはハンドリングが楽なデータであるといえます。US限定ですが。


4. ハンズオンの内容

4-1. OpenEVでWELDデータの可視化
gdal_mergeによりコンポジット画像を作成しOpenEVでヒストグラムの調整などをして可視化。その後はgdal_calcでログスケール変換した各バンドを同様にコンポジット→可視化。

gdal_calcはほんとうに便利で簡単なバンド演算はほとんどこれで可能。自分が従来GRASSで行なっていた単純なバンド演算のいくつかも既にgdal_calcに移行しています。

4-2. WELDデータのグリッド単位での時系列変化について
WELDではサイト上で緯度経度で指定したグリッドにおける放射量の時系列情報を得ることもできます。ハンズオンでは、CSVデータとして用意されたあるグリッドにおける時系列データをheadコマンドで確認をしつつ、awkによるデータ整形、Rによる時系列プロットを行なっていました。自分はRは普段良く使っているので問題なかったのですが、Rの未経験者にとってはなかなかしんどそうなところでした。一回目の山場です。

4-3. WELDデータの森林火災エリアのマッピング:自動分類や変化の検出
このステップは今回の二回目の山場で、脱落者も続出していました。隣席のお姉さんも途中で「オー・・・。」とかいって大きなため息をついて脱落してました。自分もついていくのが精一杯で助ける余裕はありませんでした。

まずはPKTOOLSを使って多バンドのコンポジット作成(pkcropコマンド)とマスク画像の作成(pkgetmask)。gdal_rasterizeで教師データとしての森林火災域(とおもわれる領域)のポリゴンをラスタライズを前準備として行いました。ちなみにPTOOLSはGDALライブラリを使ってC++で実装されたコマンド群です。PKTOOLSコマンドの多くはすでにあるような処理ばかりでしたが…、SVM分類などが簡単にできるのは便利だと思いました。

その後1つ目の教師付き分類(SVM)の処理に入ります。ラスタライズした教師領域の0.01%を処理時間節約のためにランダム抽出して教師データとして利用してpkclassify_svmコマンドでSVM分類。できたものはOpenEVで可視化して確認を行いました。

続いて2つ目の教師付き分類(Randomforest)の処理に入ります。ここが非常に長くて複雑。流れとしてはまずOFTを使ってK-means NN教師なし分類でセグメンテーション分類したうえで、教師データを用いてR上でRandomforestの教師付き分類を行います。ただひとつ残念なのがRでのRandomforest分類の結果は直接ラスターイメージとして出力されるのではなく、初めのセグメンテーション分類IDとRandomforest分類の対応を示すテキストファイルとして出力されるため、そこからラスターイメージに持っていくために処理が複雑化しているという部分です。この部分の個別の処理を知りたい方はリンク先の説明を参照下さい。

4-4. WELDデータの大気補正前と後の比較
ここではPKTOOLSのpkinfoコマンドで待機補正前と後のRGBバンドの画素値を、それぞれ個別のテキストファイルに出力して、それをgnuplotにて可視化して比較を行いました。

4-5. オープンディスカッション
最後にこのハンズオンセッションについて、難易度や作業量はどうだったか参加者の意見を聞いていました。自分が聞き取れた範囲では、少し作業量が多杉という意見がちらほら出ていました。


5. おわりに

このハンズオンでは多くのオープンソースソフトウェアを使用するため、VirtualBoxを使ってあらかじめ主催者側が環境設定したイメージを利用しました。おかげで複数のソフトウェアを組み合わせて一連の解析を行うという実践的な処理が可能となり、個人的には非常に勉強になるハンズオンだと感じました。

ちなみに今週土曜日に中部大学で第3回FOSS4Gツール勉強会@名古屋を行います。こちらは全然スパルタではないので安心して気軽にご参加ください。
https://sites.google.com/site/foss4gnagoya/foss4gnagoya20121215

あとAGUの通常のセッションの参加報告を研究室のブログにあげていますのでこちらも興味がありましたらどうぞ。
http://hello.ap.teacup.com/snowman/1602.html

2012年10月31日水曜日

FOSS4G2012(Tokyo/Osaka)のお知らせ

オープンソースGISソフトウェアの祭典、FOSS4G2012が11月に東京と大阪で行われます。 昨年はいち参加者として、初めてFOSS4G2011 Tokyoに行きましたが、今年は運営側として初めての参加となります。

FOSS4G Tokyo (11/3, 4, 5)

FOSS4G Osaka (11/7, 8)


2010年度までは自分の博士論文の執筆に追われて余裕がなく、このようなコミュニティの活動にも興味を持ってはいたのですが、まずはD論書かなきゃヤバイ(汗)ということで、密かに?FOSS4G活動(といってもGRASS GISやGDAL、QGIS、Rで解析処理をするだけでアウトプットは皆無)をしていました。


そんなこんなで蓄積した自分の知識・経験をコミュニティに還元しつつ、それを武器として人脈を広げてみようと思い、Nagoya.RというR言語の勉強会が名古屋大学で不定期に行われていたので、そちらでRによるGISやリモートセンシング解析の話をしたりWebに発表資料の公開などを始めました。するとそのWebの資料を見つけて、以前衛星ワークショップや学内勉強会などで少しだけ知り合いだった平松さんにRをつかったGISとリモセン解析に興味を持って連絡をいただき、各種FOSS4Gツールについての勉強会面白そうだねということでFOSS4Gツール勉強会@名古屋の話が立ち上がりました。


そして立ち上げにあたって、FOSS4Gの本丸にていろんな方に宣伝したりご助言などを頂けたらと思い、FOSS4G2011 Tokyoへと参加して、気がついたら中の人になってました。FOSS4G2011 Tokyoに参加して以来、多くの方々と知り合うことができたため、この1年間は自分にとっては実りのあるとしでした(プライベートでも節目の年でした)。


最後に宣伝を、11/4と11/8には様々なFOSS4Gソフトウェアを使った地形解析処理の紹介の話をします。また、11/5にはQGISのハンズオンの講師を、11/8にはRでリモートセンシングハンズオンの講師をつとめます。これを機会に研究&業務利用のGIS及びリモートセンシング環境をオープンソース環境への切り替えを!

2012年6月22日金曜日

Octaveを使った地上開度および地下開度の計算

地上開度/地下開度という地形量

DEM(Digital Elevation Model:デジタル標高モデル)から計算できる地形量には様々なものがあります。有名なものとしては「傾斜」、「斜面方位」などが挙げられ、他にも斜面が凸か凹かの具合を示す「曲率」などもあります。
今回とりあげる地上開度というのは、ある地点で空を見上げたとき、どれくらい開けているかを示すものです。谷底では地上開度は小さくなり、尾根上では大き くなります。また、同じ谷底地形の中でも、V字谷のようにすぐそばに斜面が迫っている谷では小さく、U字谷のように谷底に広い平坦地が広がっている谷では V字谷と比べて大きくなります。

地上開度/地下開度のアルゴリズム

Prima et al. (2006)で地上開度の計算アルゴリズムが説明されています。細かい説明はこの論文をご覧いただくとして、ざっくり説明すると以下のとおりです。
まず計算するにあたって探索半径というものを決めます。そのあとDEMの各グリッドにおいて、東西南北とそれらの中間の計8方位の探索半径内のグリッドの 標高値と中心グリッドの標高値から仰角を計算して8方位それぞれにおける最大仰角を計算します。そしてその8つの値から平均最大仰角を求めて、最終的に天 頂から平均最大仰角を引いたものが地上開度となります。
ちなみに地下開度の場合は仰角を俯角に、天頂を天底(←この日本語初めて知りました、英語ではNadir)に置き換えます。

Octaveで計算してみよう

前置きが長くなりましたが、Octaveで地上開度および地下開度を計算するOctaveプログラムをこちらに公開します。
  • openness_calc.m : 地上開度/地下開度のメインの計算部分
  • gdalread.m : GeoTIFFなどのGDALライブラリで対応しているラスターデータをOctaveに読み込む関数プログラム
  • gtiffwrite.m : 計算結果をGeoTIFFに書きだす関数プログラム
※2012-6-22 17:10 追記
gdalreadとgtiffwriteはGDAL/OGRのgdal_translateというプログラムを利用しているため、OSGeo4Wというインストーラを利用してGDAL/OGRをインストールしておく必要があります。
OSGeo4W (http://trac.osgeo.org/osgeo4w/)


使用方法:
  1. 上記3つのプログラムと計算したいDEM(例. aster2003_dem30.tif)ファイルを同じディレクトリにおきます。※DEMは投影座標じゃなきゃダメ。
  2. Octaveを起動して、そのディレクトリに移動します。
  3. Octaveコンソールで、
openness_calc('aster2003_dem30.tif', 300, 'up')

というようにコマンドを打ちます。上のオプションは1つめがファイル名、2つめが任意の探索半径、3つめが地上開度(up)か地下開度(down)かの選択フラグです。

結構計算時間はかかります。フリーソフトのSAGA GISでも開度計算はできるらしいのでそちらの方が楽かも(汗)。上記ソースコードをもとに計算アルゴリズムをチューニングして遊びたいという方はこのOctaveプログラムは参考になるかもしれません。


引用文献
Prima, O. D. A., Echigo, A., Yokoyama, R., & Yoshida, T. (2006). Supervised landform classification of Northeast Honshu from DEM-derived thematic maps. Geomorphology, 78(3-4), 373-386. Elsevier Science Bv. doi:10.1016/j.geomorph.2006.02.005

2012年3月29日木曜日

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

2016/12/22現在の最新バージョンはGMT 5.3.1です。そちらに対応させた記事を書いていますので、この古い記事ではなくそちらを参考にすることをおすすめします。

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



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


1. はじめに


WindowsでGMTを使う場合は、GMTそのもの以外にも以下のツールが必要となってきます。


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

2. インストール

GMTサイトのwindows用ダウンロードのページから、

  • Windows (64 bit)の人は、 gmt-4.5.7_install64.exe
  • Windows (32 bit)の人は、 gmt-4.5.7_install32.exe

をダウンロードして、インストールを行います。

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

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

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

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


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

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

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

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

すでにGsviewをインストールしていれば、EPSファイルに関連付けられているはずなのでダブルクリックで作成した図を見ることができると思います。 Adobe Illustratorをインストールしている場合はそちらに関連付けられているかもしれません。

その場合は図を見るだけならGSview、編集する場合はIllustratorと使いわけると良いと思います(Illustratorの起動はGSviewより重いため)。

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