TOPEXをGRASSで算出してみる

もくもく単純データ処理に嫌気がさしてきたので
気分転換を兼ねて,
以前相談を受けていて放置していた
TOPEX(Topographic Exposure)の算出を
Grassをつかってやってみました。

GrassでTOPEXを計算,ていうのはこの方がやっていたので参考にしました。
Calculating Togographic Exposure With Grass

もともと俯角はすべて0°とするものだったらしいけど
負の値として計算することもあるようですね。
たとえばこれ(P47-48)
EMPIRICAL MODELLING OF WINDTHROW RISK USING GEOGRAPHIC INFORMATION SYSTEMS

DEMを基にしたTOPEXは
こんな感じで計算するようです。

1) 対象セルA00から8方位の計測線を延ばす
2) その延長上かつ特定の計測範囲内において
一定の間隔でA00からの高低角を算出する
(下図では　A00-A33 , A00-A55
およびA00-A77 間の高低角を算出する)
3) 求めた高低角の最大値を
その方位の高低角と定義する
4) 各方位で同様の計算を行い
8方位すべての値を積算したものをTOPEXとする

↓かなりアレな図ですが(平面図です)背景のメッシュはDEMだと思ってください…

どうやら地上開度とアイディアはほとんど同じ。
開度による地形特徴の表示

とりあえず

1)Grass を外部からうごかす環境設定済み
2)DEMはGrassにインポート済み
3)俯角もそのまま計算する(負の値として処理)

という前提で進みます。

引数多すぎとか8方位なのにいらん計算してるとか
そもそもGeotiffのインポートからやれよとか
いろいろあるけど気分転換なので　気にしない　気にしない

ここからダウンロードして解凍
topog フォルダをフォルダごと
windowsであれば
C:\Python27\Lib\
にコピーしてください。

それで

python (command line)
または
コマンドプロンプト(最初にpythonと入力しENTERする)

などで以下のようにコマンドを入力して実行してみてください。

# モジュールのインポート
from topog import topex as t
# t.topex("gisdataディレクトリ","ロケーション","マップセット",
#　　　　　　"入力DEM名",DEMの解像度(m),計測間隔(m),計測範囲(m),"出力ラスタ名")
# 入力例
t.topex("C:\\temp\\grass","hokk","PERMANENT","dem",10,100,2000,"topex_raster")

# モジュールのインポート

from topog import topex as t

# t.topex("gisdataディレクトリ","ロケーション","マップセット",

#　　　　　　"入力DEM名",DEMの解像度(m),計測間隔(m),計測範囲(m),"出力ラスタ名")

# 入力例

t.topex("C:\\temp\\grass","hokk","PERMANENT","dem",10,100,2000,"topex_raster")

topex.pyの中身はこうなっています。

## import module ##
import grass.script as grass
import grass.script.setup as gsetup
import math
import os

def topex(gisdb,location,mapset,dem,res,interval,range_max,output):

    # set up grass path #
    gisbase = os.environ['GISBASE']
    # set GRASS gisdatabase, location, mapset #
    gsetup.init(gisbase, gisdb, location, mapset)
    # set region as same as input DEM (extent and resolution)
    grass.run_command("g.region", rast=dem)
    # calculate pixel based interval
    interval_p = int(round(interval / res))
    # calculate pixel based maximum range
    range_max_p = int(round(range_max / res))
    
    ## set formula for calculation ##
    # set initial part of formula
    formula_t1 = output + " = "
    # 8 direction as azimuth angle
    for azi in (0,45,90,135,180,225,270,315):
        # initial cell number
        num = 0
        formula = "max("
        while num &lt; range_max_p:
            num = num + interval_p
            x = int(round(math.sin(math.radians(azi))*num))
            y = int(round(math.cos(math.radians(azi))*num))
            formula = formula + "atan((" + dem + "[" + str(y) + "," + str(x) + "] - " + dem + ") / (" + str(res) + " * " + str(num) + ")),"
        # build formula for calculating max elevation angle at each direction    
        formula_t = formula[0:-1] + ")"
        # build formula for summing up value at all direction 
        formula_t1 = formula_t1 + formula_t + " + "
    # build formula for r.mapcalc
    formula_f = formula_t1[0:-3]    
    print formula_f    
    ## run r.mapcalc command ##
    grass.run_command("r.mapcalc_wrapper",expression = formula_f)
    
    print "finish"
    
if __name__ == '__main__':
    print "this is code block"

## import module ##

import grass.script as grass

import grass.script.setup as gsetup

import math

import os

def topex(gisdb,location,mapset,dem,res,interval,range_max,output):

# set up grass path #

gisbase = os.environ['GISBASE']

# set GRASS gisdatabase, location, mapset #

gsetup.init(gisbase, gisdb, location, mapset)

# set region as same as input DEM (extent and resolution)

grass.run_command("g.region", rast=dem)

# calculate pixel based interval

interval_p = int(round(interval / res))

# calculate pixel based maximum range

range_max_p = int(round(range_max / res))

## set formula for calculation ##

# set initial part of formula

formula_t1 = output + " = "

# 8 direction as azimuth angle

for azi in (0,45,90,135,180,225,270,315):

# initial cell number

num = 0

formula = "max("

while num < range_max_p:

num = num + interval_p

x = int(round(math.sin(math.radians(azi))*num))

y = int(round(math.cos(math.radians(azi))*num))

formula = formula + "atan((" + dem + "[" + str(y) + "," + str(x) + "] - " + dem + ") / (" + str(res) + " * " + str(num) + ")),"

# build formula for calculating max elevation angle at each direction

formula_t = formula[0:-1] + ")"

# build formula for summing up value at all direction

formula_t1 = formula_t1 + formula_t + " + "

# build formula for r.mapcalc

formula_f = formula_t1[0:-3]

print formula_f

## run r.mapcalc command ##

grass.run_command("r.mapcalc_wrapper",expression = formula_f)

print "finish"

if __name__ == '__main__':

print "this is code block"

環境変数の設定とかめんどくせ　という場合は
同梱してあるtopexf.pyを使って

1)出力コマンドをテキストファイルに出力
2)grassのレイヤマネージャにあるコマンドコンソールにコピペして実行

てのがてっとりばやくて楽かもしれないですね。
これはpython2xがあれば動きます。

# モジュールのインポート
from topog import topexf as t
# t.topexf("入力DEM名",DEMの解像度(m),計測間隔(m),計測範囲(m),"出力ラスタ名","出力テキストファイル(フルパス)")
# 入力例
t.topex("dem",10,100,2000,"topex_raster01","C:\\temp\\topex_command.txt")

# モジュールのインポート

from topog import topexf as t

# t.topexf("入力DEM名",DEMの解像度(m),計測間隔(m),計測範囲(m),"出力ラスタ名","出力テキストファイル(フルパス)")

# 入力例

t.topex("dem",10,100,2000,"topex_raster01","C:\\temp\\topex_command.txt")

topexf.pyの中身はこうなっています。

## import math module ##
import math

def topexf(dem,res,interval,range_max,output,out_formula):

    # calculate pixel based interval
    interval_p = int(round(interval / res))
    # calculate pixel based maximum range
    range_max_p = int(round(range_max / res))
    
    ## set formula for calculation ##
    # set initial formula
    formula_t1 = output + " = "
    
    for azi in (0,45,90,135,180,225,270,315):
        # initial cell number
        num = 0
        formula = "max("
        while num &lt; range_max_p:
            num = num + interval_p
            x = int(round(math.sin(math.radians(azi))*num))
            y = int(round(math.cos(math.radians(azi))*num))
            formula = formula + "atan((" + dem + "[" + str(y) + "," + str(x) + "] - " + dem + ") / (" + str(res) + " * " + str(num) + ")),"
        # build formula for calculating max elevation angle at single direction    
        formula_t = formula[0:-1] + ")"
        # build formula for summing up value at all direction 
        formula_t1 = formula_t1 + formula_t + " + "
    # build formula for r.mapcalc
    formula_f = formula_t1[0:-3]    
    print formula_f    

    f = open(out_formula, 'w')
    f.write("r.mapcalc '" + formula_f + "'")
    f.close()
    
    print "finish"
    
if __name__ == '__main__':
    print "this is code block"

## import math module ##