記事作成時から、Tellusからデータを検索・取得するAPIが変更になっております。該当箇所のコードについては、以下のリンクをご参照ください。
https://www.tellusxdp.com/ja/howtouse/access/traveler_api_20220310_
firstpart.html

2022年8月31日以降、Tellus OSでのデータの閲覧方法など使い方が一部変更になっております。新しいTellus OSの基本操作は以下のリンクをご参照ください。
https://www.tellusxdp.com/ja/howtouse/tellus_os/start_tellus_os.html

宙畑で何度か紹介している衛星データのひとつ「SARデータ」について、簡単に目で見て分かる、時期の異なる衛星画像の差分抽出にチャレンジしました！

SARデータのおさらいから、実際のコード解説と複数地点での画像をご紹介します。

※本記事は宙畑メンバーが気になったヒト・モノ・コトを衛星画像から探す不定期連載「宇宙データ使ってみた-Space Data Utilization-」の第16弾です。まだまだ修行中の身のため至らない点があるかと思いますがご容赦・アドバイスいただけますと幸いです！

(1)SARの仕組み

まずは簡単にSARデータについておさらいします。

一般的に我々が知っているカメラで撮影する画像は、詳細には光学画像と呼ばれ、視覚情報が記録されたデータとなっています。上の絵であるように、物体の色や形が分かり、我々は写っているものが何であるかを判別します。

太陽やライトなど光源からの光の反射を見ているため、暗い部屋では様子が分からないのと同様に、光学画像では夜や雲がかかると様子が分からなくなってしまいます。

一方のSARデータは、触覚に例えられ物体の形状や表面の粗さが分かります。衛星から電波を反射してその跳ね返りを見ているので、夜でも雲がかかっていても観測することができます。

詳細は以下の記事をご覧ください。

(2)時期の異なる2枚の衛星画像の差分抽出

SARデータの活用事例のひとつして、時期の異なる衛星画像の差分抽出があります。

前の画像を青く、後の画像を赤く色づけて画像を合成するとどちらかの画像にしかなかったものが色づいて見えます。

この例では、リンゴはBeforeの画像にはなくAfterの画像で増えているので合成した画像では「赤」で表現され、反対にBeforeの画像にあったフォークはAfterの画像ではなくなっているので合成した画像では「青」に見えます。どちらの画像にもあるお皿は合成された画像では「白」く見えます。逆にどちらの画像でも何もない（電波が返ってこない）部分は「黒」に見えます。

このように、時刻の異なる2枚の画像をそれぞれ青と赤に着色して合成することで、その時間で何が変わったのかを知ることができるのです。

(3)APIでSAR画像を取得する

では、実際にSAR画像を使った差分抽出をTellusを使って実装してみたいと思います。

TellusでJupyter Lab(Jupyter Notebook)を使う方法は「Tellusの開発環境でAPI引っ張ってみた」をご覧ください。

現在Tellusで触れるSAR画像はPALSAR-2とASNARO-2の2種類です。今回は時期の異なる衛星画像を比較したいので、観測期間の長いPALSAR-2を取得することとします。

TellusのAPIによりALOS-2の光学画像を取得することができます。
APIの使い方はリファレンスにまとめられています。
APIリファレンス

※画像を利用する際は、データ詳細のデータポリシーを確認して、規約を違反しないように注意してください。

https://gisapi.tellusxdp.com/PALSAR-2/{scene_id}/{z}/{x}/{y}.png

※PALSAR-2の取得可能期間は2014年～2018年ですが、現在は2018年のみ公開しています。

zとxとyはTellus OSの[データ選択]>[タイル座標]（下から2番目）のチェックをいれると表示されます。

#必要なライブラリのインポート
import os, json, requests, math, cv2
from PIL import Image
from skimage import io
from io import BytesIO
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

TOKEN = 'ここには自分のアカウントのトークンを貼り付ける'
DOMAIN = 'tellusxdp.com'

# PALSAR2イメージ取得
def get_PALSAR2_img(sceneid, z, x, y):
    url = 'https://gisapi.tellusxdp.com/PALSAR-2/{0}/{1}/{2}/{3}.png'.format(sceneid,z,x,y)
    headers = {
        'Authorization': 'Bearer ' + TOKEN
    }
    
    r = requests.get(url, headers=headers)
    if r.status_code is not 200:
        raise ValueError('status error({}).'.format(r.status_code))
    return io.imread(BytesIO(r.content))
#イメージ表示
img = get_PALSAR2_img('ALOS2221192900', 13, 7277, 3227)
io.imshow(img)

トークンはマイページのＡＰＩアクセス設定（要ログイン）から取得してください。

実行すると画像が表示されます。

東京湾周辺の様子が見えていますね。

(4)時期の異なる衛星画像を比較する

これを応用して時期の異なる衛星画像の重ね合わせにチャレンジしてみます。

まずは、該当する画像のscene_idを探します。
今回はTellus OSのメタデータ検索を使います。

使い方はこちらの動画で紹介しています。

今回選んだのは以下の2つです。

①Before画像

scene_id : ALOS2221192900
撮影日：Thu, 28 Jun 2018 02:42:58 GMT

②After画像

scene_id : ALOS2243962900
撮影日：Thu, 29 Nov 2018 02:43:00 GMT

この画像をそれぞれ、青と赤に変換します。

# 色変換
TYPE_BLUE = 1
TYPE_RED = 2

def img_change_color(img,colortype):
    
    temp_img = cv2.split(img)
    
    # ゼロ埋めの画像配列
    if len(img.shape) == 3:
        height, width, channels = img.shape[:3]
    else:
        height, width = img.shape[:2]
        channels = 1
    zeros = np.zeros((height, width), img.dtype)
    #TYPE_BLUEだったら青だけ、それ以外（TYPE_RED）だったら赤だけ残し、他をゼロで埋める
    if colortype == TYPE_RED:
        return cv2.merge((temp_img[1],zeros,zeros))
    else:
        return cv2.merge((zeros,zeros,temp_img[1]))

これを使って、それぞれの画像の色を変更します。

この2枚の画像を一つの画像に重ね合わせます。

# 重ね合わせ    
def img_overray(img1, img2, accuracy):
    img_merge = cv2.addWeighted(img1, 1, img2, 1, 0)
    
    height, width, channels = img_merge.shape[:3]

    saveimg = np.zeros((height, width, 3), np.uint8)
    
    for x in range(height):
        for y in range(width):
            b = img_merge[x,y,0]
            g = img_merge[x,y,1]
            r = img_merge[x,y,2]
            
            if (r > b-accuracy and r < b+accuracy): 
                saveimg[x,y] = [b,b,b]
            else :
                saveimg[x,y] = [r,g,b]
                
    return saveimg

ここで、accuracyは2枚の画像を重ね合わせる際の調整の役割を果たします。

実行すると以下のような結果が出てきます。

湾の船の出入りが青と赤の点で見えていますね。

最後にコード全体を貼り付けておきます。

import os, json, requests, math, cv2
from PIL import Image
from skimage import io
from io import BytesIO
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

TYPE_BLUE = 1
TYPE_RED = 2

TOKEN = '自分のトークンを貼り付ける'
DOMAIN = 'tellusxdp.com'

#opncv -> pil 変換
def cv2pil(image):
    new_image = image.copy()
    if new_image.ndim == 2:  # モノクロ
        pass
    elif new_image.shape[2] == 3:  # カラー
        new_image = cv2.cvtColor(new_image, cv2.COLOR_BGR2RGB)
    elif new_image.shape[2] == 4:  # 透過
        new_image = cv2.cvtColor(new_image, cv2.COLOR_BGRA2RGBA)
    new_image = Image.fromarray(new_image)
    return new_image

# PALSAR2イメージ取得
def get_PALSAR2_img(sceneid, z, x, y):
    url = 'https://gisapi.tellusxdp.com/PALSAR-2/{0}/{1}/{2}/{3}.png'.format(sceneid,z,x,y)
    headers = {
        'Authorization': 'Bearer ' + TOKEN
    }
    
    r = requests.get(url, headers=headers)
    if r.status_code is not 200:
        raise ValueError('status error({}).'.format(r.status_code))
    return io.imread(BytesIO(r.content))

# 色変換
def img_change_color(img,colortype):
    
    temp_img = cv2.split(img)
    
    # ゼロ埋めの画像配列
    if len(img.shape) == 3:
        height, width, channels = img.shape[:3]
    else:
        height, width = img.shape[:2]
        channels = 1
    zeros = np.zeros((height, width), img.dtype)
    
    if colortype == TYPE_RED:
        return cv2.merge((temp_img[1],zeros,zeros))
    else:
        return cv2.merge((zeros,zeros,temp_img[1]))
    
# 重ね合わせ    
def img_overray(img1, img2, accuracy):
    img_merge = cv2.addWeighted(img1, 1, img2, 1, 0)
    
    height, width, channels = img_merge.shape[:3]

    saveimg = np.zeros((height, width, 3), np.uint8)
    
    for x in range(height):
        for y in range(width):
            b = img_merge[x,y,0]
            g = img_merge[x,y,1]
            r = img_merge[x,y,2]
            
            if (r > b-accuracy and r < b+accuracy): 
                saveimg[x,y] = [b,b,b]
            else :
                saveimg[x,y] = [r,g,b]
                
    return saveimg


# 比較するPALSAR-2のsceneid
#  赤くする画像
scene1 = 'ALOS2241452930'

#  青くする画像
scene2 = 'ALOS2218682930'

# タイル座標/精度
dstx = 3539
dsty = 1635
zoom = 12
aaccuracy = 30

# タイル取得範囲
tilenum = 3

# タイル4枚つなげて一つにする
img_blue = Image.new('RGB',(256 * tilenum, 256 * tilenum)) 
img_red = Image.new('RGB',(256 * tilenum, 256 * tilenum)) 
img_merge = Image.new('RGB',(256 * tilenum, 256 * tilenum)) 

for x in range(3):
    for y in range(3):
        print(dstx + x, dsty + y)
        img1 = get_PALSAR2_img(scene1, zoom, dstx + x, dsty + y)
        img2 = get_PALSAR2_img(scene2, zoom, dstx + x, dsty + y)

        img1 = img_change_color(img1, TYPE_BLUE)
        img2 = img_change_color(img2, TYPE_RED)
        
        img3 = img_overray(img1,img2,aaccuracy)
        
        img_merge.paste(cv2pil(img3), (x*256, y*256))        
        
        
#表示用
temp_img = Image.new('RGB',(256 * tilenum,256 * tilenum))
temp_img.paste(img_merge,(0,0))

im_list = np.asarray(temp_img)
plt.imshow(im_list)
plt.show()    

(5)色々な場所の変化を見てみる

同じ要領で他の場所も見てみます。
現在搭載されているのは、2018年の後半のデータのため、その時期に起こった出来事に絡めて考えてみます。

◆田んぼの収穫の様子

畑が収穫されている様子（青い箇所箇所）が見えます。田んぼの場合、水が張ってある時は暗く見え、収穫し裸地が見えると明るく見えるようになるため、After画像（青）で明るく見えているようです。

何かは分かりませんが、川沿いに3点ほどなにか明るく反射するものが減っているのも興味深いですね。

# 比較するPALSAR-2のsceneid
# （赤） 2018/06/31 青森
scene1 = ALOS2226072790'

# （青）2018/10/23 青森
scene2 = 'ALOS2238492790'

#位置情報
dstx = 7290
dsty = 3073
zoom = 13

◆土砂崩れの発生状況

2018年9月6日に発生した北海道胆振東部地震では、大規模な土砂災害被害が発生しました。
衛星でもその様子が見えています。

# 比較するPALSAR-2のsceneid
# （赤）2018/7/14 北海道
scene1 = 'ALOS2219122750'

# （青）2019/9/20 北海道
scene2 = 'ALOS2233612750'

# タイル座標
dstx = 7324
dsty = 3017
zoom = 13

◆高速道路建設の進捗状況

大きなインフラの建設状況も衛星で追うことができます。2019年度完成予定で進められている東海環状自動車道（大野神戸－大垣西）の進捗状況を衛星で見てみると、道路の建設が進められている様子が見えます。

※C3 東海環状自動車道（大野神戸－大垣西）の進捗状況／NEXCO中日本

# 比較するPALSAR-2のsceneid
# （赤）2018/11/11 岐阜
scene1 = 'ALOS2241302900'

# （青）2018/06/10 岐阜
scene2 = 'ALOS2218532900'

# タイル座標
dstx = 14408
dsty = 6466
zoom = 14

◆船舶の往来

SAR画像では、海面の上を走行する船の様子が良く見えます。
東京湾を見てみると、多くの船が日々行きかっている様子がよく分かりますね。

機械学習で船の検出ができれば、港を日々往来する船の数の検出ができ、経済状況などの推定ができるかもしれません。

# 比較するPALSAR-2のsceneid
# （赤）2018/6/28 東京湾
scene1 = 'ALOS2221192900'

# （青） 2019/11/29 東京湾
scene2 = 'ALOS2243962900'

# タイル座標
dstx = 3638
dsty = 1614
zoom = 12

(6)まとめ

今回は、PALSAR-2のSAR画像を使って時期の異なる2枚の画像の差分抽出にチャレンジしました。

今回記載した内容はあくまで超ざっくり2枚の画像の差分を見るものではありますが、それでも様々な場所の変化を見ることができました。

みなさんも気になるエリアでぜひ試してみてください！

また、観測画像が多いため今回はPALSAR-2で試行しましたが、Tellusに搭載されているもう一つのSAR画像ASNARO-2でも同様のことができるはずです。

ASNARO-2では解像度が良い半面、わずかなズレも見えてしまうので、どこまで差分を見ることができるのか、こちらもぜひ試してみていただけると幸いです。

「Tellus」で衛星データを触ってみよう！

日本発のオープン＆フリーなデータプラットフォーム「Tellus」で、まずは衛星データを見て、触ってみませんか？
★Tellusの利用登録はこちらから