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

記事作成時から、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

私たちが普段見ているテレビの映像や写真の画像は、赤・青・緑の3色を用いて表現されています。
これに対して、衛星が撮影している画像は、赤・青・緑以外にも短波長赤外線や近赤外線など、より多くの種類の波長で構成されています。
どのような波長のデータを持っているかは衛星ごとに異なります。

この波長の組み合わせを変えてブレンドした画像を作ることで、いつも通りの衛星データを見るのとはまた異なる情報を読み取れるようになります。

たとえば、正規化水指標（NDWI）を使えば画像の中で水がある場所を濃い色で、正規化積雪指標（NDSI）を使えば画像の中で雪のある場所を濃い色で強調できます。

このように、目的とする対象の度合いを見えるような組み合わせを用いて比較できるよう正規化した指標のことを「正規化指標」と呼びます。

正規化指標とその種類については、以下の記事でわかりやすく紹介されています。

今回はLandsat-8の画像データを用いて、任意の波長を組み合わせ、様々な正規化指標に基づくデータを取得する方法を紹介します。

Landsat-8の画像データを取得し、正規化植生指標（NDVI）で画像中の植物を目立たせる

Tellusのバンド別画像合成APIでは、正規化植生指標（NDVI）を計算するためのプリセットが用意されています。
Landsat-8 API（APIリファレンス）
presetの値にndviを指定することで、NDVIを計算することができます。

ここからは、Landsat-8の衛星データを用いて、NDVIのプリセットを使ったバンド別画像合成APIのサンプルコードを見ていきましょう。

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

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

# Landsat-8のシーン情報を取得するメソッド
def get_landsat8_scene(min_lat, min_lon, max_lat, max_lon):
    url = "https://gisapi.tellusxdp.com/api/v1/landsat8/scene" \
        + "?min_lat={}&min_lon={}&max_lat={}&max_lon={}".format(min_lat, min_lon, max_lat, max_lon)
    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    
    r = requests.get(url, headers=headers)
    return r.json()

# Lantsat-8のバンド別画像合成API実行メソッド
def get_landsat8_image(z, x, y, path, row, productId, params={}):
    url = "https://gisapi.tellusxdp.com/blend/landsat8/{}/{}/{}/{}/{}/{}.png".format(path, row, productId, z, x, y)
    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    
    r = requests.get(url, headers=headers, params=params)
    return io.imread(BytesIO(r.content))

パラメータを指定し、APIを呼び出すためのメソッドを作成します。

まずは座標を指定してLandsat-8のシーン情報を取得し、取得件数を表示してみましょう。

landsat8_scenes = get_landsat8_scene(35.530107, 138.693537, 35.460226, 138.777995)
len(landsat8_scenes)

実行結果

62

シーン情報が62件取得できました。
今回は50番目の画像を利用します。

# 取得したシーンの何番目を使うか指定する
img_num = 50

print(landsat8_scenes[img_num])

# サムネイルを確認
io.imshow(io.imread(landsat8_scenes[img_num]["thumbs_url"]))

{'acquisitionDate': 'Wed, 29 Oct 2014 01:22:10 GMT', 'c_lat': 36.006600000000006, 'c_lon': 138.576, 'cloudCover': 0.26, 'date_added': 'Tue, 25 Jun 2019 00:00:00 GMT', 'entityId': 'LC81080352014302LGN01', 'max_lat': 37.0795, 'max_lon': 139.91, 'min_lat': 34.9337, 'min_lon': 137.242, 'path': 108, 'productId': 'LC08_L1TP_108035_20141029_20170418_01_T1', 'row': 35, 'thumbs_url': 'https://s3-us-west-2.amazonaws.com/landsat-pds/c1/L8/108/035/LC08_L1TP_108035_20141029_20170418_01_T1/LC08_L1TP_108035_20141029_20170418_01_T1_thumb_small.jpg', 'tile_path': 'landsat8/108/035/LC08_L1TP_108035_20141029_20170418_01_T1/{z}/{x}/{y}.png'}

続いて、バンド別画像合成APIの呼び出しに必要な情報と、画像のうち見たい箇所のズーム率とタイル番号を指定します（タイル番号については こちらの記事 をご参考ください）。

# ズーム率
zoom = 12

# タイル番号
tile_num = (3626,1615)

# API呼び出しパラメータを設定
path = landsat8_scenes[img_num]['path']
row = landsat8_scenes[img_num]['row']
productId = landsat8_scenes[img_num]['productId']

指定したタイルには、河口湖が写っています。
画像右側には湖と市街地があり、湖の周囲には植物があります。

最後に、プリセットをNDVIに設定してバンド別画像合成APIを実行しましょう。

# presetをndviに設定
params = {
    'preset' : 'ndvi'
}
img = get_landsat8_image(zoom, tile_num[0], tile_num[1], path, row, productId, params)
io.imshow(img)

こちらのコードを実行すると、以下のような画像が出力されます。

先ほど見た画像の中で、湖と市街地で示した部分が赤っぽく表示される一方、その周囲は緑色に表示されています。
この緑色の部分が、正規化指標のNDVIで強調表示された植物です。

このように、NDVIはAPIでプリセットが用意されているため、簡単にデータを分析することができます。
現在、プリセットが用意されている正規化指標はNDVIのみですが、バンド別画像合成APIでは、クエリのr,g,bの値を設定することで、任意の波長の組み合わせを取得できます。

波長を指定してLandsat-8の画像データを取得する

正規化指標を計算するためのメソッドを新しく作ります。
このメソッドは、band1とband2の値を受け取り、(band1 – band2 / band1 + band2)の式を実行した結果を返します。

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)

次に、必要な波長を指定してデータを取得しましょう。
先程のサンプルコードで作成した、Landsat-8のバンド別画像合成APIを実行するメソッドを以下に再掲します。

def get_landsat8_image(z, x, y, path, row, productId, params={}):
    url = "https://gisapi.tellusxdp.com/blend/landsat8/{}/{}/{}/{}/{}/{}.png".format(path, row, productId, z, x, y)
    headers = {
        "Authorization": "Bearer " + TOKEN
    }
    
    r = requests.get(url, headers=headers, params=params)
    return io.imread(BytesIO(r.content))

params={}という引数に、先程はNDVIのプリセットを指定しました。
プリセットを使わずに任意の波長を取得する場合は、paramsの引数に、r,g,bの各パラメータに波長のバンドを設定することができます。

ここが今回のポイントです。
利用したい正規化指標に合わせて、APIで取得する波長を設定します。

正規化植生指標（NDVI）は、paramsの値に近赤外(NIR)と赤(Red)を設定することでも計算することができます。
Landsat-8では赤は4、近赤外は5にマッピングされているので、rの値に4、gの値に5をセットします。

# NDVI
# Red(4)とNIR(5)を設定する
params = {
    'r': 4,
    'g': 5
}

どの波長がどのバンドに入ってるかは、Landsat-8のデータ詳細ページに記載されています。
Landsat-8のデータ詳細

最後に、作成した関数を実行し、NDVIの画像データを表示します。

# Landsat-8のバンド別画像合成APIを実行する
img = get_landsat8_image(zoom, tile_num[0], tile_num[1], path, row, productId, params)

# 正規化指標の計算をする
calc_img = calc_normalized_index(img[:,:,1], img[:,:,0])

io.imshow(calc_img)

画像中の植物が強調されて青く表示されています。
一方、植物がない湖や市街地の部分は濃い赤色になっています。

プリセットを用いて取得した画像と比較してみましょう。
正規化指標により植物が強調表示されており、その箇所が一致していることがわかります。

このように、APIで取得した結果を用いた計算でも、プリセットを利用した場合と同様の結果を求めることができます。

さまざまな指標を計算する（NDWI,NDSI）

波長のバンド値を変えることで、他の正規化指標も計算することができます。

例えば、正規化水指標（NDWI）で地表面における水域を強調表示することができます。
この指標は水を目立たせる指標なので、海と川が写っている場所を見てみましょう。

# 東京湾と荒川
zoom = 11
tile_num = (1819,806)

path = 107
row = 35
productId = 'LC08_L1TP_107035_20170320_20170328_01_T1'

正規化水指標（NDWI）では赤(R)と短赤外(SWIR)の波長を使います。
NDWIを求める際の短赤外の反射率は1560-1660nm(SWIR2)が適当なため、短赤外にはBand6を利用します。
rの値に4、gの値に6をセットして、NDWIを計算してみましょう。

# NDWI 
# Red(4)とSWIR2(6)を設定する
params = {
    'r': 6,
    'g': 4
}

最後に作成した関数を実行し、NDWIの画像データを表示します。

img = get_landsat8_image(zoom, tile_num[0], tile_num[1], path, row, productId, params)
calc_img = calc_normalized_index(img[:,:,1], img[:,:,0])
io.imshow(calc_img)

NDWIを利用した画像が取得できました。
この画像では、水域である東京湾と荒川の部分が青色に強調されています。

同様に、正規化土壌指標（NDSI）を調べてみましょう。
地表面における砂やコンクリートなどの分布を示す指標のため、砂が多い土地、鳥取砂丘を見てみます。

# 鳥取砂丘が含まれるタイル
zoom = 12
tile_num = (3575,1614)

path = 111
row = 35
productId = 'LC08_L1TP_111035_20180607_20180615_01_T1'

正規化土壌指標（NDSI）を調べるためには、短赤外(NRI)と近赤外(SWIR)波長を使います。
NDSIを求める際の短赤外にはBand6を利用します。paramsの値を書き換えて実行してみましょう。

# NDSI 
# NIR(5)とSWIR2(6)を設定する
params = {
    'r': 5,
    'g': 6
}
img = get_landsat8_image(zoom, tile_num[0], tile_num[1], path, row, productId, params)
calc_img = calc_normalized_index(img[:,:,1], img[:,:,0])
io.imshow(calc_img)

こちらも、鳥取砂丘が濃い青色で強調されました。
砂浜がある海岸線も青く表示されています。

今回はLandsat-8のデータを用いて、正規化指標の計算方法を紹介しました。

これ以外にも様々な正規化指標があり、今回紹介したバンド別合成APIのパラメータを書き換えることで計算することができます。
正規化指標の使い方を覚えて、衛星データをより詳しく分析してみましょう！

参考資料
リモートセンシングによる土地被覆分類へのSOM適応