記事作成時から、Tellusからデータを検索・取得するAPIが変更になっております。該当箇所のコードについては、以下のリンクをご参照ください。
https://www.tellusxdp.com/ja/howtouse/access/traveler_api_20220310_
firstpart.html

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

本記事ではTellusのJupyter Notebookを使ってAVNIR-2/ALOSによる光学画像から緑地を切り出す方法を紹介します。

TellusでJupyter Notebookを使う方法は、「Tellusの開発環境でAPI引っ張ってみた」をご覧ください。

(1)AVNIR-2のAPIでNDVI画像を取得する

以前の記事「【ゼロからのTellusの使い方】Jupyter NotebookでAVNIR-2/ALOSによる光学画像を取得する」では、AVNIR-2/ALOSが取得した光学画像をTellusのJupyter Notebook上で利用する方法を紹介しました。

このときはTrue ColorとNatural Colorの画像合成を行いましたが、AVNIR-2のデータはこれ以外にも様々な方法で活用することができます。

今回は地上の植物の繁殖度合いを示す植生活性度(NDVI)を色で表現した画像を取得してみましょう。

一般的に植物は、赤の波長の光ほどよく吸収し、近赤外の波長の光はよく反射する性質があります。

この性質を利用したのがNDVIです。

NDVIを自分で計算して画像を作ることもできますが、APIにpreset=ndviと与えるだけでNDVI合成の画像を簡単に取得することができます。

この解析ではopenCVというライブラリを利用します。

TellusのJupyter Notebook上でopenCVを使うにはインストール作業が必要です。

画面の右上の「New」から「Terminal」をクリックして黒い画面（ターミナル）を開きます。

そこで「pip install opencv-python」を実行します。

「Successfully installed opencv-python-バージョン名」と出れば成功です。

準備が完了したので、Pythonのサンプルコードを実行してみましょう。

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

TOKEN = "ここには自分のアカウントのトークンを貼り付ける"

def get_AVNIR_image(zoom, xtile, ytile, opacity=1, r=3, g=2, b=1, rdepth=1, gdepth=1, bdepth=1, preset=None):
    url = "https://gisapi.tellusxdp.com/blend/{}/{}/{}.png?opacity={}&r={}&g={}&b={}&rdepth={}&gdepth={}&bdepth={}".format(zoom, xtile, ytile, opacity, r, g, b, rdepth, gdepth, bdepth) 
    if preset is not None:
        url += '&preset='+preset
    headers = {
        "content-type": "application/json",
        "Authorization": "Bearer " + TOKEN
    }
    r = requests.get(url, headers=headers)
    return io.imread(BytesIO(r.content))

true_img1 = get_AVNIR_image(zoom=13,xtile=7276,ytile=3226)
ndvi_img1 = get_AVNIR_image(zoom=13,xtile=7276,ytile=3226,preset='ndvi')
io.imshow(np.hstack((true_img1, ndvi_img1)))

「Run」をクリックすると、東京の築地周辺の光学画像が表示されます。

NDVI合成では植物の繁殖度合いが高いほど緑に、低いほど赤に描画されます。

そのため、画像中央の浜離宮恩賜庭園や左端の芝公園、左上の日比谷公園の辺りが、NDVI合成では黄緑色に光って見えています（※少し色の変わり目がわかりにくいかもしれません）。

(2)HSV変換

緑地を切り出したければ、NDVI合成の黄色や緑の部分だけを抽出すればOK。色の抽出によく使われるのがHSV変換です。

RGBが色を赤、緑、青の3つの情報で表現するのに対し、HSVは色相（色の種類）、彩度（鮮やかさ）、明度（明るさ）の3つで色を表現します。

※より詳しくは「HSV色空間」で検索してみてください。

色の変換や抽出にはOpenCVという画像処理ライブラリを使ってみましょう。OpenCVでは色相は0 – 180の範囲で処理されます。

また彩度と明度は0 – 255の範囲になります。

def masking(img, hsv_min, hsv_max):
    # OpenCVで扱うために色順をRGB(RGBA)からBGRに入れ替える
    cv_bgr = cv2.cvtColor(img, cv2.COLOR_RGBA2BGR) 
    # BGRからHSVに変換
    cv_hsv = cv2.cvtColor(cv_bgr, cv2.COLOR_BGR2HSV)
    
    # 指定した範囲の色以外を黒に変換する
    mask = cv2.inRange(cv_hsv, hsv_min, hsv_max)
    masked_img = cv2.bitwise_and(cv_bgr, cv_bgr, mask=mask)

    # 指定した範囲の色の割合
    rate = len(np.where(mask.ravel()!=0)[0])/len(mask.ravel())

    # 得られた画像の色順をBGRからRGBに戻す
    return cv2.cvtColor(masked_img, cv2.COLOR_BGR2RGB), rate


green_min = np.array([25,0,0]) 
green_max = np.array([75,255,255]) 

これが黄色や緑の部分だけを抽出するための関数です。今回は色相が25から75の範囲を抽出します。

では、この関数にNDVI合成した画像を渡し、結果を表示しましょう。

masked_img1, rate = masking(ndvi_img1, green_min, green_max)
io.imshow(masked_img1)
print(rate)

黄緑色の部分だけが抽出され、浜離宮恩賜庭園や芝公園、日比谷公園、また埋立地の護岸の植物がくっきりとわかるようになりました。

また、この画像内の植物の割合は約9%(0.085273……)と表示されていることが分かります。

(3)郊外の画像を見てみよう

都心部と比較して郊外はどのような色になるかも見てみましょう。

true_img2 = get_AVNIR_image(zoom=13,xtile=7253,ytile=3230)
io.imshow(true_img2)

これは同じズーム率での河口湖周辺の座標です。

中央に位置する湖と、その周りの居住地や木々の生い茂った山間部の境界が、True Color合成でもはっきりと分かります。

ここを先ほどと同様にNDVI合成で取得して、植物の生えている場所を抽出してみましょう。

ndvi_img2 = get_AVNIR_image(zoom=13,xtile=7253,ytile=3230,preset='ndvi')
io.imshow(masking(ndvi_img2, green_min, green_max))

浜離宮恩賜庭園の周辺と比べて濃い緑の部分が多いですね。やはり郊外のほうが植生が活性化している、つまり植物の繁殖が活発なことがひと目でわかります。

居住地の辺りは黄緑がまだら模様になっています。

これは居住地に点在する畑や屋敷林が抽出されているためで、公園以外の土地のほとんどが黒くフィルターされていた都心部との大きな違いです。

画像内の植物の割合は約58%で、都心部の約9%と比べて自然が豊かなことが感覚だけでなく数字からもわかります。

以上、Jupyter NotebookでAVNIR-2/ALOSの光学画像から緑地を切り出し、画像内の植物の割合を算出する方法でした。

今回の抽出はあくまで一例です。

「NDVI合成の黄色や緑の部分が緑地である」と仮定して、

green_min = np.array([25,0,0]) 
green_max = np.array([75,255,255]) 

色相の下限を25、上限を75と設定して抽出を行いましたが、精度をより高めるにはこの値の調整が必要となります。

また、NDVI合成とHSV変換を用いるよりももっと効果的な抽出条件があるかもしれません。

見たいモノや目的、求める精度が変われば、最適な抽出条件も変わります。それを見つけるのが画像処理の醍醐味でもあり大変なところです。

機械学習で抽出条件を調べるのも良いでしょう。緑地の切り出しを足がかりに、より高度な衛星データの利用法を模索してみてください。