人工衛星画像の無料環境なプラットフォーム「Tellus」にてデータ解析． -PALSAR-2編-

はじめに．

何度かに分けて「Tellus」に提供されている衛星画像の解析処理方法にいて紹介してきました．
ここで，PALSAR-2の衛星画像を例として，衛星画像の取得，複数画像の組み合わせ処理，Google Earthでの表示するためのフォーマット処理のフローをまとめてみました．これにより，以下のスクリーンのようにGoogle EarthにてTellusで生成した画像を表示・利用することができるようになります．
(参考ページの情報を重複しているところが多々ありますが，全体のフローがわかるのが有用かと思いまとめてみました．)

(追加更新 2019/09/29)

作成するプロダクトの例として，ここで作成したPALSAR-2の東京の観測データを共有します．Google EarthをインストールしたPCにダウンロードし，ダブルクリックで開いてください．

PALSAR-2の観測画像の解析例：東京
（注:PALSAR-2のデータポリシーにより解像度を修正しています．)

参考ページ

人工衛星画像のオープンフリーなプラットフォーム「Tellus」にてデータ解析．
人工衛星画像の無料環境なプラットフォーム「Tellus」にてデータ解析． -ASNARO-1編-
人工衛星画像の無料環境なプラットフォーム「Tellus」にてデータ解析． -Google Earth編-

PALSAR-2とは．

PALSAR-2はJAXAが開発・運用している「だいち２号」（ALOS-2)に搭載している電波による地表面の状態を撮影するセンサになります．（説明が難しい．．．地表を電波のメガネでみている，という表現がわかりやすいのような，わかり難いような．）

陸域観測技術衛星2号「だいち2号」（ALOS-2）（クレジット：宇宙航空研究開発機構）

PALSAR-2の特徴は，地上に対して斜め方向に電波を放射し，地表から反射してくる電波を受け取って，その強弱を画像として表示します．

参照：合成開口レーダ（SAR）のキホン～事例、分かること、センサ、衛星、波長～
何度も書きますが，空畑の記事はわかりやすいです．
一言でいうと，川や池，道路などの構造物がないもの反対方向に反射するため信号が弱いため黒くなり，建物や森などの構造物は衛星方向にも電波が反射されるため信号が強く白く写ります．
PALSAR-2の画像の一度の撮影範囲やその空間分解能は以下となります．

[ALOS-2プロジェクト / PALSAR-2]
(https://www.eorc.jaxa.jp/ALOS-2/about/jpalsar2.htm)　（クレジット：宇宙航空研究開発機構）

いろいろな観測モードがあって分かり難いですが，Tellusに提供されているのは，この中の"高分解能 3m"になります．一度の撮影領域は５０km四方になります．（2019/9/23時点）

関心のある領域のPALSAR-2の観測画像をTellusにてチェックし，その画像をGoogle Earthで表示できるよう処理していきます．

PALSAR-2の撮影画像の解析．

Tellusへの登録，および総合開発環境の立ち上げについては，こちらのページをご参考にしてください．
人工衛星画像のオープンフリーなプラットフォーム「Tellus」にてデータ解析．

今回の解析に必要なモジュールを実行します．

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

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

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

import simplekml
import zipfile

TOKEN = '********-****-****-****-************'
DOMAIN = 'tellusxdp.com'

このあとGoogle Earthに解析したプロダクトを表示するための地理情報を処理するための「simplekeml」というモジュールが必要です．インストールされていない方は， pip installなど，それぞれの環境に応じてインストールしてください．
また，TOKENは，それぞれにTellusのアカウント取得時に提供される個々の情報を設定してください．

次に，TellusでのPALSAR-２の画像取得，取得した複数の画像の組み合わせ処理，画像の地理情報の算出するための各関数を実行します．

# 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()

＃PALSAR2のSAR画像を取得する．
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))

# PALSAR2の複数画像の組み合わせ処理
def getPalsar2Img0(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[0]['entityId'], zoom, item_h, item_v)
            else:
                img_h=np.concatenate((img_h, get_Palsar2_scene_tile0(res[0]['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[0]['acquisitionDate'])
    print(img.shape)
    return img

# タイルの北西端の緯度・経度を算出する．
def tilelatlon(x, y, z):
    lon = (x / 2.0**z) * 360 - 180 # 経度（東経）
    mapy = (y / 2.0**z) * 2 * pi - pi
    lat = 2 * atan(e ** (- mapy)) * 180 / pi - 90 # 緯度（北緯）
    #print (lat,lon)
    return lat,lon

次に，以下を実行し　関心領域の緯度・経度情報や，その地域の地理院タイルの座標情報を取得します．

# 関心領域の国土地理院タイル座標の取得．
from IPython.display import HTML
HTML(r'<iframe width="960" height="480" src="https://tellus-map-test.herokuapp.com/" frameborder="0"></iframe>')

Tellusが提供している衛星画像は，この地理院座標のタイル座標で分割されています．
このページで画像の座標情報をメモし，それをTellusのAPIで実行し，，，というのは手間です．　そのため，上記のアプリケーションを実行し，これらの情報を取得・コピペしています．このアプリケーションは以下をご参考ください．（便利なアプリケーションを提供いただきありがとうございます．）
参考情報：Tellusで超低高度衛星つばめの衛星画像タイルを取得して連結してみた

こちらを実行し，関心領域の地理院座標を取得します．使い方のコツは，関心領域まで地図を拡大した後，範囲指定の横のZoomを希望の数字にしてから，”範囲指定”をクリックし範囲をマウスで示します．PALSAR-2の衛星画像のZoomは９から１６になります．できるだけ高分解能の画像を扱いたいときは１６にします．このとき，選択する範囲が広いと多くの画像を扱うことになりますので，適切な値になるよう気をつけてください．

取得した地理情報を以下にコピペします．

# 関心領域の国土地理院タイル座標をコピペ．
AREA = {"min_lat":35.642,"min_lon":139.675,"max_lat":35.777,"max_lon":139.837}
lat= {"lat":[25785,25815]}
lon = {"lon":[58195,58224]}
zoom = {"zoom":16}

次に，該当する領域のTellusのPALSAR-2の画像情報を取得します．

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

こちらを実行したアウトプットが以下となります．

number of scenes : 7
Thu, 29 Nov 2018 02:43:05 GMT
Thu, 29 Nov 2018 02:42:57 GMT
Thu, 18 Oct 2018 02:43:04 GMT
Thu, 18 Oct 2018 02:42:56 GMT
Thu, 06 Sep 2018 02:42:56 GMT
Thu, 28 Jun 2018 02:43:03 GMT
Thu, 28 Jun 2018 02:42:55 GMT

関心領域のPALSAR-2の衛星画像は７シーンあることがわかりました．最新の画像は２０１８年１１月２９日のものになります．
この最新画像の全体画像（サムネイル）を見てみます．

thumbs = io.imread(res[0]['thumbs_url'])
io.imshow(thumbs)

実行結果が以下となります．

東京と神奈川の５０km四方の撮影画像であることがわかりました．
それでは，この画像のより詳細なメタ情報を取得します．

res[0]

実行結果が以下となります．

{'acquisitionDate': 'Thu, 29 Nov 2018 02:43:05 GMT',
 'censor': 'HH',
 'clat': 35.617,
 'clon': 139.659,
 'cloudCover': 'N/A',
 'entityId': 'ALOS2243962900',
 'max_lat': 35.963,
 'max_lon': 140.034,
 'min_lat': 35.268,
 'min_lon': 139.288,
 'path': 0,
 'productId': 'SARD000000232681-00007-066-011',
 'row': 0,
 'sceneId': 'ALOS2243962900-181129',
 'thumbs_url': 'https://tile.tellusxdp.com/thums/PALSAR2/ALOS2243962900.PNG'}

先程は撮影日の情報だけでしたが，それ以外の情報も取得することができました．
ここで重要なのは，'entityId'になります．こちらから，希望するタイル座標の画像を取得します．

img_PAL_tokyo = getPalsar2Img0(zoom["zoom"],lat["lat"][0], lat["lat"][1], lon["lon"][0], lon["lon"][1])]

こちらを実行すると，以下の画像を合成・取得することができました．

次に，次のコードを実行し画像を一時保存した後，作成したプロダクトをローカルで扱うことができるよう加工（二次生成物）であるjpeg画像に変換し保存します．

Image.fromarray(img_PAL_tokyo).save('Palsar2_Tokyo_16.png')
path_png = 'Palsar2_Tokyo_16.png'
rgb_im = Image.open(path_png).convert('RGB')
rgb_im.save('Palsar2_Tokyo_16.jpg')

以上で，Tellusにて関心領域のPALSAR−２の観測画像を取得することができました．
次に，この画像をGoogle Earthにて表示できるよう地理情報を取得します．

Google Earthで表示．

Google EarthにTellusにて生成した衛星画像を表示させるには，この画像の地理情報が必要となります．
ここでは，地理情報として，同地域のkmlフォーマットのファイルを生成します．ｋｍｌなどの地理情報については，以下をご参考ください．

人工衛星画像の無料環境なプラットフォーム「Tellus」にてデータ解析． -Google Earth編-

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='Palsar2_Tokyo_16')
ground.icon.href = 'Palsar2_Tokyo_16.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("Palsar2_Tokyo_16.kml")

作成したkmlファイルは以下になります．

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2">
    <Document id="6">
        <GroundOverlay id="7">
            <name>Palsar2_Tokyo_16</name>
            <Icon id="8">
                <href>Palsar2_Tokyo_16.jpg</href>
            </Icon>
            <LatLonBox>
                <north>35.77771427205079</north>
                <south>35.63944106897392</south>
                <east>139.8394775390625</east>
                <west>139.6746826171875</west>
                <rotation>0</rotation>
            </LatLonBox>
        </GroundOverlay>
    </Document>
</kml>

このkmlファイルと，先に作成したPALSAR-2の衛星画像を一つのファイルに格納・圧縮します．

with zipfile.ZipFile('Palsar2_Tokyo_16.kmz', "w", zipfile.ZIP_DEFLATED) as zf:
    zf.write('Palsar2_Tokyo_16.kml', arcname='Palsar2_Tokyo_16.kml')
    zf.write('Palsar2_Tokyo_16.jpg', arcname='Palsar2_Tokyo_16.jpg')

こちらを実行すると，左のウィンドウに"Palsar2_Tokyo_16.kmz"というファイルが作成されたのがわかります．このファイルを右クリックでDownloadし，Google EarthがインストールされたPCにてダブルクリックして起動してみてください．このページのトップの画像のように，PALSAR-2の撮影画像がGoogle Earthに表示されます．

PALSAR-2の画像のときに，地表に構造物が少ないときは暗く（黒く），構造物があるところは明るく（白く）写ると説明しましたが，同じ建物が密集しているところでも，明るいところと暗いところがあるのがわかるかと思います．なぜ，このような明暗ができるのか，PALSAR-2の撮影方法からその特徴を考えてみるなど，いろいろ推察してみるのは面白いかと思います．

おわりに．

Tellusにて提供されているPALSAR-2の撮影画像を参考に，データの取得・変換方法と，ローカルでGoogle Earthに表示するための処理について紹介しました．

PALSAR-2などのSARセンサの撮影画像は，ASNARO-1のような光学センサの画像とことなり，白黒の画像となります．これらは画像そのものを使うのではなく，画像を数値処理したり，他の情報と組み合わせることで，新しい情報に変換し利用されています．
利用例として，空畑で紹介されていましたので，ご参考ください．

海外に学ぶ、ビジネス効率を加速させる衛星データ活用術！@空畑

全体コード．

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

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

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

import simplekml
import zipfile

TOKEN = '********-****-****-****-************'
DOMAIN = 'tellusxdp.com'

# 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()

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))

# PALSAR2の最新の画像のタイル処理
def getPalsar2Img0(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[0]['entityId'], zoom, item_h, item_v)
            else:
                img_h=np.concatenate((img_h, get_Palsar2_scene_tile0(res[0]['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[0]['acquisitionDate'])
    print(img.shape)
    return img

# タイルの北西端の緯度・経度を算出する．
def tilelatlon(x, y, z):
    lon = (x / 2.0**z) * 360 - 180 # 経度（東経）
    mapy = (y / 2.0**z) * 2 * pi - pi
    lat = 2 * atan(e ** (- mapy)) * 180 / pi - 90 # 緯度（北緯）
    #print (lat,lon)
    return lat,lon

# 関心領域の国土地理院タイル座標の取得．
from IPython.display import HTML
HTML(r'<iframe width="960" height="480" src="https://tellus-map-test.herokuapp.com/" frameborder="0"></iframe>')

# 上記にて取得した地理情報をコピペ．
AREA = {"min_lat":35.642,"min_lon":139.675,"max_lat":35.777,"max_lon":139.837} 
lat= {"lat":[25785,25815]}
lon = {"lon":[58195,58224]}
zoom = {"zoom":16}

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

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

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

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

# 取得した画像の保存．
Image.fromarray(img_PAL_tokyo).save('Palsar2_Tokyo_16.png')
path_png = 'Palsar2_Tokyo_16.png'
rgb_im = Image.open(path_png).convert('RGB')
rgb_im.save('Palsar2_Tokyo_16.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='Palsar2_Tokyo_16')
ground.icon.href = 'Palsar2_Tokyo_16.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("Palsar2_Tokyo_16.kml")

# 地理情報と画像とをあわせたkmzを出力．
with zipfile.ZipFile('Palsar2_Tokyo_16.kmz', "w", zipfile.ZIP_DEFLATED) as zf:
    zf.write('Palsar2_Tokyo_16.kml', arcname='Palsar2_Tokyo_16.kml')
    zf.write('Palsar2_Tokyo_16.jpg', arcname='Palsar2_Tokyo_16.jpg')