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

森林限界と呼ばれる、ある高度以上で樹木が全く見られなくなる地帯が山々には見られます。今回は標高データを取得できるASTER GDEM 2.0と、植生の分布を可視化できるLandsat-8のデータを使い、各標高ごとの植生分布を可視化してみたいと思います。可視化する対象として、昨年の夏に北アルプスに位置する穂高三山（北穂高、奥穂高、前穂高）に登ってきたので、北アルプスの山々の標高と植生の関係を実際の写真も交えながら可視化して見たいと思います。

1.森林限界とは

森林限界とは高緯度地域や高山などにおいて、高木が育たなくなる限界高度のことを差す用語です。日本アルプスなどの高山を登る登山家の間では登った高さを表すために度々耳にする言葉です。

樹木が育つためには適度な気温と湿度、水分、そして日光が必要となりますが、これらの条件のいずれかが欠ける場所では、樹木が全く育っていない景観が発生します。

参考情報：森林限界（ウィキペディア）

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

はじめにASTER GDEM 2.0の標高データを取得します。取得方法は、富士山の標高をプロットした下記の宙畑の記事を参考にしています。

また、植生の分布を見るためにLandsat-8のバンド（波長）別合成画像を取得します。植生の分布を見るために利用するプリセットとして、今回はNatural Colorを指定します。

Natural Colorによる合成画像の説明は下記の記事も参考にしてご覧ください。

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

それではまず標高データを取得してみたいと思います。

# 必要なライブラリを読み込む
import requests, json
from skimage import io
from io import BytesIO
import numpy as np
import math
import matplotlib.pyplot as plt
%matplotlib inline

TOKEN = 'ここにご自身のトークンを貼り付けてください'

def get_astergdem2_dsm(zoom, xtile, ytile):
    """
    ASTER GDEM2 の標高タイルを取得

    Parameters
    ----------
    zoom : str 
        ズーム率 
    xtile : str 
        座標X 
    ytile : number 
        座標Y
    Returns
    -------
    img_array : ndarray
        標高タイル(PNG)
    """
    url = " https://gisapi.tellusxdp.com/astergdem2/dsm/{}/{}/{}.png".format(zoom, xtile, ytile)
    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    
    r = requests.get(url, headers=headers)
    return io.imread(BytesIO(r.content))

それでは実際に標高データを表示してみます。今回は奥穂高岳に焦点を当てて、そのタイル座標を指定しています。

# 穂高山が存在しているタイル座標を指定
z = 10
x = 903
y = 401
gdem_img = get_astergdem2_dsm(z, x, y)
io.imshow(gdem_img)

このような画像が取得できたかと思います。

次にLandsat-8のシーンIDからバンド別画像合成APIで使用するパスを取得する関数と、それを使って合成画像を取得する関数を定義します。

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

次に、合成画像も表示してみましょう。
まずはTellusOS上で利用可能なLandsat-8のシーンを探します。

利用するシーンが決まったら先ほど定義した関数を呼び出し、合成画像を表示してみます。
今回はpresetとしてnaturalを指定します。

scene_id = "LC81090352018176LGN00"
natural_img = get_landsat_image(scene_id, z, x, y, preset="natural")

io.imshow(natural_img)

下記のような画像が表示されたかと思います。
緑色のところは植物が生育しているエリアですが、赤紫色になるにつれて植物が生育していないエリアとなっています。

ちなみに、元のシーンを確認すると赤紫色がかっている箇所のあたりに奥穂高岳（標高3190メートル）があり、左下にかけて広範囲で赤紫色になっている箇所は西穂高岳（標高2909メートル）があります。右上の白くモコモコとしたものは雲が写り込んでおり山ではないのでご注意ください。

それではこれらの画像のピクセルデータから、指定した列または行の植生分布のRGB値と標高のデータを取得する関数、およびそれらを使って標高ごとの植生分布をプロットする関数を定義します。

def get_colors(tile, index, direction, z, x, y, xtile_size=1, ytile_size=1, u=1):
    """
    タイル画像から植生に応じたRGB値を取得する

    Parameters
    ----------
    tile : ndarray 
        植生タイル
    index : number 
        走査するピクセルのインデックス 
    direction : number 
        走査方向(0で上から下、1で左から右)
    zoom : number 
        ズーム率 
    xtile : number 
        座標X 
    ytile : number 
        座標Y
    xtile_size : number 
        X方向のタイル数 
    ytile_size : number  
        Y方向のタイル数
    u : number 
        標高分解能[m]（デフォルト 1.0[m]）
    """
    if direction == 0: 
        line = tile[:,index]
    else:
        line = tile[index,:]

    colors = []
    colors_err = []
    for k in range(len(line)):
        R = line[k,0]
        G = line[k,1]
        B = line[k,2]
        try:
            colors.append([R, G, B, u])
        except ValueError as e:
            colors.append([0,0,0,0])
            colors_err.append((i,j))
    return colors

def get_tile_bbox(zoom, topleft_x, topleft_y, size_x=1, size_y=1):
    """
    タイル座標からバウンダリボックスを取得
    https://tools.ietf.org/html/rfc7946#section-5
    """
    def num2deg(xtile, ytile, zoom):
        # https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
        n = 2.0 ** zoom
        lon_deg = xtile / n * 360.0 - 180.0
        lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
        lat_deg = math.degrees(lat_rad)
        return (lon_deg, lat_deg)
    
    right_top = num2deg(topleft_x + size_x , topleft_y, zoom)
    left_bottom = num2deg(topleft_x, topleft_y + size_y, zoom)
    return (left_bottom[0], left_bottom[1], right_top[0], right_top[1])

def calc_height_chiriin_style(R, G, B, u=100):
    """
    標高タイルのRGB値から標高を計算する
    """
    hyoko = int(R*256*256 + G * 256 + B)
    if hyoko == 8388608:
        raise ValueError('N/A')
    if hyoko > 8388608:
        hyoko = (hyoko - 16777216)/u
    if hyoko < 8388608:
        hyoko = hyoko/u
    return hyoko


def plot_height(tile_gdem, tile_natural, index, direction, z, x, y, xtile_size=1, ytile_size=1, u=1):
    """
    高度のグラフをプロットする

    Parameters
    ----------
    tile_gdem : ndarray 
        標高タイル
    tile_natural : ndarray 
        natural画像タイル
    bbox : object 
        標高タイルのバウンディングボックス 
    index : number 
        走査するピクセルのインデックス 
    direction : number 
        走査方向(0で上から下、1で左から右)
    zoom : number 
        ズーム率 
    xtile : number 
        座標X 
    ytile : number 
        座標Y
    xtile_size : number 
        X方向のタイル数 
    ytile_size : number  
        Y方向のタイル数
    u : number 
        標高分解能[m]（デフォルト 1.0[m]）
    """

    bbox = get_tile_bbox(z, x, y, xtile_size, ytile_size) 
    width = abs(bbox[2] - bbox[0])
    height = abs(bbox[3] - bbox[1])
    if direction == 0: 
        line = tile_gdem[:,index]
        per_deg = -1 * height/tile_gdem.shape[0]
        fixed_deg = width/tile_gdem.shape[1] * index + bbox[0]
        deg = bbox[3]
        title = 'height at {} degrees longitude'.format(fixed_deg)
    else:
        line = tile_gdem[index,:]
        per_deg = width/tile_gdem.shape[1]
        fixed_deg = -1 * height/tile_gdem.shape[0] * index + bbox[3]
        deg = bbox[0]
        title = 'height at {} degrees latitude'.format(fixed_deg)
        
    heights = []
    heights_err = []
    degrees = []
    for k in range(len(line)):
        R = line[k,0]
        G = line[k,1]
        B = line[k,2]
        try:
            heights.append(calc_height_chiriin_style(R, G, B, u))
        except ValueError as e:
            heights.append(0)
            heights_err.append((i,j))
        degrees.append(deg)
        deg += per_deg

    fig, ax= plt.subplots(figsize=(8, 4))

    colors = np.array(get_colors(tile_natural, index, direction, z, x, y)) / np.array([255, 255, 255, 1])
    plt.title(title)
    plt.scatter(degrees, heights, s=20,c=colors)
    ax.text(0.01, 0.95, 'Errors: {}'.format(len(heights_err)), transform=ax.transAxes)
    plt.show()

ようやく準備が整いました。それでは奥穂高岳を通るように植生分布をプロットしてみたいと思います。

奥穂高岳のあたりを通りつつ、うまく雲を避けながらプロットしてみます。画像でいうと赤い線のあたりが対象とお考えください。

plot_height(gdem_img, natural_img, 10, 1, z, x, y)

下記のようなプロットが表示されたかと思います。

これをみるとおよそ2000メートルあたりから植生が変わり始め、2500メートルを超えると全く植物がみられなくなるというように読み取ることができそうです。

ちなみに、北アルプスには涸沢ヒュッテ（標高2309メートル）と呼ばれるオンシーズンにはテントで埋め尽くされるぐらい有名なキャンプ地があります。

そのキャンプ地の様子が下記になります。

まだ低木による緑があることが見受けられますが、だんだんと高木が少なくなっている様子がみて取れるかと思います。

次に別方向にもプロットして標高と植生の分布の変化をみてみたいと思います。

奥穂高のあたりを通るように、今度は列でプロットしてみます。

plot_height(gdem_img,natural_img, 150, 0, z, x, y)

やはり北アルプスにおいては2000〜2500メートルあたりからが、植生が変化する境界となっていそうですね。

ちなみに、3000メートルを超える奥穂高岳からみた景色が下記になります。

全く植物が生えておらず、荒涼としている様子がお分りいただけるかと思います。

4.その他の地域(北海道)の森林限界も可視化してみる

森林限界は緯度によっても大きく変わります。せっかくなので比較として、北海道にある大雪山系のあたりも標高と植生の関係をプロットしてみたいと思います。

先に使用するシーンを探しておきます。北海道の大雪山系では7月中旬ごろまで残雪が残っているため、7月下旬から8月のシーンを利用するように探します。

長くなってしまったのでプロットの結果のみ掲載させていただきます。

# 大雪山系の中でも最高峰である旭岳周辺の標高データを取得
z2 = 10
x2 = 918
y2 = 373
gdem_img2 = get_astergdem2_dsm(z2, x2, y2)
io.imshow(gdem_img2)

scene_id2 = "LC81070302017207LGN00"
natural_img2 = get_landsat_image(scene_id2, z2, x2, y2, preset="natural")

plot_height(gdem_img2, natural_img2, 90, 0, z2, x2, y2)

このプロット画像をみると、北海道の大雪山系においても2000メートル前後から植生の変化が始まっているように読み取れそうです。しかし、実際には北海道の大雪山系では1000〜1500メートルあたりから森林限界が生じると言われていることと一見矛盾しているように思えます。

しかし、大雪山系の旭岳においては、高木の代わりに低木のハイマツなどの高山植物がよく生育していることが知られており、そのハイマツの存在が1500メートル以上の植生に影響を及ぼしているものと考えられそうです。

これらの高木と低木も区別できると、より森林限界をはっきりと可視化できるかもしれません。

参考：大雪山の大自然～高山植物

5.まとめ

今回はASTER GDEM 2.0の標高データと、Landsat-8の植生分布合成画像を利用してどの程度の標高から森林限界が生じるのかを可視化して見ました。その結果、本州の北アルプスにおいては約2000メートルほどの高度から樹木が見られなっていました。また、より緯度の高い北海道の山脈においてはそれよりも低い高度においても、森林限界が生じることが読み取れました。

今回は私も何度か登ったことのある北アルプスと、比較対象として高緯度地域にある大雪山を取り上げて見ましたが、逆に九州などの南の方の山や、さらに国外の山などとも比較できたりすると面白いかもしれません。Wikipediaの森林限界の説明によれば、回帰線を挟む低緯度の地域では、季節の影響を受けにくいためおよそ3000メートルという高さから森林限界が生じるそうです。

また、地球温暖化など環境の変化によって、森林限界の境界線は移動することが知られています。遥か昔の氷河時代においては大規模な森林限界の垂直移動・水平移動が生じたそうです。衛星データを時系列に沿って比較してみることで、森林限界の移動を可視化することができると、生態系の保全活動などにも役立つかもしれません。

「Tellus」で衛星データを触ってみよう！

日本発のオープン＆フリーなデータプラットフォーム「Tellus」で、まずは衛星データを見て、触ってみませんか？
★Tellusの利用登録はこちらから