雪崩危険箇所を衛星から取得した標高データと重ねて分析してみた
国土地理院が提供している標高タイルと、国土数値情報が提供する土砂災害・雪崩メッシュデータ、および雪崩危険箇所を組み合わせて、雪崩に対して注意が必要な場所について、地理情報を可視化して考察してみました。
2022年8月31日以降、Tellus OSでのデータの閲覧方法など使い方が一部変更になっております。新しいTellus OSの基本操作は以下のリンクをご参照ください。
https://www.tellusxdp.com/ja/howtouse/tellus_os/start_tellus_os.html
国土地理院が提供している標高タイルと、国土数値情報が提供する土砂災害・雪崩メッシュデータ、および雪崩危険箇所を組み合わせて、雪崩に対して注意が必要な場所について、地理情報を可視化して考察してみました。
(1)はじめに
日本は世界でもまれに見る雪国で、積雪が1mを超える場所にこれだけの人口が住んでいる国はそう多くありません。
そして、雪は時に様々な災害を引き起こします。特に雪崩は積雪山間部では毎年のように発生しています。しかしながら、雪崩は実際に目にする機会が少なく、どのような状況で起こりやすいのか、今一つイメージできない人が多いのではないでしょうか。
そこで、本記事では衛星データと雪崩の発生箇所、雪崩危険箇所のオープンデータを組み合わせて、雪崩が起こりやすい地形の条件について地図上に可視化して考察してみたいと思います。
(2)雪崩発生現場データの説明(取得方法、データポリシー)
使用するデータは下記の通りです。
- ●土砂災害・雪崩メッシュデータ（国土数値情報）
- ●標高タイル（国土地理院）
- ●雪崩危険箇所（ハザードマップポータルサイト）
まず、国土数値情報から土砂災害・雪崩メッシュデータを全国分取得します。こちらは平成18年度〜21年度に起こった土砂災害をメッシュデータとしてまとめたものになります。ダウンロードすると下記のようにshapeファイル形式のデータが格納されています。
こちらのshapeファイルをgeojson形式に変換したのち、全国分のデータを一つにマージしました。さらに、その中から災害タイプが「雪崩」となっているデータだけを抽出したgeojsonファイルを作成しました。それぞれ作成したデータは下記に置いてあります。
国土地理院の標高タイルは航空レーザー測量により、計測された標高のデータになります。APIが公開されているのでそちらを利用することで標高のタイル画像をpng形式で取得することができます。
ハザードマップポータルサイトでは雪崩危険箇所のタイルデータをAPI形式で公開しています。こちらは国土数値情報にある平成22年度土砂災害危険箇所をもとに国土地理院が加工したデータになります。
上記のハザードマップポータルサイトのAPIを使うことでこの雪崩危険箇所のタイル画像をpng形式で取得できますので、先ほどの標高タイルと組み合わせて、雪崩による被害が起こりやすい箇所と地形の関係をみていく際に利用します。
(3)標高データと雪崩危険箇所のデータを重ね合わせて確認してみた
それでは先ほどの土砂災害・雪崩メッシュデータをTellusOS上にアップロードして中身を確認してみます。
災害タイプを雪崩だけに絞るとあまりデータ量が多くないことがわかります。もう少し拡大してみましょう。
雪崩メッシュデータは主に北陸方面にデータが偏っていることがわかります。
次にTellusOS上でタイル座標を表示させ、雪崩メッシュデータが存在しているタイルを確認してみます。
上記で確認できたタイル座標を使って、国土地理院の標高データと、ハザードマップポータルサイトから取得した雪崩危険箇所を重ねて表示してみました。
上記の図の中のz=10, x=906, y=398の地点に関して、国土地理院の標高タイルから取得した画像がこちらになります。
一部、完全に黒く塗りつぶされてしまっているエリアが見受けられますが、こちらはそもそも測量されていない地点と思われます。後ほどズームレベルを調整して、雲の影響がないエリアに絞って確認していくこととします。
次に、ハザードマップポータルサイトから、z=10, x=906, y=398地点の雪崩危険箇所のPNG画像を取得してみます。
これらの画像を重ね合わせたものが次の画像になります。
無事に重ね合わせてみることができました。なんとなく標高と雪崩危険箇所に規則性がありそうに見えますが、このズームレベルではまだよくわかりませんね。もう少し拡大して確認してみることにします。
上記のz=11, x=1813, y=796地点の標高データおよび雪崩危険箇所のデータを重ね合わせてみてみることにします。
さらにもう１段階近づいた画像も取得してみます。
ズームレベルが１２になると、標高の高低と雪崩危険箇所の分布の傾向がはっきりとみて取れますね。
(4)検証結果
次に、雪崩の危険箇所と標高にどんな関係があるのか考察してみます。
下記はズームレベル１２における複数箇所の標高と雪崩危険箇所を重ね合わせた画像になります。
これらの画像を見ると、標高画像の脈のように黒くなっている標高の低い地点、つまり谷に沿って雪崩危険箇所が集中していることがわかります。特に谷と谷とが合流する地点「出合い」に雪崩が起こりやすい場所が点在しているように見えます。
出合いとは、
複数の地形が合流する地点の呼称。
例えば、二つの沢(谷、川)が合流する地点 など。
南アルプスの”野呂川出合”など、そのまま地名になっている場合もある。
逆に、谷に対して雪崩が起きやすい地点は、一般的に冬に西風または西北風が吹くことから雪の吹き溜まりのできやすい風下斜面に当たる東または東南に向いている方角と言われています。
一方で、雪崩危険箇所の分布を見ると、谷に対して東西均等に雪崩の起きやすいとみなされている地点が多く存在しているように見えます。
これは、雪崩危険箇所を各自治体が設定する際、雪崩の起きやすさには傾斜や樹木の多少などの複合的な要因が絡んでおり、一概に方位だけでその発生度合いが決まるわけではないことを考慮しているのだと思われます。
また、高さのある木々が少ない、植生がまばらな場所では雪崩が起きやすいと言われています。
中高木が密に生えている斜面では雪崩が発生しにくい一方、低木林やまばらな植生の斜面では雪崩発生の危険が高くなります。笹や草に覆われた斜面などは裸地よりも発生しやすい地形です。
最大で時速200㎞ものスピードに！ 雪崩（なだれ）から身を守るために
そこでLandsat-8のナチュラルカラー合成画像と雪崩危険箇所を重ねて、実際に植生がまばらな土地で雪崩の危険性が高いのかを確認してみました。その結果が下記の画像になります。
植生が分布している箇所が緑色で強調されています。これを見ると確かに雪崩が起きやすい危険箇所は植生がまばらな箇所との境界に多く分布しているように読み取れます。
以上のことから、谷と谷の出合、特に植生がまばらになる箇所では雪崩が発生する危険性が高いと言えそうです。衛星画像から植生や地形を確認することで、積雪期の雪崩危険箇所を予測することができそうです。
また、西高東低の気圧配置から崩れやすい雪庇（せっぴ）などが南東方向に生じやすいと言われていますが、雪崩危険箇所を見る限りでは必ずしも一方向だけに注意していれば良いというわけでもなさそうです。
(5)まとめ
衛星画像と雪崩危険箇所の分布を重ね合わせることで、地形と雪崩の危険性との関連を画像上で確認することができました。衛星データを利用することで自然現象を可視化し、視覚的に把握しやすくすることができます。これから冬にかけて雪山登山を楽しまれる方も多いと思います。山行に出かける前にぜひ、ルート上に雪崩が発生しやすい箇所がないか、衛星画像で確認してみてはいかがでしょうか。今回使用したコードを下記にまとめましたので、適宜参考にしていただけたらと思います。今回の記事が雪崩による被害の防災に役立てていただけたら幸いです。
（ソースコードまとめ）
標高データと雪崩危険箇所を重ねるコード
import os, requests, json
from skimage import io
from io import StringIO
from io import BytesIO
import string
from time import sleep
import pandas as pd
import numpy as np
from PIL import Image
import cv2
from colour import Color
from matplotlib.colors import LinearSegmentedColormap
import matplotlib.pyplot as plt
from matplotlib import cm
import math
img_dir_path = "imgs/"
if not os.path.exists(img_dir_path):
os.makedirs(img_dir_path)
def get_hyoko_tile(z, x, y, params={}):
"""
国土地理院標高タイル情報を取得する
Parameters
----------
z : int
タイルのズーム率
x : int
タイルのX座標
y : int
タイルのY座標
option : dict
APIパスのオプション(entityId)
params : dict
クエリパラメータ
Returns
-------
img: ndarray
タイル画像
http://cyberjapandata.gsi.go.jp/xyz/dem5a/15/29100/12899.txt
"""
url = 'https://cyberjapandata.gsi.go.jp/xyz/dem5a/{}/{}/{}.txt'.format(z, x, y)
r = requests.get(url, params=params)
if r.status_code is not 200:
raise ValueError('status error({}).'.format(r.status_code))
#標高値がない区画は0mに置換する
maptxt = str.replace(r.text, u'e', u'-0.0')
Z = pd.read_csv(StringIO(maptxt), header=None)
Z_array = np.array(Z)
return Z_array
# reference
# http://memomemokun.hateblo.jp/entry/2018/11/01/221545
def get_avalanche_tile(z, x, y, params={}):
"""
雪崩危険箇所タイル情報を取得する
Parameters
----------
z : int
タイルのズーム率
x : int
タイルのX座標
y : int
タイルのY座標
option : dict
APIパスのオプション(entityId)
params : dict
クエリパラメータ
Returns
-------
img: ndarray
タイル画像
https://disaportaldata.gsi.go.jp/raster/05_nadarekikenkasyo/{}/{}/{}.png
"""
url = 'https://disaportaldata.gsi.go.jp/raster/05_nadarekikenkasyo/{}/{}/{}.png'.format(z, x, y)
r = requests.get(url, params=params)
if r.status_code is not 200:
raise ValueError('status error({}).'.format(r.status_code))
return io.imread(BytesIO(r.content))
#指定するタイル座標は適宜必要なものに置き換えてください。
(z, x, y) = (12, 3616, 1594)
hyoko = get_hyoko_tile(z,x,y)
fig1 = plt.figure()
plt.imshow(hyoko,cmap="gray")
plt.colorbar()
fig1.savefig('imgs/gdem.png')
#雪崩危険箇所データを読み込む
avalanche = get_avalanche_tile(z,x,y)
fig2 = plt.figure()
plt.imshow(avalanche[:,:,0], cmap="Reds")
plt.colorbar()
fig2.savefig('imgs/avalanche.png')
img1 = cv2.imread('imgs/gdem.png')
img2 = cv2.imread('imgs/avalanche.png')
# RGB形式に変換
img1 = cv2.cvtColor(img1, cv2.COLOR_BGR2RGB)
img2 = cv2.cvtColor(img2, cv2.COLOR_BGR2RGB)
fig3 = plt.figure()
blended = cv2.addWeighted(src1=img2,alpha=0.7,src2=img1,beta=0.3,gamma=0)
plt.imshow(blended)
fig3.savefig("imgs/belended_{}_{}_{}.png".format(z, x, y))
Landsat-8の衛星画像と雪崩危険箇所をプロットするコード
TOKEN = 'ここにご自身のトークンを貼り付けてください'
def get_landsat_tile_path_from_scene_id(id):
"""
Landsat-8のシーンIDから、バンド別合成画像APIで使用するパスを取得する
Parameters
----------
id : str
シーンID
Returns
-------
path : str
APIパス文字列
"""
url = "https://gisapi.tellusxdp.com/api/v1/landsat8/get_scene/{}".format(id)
headers = {
"content-type": "application/json",
"Authorization": "Bearer " + TOKEN
}
r = requests.get(url, headers = headers)
data = json.loads(r.content.decode('utf-8'))
return data["tile_path"]
def get_landsat_image(scene_id, zoom, xtile, ytile, opacity=1, r=3, g=2, b=1, preset=None):
"""
Landsat-8のバンド別合成画像を取得する
Parameters
----------
scene_id : str
シーンID
zoom : str
ズーム率
xtile : str
座標X
ytile : number
座標Y
opacity : number
透過率
r : number
Redに割り当てるバンド
g : number
Greenに割り当てるバンド
b : number
Blueに割り当てるバンド
preset : str
使用するプリセット
Returns
-------
img_array : ndarray
バンド別合成画像タイル(PNG)
"""
host = "https://gisapi.tellusxdp.com"
path = "/blend/" + get_landsat_tile_path_from_scene_id(scene_id).format(z=zoom, x=xtile, y=ytile)
queries = "?opacity={}&r={}&g={}&b={}".format(opacity, r, g, b)
if preset is not None:
queries += "&preset=" + preset
url = host + path + queries
headers = {
"content-type": "application/json",
"Authorization": "Bearer " + TOKEN
}
r = requests.get(url, headers=headers)
return io.imread(BytesIO(r.content))
scene_id = "LC81080342013155LGN01"
# use composition preset nvdi
landsat_img = get_landsat_image(scene_id, z, x, y, preset="natural")
avalanche_img = get_avalanche_tile(z,x,y)
fig4 = plt.figure()
landsat_blended = cv2.addWeighted(src1=avalanche_img,alpha=0.5,src2=landsat_img,beta=0.7,gamma=0)
plt.imshow(landsat_blended)
fig4.savefig("imgs/landsat_belended_{}_{}_{}.png".format(z, x, y))
