はじめに

この記事は、シスコシステムズ合同会社の同士による Cisco Advent Calendar 2019 の 21 日目として投稿しています。
そもそも Cisco と衛星データってどいういう繋がりがあるの？という疑問もあると思いますのでそのあたりも含めてご紹介しようと思います。

この記事でできるようになること

衛星データとは何かがわかる
衛星データプラットフォーム Tellus とは何かがわかる
Tellus 開発環境の利用方法がわかる
Tellus 開発環境を利用して衛星データの解析ができる
NDVI, NDWIとは何かがわかる

最終的にこんなことができます

高精細な衛星画像を日付,位置を選択して表示
Tue, 31 Jul 2018 12:42:49 GMT, 北緯35.675 - 35.605 : 東経139.745 - 139.835
衛星画像から植生状況（緑地部分）を可視化
植生状況を数値化(NDVI)

衛星データプラットフォーム Tellus とは？

公式ページから

Tellus（テルース）は、政府衛星データを利用した新たなビジネスマーケットプレイスを創出することを目的とした、日本初のオープン＆フリーな衛星データプラットフォームです。複数のデータをかけ合わせ、新たなビジネス創出を促進するためのあらゆるファンクションを提供します。
https://www.tellusxdp.com/ja/about

これまでの衛星データ活用では1枚あたり数十万円するような画像データを特殊なデータフォーマットで扱わなければなりませんでした。価格・ストレージ・コンピューティングリソースの観点で様々な課題があり、これまで豊富な衛星データがありながら、活用はとても限定的でした。これらの課題を解決するために経済産業省が2020年までの3年間のプロジェクトとして立ち上げたのが"政府衛星データのオープン&フリー化及びデータ利活用促進事業"その中止的存在となるものがTellusとなります。

特徴

すべてのデータが無償で利用可能（上限あり）
衛星データだけでなく地上データ（人流など）も利用可能
データだけでなく解析リソース（CPU/GPU/Storage）も利用可能
マーケットプレイスによりデータやアプリケーション、アルゴリズムの取引が可能
活用支援のためのトレーニングプログラムが充実

Ciscoとの関連性

Ciscoと宇宙ビジネスについての関連性ですが、Ciscoは2003年（実験衛星）、2009年（商用衛星）に宇宙空間へ度々IPルータを投入しています。衛星には、気象衛星、観測衛星、測位衛星、通信衛星の4種類が存在します。この内Ciscoが過去に取り組んだのは通信衛星となります。
https://www.cisco.com/c/ja_jp/about/newsroom/archives/us/2010/prod-072010b.html

TellusではxData Allianceとして様々な領域で様々な事業者・研究機関・団体が参画しており、Ciscoも2019年2月21日から参画しています。今回は上記の種類でいうと観測衛星を使用した新たなビジネス創出を目的としています。

https://www.xdataalliance.com

また、最近ではTokyo Moonshot Challengeという新しい都市体験を提供するサービスデザインのアイデア創出・具現化コンテストも行っていますので、ご興味のある方はご参加いただければと思います。

Tellus OS を利用してみよう

それでは実際にTellusを使用して衛星データを活用してみましょう。
利用の手順は以下の通りとなります。
1. Tellusへのユーザー登録
2. Tellus OSへのログイン
3. Tellus 開発環境へのアクセス
4. 使用する言語を選択
- R及びPythonが使用できますが、今回はPythonを使用します
5. Jupyter NotebookでCodeを記述
6. 実行及び結果の確認

Tellusにユーザー登録をしてみよう

Tellusにアクセスし、必要情報を入力してユーザー登録を行います。

若干入力箇所が多いかもしれませんが、めげないで記入しましょう。

Tellus 開発環境を利用してみよう

ユーザー登録が終了したら、実際にTellusを利用してみましょう。Tellusには2種類の利用環境が存在します。
まずはTellus Webページの右上かその下のボタンからTellus OSにログインしましょう。

1. Tellus OS

GUIを利用して直感的な操作が可能な環境が用意されております。
プリセットでデータが用意されていますので、大まかであればこちらを利用しての解析が可能です。こちらである程度解析対象を絞り込んでから、もし必要であれば下記の開発環境で詳細な解析を実施するという流れです。

2. Tellus 開発環境

PythonやRから利用できるAPIが用意されております。APIだけでなくComputingやStorage,GPUなどの環境も提供されております。開発環境へはTellus OS上からアクセス可能です。Tellus OS上の以下のボタンをクリックして開発環境にアクセスできます。

なお、開発環境へのアクセスは申請が必要です。

開発環境申請方法

マイページにアクセスします
マイページの開発環境申請で必要な開発環境を選択して申請します。
- 現在申請が集中しておりますためお時間がかかる場合があります

API TOKENの登録

開発環境が利用可能になりましたらAPI TOKENの登録をします。マイページから登録が可能です。
Keyに好きな名前をつけましょう。

Jupyter Note Book の起動

Tellus OSにログインして、開発環境に入りましょう。

実際に Code を書いてみよう

ログインしたらWorkディレクトリの下に移動します。
Python3 Notebookを新たに作りましょう。

Jupyter Notebookの説明はこちらでは省略しますが、使用されたことのある方にはおなじみの環境ですね。
今回使用するモジュールは以下のとおりです。はじめにこれらをImportします。

名称	概要
Numpy	Pythonの数値計算を行うライブラリ
Scikit-image	Pythonの画像解析ライブラリ
Requests	Python用Apache2 Licensed ベースのHTTP通信ライブラリ
Io	I/Oを扱うPythonの標準モジュール
PIL	Pythonの画像処理ライブラリ
Matplotlib	Pythonのグラフ描画ライブラリ
math	Pythonの計算用ライブラリ

# 必要なモジュール、ライブラリをインポート
import numpy as np 
from skimage import io, color, img_as_ubyte, filters, transform
import requests
from io import BytesIO
from PIL import Image
import matplotlib.pyplot as plt
%matplotlib inline
import math

ここで、今回使用する衛星のご紹介をいたします。

衛星名称	ASNARO-1
データプロバイダ	PASCO CORPORATION
搭載期間	2019年〜2019年
地表面解像度	0.5m(パンクロ)、2m(マルチ)、0.5m(パンシャープン)
利用例	高分解能光学衛星データを用いた水域における船舶検出
波長	パンクロ（白黒画像）:450-860 nm マルチ（カラー画像）: Band1：400-450 nm (Ocean Blue) Band2：450-520 nm (Blue) Band3：520-600 nm (Green) Band4：630-690 nm (Red) Band5：705-745 nm (Red Edge) Band6：760-860 nm (NIR)

アクセストークンを指定します。

# APIアクセストークンをAPI_TOKEN に格納
API_TOKEN = 'YOUR TOKEN HERE'

今回は、東京湾周辺の画像を取得してみたいと思いますので、以下の通り緯度・経度を入力します。

# 緯度の最小値
min_lat = 35.675
# 経度の最小値
min_lon = 139.745
# 緯度の最大値
max_lat = 35.605
# 経度の最大値
max_lon = 139.835

衛星データをAPIで取得するためのURL、ヘッダ、パラメータを指定します。

# 取得する衛星画像のサーバURL、ヘッダ、パラメータを指定  
# ASNARO-1のURL
urls_ASNARO1 = 'https://gisapi.tellusxdp.com/api/v1/asnaro1/scene'
# ヘッダ
headers = {'Authorization': 'Bearer ' + API_TOKEN}
# パラメータ（緯度、経度の最大値、最小値）を指定
parameter = { 'min_lat': min_lat, 'min_lon': min_lon, 'max_lat': max_lat, 'max_lon': max_lon}

サーバURLとパラメータはTellusのAPI Referenceで確認可能です。
APIリファレンス : https://www.tellusxdp.com/ja/dev/api

衛星データのAPIには主に2種類あります。
1. 各衛星がどのようなシーンを持っているか
2. 各シーンごとの画像データ・波長別画像データ・SARデータ

まずは、ASNARO-1がどのようなシーンを持っているか確認してみましょう。

# requestsモジュールで指定したエリアを観測したシーンを検索し、リストをJSON形式で取得
r_ASNARO1 = requests.get(urls_ASNARO1, params=parameter, headers=headers).json()

# 中身を表示してみましょう
r_ASNARO1[0]
{'acquisitionDate': 'Wed, 31 Oct 2018 12:31:16 GMT',
 'clat': 35.6602740014593,
 'clon': 139.701407725281,
 'cloudCover': 20.0,
 'entityId': '20190115064426419_AS1',
 'max_lat': 35.7191765745438,
 'max_lon': 139.786374532312,
 'min_lat': 35.6013160074949,
 'min_lon': 139.61632039109,
 'path': 0,
 'productId': '20190115064426419_AS1_D01_L1B',
 'row': 0,
 'thumbs_url': 'https://tile.tellusxdp.com/thums/ASNARO-1/20190115064426419_AS1.PNG'}

この緯度経度の範囲に何枚の撮影データがあるか確認してみましょう

len(r_ASNARO1)
28

28枚あることがわかりました。
それではその中身を表示してみましょう。

# ASNARO-1の検索結果を表示
print("ASNARO-1のリスト [撮影日時, entityId, productId, cloudCover]")
for res in r_ASNARO1:
    print(", ".join([res['acquisitionDate'],
        str(res['entityId']), str(res['productId']), str(res['cloudCover'])]))

ASNARO-1のリスト [撮影日時, entityId, productId, cloudCover]
Wed, 31 Oct 2018 12:31:16 GMT, 20190115064426419_AS1, 20190115064426419_AS1_D01_L1B, 20.0
Sat, 20 Oct 2018 12:52:04 GMT, 20181224064142828_AS1, 20181224064142828_AS1_D01_L1B, 19.0
Tue, 31 Jul 2018 12:42:49 GMT, 20181224064003777_AS1, 20181224064003777_AS1_D01_L1B, 0.0
Sat, 30 Jun 2018 12:17:46 GMT, 20181224062516884_AS1, 20181224062516884_AS1_D01_L1B, 13.0
Fri, 08 Jun 2018 12:25:21 GMT, 20190115062640786_AS1, 20190115062640786_AS1_D01_L1B, 20.0
Sat, 02 Jun 2018 12:36:46 GMT, 20181224062257563_AS1, 20181224062257563_AS1_D01_L1B, 18.0
Sat, 21 Apr 2018 12:38:02 GMT, 20181225033937173_AS1, 20181225033937173_AS1_D01_L1B, 0.0
Sun, 25 Mar 2018 12:27:00 GMT, 20181224062059947_AS1, 20181224062059947_AS1_D01_L1B, 0.0
Wed, 03 Jan 2018 15:30:59 GMT, 20181224061906786_AS1, 20181224061906786_AS1_D01_L1B, 0.0
Tue, 02 Jan 2018 02:16:36 GMT, 20181224061858253_AS1, 20181224061858253_AS1_D01_L1B, 0.0
Fri, 29 Dec 2017 15:24:58 GMT, 20181224061846634_AS1, 20181224061846634_AS1_D01_L1B, 1.0
Sun, 05 Nov 2017 13:55:39 GMT, 20181224061816204_AS1, 20181224061816204_AS1_D01_L1B, 0.0
Sun, 09 Jul 2017 15:40:45 GMT, 20181224061651885_AS1, 20181224061651885_AS1_D01_L1B, 0.0
Mon, 22 May 2017 14:06:36 GMT, 20181224061523296_AS1, 20181224061523296_AS1_D01_L1B, 1.0
Sat, 06 May 2017 15:29:27 GMT, 20181224061506752_AS1, 20181224061506752_AS1_D01_L1B, 8.0
Tue, 04 Apr 2017 15:27:01 GMT, 20181224061413592_AS1, 20181224061413592_AS1_D01_L1B, 1.0
Fri, 24 Mar 2017 13:57:55 GMT, 20181224061303679_AS1, 20181224061303679_AS1_D01_L1B, 0.0
Sat, 18 Mar 2017 15:39:14 GMT, 20181224061254113_AS1, 20181224061254113_AS1_D01_L1B, 2.0
Mon, 24 Oct 2016 14:10:19 GMT, 20181224060959386_AS1, 20181224060959386_AS1_D01_L1B, 0.0
Fri, 17 Jun 2016 15:39:32 GMT, 20181224045149110_AS1, 20181224045149110_AS1_D01_L1B, 29.0
Thu, 05 May 2016 14:08:24 GMT, 20181223152759390_AS1, 20181223152759390_AS1_D01_L1B, 0.0
Fri, 18 Mar 2016 15:31:58 GMT, 20181223152702426_AS1, 20181223152702426_AS1_D01_L1B, 0.0
Mon, 30 Nov 2015 14:31:57 GMT, 20190117124820618_AS1, 20190117124820618_AS1_D01_L1B, 33.0
Tue, 03 Nov 2015 14:32:11 GMT, 20181225123811442_AS1, 20181225123811442_AS1_D01_L1B, 0.0
Sun, 25 Oct 2015 03:35:55 GMT, 20181223152305607_AS1, 20181223152305607_AS1_D01_L1B, 0.0
Thu, 06 Aug 2015 03:12:35 GMT, 20181223060118912_AS1, 20181223060118912_AS1_D01_L1B, 13.0
Wed, 05 Aug 2015 14:11:45 GMT, 20181223060109423_AS1, 20181223060109423_AS1_D01_L1B, 0.0
Sun, 28 Jun 2015 12:41:27 GMT, 20181225122707116_AS1, 20181225122707116_AS1_D01_L1B, 10.0

ここで、"cloud Cover"というのは被雲率と呼ばれるもので、画像の中にどのくらい雲が含まれるかを数値化したものです。今回はSARという雲に影響されない合成開口レーダーではなく、通常の光学衛星画像ですので被雲率の無いなかで最も新しい以下のデータを選択することにします。

Tue, 31 Jul 2018 12:42:49 GMT, 20181224064003777_AS1, 20181224064003777_AS1_D01_L1B, 0.0

次は、2つ目のAPIを使用してシーンごとの画像データを取得してみましょう。
ASNARO-1では画像の取得にproductIdを使用します。ここではentityIdも格納しておきましょう。

# 取得したいシーンを選択
# 選択したASNARO-1衛星画像の productId、entityIdを格納
ASNARO1_entityId = "20181224064003777_AS1"
ASNARO1_productId = "20181224064003777_AS1_D01_L1B"

ここで、画像タイルと拡大率の説明をします。(z = Zoom Level, x = tile_x, y = tile_y)と表現されますが、TellusではZ = 0 - 18の間で各Zoom Levelごとに256 x 256のタイル画像をAPIで取得することが可能です。これはWebメルカトル図法をもちいたもので、メルカトル図法で投影された世界地図を、東経180°から西経180°、北緯85.0511°から南緯85.0511°の範囲で切り取ったものです。

Z=1の場合

世界全体が1枚のタイルで表現されます。

Z=1の場合

タイルが4つに分割されます。

Z=2の場合

まだまだ日本全体が1枚のタイルで収まります。

Z=7の場合

だんだん東京湾が見えてきました。

Z=12の場合

今回の表示範囲であるZ=12, X=3638, Y=1613が見えてきました。

さて、今回表示しようと思った範囲をこのようにTellus OS上からタイルで見つけることも可能ですが、この場合は以下のように"データ選択"→"タイル画像"と選択します。

ただ、毎回このようなことを実施することが面倒な場合、緯度、経度をタイル座標にする手法をPython実装する方法を考えます。
詳しいアルゴリズムは長くなってしまうので差し控えますが、以下のコードで変換が可能です。

# 指定した緯度、経度を含むタイルの座標（x, y, z）を算出する。
# 緯度経度とズームレベル zを入力として対応するタイル座標を算出する関数
# 指定可能なズームレベル zについてはTellus 開発者向けページの「APIリファレンス」で確認できます。
# https://www.tellusxdp.com/ja/dev/api
def LatLon2Tile(lat,lon, z):
    L = 85.05112878
    tileX = int( (2**(z+7)) * ((lon/180) + 1) )
    tileY = int( (2**(z+7)/math.pi) * (-1 * np.arctanh(math.sin(lat * math.pi/180)) + np.arctanh(math.sin(L * math.pi/180))) )
    return int(tileX/256), int(tileY/256)

# 今回の指定エリアの中心座標(緯度、経度)
center_lat = (min_lat + max_lat) / 2
center_lon = (min_lon + max_lon) / 2
# 今回はズームレベル zを 12とします。
z = 12
# タイル座標を取得
x, y = LatLon2Tile(center_lat,center_lon, z)

この関数を用いて実際に衛星画像を取得してきましょう。
今回はRGBとNIR（近赤外）のバンド別に画像を取得していきます。

# ASNARO-1の衛星画像を取得
# サーバURLを指定
url_ASNARO1 =  'https://gisapi.tellusxdp.com/blend/asnaro1'
# バンドの割り当てを決定
query_ASNARO1 = "r=4&g=3&b=2" 
res_ASNARO1 = requests.get("{}/{}/{}/{}/{}.png?{}".format(
                url_ASNARO1,
                ASNARO1_productId,
                z, x, y,
                query_ASNARO1),
                headers={"Authorization": "Bearer " + API_TOKEN})
RGB_ASNARO1 = np.asarray(Image.open((BytesIO(res_ASNARO1.content))))
# データ型を整数に変更
RGB_ASNARO1 = RGB_ASNARO1.astype(np.uint8)

バンド割り当てでR=Band4, G=Band3, B=Band2として取り込みます。
取得した画像にはnumpy配列でRGBA(A=Alpha Chennel)という4次元で取り込まれます。

# 7. ASNARO-1の単バンド画像を取得(赤、青、緑、近赤外)
# 赤、青、緑の単バンド 2次元配列に変換
R_ASNARO1 = RGB_ASNARO1[:,:,0]
G_ASNARO1 = RGB_ASNARO1[:,:,1]
B_ASNARO1 = RGB_ASNARO1[:,:,2]

# ASNARO-1のバンド割り当てを変更（赤：バンド4、緑：バンド4、青：バンド4）
query_ASNARO1 = "r=6&g=6&b=6" 
# ASNARO-1衛星画像を取得
res_ASNARO1 = requests.get("{}/{}/{}/{}/{}.png?{}".format(
                url_ASNARO1, 
                ASNARO1_productId, 
                z, x, y, 
                query_ASNARO1),
                headers=headers)
# 近赤外単バンド画像をNumPy２次元配列
NIR_ASNARO1 = np.asarray(Image.open((BytesIO(res_ASNARO1.content))))
NIR_ASNARO1 = NIR_ASNARO1[:,:,0].astype(np.uint8)

次に近赤外NIR=Band6を取り込みます。

取り込んだ画像を人間の目に近い合成方式であるトゥルーカラー合成で表示します。
以下はトゥルーカラー合成のイメージです。

# 9. トゥルーカラー合成
# R：赤バンド、G：緑バンド、B：青バンドとして合成する
# ASNARO-1
true_ASNARO1 = np.dstack((R_ASNARO1, G_ASNARO1, B_ASNARO1))
# ASNARO-1衛星画像の表示
plt.figure(figsize=(20,20)) # 表示サイズをインチで指定 20×20インチ
plt.imshow(true_ASNARO1)
plt.show()

256 x 256ピクセルで取得されるため、この程度引き伸ばしてしまうとかなりぼやけしまいますが、浜離宮恩賜庭園や河川からの堆積物の流れ込み状況などある程度の情報は取得できそうです。

Google Earthではだめなの？

ここまで、なかなか面倒なCodeを書いてきたと思いますが、こんなことしなくてもGoogle EarthでGCP上からもっと便利に使えるのではないか？ということを考える方もいらっしゃると思います。私も実際そうでした。
実際何を目的とするかですので、Google Earthやその他のデータソースでも目的を満たせれば全く問題ありません。

しかしながら、衛星データには他のデータには存在しない特徴がいくつもあります。
独自調べですが、以下にそれぞれの特徴をまとめてみました。

	人工衛星	Google Earth	航空ドローン	航空機	Tellus
価格	✕✕✕	◎	△	✕	◎
過去データ	◎	△	✕	◯	◎¹
頻度	◎²	✕	△	✕	◎
解像度	◎	◎	◎	◎	◎
可視波長センサ（受動）	◎	◎	◎	◎	◎
その他のセンサ（受動）	◎	✕	✕	✕	◎
SAR（レーダー）	◎	✕	✕	✕	◎
API連携	✕	◎	✕	✕	◎
衛星以外のデータ	✕	✕	✕	✕	◎

衛星データの有利なところは上記にもある通り、過去データの豊富さ、様々な波長のセンサ、SARという合成開口レーダーによるアクティブセンシングという点が挙げられます。これらを用いることにより可視光では得られなかった情報を得ることができるようになります。

また、Tellusでは気温や人工流動データなど、衛星以外から得られるデータも提供していますので、そういった点でも有利です。

NDVI, NDWIってなんだろう

では可視光で得られなかった情報とはどういったものなのでしょうか？

https://sorabatake.jp/5192/

観測衛星は様々な波長の光学センサを搭載しているため、衛星ごとに取得できる波長が異なります。
上記のように対象物によって吸収・反射する電磁波の波長がことなるためこれらを組み合わせることにより様々な情報を得ることが出来ます。以下のように代表的な組み合わせが知られています。

https://sorabatake.jp/5192/

その中でもよく用いられる指標として植物がどれくらい生い茂っているかの指標であるNDVI、
水がどれくらい存在しているかという指標であるNDWIというものが存在します。
今回各画像のピクセルはuint8$(2^8=256)$で示されるため0−255までの数値で示され、
NDVI及びNDWIは(-1 → 1)の間で正規化して示されるようになります。

正規化植生指標NDVI（赤(R)、近赤外(NIR)の組み合わせ）

NDVI = \frac{R - NIR}{R + NIR}

正規化水指標NDWI（赤(R)、短波長赤外(SWIR)の組み合わせ）

NDWI = \frac{R - SWIR}{R + SWIR}

次に、ナチュラルカラー合成という合成方法を行います。
これは、R=Rの代わりにR=NIR（近赤外）としたカラー合成方法です。植物を緑で表示することが出来ます。
具体的イメージは以下のとおりです。

# 10. ナチュラルカラー合成
# R：赤バンド、G：近赤外バンド、B：緑バンドとして合成する
# ASNARO-1
natural_ASNARO1 = np.dstack((R_ASNARO1, NIR_ASNARO1, G_ASNARO1))

# ASNARO-1衛星画像の表示
plt.figure(figsize=(20,20)) # 表示サイズをインチで指定 20×20インチ
plt.imshow(natural_ASNARO1)
plt.show()

植生部分が緑色に可視化されていて、トゥルーカラー合成で見るよりも植生領域がわかりやすく表示できます。

次に、上記でご説明したNDVIを可視化してみましょう。
RとNIRを計算式に入れていきます。

# 1. NDVI
# データ型をuint8からfloat32に変更
NIR_ASNARO1 = NIR_ASNARO1.astype(np.float32)
R_ASNARO1 = R_ASNARO1.astype(np.float32)

# NDVIを算出
NDVI_ASNARO1 = (NIR_ASNARO1 - R_ASNARO1) / (NIR_ASNARO1 + R_ASNARO1)

plt.figure(figsize=(20,20)) # 表示サイズをインチで指定 6インチ × 6インチ
plt.imshow(NDVI_ASNARO1, cmap='Reds')
plt.colorbar()
plt.show()

NDVIを計算し、植生指数を数値化することが出来ました。残念ながらASNARO-1にはSWIR(短波長赤外)センサは搭載されておりませんため、NDWIを算出することは出来ません。しかしながら、こういった場合であっても、SWIRを搭載しているLandsat8のデータを用いることによりNDWIを算出することが可能です。ただし、Landsat8の場合はZ=12の拡大率までしか対応しておりませんので、この画像の粗さが限界となります。このように、Tellusでは目的により様々な衛星のデータを切り替えて使用することが可能です。

でもなんかやっぱり画像粗くない？

記事の最初の方でお見せした画像に比べて解像度が低いと思いますが、Webブラウザの読み込みが失敗しているわけではありません。各ズームレベルでは256x256の解像度での画像が取得されるため、例えばより高精細な画像を取得したい場合はズームレベルを上げる必要があります。ズームレベルは衛星ごとにも異なりますが、同一衛星内でもセンサにより異なります。ASNARO-1ではRGBを合成した光学センサでは最大ズームレベルが18となっておりますが、今回のようにバンド別の光学センサの場合は最大ズームレベルは16までとなります。各衛星ごとのSpec詳細は以下のページに詳しく書かれておりますので、必要に応じてご参照ください。
https://www.tellusxdp.com/ja/dev/data

さて、解像度の高い画像(Z=16)を取得しようとすると256x256ピクセルのため範囲が狭まってしまいます。しかしながら、範囲を広げよう(Z=12)とすると256x256ピクセルのため今度は解像度が下がってしまいます。このような場合はZ=16の画像を連結して表示することでZ=12の範囲でZ=16レベルの解像度で画像を取得することが可能です。

def get_combined_image(get_image_func, z, topleft_x, topleft_y, size_x=1, size_y=1, option={}, params={}):
    """
    結合されたタイル画像を取得する
    Parameters
    ----------
    get_image_func : function 
        タイル画像取得メソッド
        引数はz, x, y, option, params
    z : int 
        タイルのズーム率 
    topleft_x : int 
        左上のタイルのX座標 
    topleft_y : int 
        左上のタイルのY座標 
    size_x : int 
         タイルを経度方向につなぎわせる枚数 
    size_y : int 
        タイルを緯度方向につなぎわせる枚数
    option : dict 
        APIパスのオプション(z,x,y除く)
    params : dict 
        クエリパラメータ
    Returns
    -------
    combined_image: ndarray
        結合されたタイル画像
    """
    rows = []
    blank = np.zeros((256, 256, 4), dtype=np.uint8)

    for y in range(size_y):
        row = []
        for x in range(size_x):
            try:
                img = get_image_func(z, topleft_x + x, topleft_y + y, option)
            except Exception as e:
                img = blank

            row.append(img)
        rows.append(np.hstack(row))
    return  np.vstack(rows)

ここでは、まず空のnumpy行列を作成したあとに一つずつ画像を入れていくアルゴリズムをご紹介します。
Zが1増えるごとに画像は4倍、4倍となっていくため、Z=12からZ=16にすると、$4^{4}=256$タイルが並ぶことになります。これを図示すると下記のような感じです。

これを実際に実装したものが上記のCodeになります。

def get_asnaro1_image(z, x, y, option={}, params={}, query=''):
    """
    ASNARO-1のタイル画像を取得する
    Parameters
    ----------
    z : int 
        タイルのズーム率 
    x : int 
        タイルのX座標 
    y : int 
        タイルのY座標 
    option : dict 
        APIパスのオプション(productId)
    params : dict 
        クエリパラメータ
    Returns
    -------
    img: ndarray
        タイル画像
    """
    url = 'https://gisapi.tellusxdp.com/blend/asnaro1/{}/{}/{}/{}.png?{}'.format(option['productId'], z, x, y, query)
    headers = {
        'Authorization': 'Bearer' + API_TOKEN
    }
    r = requests.get(url, headers=headers, params=params)
    if not r.status_code == requests.codes.ok:
        r.raise_for_status()
    return io.imread(BytesIO(r.content))

つなぎ合わせる元の画像である256x256の画像を取得するアルゴリズムを今回は上記のように書いています。
下記では、トゥルーカラー合成のためにR=Band4, G=Band3, B=Band2として画像を取得していきます。

# ASNARO-1の衛星画像を取得
def rgb_asnaro1(z, x, y, option):
    # サーバURLを指定
    url_ASNARO1 =  'https://gisapi.tellusxdp.com/blend/asnaro1'
    # バンドの割り当てを決定
    query_ASNARO1 = "r=4&g=3&b=2" 
    res_ASNARO1 = requests.get("{}/{}/{}/{}/{}.png?{}".format(
                url_ASNARO1,
                option['productId'],
                z, x, y,
                query_ASNARO1),
                headers={"Authorization": "Bearer " + API_TOKEN})
    RGB_ASNARO1 = np.asarray(Image.open((BytesIO(res_ASNARO1.content))))
    # データ型を整数に変更
    #RGB_ASNARO1 = RGB_ASNARO1.astype(np.uint8)
    return RGB_ASNARO1.astype(np.uint8)

下記では、ナチュラルカラー合成のためにR=Band6, G=Band6, B=Band6として画像を取得していきます。

# ASNARO-1の衛星画像を取得
def nir_asnaro1(z, x, y, option):
    # サーバURLを指定
    url_ASNARO1 =  'https://gisapi.tellusxdp.com/blend/asnaro1'
    # バンドの割り当てを決定
    query_ASNARO1 = "r=6&g=6&b=6" 
    res_ASNARO1 = requests.get("{}/{}/{}/{}/{}.png?{}".format(
                url_ASNARO1,
                option['productId'],
                z, x, y,
                query_ASNARO1),
                headers={"Authorization": "Bearer " + API_TOKEN})
    RGB_ASNARO1 = np.asarray(Image.open((BytesIO(res_ASNARO1.content))))
    # データ型を整数に変更
    #RGB_ASNARO1 = RGB_ASNARO1.astype(np.uint8)
    return RGB_ASNARO1.astype(np.uint8)

下記では中心座標からTop Left XとTop Left Yの座標を算出しています。
+1しているのは中心座標から1/2を引き算してしまうと座標が戻りすぎてしまうために行っている補正です。

z = 16
size_x = 16
size_y = 16
x, y = LatLon2Tile(center_lat,center_lon, z)
topleft_x = int(x - size_x / 2 + 1)
topleft_y = int(y - size_y / 2 + 1)
option = {
    'productId': ASNARO1_productId
}
query_ASNARO1 = "r=4&g=3&b=2"

トゥルーカラー合成画像を実際につなぎ合わせていく関数に対して上記のパラメータを入れていきます。

combined = get_combined_image(rgb_asnaro1, z, topleft_x, topleft_y, size_x, size_y, option, params_ASNARO1)

ナチュラルカラー合成画像を実際につなぎ合わせていく関数に対しても上記パラメーを入れていきます。

combined_nir = get_combined_image(nir_asnaro1, z, topleft_x, topleft_y, size_x, size_y, option, query_ASNARO1)

各バンド別に画像を配列に格納していきます。

# ASNARO-1の単バンド画像を取得(赤、青、緑、近赤外)
# 赤、青、緑の単バンド 2次元配列に変換
COMR_ASNARO1 = combined[:,:,0].astype(np.uint8)
COMG_ASNARO1 = combined[:,:,1].astype(np.uint8)
COMB_ASNARO1 = combined[:,:,2].astype(np.uint8)

# 近赤外単バンド画像をNumPy２次元配列
COMNIR_ASNARO1 = combined_nir[:,:,0].astype(np.uint8)

まずはじめに、トゥルーカラー合成です。

# トゥルーカラー合成
# R：赤バンド、G：緑バンド、B：青バンドとして合成する
# ASNARO-1
true_COMASNARO1 = np.dstack((COMR_ASNARO1, COMG_ASNARO1, COMB_ASNARO1))
# ASNARO-1衛星画像の表示
plt.figure(figsize=(20,20)) # 表示サイズをインチで指定 20×20インチ
plt.imshow(true_COMASNARO1)
plt.show()

かなり高精細な画像が取得できましたね。
ピクセルでいうと256 x 16 = 4096ピクセルの画像です。

次に、ナチュラルカラー合成です。

# ナチュラルカラー合成
# R：赤バンド、G：近赤外バンド、B：緑バンドとして合成する
# ASNARO-1
natural_COMASNARO1 = np.dstack((COMR_ASNARO1, COMNIR_ASNARO1, COMG_ASNARO1))

# ASNARO-1衛星画像の表示
plt.figure(figsize=(20,20)) # 表示サイズをインチで指定 20×20インチ
plt.imshow(natural_COMASNARO1)
plt.show()

こちらのかなり高精細な画像が取得できました。

# NDVI
# データ型をuint8からfloat32に変更
COMNIR_ASNARO1 = COMNIR_ASNARO1.astype(np.float32)
COMR_ASNARO1 = COMR_ASNARO1.astype(np.float32)

# NDVIを算出
COMNDVI_ASNARO1 = (COMNIR_ASNARO1 - COMR_ASNARO1) / (COMNIR_ASNARO1 + COMR_ASNARO1)

plt.figure(figsize=(20,20)) # 表示サイズをインチで指定 6インチ × 6インチ
plt.imshow(COMNDVI_ASNARO1, cmap='Reds')
plt.colorbar()
plt.show()

NVDIについてもかなり高精細な画像が出来ています。

おまけ

本記事では主に衛星データの利用方法について記載していますが、画像といえば画像認識をやってみたくなリますよね。
今回の記事では触れませんが、もちろんTellus上で画像認識も可能です。
停泊中の船舶を機械学習モデル³を用いて検出した結果が以下になります。

確信度0.8以上のものを抜き出していますが、全体で326隻の停泊中船舶を発見することが出来ました。
陸地の方の川に停泊中の船も検出しています。
航行中の船舶は検出対象から除外していますので、波を引いている船舶は検出されていないことがわかります。

Reference

空畑
- Tellusの使用方法や宇宙ビジネスについてのニュースを取り扱っています
Tellus
- Tellus OSや開発環境へのアクセスはこちらから行います

免責事項

本サイトおよび対応するコメントにおいて表明される意見は、投稿者本人の個人的意見であり、シスコの意見ではありません。本サイトの内容は、情報の提供のみを目的として掲載されており、シスコや他の関係者による推奨や表明を目的としたものではありません。各利用者は、本 Web サイトへの掲載により、投稿、リンクその他の方法でアップロードした全ての情報の内容に対して全責任を負い、本 Web サイトの利用に関するあらゆる責任からシスコを免責することに同意したものとします。

今後の搭載衛星データによる ↩
衛星コンステレーションなど頻度の高い手法がある ↩
SSD512 Sliding Window(Stride 256, Window Size 512/512) ↩

Tellusによる初めての衛星データ解析

はじめに

この記事でできるようになること

最終的にこんなことができます

衛星データプラットフォーム Tellus とは？

公式ページから

特徴

Ciscoとの関連性

Tellus OS を利用してみよう

Tellusにユーザー登録をしてみよう

Tellus 開発環境を利用してみよう

1. Tellus OS

2. Tellus 開発環境

開発環境申請方法

API TOKENの登録

Jupyter Note Book の起動

実際に Code を書いてみよう

Z=1の場合

Z=1の場合

Z=2の場合

Z=7の場合

Z=12の場合

Google Earthではだめなの？

NDVI, NDWIってなんだろう

正規化植生指標NDVI（赤(R)、近赤外(NIR)の組み合わせ）

正規化水指標NDWI（赤(R)、短波長赤外(SWIR)の組み合わせ）

でもなんかやっぱり画像粗くない？

おまけ

Reference

免責事項