はじめに．

　人工衛星画像の無料環境なプラットフォームである「Tellus」を用いて，PALSAR-2とASNARO-1の人工衛星画像の解析を行ってきました．

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

次に，作成したプロダクト（生成物）をローカルでも評価できるようにしてみます．
QGISやArcGISなどのGISソフトにて表示・評価することが考えられますが，専用のGISソフトに慣れていない方も多いため，多くの方に馴染みのあるGoogl Earth(pro)による表示を試します．　
QGISについては，以下のサイトをご参考ください．

QGISによる衛星画像のカラー画像化．

kmlとkmz．

Tellusにて作成したプロダクト（今回は，jpeg画像）をGoogle Earthに表示するためには，その画像と一緒に地理情報が必要になります．
同じ画像ファイルでも，Geotiffであれば，画像のメタ情報として地理情報が含まれていますが，jpegなどの一般的な画像情報には地理情報は含まれていません．そこで，対象とする画像の地理情報のファイルを作成し，そのファイルの画像ファイルのメタ情報として利用します．
ここで利用する地理情報のメタファイルのフォーマットとして利用するのがkmlになります．

kmlは，Keyhole Markup Languageの頭文字で，地理情報の世界的な規格団体であるOGCが定めた規格になります．と，細かな話はWikipediaを参考にしてください．

たとえば，kmlファイルは以下のような情報があります．

botanicalgarden.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="17">
        <Placemark id="19">
            <name>Kirstenbosch</name>
            <Point id="18">
                <coordinates>18.432314,-33.988862,0.0</coordinates>
            </Point>
        </Placemark>
    </Document>
</kml>

このkmlファイルは，地図上のあるポイントを示しています．上記のファイルをメモ帳などの適当なテキストソフトに貼り付けて，kmlの拡張子で保存してください．その後，Googl Earthがインストールされていれば，そのファイルをダブルクリックすると，Googl Earthのある位置に移動するのががわかります．

このkmlファイルに対象とする画像の情報を埋め込み，このkmlと画像を圧縮して一つにしたファイルがkmzになります．

Google Earthに表示する．

みなさんご存知かと思いますので，Google Earthの説明はここでは省きます．以下の，Google Earthのサイトをご参考にしてください．

Google Earth

私は自分がいつも使っているMacbook AirにGoogle Earth Proをインストールし使っています．

pythonにてkmlファイルを作成するためにsimplekmlというモジュールを利用しました．こういうもの作ってオープンにしてくることに感謝します．大げさですが，こういった活動は人類の財産になっていると考えています．
利用に関するドキュメントは以下にありますので，辞書代わりに使ってください．

simplekml Documentation Release 1.3.1

このドキュメントを参考に，これまでのTellusの総合開発環境にて作成した茨城県つくば市の衛星画像をローカルのGoogle Earthで利用できるようにします．

import simplekml

まず，こちらを実行しsimplekmlのモジュールをimportします．simplekmlのモジュールがない場合は，jupyter labのotherからterminalを実行し，simplekmlをインストールしてください．

今回は，Tellusにて作成したASNARO-1の茨城県つくば市の画像を利用します．
この画像を作成するときに用いた，画像の座標情報は以下になります．（座標の求め方は，人工衛星画像の無料環境なプラットフォーム「Tellus」にてデータ解析． -ASNARO-1編-をご参考ください．）

lat= {"lat":[102863,102881]}
lon = {"lon":[233091,233111]}
zoom = {"zoom":18}

この地理院座標から緯度（latitude)と経度(longtitude)を求めます．ここでは，以下の関数を用いて求めました．

#タイルの北西端の緯度・経度を算出する．
# -*- coding:utf-8 -*-
from math import pi, e, atan

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

例えば，複数の画像ファイルを合成した画像の北西端の緯度・経度は以下を実行し求めます．

tilelatlon(lon["lon"][0], lat["lat"][0], zoom["zoom"])

この結果が以下となります．

(36.08573111062387, 140.10177612304688)

この地理院座標から緯度・経度を求める関数を用いて，合成画像の四隅の緯度・経度を求めます．

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='Asnaro1_Tsukuba_18')
ground.icon.href = 'Asnaro1_Tsukuba_18.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("tsukuba.kml") #kmlのファイル名(任意)

こちらが作成したkmlファイルになります．

tsukuba.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="81">
        <GroundOverlay id="82">
            <name>Asnaro1_Tsukuba_18</name>
            <Icon id="83">
                <href>Asnaro1_Tsukuba_18.jpg</href>
            </Icon>
            <LatLonBox>
                <north>36.08573111062387</north>
                <south>36.064641955171425</south>
                <east>140.130615234375</east>
                <west>140.10177612304688</west>
                <rotation>0</rotation>
            </LatLonBox>
        </GroundOverlay>
    </Document>
</kml>

このkmlファイルと対象の画像ファイルをkmzのファイルとして圧縮します．

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

これで，Asnaro1_Tsukuba_18.kmzという圧縮ファイルが作成されます．　こちらのファイルをローカルにダウンロードし実行してください．　Google Earthが立ち上がり，対象の画像がGoogle Earthに貼り付けられます．　以下が表示例です．

衛星画像の上に衛星画像が貼り付けられているので，これでは分かり難いですね．
近づけて見てみます．

左図：Tellusにて作成したASNARO-1の衛星画像　右図：Google Earthのベースマップ

Google Earthの左ウィンドウのチェック欄から対象とする画像のチェックを外すと，Google Earthのベースマップが表示されます．衛星画像はある角度から撮像されているため，左のASNARO-1の画像はビルの影が大きく伸びているのがわかります．

おわりに．

　Tellusにて作成したプロダクトがTellus上で使えないのはもったいなく，　また特別なGISソフトがなくても表示できたほうが便利だと考え，Google Earthでの表示方法を作りました．
　今回は，ASNARO-1の画像になりますが，それぞれの衛星から求めた価値の高い情報を汎用的な地図アプリケーションに表示させることができれば，利用の拡大につながると考えています．
　この情報がみなさんのお役に立てば幸いです．コメントやフィードバックがいただければ嬉しいです．励みになります．

最後に，今回作成したコードをまとめます．

import simplekml
import zipfile

from math import pi, e, atan

#タイルの北西端の緯度・経度を算出する．
# -*- coding:utf-8 -*-
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

#ASNARO-1の衛星画像より求めた関心域（茨城県つくば市）の地理院座標の情報
lat= {"lat":[102863,102881]}
lon = {"lon":[233091,233111]}
zoom = {"zoom":18}

#合成画像の右上・左上・右下・左下の四隅の地理情報
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='Asnaro1_Tsukuba_18')
ground.icon.href = 'Asnaro1_Tsukuba_18.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("tsukuba.kml")

#kmzファイルの作成
with zipfile.ZipFile('Asnaro1_Tsukuba_18.kmz', "w", zipfile.ZIP_DEFLATED) as zf:
    zf.write('tsukuba.kml', arcname='tsukuba.kml')
    zf.write('Asnaro1_Tsukuba_18.jpg', arcname='Asnaro1_Tsukuba_18.jpg')