最大傾斜量マップ（雪崩関連）

 

 

将来、また引っ張り出して計算するかもしれないので、そのための備忘。

 

 

最終的にはこうなった

結論から先に。斜度を国土地理院がやっている方法よりさらに良い方法を思いついたので改めて斜度計算をして地図に落とした。その結果、危険斜度があらたに浮かび上がり、より安全なトレースラインを引けた。

 

国土地理院が提供している斜度マップ。赤いのは雪崩危険地帯。

 

 

この記事の方法で得られた斜度マップ。まぁぱっとみても等高線の按配からいってこのほうが正確に見える。青い部分（例）が新たに現れた雪崩危険斜度。探せば他にもっとある。

 

 

 

 

この斜度を考慮して新たに引いたトレースライン。

 

 

 

 

 

 

何が問題で、どうやったのか

 

 

あらまし

国土地理院は既に等高線マップ（標高データ）から、斜度計算を既にしている。さらに、地理院地図に重ねるレイヤーとして提供している。しかし、国土地理院の方法は、二点ばかり欠点がある。国土地理院が算出する斜度は、提供している二次元標高分布データにおいて、ある格子点を中心に、上下６点、左右６点を使って斜度を求めている。直感的な説明をすればこうだ。南北方向、東西方向の２つの斜度を、それぞれ６点の標高データを平均することで算出。それをさらに平均した斜度を求めている。やり方としては、二乗平均の平方根を取ることで、ある中心点における代表的な斜度を定めていることとなる。各格子が、それぞれ一定（平滑）の斜面を持つとみなせる程度に格子間隔が十分に細かければそれでよいのだが、１０ｍ×１０ｍの格子ではそうでない場所が多く見受けられる。たとえば、北東、北西、南東、南西に斜面が走っている時に、その斜面が本来持つ斜度が見えずらくなるので、危険斜度部分で鈍されて緩い斜度とみなされている部分が出てくる。また平滑した斜面を持つ１つの格子とみなして計算しているから、格子周りで斜度の向きや大きさが異なるような時、最大斜度はどっちの方向でその斜度は？ という探索が重要なのだが、その点も無視されている。

 

ということで、この記事で求める斜度計算は４方向にそれぞれ正確な斜度を幾何学的に求めることにする。かつ、そのうちの最大斜度を採用することで斜度分布を得る。適切には、北方向、南、西、東、北東、北西、南東、南西方向の合計８方向をきちんと計算して最大斜度を求める必要がある。今回のサンプルでは東西南北に対して斜めに走った斜面で検証しているので結果は変わらない。

 

 

 

 

雪崩と傾斜度

雪崩発生のメカニズムと傾斜度の関係は比較的よく調査されている。ただし、植生の状況や、山の斜面が東西のどっちを向いているか（積雪量が全然違う）、雪崩を誘引するような特異点（突起、雪庇ができる・崩れる、落雪がおきやすい）があるか（持論）などで、事情は異なる。しかし、そういっていては定量的に推し量れない。一般的には斜度と雪崩の発生可能性の高さについてはこうされている。

 

 

 

 

 

国土地理院の標高データ

日本では国土地理院が、基盤地図情報数理標高モデルという標高タイルを公表しており、標高データを自由にダウンロード可能である。タイル幅は、測量方法により５ｍ×５ｍ（平野部、特定の火山地帯）、または、１０ｍ×１０ｍ（その他の地域）となっている。５ｍ幅といえば、歩幅で１０歩弱程度であるので、相当に精度は高い。前者は写真測量と航空レーザー測量によるもので全国は網羅していないものの、最大粗さでも１０ｍ×１０ｍ幅で国土を網羅しているということである。

 

 

国土地理院による斜度量区分マップ

国土地理院でも、上記の標高タイル（標高）を用いて全国の斜度量を算出して地図として公開している。説明資料はコチラ。

 

地図はコチラ。

 

 

国土地理院が採用している斜度量の欠点？

標高タイルにおける格子点８点をもちいている。直感的な説明としては、h₂₂を中心として、東西方向のそれぞれ二方向について、南北、東西側のそれぞれ３点、合計６点を用いて、Sx、Syを算出している。これをもとにSを、二乗平均平方根S、として斜度量としている。このSの逆正接を求めることに寄って、斜度D（角度）を求めている。

 

 

これにより上の地図を得ているのだが、２点ほど問題があるようにみえる。欠点１：h₂₂を中心に合計８点を用いて、東西、南北というｘｙ平面の２方向に分解してそれぞれの斜度を出しているため（ことになるので）、北西、北東、南西、南東方向に斜面が走っている場合にずれが大きいだろう（鈍る）。欠点２：そのメッシュにおける最大斜度こそが重要であるにも関わらず、合計８点に捜索範囲を広げ、二乗平均の平方根を取ってしまっているため、より簡便に滑らかに表現は出来るものの、結果的には、特異的な場所について斜度を丸めて捨てることとなる。

 

 

 

格子点における最大斜度を活かす方法

ぱっと思いつく単純な方法は、h₂₂を中心として、東西南北のそれぞれ４点を用いて、三角関数（逆正接）でそれぞれの方向の斜度が分かるので、４方向のうち最大の斜度のものを使用すればよい。これにより欠点２：は克服できるが、欠点１：は残る。斜度を東西南北の２方向で定義してしまっているので、下のような地形に格子点があるときに斜度を正確に表現できていない。実際、これにより計算した場合には、本来の斜度よりゆるい斜度が算出されてしまう。たとえば、その地点における最大斜度は、北西方向に３０°であったとしても、この計算では西方向に２６°、といった具合となる。国土地理院が採用している、格子８点を用いた二乗平均平方根の手法も同様の傾向がでる。

 

 

上の橙の四角のブロックは画像の関係で少し滲んでいるが、１つ１つが、まさに１０ｍ幅の格子である。等高線の間隔に比べて十分に小さいとは言えない。よって上に示したような格子に置いて斜度を求める時に、国土地理院が採用している代表斜度を求める手法では、正確に斜度を得られないことは視覚的に分かろう。

 

 

 

幾何学的に角度を求める

空間上の３点の格子点を通る平面と地平との成す角度を求めることで欠点１はクリアできる。たとえば、h₂₂を中心として h₁₂, h₂₁を通るXYZ平面と地平面（水平面）との成す角度（斜度）を求めることとなる。h₂₂を中心にこのようにして、３点をとっていくと、h₂₂を中心に格子点を合計９点つかうならば、８通り（変則的な取り方をすればもう少し可能。最大８４通り）の平面を規定できるが、ここでは簡易的に国土地理院との違いを明確に出すためにも、上記例の取り方で４通りで行う、つまり斜め。正しくは最低でも８通りで行うべし。すなわち、xの正負の方向, yの正負の方向で４つ、x軸, y軸と４５度の方向に４方向。すなわち８方向となる。

 

 

 

ところで、こんな平面を使わずともh₂₂の左上h₁₁の格子点を使えば、二点間の格子点により逆正接で簡単に左上方向の斜度は出せる。実際やってみたが、これで斜度を出すと極端なケースが多すぎて斜度の分散が大きくなりすぎるようだ。よって上の方法を採用するに至った。ということである程度はこの方法でも鈍らせている。けれど上図を見ても視覚的にはその方が正しかろうことは分かる。二乗平均平方根と二点間斜度の間に位置する方法と考えてもらえばいいのではないかと。

 

 

平面と平面の成す角

(x , y, z) 空間において、平面a₁x + b₁y + c₁z +d₁ =0、a₂x + b₂y + c₂z +d₂ =0 があるとき、この２平面のの成す角は、それぞれの平面の法線ベクトルは (a₁, b₁, c₁)、(a₂, b₂, c₂)であるから、cosθ = | a₁a₂ + b₁b₂ + c₁c₂ | / [√(a₁² + b₁² + c₁²) √(a₂² + b₂² + c₂²) ]　で定義される。よって θ = cos^-1 [ | a₁a₂ + b₁b₂ + c₁c₂ | / [√(a₁² + b₁² + c₁²) √(a₂² + b₂² + c₂²) ]　 ]　（式１）

 

現在考えている一方の平面は、地平であるため、この片方の平面は z=0 となる。つまり、a₁=b₁=d₁=0, c₁=1（任意）

 

もう一方の平面は、h₂₂を中心に、あと２点を選ぶこととなるが、h₂₂は便宜的に 原点(0,0)とするとプログラム上のコーディングと計算が容易となる。他二点を選んだ時に、便宜的に中心点を0, 他を1, 2と名前をつけ、法線ベクトル (a₂/d₂, b₂/d₂, c₂/d₂)　なる平面２を、改めて、(a₂, b₂, c₂)により定義してから説明する。また、選んだ２点、それぞれのz値を, z₀ , z₁ , z₂ とすれば、原点を通ることから、c₂ z₀ + 1 = 0となる。整理すると,  z₀ a₂x+ z₀ b₂y + -z + z₀= 0となる。格子間隔10(m)と、 z₁ , z₂ を代入して、a₂, = (z₁ – z₀)/10z₀、b₂ = (z₂ – z₀)/10z₀ と求まる。ここで求まった法線ベクトルはすべて分母にz₀を有するので、全てにz₀を掛けて、あらためてこれを法線ベクトル(a₂, b₂, c₂)とする。z₀、z₁、z₂は、地理院地図におけるそれぞれの標高データであるから、これを（式１）に代入することで斜度θが求まる。この計算をh₂₂を中心に、n 通り (n=最小４, 適当なのは８, 最大で８４=₉C₃=₉P₃/3!)行い、Max ( θn )が格子点h₂₂における最大の斜度、最大の傾斜角となる。

 

 

 

 

傾斜度マップ（グラデーション、原図）

上記計算により求まった傾斜度マップ。

 

 

段階的斜度区分で示した雪崩危険斜度マップ

注）４５°～５５°も橙色としている。新たに出現した危険斜度でパッと目についた場所だけ、青く丸印をつけた。

 

 

 

 

国土地理院による同マップは下図。地理院地図が提供している斜度マップでは現れていない雪崩危険斜度が上図で現れている。大東岳山頂周りでも簡単に違いが分かる。

 

 

 

 

 

 

新たに引かれる予定トラック

 

数十ｍほど東へ寄ったほうがよさそうだ

 

 

---

VBA ソースコード

※色塗りがそのまま可能なのでVBAを使っているが計算は遅い。メモリ１６ＧＢ、Core i7 （二年前のレベル）で数分かかる、描画が同時におこると極めて遅くなるため、Run前はすべてのセルを無色にして、条件付き書式等は一旦すべて解除するべし。計算するだけならfortranがベター。ただし、格子点、標高をきちんと別ファイルに作る必要があるのでその手間はかかる。日本全土を計算するならば当然その方法で。ちゃんとした解析屋さんがみればまず倍精度じゃないからだめだ！とか叱られそうだが！

 

Sub eve()

ii = 3

For j = 0 To 749
For i = 0 To 1124
ii = ii + 1
Cells(j + 3, i + 3) = Cells(ii, 1)
Next
Next

End Sub

Sub gradient2()
‘3次元平面と水平面の成す角を求め、最大斜度を決定する。（4通り）
‘必要に応じて１２通り

For j = 4 To 1126 ‘1126 ‘1126 ‘列
For i = 5 To 751 ‘751 ‘行
g5 = 0
‘右上
‘ (0,0,1) ( cells(i-1,j)/10, cells(i,j+1)/10, -1)
cos1 = 1 / Sqr(((Cells(i – 1, j) – Cells(i, j)) / 10) ^ 2 + ((Cells(i, j + 1) – Cells(i, j)) / 10) ^ 2 + 1)
rad = ACN(cos1)
g1 = GetAngle(rad)

cos2 = 1 / Sqr(((Cells(i – 1, j) – Cells(i, j)) / 10) ^ 2 + ((Cells(i, j – 1) – Cells(i, j)) / 10) ^ 2 + 1)
rad = ACN(cos2)
g2 = GetAngle(rad)

cos3 = 1 / Sqr(((Cells(i + 1, j) – Cells(i, j)) / 10) ^ 2 + ((Cells(i, j + 1) – Cells(i, j)) / 10) ^ 2 + 1)
rad = ACN(cos3)
g3 = GetAngle(rad)

cos4 = 1 / Sqr(((Cells(i + 1, j) – Cells(i, j)) / 10) ^ 2 + ((Cells(i, j – 1) – Cells(i, j)) / 10) ^ 2 + 1)
rad = ACN(cos4)
g4 = GetAngle(rad)

h = Application.WorksheetFunction.Max(g1, g2, g3, g4)

If h < 30 Then
‘h = 0
Cells(755 + i, j).Interior.Color = RGB(255, 255, 255)
End If
If h >= 30 And h < 35 Then
‘ h = 1
Cells(755 + i, j).Interior.Color = RGB(255, 153, 0)
End If
If h >= 35 And h < 45 Then
‘ h = 2
Cells(755 + i, j).Interior.Color = RGB(255, 51, 0)
End If
If h >= 45 And h < 55 Then
‘h = 3
Cells(755 + i, j).Interior.Color = RGB(255, 153, 0)
End If
If h >= 55 And h < 60 Then
‘h = 3
Cells(755 + i, j).Interior.Color = RGB(255, 255, 0)
End If
If h >= 60 Then
‘h = 4
Cells(755 + i, j).Interior.Color = RGB(255, 255, 255)
End If
‘Cells(755 + i, j) = Round(h, 1)

‘Debug.Print cos1, cos2, cos3, cos4
‘If g6 < h Then
‘ g6 = g5
‘End If

Next
Next
‘Debug.Print g6
‘125% 162% on fireworks
End Sub
Function ACN(x) ‘逆余弦 arccos
Select Case x
Case Is = -1
ACN = pi
Case Is = 1
ACN = 0
Case Else
ACN = pi / 2 – Atn(x / Sqr(1 – x ^ 2))
End Select

End Function

Function pi()
pi = Atn(1) * 4
End Function

Function GetRadian(ByVal angle)
GetRadian = angle * (pi / 180)
End Function

Function GetAngle(ByVal radian)
GetAngle = radian / (pi / 180)
End Function