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

2 件のコメント:

  1. 参考にさせてもらってます。
    python用に移植できないかとソースコードを読ませてもらいましたが、upのコードが違う気がします。四方位?

    返信削除
    返信
    1. 読んでくださってありがとうございます!もうだいぶ昔に書いたコードなのでよく覚えていないところもありソースコードも汚くて読みにくいかもしれません(^_^;)。upについて、downとはコードの記述の順番が異なり分かりにくいのですが、8方位を計算してあります(for ループの中)。また何かわからないことがあればお気軽に聞いて下さい!

      削除