エクセルで作る地形解析図（1）
１．地形解析図の例そもそも、地形解析とは何か、というのは気が向いたらまた後日。このページはけっこう適当なので、国総研資料、第204号の資料１を先に見ておいた方が間違いないかもしれないことをお断りしておく。 ＜追記 2022.08.28＞ だいぶあとになってから、ここに書いていることがまとめられている論文を見つけた。相当検索していたはずなのに、どうして今まで引っかからなかったのだろう、、、 * 太田、八戸（2006）：「数値標高モデルによる地形計測の現状と応用例」、応用地質、Vol.46、No.6、pp.347-360 下図は利尻島。左から陰影図＋段彩図、傾斜量図、地上開度図。 
国土地理院基盤地図情報（10mDEM）および同閲覧コンバートソフトを使用
以下はいずれもエクセルで作れる図。図によって、例えば計算をする際に傾斜量で緩傾斜部に重みをつけるなどのことをする方法があるが、その中には計算をしなくても閲覧コンバートソフトでわりと自在に変化させられるものもある。 傾斜量図、 斜面勾配図、 斜面方位図、 デコボコ図、 斜面曲率、 地形変換線、 尾根・谷、水系密度、 texture図、 陰影図、 地上開度図、地下開度図、 尾根谷度図、 Convergence Index、 接峰面図、接谷面図、起伏量図、 起伏指標、 ラプラシアン、エッジ、 高度分散量図、高度分散異常図 
【 傾斜量図 】地形の傾斜の程度を隣接セルとから ”傾斜量”という値として算出したもの。国土地理院のこちら参照。計算式についてはいくつか提案されている。 Sx＝{(H11-H13)/(2*Dx) + (H21-H23)/(2*Dx) + (H31-H33)/(2*Dx)}／3 ＝{ H11+H21+H31−(H13+H23+H33)}／(6×Dx) Sy＝{ H11+H12+H13−(H31+H32+H33)}／(6×Dy) S＝sqrt（Sx^2+Sy^2） なお、次の斜面勾配との違いも含めて、作図に関して後述。 * 神谷泉、田中耕平、長谷川裕之、黒木貴一、早田靖博、小田切聡子、政春尋志（1999）：「傾斜量図の作成とその応用」、情報地質、vol.10、No.2、pp.76-79
【 斜面勾配図 】傾斜量と同義であるが、XY方向に隣接する4近傍に重みをつけ、傾斜角を斜面勾配I（度）として求めたもの。この算出法はGISソフトとしておそらく最も使われているArcGISで用いられている（ArcGISのヘルプより）。Sx＝{(H11＋2×H21＋H31)−(H13＋2×H23＋H33)}／(8×Dx) Sy＝{(H11＋2×H12＋H13)−(H31＋2×H32＋H33)}／(8×Dy) I ＝arctan { (Sx)^2＋(Sy)^2 }^0.5×180／π 記号は傾斜量図と同じ 【斜面方位図】
Hx＝{ H11+H21+H31−(H13+H23+H33)}／(3×Dx) Hy＝{ H31+H32+H33−(H11+H12+H13)}／(3×Dy) θ＝arctan（Hx / Hy）×180／π ※この計算式の場合、角度と方角との関係は 右表のようになる。 ちなみに、ArcGISの数式はこれが使われている（ArcGISのヘルプより）。 Hx＝{ H13+2*H23+H33−(H11+2*H21+H31)}／8 Hy＝{ H31+2*H32+H33−(H11+2*H12+H13)}／8 θ＝arctan（-Hx / Hy）×180／π 同じ箇所の段彩陰影図と斜面方位図を並べて示す。この斜面方位図に陰影はつけていないが、この例では北を白、南を黒としている。このため北から光を当てている陰影図のようにも見える。カラー設定に示すように西を茶、東を緑としているが、南北の白黒を逆に設定すると、尾根谷が逆転して見える。方向依存性のためと思われる。 
国土地理院基盤地図情報（10mDEM）および同閲覧コンバートソフトを使用
また、上図の左端上半は河床幅が広くなりかなりの緩勾配と思われるが、斜面方位図だと勾配にかかわらず方位がしっかり出ている。右端上半の緩斜面の微地形もよくわかる。等高線図や陰影図からだと読み取りにくい極緩勾配斜面の傾斜方向がはっきりわかるのもこの図の特徴と思う。 
【 デコボコ図 】地形の凹凸を表現するものとして作成してみた。最初は斜面の曲率計算をしたくて計算式を探していたのだが、難しい式ばかり出てきて私には手におえそうもなかった。それで、斜面形状の凹凸を表現するんだったら難しい計算をしなくたって、単純に両隣セルとの標高差を表現したっていいじゃないの？ と思って私が勝手に作成したものだ。このときは、ラプラシアンと曲率がほぼ同義だということを知らなかった。 まだそんな何箇所も作成していないが、どうやら緩斜面の微地形を判読するのに使えそうだ。右図は傾斜量図にデコボコ図を40％透過で重ねてみた。尾根のように上に凸となる箇所を茶色、谷のように下にへこむ箇所を黒で示している。右図のクリックで、同じ地区を比較のために並べた図を示す。 同じデータを閲覧コンバートソフトで色指定することにより、尾根と遷急ゾーンだけを図示したり、水系と遷緩ゾーンだけを図示したりもできる。 ラプラシアンとの違いについては別途考察。 
曲率とラプラシアンの特徴について、下記曲率の項に示した文献に書いてあった。『凹凸の指標として, ラプラシアンが使われることが多い。ところがラプラシアンには, 値が地表面の傾斜に依存するという欠点がある。』『ラプラシアンが特にある物理モデルを表現するような場合以外は, 地表面の曲面的な凹凸の指標としてラプラシアンを使うのは不適切である。』 ・・・ 西田ほか（1997）より抜粋 確かに、垂直に近い面のラプラシアンを計算すると、水平距離が0に近づき、水平距離で微分するラプラシアンは無限大に近づくんだろうが、実際の斜面を考えると、勾配45度以上の斜面の凹凸が問題になることは少ないんじゃないだろうかと思ったのだが、だめらしい。 【 斜面曲率 】
（1） 西田ほか（1997）の式曲率の算出方法を探したところ、西田ほか（1997）を見つけた。計算式は右の偏微分方程式であるが、数値計算する際には右下図のaとbの両方のケースを計算し、それを平均して求めている。単位は、通常、標高もグリッド間隔（セルサイズ）も m なので、 m^-1 となる。 K＝{Kxx (1+Ky^2) + Kyy (1+Kx^2) - 2 Kx Ky Kxy} ／{2 (1+Kx^2+Ky^2)^(3/2)} 
右図aのケースKx＝(H23-H21)／(2*Dx) Ky＝(H32-H12)／(2*Dy) Kxx＝(H21+H23 - 2*H22)／Dx^2 Kyy＝(H12+H32 - 2*H22)／Dy^2 Kxy＝(H31-H13-H33+H11)／(2*(Dx^2+Dy^2)) 右図bのケース Kx＝(H33-H11)／(2*(Dx^2+Dy^2)^0.5) Ky＝(H13-H31)／(2*(Dx^2+Dy^2)^0.5) Kxx＝(H33+H11 - 2*H22)／(Dx^2+Dy^2) Kyy＝(H13+H31 - 2*H22)／(Dx^2+Dy^2) Kxy＝(H23-H21-H32+H12)／(Dx^2+Dy^2) * 西田顕郎、小橋澄治、水山高久（1997）：「数値地形モデルに基づく地震時山腹崩壊斜面の地形解析」、砂防学会誌、Vol.49、No.6、pp.9〜16 （2） ArcGISの式 ArcGISでは曲率の計算式としてこの数式が使われている。100倍しているのは、 曲率 = - 2 * (D + E) * 100 D = [(H21 + H23) /2 - H22] / Dx^2 E = [(H12 + H32) /2 - H22] / Dy^2 * ArcGIS 10 Desktop ヘルプの中で、本ページに特に関係するのは サーフェス ツールセット 【 地形変換線 】
国土地理院基盤地図情報（10mDEM）および同閲覧コンバートソフトを使用
【尾根・谷、水系密度】
この項は岩橋（1994）による。尾根・谷の検出は、メディアンフィルタを使う。DEM値と、それを中心に含む3×3セルの中間値（メディアン）を比較し、DEM値の方が大きければ尾根、小さければ谷とする。尾根を100、谷を-100、それ以外を0とし、それぞれに茶、濃青、白を割り当てて作図したのが右図（東西約3.5km）。曲率を使うよりよっぽど簡単に図化できてしまった。孤立した点（セル）をノイズとして除去フィルタをかければよりよさそう。 DEMを使って水系図を自動で抽出させるのに、0次谷をラプラシアンで排除している文献も見た。水系を抽出する計算が載っていないのでどうしているのかわからないが、メディアンフィルタと曲率を組み合わせると精度の高い水系図や地形変換線図ができるかもしれない。 ところでメディアンフィルタがなぜ使えるのか？ 同じ曲率を示す尾根上の点と遷急点があったとする。尾根上の点は計算される9点のうち上位3点に入る確率が高く、遷急点は中位3点に入る確率が高い。だから尾根上の点は抽出されるけど遷急点は抽出されにくいということなんだろう。 谷を100、それ以外を0として水系図としたのが下左図（東西約6km）。半径100mに入る谷としたセルの数を円フィルターを使ってカウントして図化したのが下右の水系密度分布図。当地の場合、10mDEMのセルサイズが南北と東西とで大きく異なるため、東西は最大23セル、南北は最大17セル、100mに入るセル数は287となっている。最大値を茶、最小値0を薄青、中間値に白を割り当て、水系図を重ねている。 なお密度分布を求めるのに、岩橋（1994）で半径10セルとしていたのでそのままやってみたが、DEMのメッシュ間隔がそもそも違うので下右図のように違和感のある図になってしまった。しかしこういう手法もあることを知った。半径250mで計算し直すとこうなった（東西約7km）。 ところで、空中写真から作成した10mDEMがたまたま手元にあるのでこれでもやってみた。山地はそこそこうまく表現できてるけど緩斜面はノイズが多そうだな、と思ってよく見たら、市街地は人工構造物、特に道路法面や河川堤防がきれいに出ていることに気づいた。さすがは航測10mDEM。どうりで水系密度がやたら高いはずだ ^^； 
【 texture図 】texture図 と呼ぶものなのかは知らない。地形の細かい凹凸（たぶん数m〜10mオーダー）のことを、”表層のきめの粗さ、細かさ” と表現するらしく、この微地形模様のことを texture と称するらしい。【 陰影図 】陰影図は段彩図と並んで最もよく見る図だと思う。これもエクセルで作れるが、実際的には少々面倒。陰影図を作るには、光源の高度と方位角が指定されていなければならない。ということは、光源位置を変えるたびに全部計算し直さなければならない。閲覧コンバートソフトがあるんだったら陰影図はそれでできるので、エクセルでいちいち計算することもないかな、と思う。 やり方としては、傾斜角度と傾斜方向を算出し、光源方向との角度を計算する。セルが光源方向を向いていれば明るく、反対方向ならば暗くする。 なお、陰影図は光源位置によって地形の見え方が大きく変わる。このため、複数の光源位置で作成した陰影図を重ねた 多光源陰影図 というものもある。多光源陰影図を作成するんだったら、マクロを作っておけば、計算時間はかかるだろうがエクセルの方が便利かもしれない。まだ試したことはないが、いずれ気が向いたらやってみたい。
【 地上開度図、地下開度図 】
一方、地下開度図は着目点における鉛直下方にある天底と地表面との角度の最低のものを8方向測定し平均したもので。すなわち、谷底であれば角度が大きく、山上であれば小さくなる。同様に地下開度を明度と比例させて作成した地下開度図では、谷筋が白く強調され、尾根筋がぼやけた表現となる。 * 横山隆三、白沢道生、菊池祐（1999）：「開度による地形特性の表示」、写真測量とリモートセンシング、vol.38、No.4 なお、オープンソースのSAGA GISというフリーウェアには、Morphometric Protectection Indexという解析項目があって、地上開度図が算出できる。計算時間も、Excelにくらべるとあっという間だ。 地下開度図の項目は見当たらないが、地形を反転させて （単純には、例えば標高値をマイナスにする） 地上開度図の計算をすると、地下開度図が求まる。 【 尾根谷度図 】地上開度は天頂と地形断面の上側接線の角度、地下開度は天底と地形断面の下側接線の角度。地上開度と地下開度との関係から傾斜の値と独立した尾根か谷かを示す指標として、着目点と地形断面上の点を結ぶ直線の角度の最大値と最低値の平均値を尾根谷度として定義された（千葉、鈴木、2004）。尾根谷度＝（地上開度−地下開度）／2 尾根谷度は、尾根部ではプラス、谷部ではマイナス、平野や広い平面状斜面では0の値をとる。ただし、単純に地上開度と地下開度の値を平均するのではなく、両者に重みをつけて合成することにより地形の半透明感や尾根と谷の強調の度合いを微調整できるそうだ（千葉ほか、2007）。 * 千葉達朗、鈴木雄介（2004）：「赤色立体地図−新しい地形表現手法」、応用測量論文集、vol.15、pp.81-89 * 千葉達朗、鈴木雄介、平松孝晋（2007）：「地形表現手法の諸問題と赤色立体地図」、地図、vol.45、No.1 【 Convergence Index 】QGISのプロセッシングにはいろんなアルゴリズムがついていて、Geoalgolithmsという項目もあった。その中にLaplacian filterもあってmethodにkernel1〜3と3種類選択できる。でも、その違いが判らない。Convergence Indexなるものもあった。わからないまま走らせてみると、尾根谷度図のようなものが出てきた。SAGA GISのhelp？を見ると、convergentとdivergentの指標だというから曲率のことらしいが、どうやって計算しているのかな？ 下左図はこれで出てきた図。右はSAGA GISで算出した地上開度図と反転地下開度図とを重ねたもの。やっぱり違うか。   【 接峰面図、接谷面図、起伏量図 】升本先生、木庭先生のページ参照接峰面図・接谷面図には方眼法と谷埋め法とがあるが、エクセルで簡単に作成できるのは方眼法。 接峰面図＝ ひとつの方眼範囲の最大値 接谷面図＝ ひとつの方眼範囲の最小値 起伏量図＝ ひとつの方眼範囲の最大値− ひとつの方眼範囲の最小値 谷埋め法に準ずると思われる接峰面図の作成法が 中山、隈元（2000）の３(1)に記述してあった。これは繰り返す手間があるが、Excelで簡単にできた。右に同じ個所の尾根谷度図と並べて示すが、尾根の形状がそのまま残り、たまたま作図したこの範囲だと地質による差がよく出ている。 * 中山大地、隈元崇（2000）：「細密DEMによる研究展望」、『デジタル観測手法を統合した里山のGIS解析 ―東京大学空間情報科学研究センター公開シンポジウム―』、杉盛啓明ほか編、地域環境GIS研究会、31-34 【 起伏指標 】QGISソフトのラスタ解析項目に「地形解析」というのがある。地形解析の項中には、この5項目の計算が用意されている。傾斜 傾斜方位 陰影解析 レリーフ 起伏指標 傾斜は斜面勾配を角度（度）で算出するもので、ArcGISと同じ計算式を使っているらしい。 傾斜方位の計算式は不明だが、これもArcGISと同様だろう。ただ、平坦面の表示をどうするのか設定するところがない。 陰影解析とレリーフというのが並んでいるのがわからない。陰影解析を実行すると、陰影図（レリーフマップ）が作成される。陰影図はHillshade mapかな？ ただ、陰影解析だと光源方向を任意に指定できるが、レリーフでは指定できない。 そして初めて知ったのが起伏指標。これを実行しようとすると、開くウィンドウには “粗度” とタイトルがでる。計算式を探してみると、どうやら TRI（Terrain Ruggedness Index）という指標で、目的のセルと周囲8近傍のそれぞれとの差の絶対値の平均 らしい。つまり、自分と周囲との平均高度差ということになる。 言い換えると、自分の隣がどのくらい高い／低いか、ということ。ところで、傾斜は自分の隣との最大傾斜方向での角度。そう考えると、起伏指標ってのは自分の隣との平均高度差なので、セル幅で除すると平均角度ということになる。 下図は、とある火山地形。左が傾斜、右が起伏指標をQGISで計算し、色合いが似るように割り当てたもの。案の定、ほとんど見分けのつかないものが作成された。  ところで、このTRIという指標は、下記が最初に提案したものらしい。 * Riley et al.(1999) ; A terrain ruggedness index that quantifies topographic heterogeneity. Intermountain Journal of Sciences 5(1-4):23-27. 検索するとこの論文のPDFが出てきた。これによると、オリジナルの計算式は差の二乗和の平方根をTRIと定義しているようだ。二乗平均平方根ともやや異なるが無視できるということか。QGISでは差の絶対値の平均になっている。計算速度を上げるためかな？ まぁ、二乗和でも平均でも、作図するときの色調補正で同じになるからどちらでもいいんだけど。 【 ラプラシアン、エッジ 】この項を最初に作成した頃はラプラシアンが何を意味しているのか全く知らなかった。その後、ものすごく重要な指標であることを知って、それまで知らなかったことが恥ずかしくなった。 （2013.9.10 追記）デジタル画像で使われるラプラシアンフィルタは、目標セルと周囲のセルとの大きな変化を抽出することを目的とするらしい。地形データにおいては、地形の凹凸および急激な傾斜変化点を表す指標となる。尾根谷度と逆にプラス側で凹地形（目標点より周りが高い）、マイナス側で凸地形（目標点より周りが低い）を表す。 デジタルラプラシアン＝H11＋H12＋H13＋H21＋H23＋H31＋H32＋H33−8×H22
エッジはラプラシアンの絶対値をとったもので、値が大きいほど大きな傾斜変換点であることを示す。画像処理においては閾値を設けて2値化して輪郭検出に使ったりするようだが、地形では値の大小の程度も重要だと思う。 【 高度分散量図、高度分散異常図 】高度分散量は単位面積中に含まれる標高値データの標準偏差と定義され（Ohmori、1978）、起伏量と類似する。同一の侵食基準面をもつ地域では、高度分散量はエリアの平均標高の関数として近似できる（Ohmori、1978）。近似曲線を地域別および地質単元別に比較すると地形は地質や構造運動を反映していることが明らかとなり、逆に近似曲線からはずれる個所は局所的に別な要因が作用していることが考えられる（例えば断層運動など）ことから、高度分散量の近似曲線からのずれ（乖離度）を高度分散量異常として定義した（三箇、藤原、2000）。 * Ohmori, H.（1978）：「Relief structure of the Japanese mountains and their stages in geomorphic development」、Bull. dept. Geogr. Univ. Tokyo、No.10、p.31-85 * 三箇智二、藤原治（2000）：「高度分散量を用いた数値地図からの地形・地質的特長の抽出」、情報地質、vol.11、No.2、pp.118-121 |
