AVNIR-2/ALOSのAPIを用いて波長別画像を取得する

Tag

はじめに

Tellusでは、ALOS(JAXAの陸域観測技術衛星だいち)に搭載された光学センサであるAVNIR-2が取得した画像を利用することができます。AVNIR-2は、2006年から可視光領域から近赤外領域にかけての波長で地上を観測していました(2011年に運用停止)。

本チュートリアルではTellusのJupyterLabを使ってAVNIR-2/ALOSによる光学画像から緑地を切り出す方法を紹介します。APIの使い方はリファレンスにまとめられています。

※画像を利用する際は、利用ポリシーを確認して、規約を違反しないように注意してください。

https://gisapi.tellusxdp.com/api/v1/blend/{z}/{x}/{y}.png?opacity={opacity}&r={r}&g={g}&b={b}&rdepth={rdepth}&gdepth={gdepth}&bdepth={bdepth}&preset={preset}
Name Input Description
z number (path) 必須 タイルマップのZ座標 必須項目 デフォルトなし ( 8 ~ 16 )
x number (path) 必須 タイルマップのX座標 必須項目 デフォルトなし
y number (path) 必須 タイルマップのY座標 必須項目 デフォルトなし
preset string (query) 任意 true / false / natural / ndvi のいずれかを指定 デフォルトなし
opacity number (query) 任意 透過度 0で完全に透過、1で不透過 デフォルト1.00 ( 0.00 ~ 1.00 )
r number (query) 任意 RGBカラーのRに入るバンドの番号を指定 デフォルトr,g,b = 4,3,2 ( 1,2,3,4,5,6 )
g number (query) 任意 RGBカラーのGに入るバンドの番号を指定 デフォルトr,g,b = 4,3,2 ( 1,2,3,4,5,6 )
b number (query) 任意 RGBカラーのBに入るバンドの番号を指定 デフォルトr,g,b = 4,3,2 ( 1,2,3,4,5,6 )
rdepth number (query) 任意 Rの濃さを指定 デフォルト1.0 ( 0.00 ~ 2.00 )
gdepth number (query) 任意 Gの濃さを指定 デフォルト1.0 ( 0.00 ~ 2.00 )
bdepth number (query) 任意 Bの濃さを指定 デフォルト1.0 ( 0.00 ~ 2.00 )

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

「AVNIR-2/ALOSの光学画像を取得」ではAPIを用いて地図タイル画像と波長合成画像を取得しました。今回は取得した波長合成画像を利用して、画像全体から見たとき、植物の被覆が全体の何割を占めるのかを見積もります。

一般的に植物は、赤の波長の光をよく吸収し、近赤外の波長の光はよく反射する性質があります。この性質を利用した指標をNDVIと呼びます。

devAvnir_api2_20200220_6.png

この値が1に近いほどその地点は植物が繁殖しています。

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

この解析ではOpenCVというライブラリを利用します。TellusのJupyterLab上でOpenCVを使うにはインストール作業が必要です。

画面の左下の「Terminal」をクリックして操作画面(ターミナル)を開きます。

devAvnir_api2_20200224_1.png

そこで「pip install opencv-python」を実行します。「Successfully installed opencv-python-バージョン名」と出れば成功です。

devAvnir_api2_20200220_7.png

準備が完了したので、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 = {
        "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)))

devAvnir_api2_20200220_3.png

ツールにある▶マークをクリック、もしくはCtrl+Enterキーを押すと、東京の築地周辺の光学画像が表示されます。

devAvnir_api2_20200220_2.png

左側がTrue color合成、右側がNDVI合成
Credit:Original data provided by JAXA

NDVI合成では植物の繁殖度合いが高いほど緑に、低いほど赤に描画されます。そのため、画像中央の浜離宮恩賜庭園や左端の芝公園、左上の日比谷公園の辺りが、NDVI合成では黄緑色に光って見えています(※少し色の変わり目がわかりにくいかもしれません)。

HSV変換

緑地を切り出したければ、NDVI合成の黄色や緑の部分だけを抽出します。色の抽出によく使われるのがHSV変換です。RGBが色を赤、緑、青の3つの情報で表現するのに対し、HSVは色相(色の種類)、彩度(鮮やかさ)、明度(明るさ)の3つで色を表現します。

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

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

水色 ピンク
0 30 60 90 120 150

また彩度と明度は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)
    # 指定した範囲の色の割合
    #.ravelのメソッドで高次元配列は一次配列に置き換えられ、値がゼロでない部分のみを抽出している
    # len()でその値の数を示し、全体のピクセルの数で割り算をし、割合とする
    rate = len(np.where(mask.ravel()!=0)[0])/len(mask.ravel())
    # 得られた画像の色順をBGRからRGBに戻す
    return cv2.cvtColor(masked_img, cv2.COLOR_BGR2RGB), rate

これが黄色や緑の部分だけを抽出するための関数です。今回は色相が25から75の範囲を抽出します。この関数にNDVI合成した画像を渡し、結果を表示しましょう。

# 色相が25から75の範囲で抽出を行う
green_min = np.array([25,0,0]) 
green_max = np.array([75,255,255])
masked_img1, rate = masking(ndvi_img1, green_min, green_max)
io.imshow(masked_img1)
print(rate)

devAvnir_api2_20200220_4.png

Credit:Original data provided by JAXA

黄緑色の部分だけが抽出され、浜離宮恩賜庭園や芝公園、日比谷公園、また埋立地の護岸の植物が明瞭になりました。また、この画像内の植物の割合は約9%(rate = 0.085273……)と求められていることが分かります。

郊外の画像を見てみよう

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

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

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

Credit:Original data provided by JAXA

中央に位置する湖と、その周りの居住地や木々の生い茂った山間部の境界が、True Color合成でもはっきりと分かります。先ほどと同様にNDVI合成で取得して、植物の生えている場所を抽出してみましょう。

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

devAvnir_api2_20200220_8.png

Credit:Original data provided by JAXA

浜離宮恩賜庭園の周辺と比べて濃い緑の部分が目立ちます。郊外のほうが植物の繁殖が活発であることが言えるでしょう。 また、居住地の辺りは黄緑がまだら模様になっています。これは居住地に点在する畑や屋敷林が抽出されているためで、公園以外の土地のほとんどが黒くフィルターされていた都心部との大きな違いです。画像内の植物の割合は約58%で、都心部の約9%と比べて自然が豊かなことが感覚だけでなく数字からもわかります。 以上、JupyterLabでAVNIR-2/ALOSの光学画像から緑地を切り出し、画像内の植物の割合を算出する方法でした。今回の抽出はあくまで一例です。NDVI合成の黄色や緑の部分が緑地である」と仮定して、

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

色相の下限を25、上限を75と設定して抽出を行いましたが、精度をより高めるにはこの値の調整が必要となります。また、NDVI合成とHSV変換を用いるよりももっと効果的な抽出条件があるかもしれません。 見たいモノや目的、求める精度が変われば、最適な抽出条件も変わります。それを見つけるのが画像処理の醍醐味でもあり大変なところです。 機械学習で抽出条件を調べるのも良いでしょう。緑地の切り出しを足がかりに、より高度な衛星データの利用法を模索してみてください。 今回使用したスクリプトです。

import os, json, requests, math
import numpy as np
import cv2
from skimage import io
from io import BytesIO
import matplotlib.pyplot as plt
get_ipython().run_line_magic('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 = {
                "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)))
# 指定した色の範囲内(hsv_min,hsv_max)でハイライトする
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())
    return cv2.cvtColor(masked_img, cv2.COLOR_BGR2RGB), rate
green_min = np.array([25,0,0])
green_max = np.array([75,255,255])
masked_img1, rate = masking(ndvi_img1, green_min, green_max)
io.imshow(masked_img1)
print(rate)
true_img2 = get_AVNIR_image(zoom=13,xtile=7253,ytile=3230)
io.imshow(true_img2)
ndvi_img2 = get_AVNIR_image(zoom=13,xtile=7253,ytile=3230,preset='ndvi')
masked_img2,rate = masking(ndvi_img2,green_min,green_max)
io.imshow(masked_img2)
print(rate)