2022年8月31日以降、Tellus OSでのデータの閲覧方法など使い方が一部変更になっております。新しいTellus OSの基本操作は以下のリンクをご参照ください。
https://www.tellusxdp.com/ja/howtouse/tellus_os/start_tellus_os.html

突然ですが、みなさんは「雪渓」という言葉を聞いたことがありますか？

高山など標高の高い場所の谷や沢では、積雪が一年中溶けずに残っている地帯があり、それが「雪渓」です。特に国内では「白馬大雪渓」「針ノ木大雪渓」「剱沢大雪渓」が有名で、「日本三大雪渓」と言います。

今回はLandsat8の衛星データから、地表面における積雪域を示す指標であるNDSI（Normalized Difference Snow Index）正規化積雪指標を用いて、日本三大雪渓として有名な白馬岳の大雪渓を衛星画像から確認してみたいと思います。

1.白馬大雪渓とは

本記事では、「日本三大雪渓」の中でも、登山者に特に人気のある「白馬大雪渓」について、衛星画像を使って確認してみたいと思います。

参考）日本三大雪渓のひとつ白馬大雪渓へ お手軽トレッキング

白馬大雪渓は白馬岳に向かう途中に位置する、幅100m、長さ3.5kmに及ぶ広大な雪渓です。白馬岳は標高2,932mの高さの山で、6〜8月の夏季においてもその平均気温は13.4℃と避暑地としても人気のコースになっています。

白馬岳の名前の由来は、春の雪解けにより露出した岩肌が「代掻き馬」のように黒い雪形が見られたことから「代馬(しろうま)岳」と呼ばれるようになったそうです。

白馬岳と大雪渓の位置関係は地理院地図で確認できます。

2.検出方法と使用するデータの説明

今回は積雪の位置を検出するために、正規化積雪指数（NDSI）を利用します。NDSIは地表面における積雪域を示す指標です。雪の反射特性が緑の波長付近で高く、短波長赤外（SWIR）付近において低いことを利用して作られています。この指標を利用すると、雲と積雪を識別することも可能です。NDSIは下記の計算式から求めることができます。

NDSI = (Green – SWIR) / (Green + SWIR)

正規化指標には様々なものがあります。その他の指標については下記の記事も参考にしてください。

緑の波長と短波長赤外を取得するために、Landsat-8の衛星データを利用します。緑および短波長赤外の波長はバンド3とバンド6から取得できるので、それらを組み合わせてNDSIを計算します。

3.コードと結果のご紹介

それでは実際にLandsat8の衛星画像を取得していきましょう。

import os, requests, json
from skimage import io
from skimage import img_as_ubyte
from skimage import filters
from io import BytesIO
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

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))

シーンIDおよびいくつかの引数を渡すことでLandsat8の衛星画像を取得する関数を定義しています。

実際に白馬岳周辺のシーンをTellusOS上で探してみます。

Tellus OSの画像で雪か雲かわからない旨を伝える構成にすると流れとして良いかもしれません。

夏がすぎても雪が残っていることを確認するために、冠雪直前の10月ごろのシーンを選択してみました

上記のシーンのIDを先ほど定義した関数に渡して実際の画像を確認してみます。

scene_id = "LC81090352015280LGN01"
z = 12
x = 3615
y = 1597
# use composition preset nvdi
true_img = get_landsat_image(scene_id, z, x, y, preset="true")

io.imshow(true_img)

上記のコードを実行すると下記のような画像が得られたかと思います。

なんとなく白っぽく残っている部分が雪渓のように見えますが、画像の中には雲も混じっていることが見えるかと思います。

それではNSDIを取得して、積雪部分を強調した画像を作成してみましょう。

def calc_normalized_index(band1, band2):
    shape = band1.shape
    assert shape == band2.shape, 'Two images are different sizes.'
    
    b1 = band1.flatten().astype(np.float32)
    b2 = band2.flatten().astype(np.float32)
    numer = b1 - b2
    denom = b1 + b2
    normalized_index = np.where(denom != 0, numer / denom, 0)
    return normalized_index.reshape(shape)

まず先に、正規化した指標を計算する関数を定義しておきます。
こちらの関数を作成するにあたって、下記の記事を参考にしています。

それでは、上記関数を使ってNSDI指標を計算してみます。

# Landsat-8のバンド別画像合成APIを実行する
img = get_landsat_image(scene_id, z, x, y, r=3, g=6)


# NDSI正規化指標の計算
band_green = img[:,:,1]
band_swir = img[:,:,0]
ndsi_img = calc_normalized_index(band_green, band_swir)

io.imshow(ndsi_img)

band3(Green), band6(SWIR)を利用するため、r=3, g=6というように引数に指定しています。

その結果下記のような画像が取得できたかと思います。

赤く強調されている箇所が積雪が残っている部分になります。冒頭の地理院地図から、おおよその白馬岳と大雪渓の位置を図示すると下記のようになります。

大雪渓のあたりに、特に濃い赤色の積雪地帯があることが確認できます。

白馬岳一帯の積雪の分布を図に示したので、せっかくなので他の時期の積雪量と比較してみたいと思います。再びTellusOS上で該当するシーンを探します。

今度は5月のシーンを選択してみます。

先ほどと同様にまずはtrueカラーの画像を表示して雪の様子を確認してみましょう。

scene_id2 = "LC81090342016139LGN01"
true_img_05 = get_landsat_image(scene_id2, z, x, y, preset="true")
io.imshow(true_img_05)

まだ溶けきっていない積雪が広範囲にわたって残っていることがわかります。

ここから再びNDSIを計算して、積雪範囲を強調してみます。

# Landsat-8のバンド別画像合成APIを実行する
img_05 = get_landsat_image(scene_id2, z, x, y, r=3, g=6)

# NDSI正規化指標の計算
band_green_05 = img_05[:,:,1]
band_swir_05 = img_05[:,:,0]
ndsi_img_05 = calc_normalized_index(band_green_05, band_swir_05)

io.imshow(ndsi_img_05)

積雪の範囲がだいぶ広範囲ですね。

先ほどの10月の画像と並べて比較してみましょう。

io.imshow(np.hstack((ndsi_img_05, ndsi_img)))

約5ヶ月の間に多くの雪が溶けて無くなっていますね。
しかし、それでも残っている部分もあり、そこに雪渓があることがわかります。

さらに、1年を通じての積雪量の変化を確認するために1月から12月までのシーンをつなぎ合わせて、GIF画像を作ってみます。

scene_ids = [
    "LC81080352016004LGN02", # 1 # 雲ががっつり被ってしまっている、が他に1月のシーンが存在しない
    "LC81090352016059LGN01", # 2
    "LC81090342016091LGN01", # 3
    "LC81090342014117LGN01", # 4
    "LC81090342016139LGN01", # 5
    "LC81090342015152LGN01", # 6
    "LC81080352015193LGN02", # 7
    "LC81090342015216LGN01", # 8
    "LC81090342013226LGN01", # 9
    "LC81090352015280LGN01", # 10
    "LC81080352018329LGN00", # 11
    "LC81090342018336LGN00", # 12
]

images = []

for scene_id in scene_ids:
    # Landsat-8のバンド別画像合成APIを実行する
    img = get_landsat_image(scene_id, z, x, y, r=3, g=6)

    # NDSI正規化指標の計算
    band_green = img[:,:,1]
    band_swir = img[:,:,0]
    ndsi_img = calc_normalized_index(band_green, band_swir)

    images.append(ndsi_img)
    
import matplotlib.animation as animation
from IPython.display import HTML

fig = plt.figure()
ims = []
for img in images:
    ims.append([plt.imshow(img)])

ani = animation.ArtistAnimation(fig, ims, interval=1000, repeat_delay=1000)
ani.save('shirouma_scenes.gif', writer='pillow')
plt.close()
HTML('')

下記のようなGIF画像が生成されます。

画像生成に使用しているライブラリを変えているので生成された画像の色合いが変わっていますが、ここでは青色の部分が雪渓の分布を示しています。

1月、2月の画像については残念ながら雲に覆われている範囲が広く、積雪の分布が不明確になってしまっています。しかし3月から10月にかけて徐々に積雪の範囲がなくなっていき、11月に入ってからまた大雪渓が姿を表すまでの様子が見て取れます。

4.他の大雪渓についても確認してみる

せっかくなので、三大雪渓の他の雪渓に関しても比較してみてみることにします。

今回使用した２つのシーンのどちらにも劔沢雪渓も写っていたので、劔沢雪渓についてもNDSI値を計算して画像を作成してみます。

劔沢雪渓の位置に関しても地理院地図でご確認ください。

# 剱岳が収まるタイルを指定
z2 = 12
x2 = 3613
y2 = 1599

# 9月ごろの画像を取得
scene_id_09 = "LC81090342013226LGN01"
img_tsurugi_09 = get_landsat_image(scene_id_09, z2, x2, y2, r=3, g=6)

# NDSI正規化指標の計算
band_green_tsurugi_09 = img_tsurugi_09[:,:,1]
band_swir_tsurugi_09 = img_tsurugi_09[:,:,0]
ndsi_img_tsurugi_09 = calc_normalized_index(band_green_tsurugi_09, band_swir_tsurugi_09)

# 5月ごろの画像を取得
scene_id_05 = "LC81090342016139LGN01"
img_tsurugi_05 = get_landsat_image(scene_id_05, z2, x2, y2, r=3, g=6)

# NDSI正規化指標の計算
band_green_tsurugi_05 = img_tsurugi_05[:,:,1]
band_swir_tsurugi_05 = img_tsurugi_05[:,:,0]
ndsi_img_tsurugi_05 = calc_normalized_index(band_green_tsurugi_05, band_swir_tsurugi_05)

io.imshow(np.hstack((ndsi_img_tsurugi_05, ndsi_img_tsurugi)))

春先の5月と、雪が溶けきった9月ごろの画像を並べてみました。

9月ごろの画像右側に、夏を越えても残っている雪渓が存在することが確認できます。

剣岳に残る雪の様子は下記のYoutube動画も参考にしてみてください。

5.まとめ

Landsat-8の衛星データからNDSIの値を計算し、積雪範囲を可視化することで大雪渓の位置を確認することができました。

また、時期の異なる衛星画像同士を比較することで、雪解けによって積雪範囲の割合の変化を可視化することができました。

今回は日本三大雪渓を対象に積雪範囲を確認していきましたが、ある程度標高の高い山であれば大小の差はあれ雪渓を見かけることはよくあります。

登山に出かける前に、同時期の衛星画像で積雪範囲を確認しておくことで、登山ルート上のどこで雪渓がみられそうか把握しておくと、登山の楽しみが一層増すかもしれません。ぜひお試しください！