はじめに

今回は、Tellusにて取得した衛星画像を自分の端末にダウンロードして利用する例をご紹介します。

衛星画像はPALSAR2を利用し、衛星画像の取得、複数画像の組み合わせ処理、および異なる時期の画像取得方法、ローカルにてGoogle EarthなどのGISソフトで閲覧するためのフォーマット処理のフローをまとめてみました。

例えば、今回のプログラムを用いることで、以下のスクリーンのようにGoogle EarthにてTellusで生成した画像を表示・利用することができるようになります。

解析環境を構築

Tellusにて取得した画像を、Google Earthにて閲覧するためには、取得した画像に位置情報をセットとしたkmz(kml)ファイルを作成する必要があります。
kml(Keyhole Markup Language)ファイルとは、Tellusで利用できるGeoJSONファイル同様に、地理空間情報を描画する方法の一種です。kmlをzipで圧縮したファイル形式は、kmzファイルと呼ばれています。

kmlファイルを作成する際に、今回はsimplekmlを利用しました。
このライブラリにはkmlファイルを作成するための関数が組み込まれており、緯度・経度情報をインプットすることで、簡単にkmlファイルを作成することができます。
simplekmlについては、以下の公式サイトやドキュメントをご参考ください。
https://pypi.org/project/simplekml/
SIMPLEKML 1.3.1 documentation
https://simplekml.readthedocs.io/en/latest/
まずは、Tellusの総合開発環境にsimplekmlをインストールします。

#kml情報取得モジュールのsimplekmlのインストール
!pip install simplekml

次に今回用いるモジュールを宣言します。

# -*- coding:utf-8 -*-
from math import pi, e, atan

import requests
import io
import numpy as np
from PIL import Image, ImageDraw, ImageFilter,ImageOps
import sys

import json
import math
from skimage import io
from io import BytesIO
%matplotlib inline

import simplekml
import zipfile

TOKEN = '自身のトークンを貼り付けてください'
DOMAIN = 'tellusxdp.com'

衛星画像の取得および組み合わせ処理の準備

以上で開発環境で衛星画像を取得し、データを処理する準備ができました。
これから、具体的に対象とする衛星画像（PARSAR-2）を取得します。

まず、PALSAR-2の観測画像を取得するための関数を準備します。
Tellusでは、地図に重畳できる、補正（幾何補正）処理をおこなった画像(L2.1)が保存されています。
TellusのPALSAR-2の観測画像の仕様はこちらをご覧ください。
https://www.tellusxdp.com/ja/dev/data/palsar-2

TellusのAPI情報、およびコードの例は以下のサイトを参考にしてください。
Tellus APIリファレンス：
https://www.tellusxdp.com/ja/dev/api

関心域の画像の有無を、以下の関数を用いて確認します。

#PALSAR2のSAR観測画像の関心領域のメタ情報を取得する。
def getPalsar2Scene(token,area):
    url = 'https://gisapi.tellusxdp.com/api/v1/palsar2/scene'
    payload = area
    headers = { 'Authorization': 'Bearer %s' % token }
    r = requests.get(url, params=payload, headers=headers)
    #print(url)
    if r.status_code is not 200:
        raise ValueError('status error({}).'.format(r.status_code))
    return r.json()

次に、関心域の中で取得したい画像の観測日を確認し、その画像を以下の関数にて取得します。
このとき画像のタイル座標はTellus OSより取得できます。
タイル座標の検索方法や、タイル画像の組み合わせ方については「【ゼロからのTellusの使い方】大きなタイル画像を取得するには？地図タイルを扱う上で知っておきたいTips4選」をご覧ください。

#PALSAR2の画像取得の関数。
def get_Palsar2_scene_tile0(scene_id, zoom, xtile, ytile):
    url = 'https://gisapi.{}/PALSAR-2/{}/{}/{}/{}.png'.format(
        DOMAIN, scene_id, zoom, xtile, ytile)
    headers = {
        'Authorization': 'Bearer ' + TOKEN
    }

    r = requests.get(url, headers=headers)
    #print(url)
    if r.status_code is not 200:
        raise ValueError('status error({}).'.format(r.status_code))
    return io.imread(BytesIO(r.content))

Tellusの画像仕様を確認すると、PALSAR-2の画像は最大ズーム（高分解能画像）は１６とあります。
このズームレベルの１枚の画像の範囲は、おおよそ500m×500m。

そのため、関心域が１枚の画像サイズ以上であれば、複数の画像を取得し、それを縦横にならべることで大きな画像を作成します。
そのための関数が以下となります。

#PALSAR2の最新の画像のタイル処理（画像番号の選択）
def getPalsar2Img_res(i,zoom,min_v,max_v,min_h,max_h):

    img = None

    for item_v in range(min_v,max_v+1):
        img_h = None
        for item_h in range(min_h,max_h+1):
            if(img_h is None):
                img_h = get_Palsar2_scene_tile0(res[i]['entityId'], zoom, item_h, item_v)
            else:
                img_h=np.concatenate((img_h, get_Palsar2_scene_tile0(res[i]['entityId'], zoom, item_h, item_v)), axis=1)
        if (img is None):
            img = img_h
        else:
            img = np.concatenate((img, img_h), axis=0)

    io.imshow(img)
    print(res[i]['acquisitionDate'])
    print(img.shape)
    return img

ここでは、画像の観測日（getPalsar2Sceneで取得した画像情報の取得番号は、0が最新画像となります）、左上、左下、右上、右上のタイル番号を入力し、組み合わせた大きな画像を作成します。

最後に、地理情報のkmlファイルを作成し、取得した画像ファイルとを組み合わせたkmzファイルを作成します。今までに準備した関数を用いて、取得したい画像の緯度・経度情報を入力し、kmlファイルを作成します。

#画像ファイルの保存およびkmzファイルの作成（PALSAR2)
def PALSAR2_kmz(day_number, image_name, ):
    A = (res[day_number]['acquisitionDate'][5:7]+res[day_number]['acquisitionDate'][8:11]+res[day_number]['acquisitionDate'][12:16]) 
    file_name = 'PALSAR2_' + A + '_res' + str(day_number)+ '_' + str(zoom["zoom"])
    #取得した画像の保存。
    Image.fromarray(image_name).save(file_name + '.png')
    path_png = file_name + '.png'
    rgb_im = Image.open(path_png).convert('RGB')
    rgb_im.save(file_name +'.jpg')
    
    #画像の地理情報（緯度・経度）の取得。
    left_top = tilelatlon(lon["lon"][0], lat["lat"][0], zoom["zoom"])
    right_top = tilelatlon(lon["lon"][1]+1, lat["lat"][0], zoom["zoom"])
    left_down =  tilelatlon(lon["lon"][0], lat["lat"][1]+1, zoom["zoom"])
    right_down = tilelatlon(lon["lon"][1]+1, lat["lat"][1]+1, zoom["zoom"])
    
    #地理情報をkmlで出力。
    kml = simplekml.Kml()
    ground = kml.newgroundoverlay(name= file_name )
    ground.icon.href = file_name +'.jpg'
    ground.latlonbox.north = left_top[0]
    ground.latlonbox.south = left_down[0]
    ground.latlonbox.east = right_top[1]
    ground.latlonbox.west = left_top[1]
    ground.latlonbox.rotation = 0
    kml.save(file_name +'.kml')
    
    #地理情報と画像とをあわせたkmzを出力。
    with zipfile.ZipFile(file_name +'.kmz', "w", zipfile.ZIP_DEFLATED) as zf:
        zf.write(file_name +'.kml', arcname=file_name +'.kml')
        zf.write(file_name +'.jpg', arcname=file_name +'.jpg')

    print(file_name +'.kmz file is completed')
    return

以上で、 TellusにてPALSAR-2およびASNARO-1の観測画像の取得、および取得画像の地理情報との組み合わせた位置・画像ファイルの作成準備ができました。

新国立競技場のPALSAR-2の画像を取得しよう

では、準備した関数を用いてTellusにてPALSAR-2の画像を取得します。
今回の対象である新国立競技場のタイル座標は以下となります。

関心域が変わるたびに、このタイル座標をTellus OSより取得することが手間な場合もあります。
そこで、任意の領域のタイル座標を所得するアプリケーションもHerokuにて公開されているサイトがありますので、こちらを参考にタイル座標を取得してみましょう。

詳細な使い方に関する解説は以下の記事をご参考にしてください。
Tellusで超低高度衛星つばめの衛星画像タイルを取得して連結してみた

#関心領域の国土地理院タイル座標の取得。
from IPython.display import HTML
HTML(r'')

ここでは、画面にてズームレベル(ALOS-2では最大16)を選択し、その後関心域をマウスで指定します。
そうすると、対象域の緯度・経度情報と、この範囲の地理院タイル座標のxおよびyの範囲が出力されます。

この情報を下記にコピペします。

#上記にて取得した地理情報をコピペ。
AREA = {"min_lat":35.672,"min_lon":139.708,"max_lat":35.684,"max_lon":139.728}
lat= {"lat":[25806,25808]}
lon = {"lon":[58201,58204]}
zoom = {"zoom":16}

次に、この範囲のPALSASR-2の画像リストを以下を実行し取得します。

#関心領域のPALSAR2のSAR観測画像のメタ情報の取得。
res=getPalsar2Scene(TOKEN,AREA)
print('number of scenes : {}'.format(len(res)))
for i in range(len(res)):
    print(res[i]['acquisitionDate'])

この結果から、同範囲で取得されている画像は、2018年6月28日、10月18日、11月29日の3候補あることがわかります。

今回は、最新な2018年11月29日の同地域を含むPALSAR-2の観測画像をみてみます。

#サムネイル画像を表示
thumbs = io.imread(res[0]['thumbs_url'])
io.imshow(thumbs)

PALSAR-2の1回（シーン）の観測範囲は観測モードによって異なりますが、Tellusに提供されている画像は50km×50kmの画像になります。
最新の観測画像のシーン情報を取得するには下記のコマンドを実行します。

#最新の衛星画像のメタ情報を取得。
res[0]

このシーン情報では、観測日だけでなく、観測画像の中心の緯度・経度や、プロダクト番号などの情報を得ることができます。

では、新国立競技場を撮影したALOS-2の最新観測画像を取得します。
この関数の”0”は、先程確認した観測画像履歴の3シーンの最新の画像を選択する、ということを意味します。

#関心領域の複数画像の合成画像を取得。
img_PAL_tokyo = getPalsar2Img_res(0,zoom["zoom"],lat["lat"][0], lat["lat"][1], lon["lon"][0], lon["lon"][1])

取得した新国立競技場のSAR画像では、競技場の形がよくわかりますね。

では、この画像のkmzファイルを作成します。

#PALSAR2 kmzファイル作成（res_number,  file_name)
PALSAR2_kmz(0, img_PAL_tokyo)

上記のファイル名のkmzファイルが作成されました。

次に、過去の新国立競技場のSAR画像を取得しました。
今回は、3番目のシーンである”Thu, 28 Jun 2018 02:42:58 GMT”の観測画像を取得します。

#関心領域の複数画像の合成画像を取得(シーン選択)。
img_PAL_tokyo_res2 = getPalsar2Img_res(2,zoom["zoom"],lat["lat"][0], lat["lat"][1], lon["lon"][0], lon["lon"][1])

2018年11月29日の観測画像と比べると、6月28日の画像は、新国立競技場の画像左部分が白い、つまり反射している信号が強いことがわかります。

この画像も同じく、位置情報を組み合わせたkmzファイルを作成します。

#PALSAR2 kmzファイル作成（res_number,  file_name)
PALSAR2_kmz(2, img_PAL_tokyo_res2)

ここで作成したkmzファイルを自分のPC（ローカル）にダウロードします。

Image("kmz_download.png")

作成したファイルをマウスで選択し右クリックするとDownloadができます。
ここでダウンロードしたkmzファイルをダブルクリックすると、Google EarthをインストールしているPCであれば、
自動的にGoogle Earthが起動し、ダウンロードしたALOS-2の画像が地図に貼り付けられます。
Google Earthをダウンロードしていない場合でも、以下のようにブラウザにて「プロジェクト」よりファイルをアップロードし、描画することもできます。

おわりに

以上、Tellusに保存されているALOS-2の任意の範囲の画像の取得方法、およびそれを自分のPC（ローカル）にて地理ソフト（今回はGoogle Earth)にて閲覧するためのkmzファイルの作成方法を紹介しました。

今回ご紹介したコードを用いて、観測日の異なるALOS-2の観測画像を比較してみましょう。SAR画像では色合いが変化した場合に、何がどのように変化したのか、理解するには難しい場合があります。そのため、光学画像を用いて比較することも有効となります。
次回の記事では、同日付近のASNARO-1の観測画像で確認し、同様にGoogle Earthに重畳する方法をご紹介します。

なお、Tellusの画像は、データポリシーによりダウンロード可能な条件が異なります。他の画像で同様の処理を行う場合には、データポリシーをご確認の上、ご利用ください。
馴染みのあるソフトで、自宅や馴染みのある場所の人工衛星の観測画像をみてみると、多くの気付きがあるかと思います。ぜひ試してみてください。