ASNARO-1の波長合成画像を取得する

Tag

はじめに

Tellusでは、ASNARO-1(一般財団法人宇宙システム開発利用推進機構とNECが開発した衛星)に搭載された詳細撮影を目的とした光学センサで、解像度0.5mの白黒画像と2mのカラー画像を提供しています。ASNARO-1は、2014年から可視光領域から近赤外領域にかけての波長で地上を観測しています(2020年時点で運用中)。

ASNARO-1の光学画像を取得するでは、シーン検索APIを利用して衛星のシーンデータを取得し、さらに地図タイル画像APIにより、地図タイル画像を取得することを学びました。このチュートリアルを通してAPIを利用してASNARO-1の波長合成画像を取得する方法を学びます。

ASNARO-1の波長合成画像をAPIを利用して取得する

波長合成画像APIで衛星データに含まれるproductIdを知る必要があります。地図タイル画像を取得した方法と同様に、シーン検索APIを利用してシーンデータを入手しましょう。 APIの使い方はリファレンスにまとめられています。 ※画像を利用する際は、利用ポリシーを確認して、規約を違反しないように注意してください。

https://gisapi.tellusxdp.com/api/v1/asnaro1/{productid}/{z}/{x}/{y}.png?opacity={opacity}&r={r}&g={g}&b={b}&rdepth={rdepth}&gdepth={gdepth}&bdepth={bdepth}&preset={preset}

devAsnaro_api1_20200220_table1.png

import os, json, requests, math
from skimage import io
from io import BytesIO
import matplotlib.pyplot as plt
%matplotlib inline
TOKEN = "ここにトークンを貼る"
def get_ASNARO_scene(min_lat, min_lon, max_lat, max_lon):
	url = "https://gisapi.tellusxdp.com/api/v1/asnaro1/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()
scenes = get_ASNARO_scene(20.425278, 122.933611, 45.557222, 153.986389)
ext_scene = scenes[7]
print(ext_scene)

定義した関数を利用して任意の範囲のシーンを取得しました。今回は8番目に含まれるシーンをext_sceneの変数に入れてあります。実行した結果は以下のようになります。

devAsnaro_api1_20200220_3.png

出力結果に含まれるproductIdを波長合成画像APIで用います。thumbs_urlを利用してサムネイル画像を表示してみましょう。

img_thumbs = io.imread(ext_scene['thumbs_url'])
io.imshow(img_thumbs)

devAsnaro_api1_20200220_4.png

Credit:Original data provided by NEC

波長画像合成APIでは、地図タイルを取得したときと同様にズーム率、タイルのX座標とY座標を知る必要があります。前回使用したものと同じ関数を定義して、シーンに含まれるタイル画像を取得しましょう。

def get_tile_num(lat_deg, lon_deg, zoom):
	# https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
	lat_rad = math.radians(lat_deg)
	n = 2.0 ** zoom
	xtile = int((lon_deg + 180.0) / 360.0 * n)
	ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
	return (xtile, ytile)
zoom = 16
(xtile, ytile) = get_tile_num(ext_scene['clat'], ext_scene['clon'], zoom)
print(xtile, ytile)

ズーム率16でシーンの中心座標に基づいたタイル番号は(56604, 26952)と求められました。 これにより波長合成画像を取得する準備が整いました。関数を定義します。

# デフォルトはTrue color画像
def get_asnaro1_blend(productId,zoom, xtile, ytile, opacity=1, r=4, g=3, b=2, rdepth=1, gdepth=1, bdepth=1, preset=None):
	url = "https://gisapi.tellusxdp.com/blend/asnaro1/{}/{}/{}/{}.png?opacity={}&r={}&g={}&b={}&rdepth={}&gdepth={}&bdepth={}".format\
	(productId, 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))

urlの変数を見ると、タイル画像を取得した場合のURLより多くの情報が与えられていることがわかります。特に重要となるのは、r、g、bの部分です。これらはRGBを示し、ここに衛星のバンドを与えることにより、どのバンドをどの色に当てはめるのかを指定することができます。ASNARO-1のバンドは以下のようになっています。

devAsnaro_api1_20200220_table2.png

関数を実行し波長合成画像を取得しましょう。

ext_scene_img_true = get_asnaro1_blend(ext_scene['productId'],16, 56604, 26952)
io.imshow(ext_scene_img_true)

devAsnaro_api1_20200220_6.png

Credit:Original data provided by NEC

True color画像を取得することができました。耕作地が森林で囲まれていることが確認できます。この画像では私たちの目で捉えた世界と同じ色合いで画像が描写されていることがわかります。 続いてNatural color画像を取得しましょう。Natural color画像では、rに赤の波長(Band4)、gに近赤外(Band6)、bに緑(Band3)を割り当てます。

ext_scene_img_natural = get_asnaro1_blend(ext_scene['productId'],16, 56604, 26952, r=4,g=6,b=3)
io.imshow(ext_scene_img_natural)

devAsnaro_api1_20200220_2.png

Credit:Original data provided by NEC

森林部分が鮮やかな緑色に見えています。このように、Natural colorでは植物が緑に強調されることがわかります。植物は近赤外領域の波長を強く反射しますので、植物は緑色を呈します。 さらにFalse color画像を取得します。今回はrに近赤外の波長(Band6)、gに赤の波長(Band4)、そしてbに緑の波長(Band3)を割り当てます。

ext_scene_img_false = get_asnaro1_blend(ext_scene['productId'],16, 56604, 26952, r=6,g=4,b=3)
io.imshow(ext_scene_img_false)

devAsnaro_api1_20200220_1.png

Credit:Original data provided by NEC

False color画像では植物が赤く表示されます。こちらでは近赤外領域の波長を赤に割り当てているためです。 最後にNDVIを用いた図を取得しましょう。今回はpresetを使っていますが、波長を指定した画像を二つ用意し、その値を直接計算することでも得ることができます。

ext_scene_img_ndvi = get_asnaro1_blend(ext_scene['productId'],16, 56604, 26952, preset = 'ndvi')
io.imshow(ext_scene_img_ndvi)

devAsnaro_api1_20200220_5.png

Credit:Original data provided by NEC

森が広がった部分が緑色で表示されていることが分かります。耕作地の一部も緑色を呈しています。このようにNDVIを用いれば、同じ耕作地であっても既に作物が育っているのか否かを判別することも可能です。 波長合成画像を得ることで、見たいもの知りたいことを強調して表示できます。強調されたピクセルは特定の範囲のピクセル値を持つため、他のピクセルと容易に区別できます。そのため、特定のピクセルを抽出すれば、画像の中に目的の対象物がどのくらい含まれているのかも知ることができます。 今回示した波長合成画像は、あくまで一部のものです。対象物の反射特性を知れば、それに当てはめる波長を変化させ、私たちの目では捉えられない、多様な世界の姿を表現することが出来ます。 波長合成APIを使って、画像の変化を楽しんで下さい。 今回使用したスクリプトです。

import os, json, requests, math
from skimage import io
from io import BytesIO
import matplotlib.pyplot as plt
%matplotlib inline
TOKEN = "ここにトークンを貼る"
def get_ASNARO_scene(min_lat, min_lon, max_lat, max_lon):
	url = "https://gisapi.tellusxdp.com/api/v1/asnaro1/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()
scenes = get_ASNARO_scene(20.425278, 122.933611, 45.557222, 153.986389)
ext_scene = scenes[7]
print(ext_scene)
img_thumbs = io.imread(ext_scene['thumbs_url'])
io.imshow(img_thumbs)
def get_tile_num(lat_deg, lon_deg, zoom):
	# https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
	lat_rad = math.radians(lat_deg)
	n = 2.0 ** zoom
	xtile = int((lon_deg + 180.0) / 360.0 * n)
	ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
	return (xtile, ytile)
zoom = 16
(xtile, ytile) = get_tile_num(ext_scene['clat'], ext_scene['clon'], zoom)
print(xtile, ytile)
# デフォルトはTrue color画像
def get_asnaro1_blend(productId,zoom, xtile, ytile, opacity=1, r=4, g=3, b=2, rdepth=1, gdepth=1, bdepth=1, preset=None):
	url = "https://gisapi.tellusxdp.com/blend/asnaro1/{}/{}/{}/{}.png?opacity={}&r={}&g={}&b={}&rdepth={}&gdepth={}&bdepth={}".format\
	(productId, 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))
ext_scene_img_true = get_asnaro1_blend(ext_scene['productId'],16, 56604, 26952)
io.imshow(ext_scene_img_true)
ext_scene_img_natural = get_asnaro1_blend(ext_scene['productId'],16, 56604, 26952, r=4,g=6,b=3)
io.imshow(ext_scene_img_natural)
ext_scene_img_false = get_asnaro1_blend(ext_scene['productId'],16, 56604, 26952, r=6,g=4,b=3)
io.imshow(ext_scene_img_false)
ext_scene_img_ndvi = get_asnaro1_blend(ext_scene['productId'],16, 56604, 26952, preset = 'ndvi')
io.imshow(ext_scene_img_ndvi)