傾斜量の算出方法

土木分野では，地形解析と呼ばれる解析あります．具体的には，標高から傾斜量や曲率などの地形的情報を算出することであり，これらの地形的情報を元に，斜面防災やハザードマップ作成などの取り組みが行われています．今回は，標高から傾斜量を算出する方法を紹介します．最後には，傾斜量を計算するためのpythonのプログラムも掲載しています．

1. 傾斜量とは

ある地点における傾斜（水平に対する角度）を示す数値です．文献によっては斜面勾配と記載されている場合もあります．傾斜量をQGISやArgGISなどのGISソフトで可視化した図を傾斜量図と呼びます．傾斜量図については，国土地理院の説明が非常に分かりやすいです．この記事では，この傾斜量図を作成ために必要な傾斜量の計算方法を紹介します．

2. 傾斜量の計算方法

傾斜量の計算方法は，これまでに様々提案されています．以下は，傾斜量の計算方法を記載している文献です．今回は佐藤らの論文に記載されている計算方法で傾斜量を算出してみます．

• 佐藤丈晴，中島翔吾，"大規模崩壊の兆候となる微地形の抽出手法－天川村における評価事例－"，日本地すべり学会誌，Vol.53，No.3, pp.29-33，2015

• 神谷泉，田中耕平，長谷川裕之，黒木貴一，早田 靖博，小田切聡子，政春尋志，"傾斜量図の作成とその応用"，情報地質，Vol.10，No.2，pp.76-79，1999

• 独立行政法人土木研究所　土砂管理研究グループ　地すべりチーム，地すべり地における航空レーザー測量データ解析マニュアル，土木研究所資料，平成21年6月

3. 傾斜量の算出

3.1. 標高データの取得

兵庫県が1mピッチの標高データを公開しているため，ここから兵庫県北部のある地域の標高データをダウンロードします．ダウンロードしたファイルは，XY座標と標高が記載されているテキストファイルなので，QGISなどでGeoTIFFファイルに変換します．下図は，QGISで，GeoTIFFファイルに変換した標高データを表示した図です．白い箇所ほど標高が高いことを示しています．

3.2. Pyhonで標高データから傾斜量を算出

佐藤らの論文に記載されている傾斜量の計算方法は，

Slop = \frac{360}{2\pi} \times \arctan \Biggl( \sqrt{\biggl(\frac{Z_{E}-Z_{W}}{2\Delta x}\biggl)^2 + \biggl(\frac{Z_{N}-Z_{S}}{2\Delta y}\biggl)^2} \Biggl)

となっています．ここで，Zは標高値であり，下付き文字のアルファベットは方角を示しています．また，$\Delta x$と$\Delta y$は，それぞれX方向とY方向のデータの物理的距離の間隔を示しています．今回利用している標高データのグリッドピッチが1mなので，$\Delta x=1$と$\Delta y=1$となります．
平面で表現されている標高データは，プログラム上で単なる2次元配列となるため，この式は計算対象のするセル（インデックス）の上下左右の標高データを参照して計算すればいいことが分かります．したがって，プログラムの流れとしては，

• 標高データのGeoTIFFファイルの読み込み
• 佐藤らの計算方法で，標高データから傾斜量を算出
• 算出した傾斜量をGeoTIFFファイルで保存

となります．以下が，この流れをPythonで記載したプログラムです．GeoTIFFファイルの読み込み・書き込みはGDALを使用しています．CUIで，実行時引数に「標高データのGeoTIFFファイルパス」と「出力する傾斜量のGeoTIFFファイルパス」を与えることで，傾斜量のGeoTIFFが出力されます．ただ，2重ループを含むため計算時間は若干長いです．C言語で書けば早いです．

import argparse, math
import numpy as np
from osgeo import gdal, gdalconst

def main():
    ExeArgParser = argparse.ArgumentParser()
    ExeArgParser.add_argument('TiffPath', type=str)
    ExeArgParser.add_argument('OutputSlopPath', type=str)
    Args = ExeArgParser.parse_args()

    # ndarray形式の標高とグリッドピッチの取得
    TiffData = gdal.Open(Args.TiffPath, gdalconst.GA_ReadOnly)
    Z = TiffData.GetRasterBand(1).ReadAsArray()
    GridPitch = abs(TiffData.GetGeoTransform()[1])
    print(Z.shape, GridPitch)

    # 傾斜量の算出
    SlopGrad = np.zeros(Z.shape)
    for i in range(1, Z.shape[0]-1):
        for j in range(1, Z.shape[1]-1):

            x = (Z[i, j+1] - Z[i, j-1]) / (2 * GridPitch)
            y = (Z[i+1, j] - Z[i-1, j]) / (2 * GridPitch)
            SlopGrad[i, j] = (360 / (2* math.pi)) * math.atan(math.sqrt(x * x + y * y))

    # 四捨五入で，小数点以下2桁に丸める
    SlopGrad = np.round(SlopGrad, decimals=2)

    # TIFFとしてデータを保存
    Dtype, NumBand = gdal.GDT_Float32, 1
    OutputTif = gdal.GetDriverByName('GTiff').Create(Args.OutputSlopPath, SlopGrad.shape[1], SlopGrad.shape[0], NumBand, Dtype)
    OutputTif.SetGeoTransform(TiffData.GetGeoTransform()) # 座標系指定
    OutputTif.SetProjection(TiffData.GetProjection()) # 空間情報を結合
    OutputTif.GetRasterBand(1).WriteArray(SlopGrad) # 傾斜量を設定
    OutputTif.FlushCache() # ディスクに書き出し
    OutputTif = None

if __name__ == "__main__":
    main()

3.3. 算出した傾斜量

下図が，傾斜量のGeoTIFFファイルをQGISで表示した図になります．黒いほど傾斜がきつい（角度が大きい）部分になります．この傾斜量図を見ると標高データの図より，山間部の部分が明確に分かるようになりますね．