はじめに．
前回，人工衛星画像を無料で解析できるオープンフリーなプラットフォームである「Tellus」でのJAXAの衛星「だいち２号」の観測画像の解析方法について紹介しました．
・人工衛星画像のオープンフリーなプラットフォーム「Tellus」にてデータ解析．
今回は，同じ「Tellus」にて利用できる衛星画像の一つであるASNARO-1の光学観測画像の解析方法について紹介します．
ASNARO-1について．
ASNAROは，Advanced Satellite with New system Architecture for Observationの略で，日本語名は”先進的宇宙システム”と呼ばれるそうです．2014年に光学衛星のASNARO-1が、2018年にレーダー衛星のASNARO-2が打ち上げられました。（Wikipedia）
ASNARO−１外観図（クレジット：NEC)
大事なところは，光学衛星とレーダ衛星の２機が運用されており，両方の衛星画像ともにTellusにて扱うことができることです．
Tellus 衛星データハンドブック
光学画像は，普段からGoogle Maｐでも使っていますので，レーダ衛星の画像よりも見慣れており，標準画像ができれば，目的に応じてどういう解析・処理をすればよいのかイメージしやすいと思います．
ASNARO-1の観測画像（データ）の主な情報は以下となります．
・観測波長帯
PAN画像（単画像，と言う意味）;450-860 nm
MS（マルチスペクトル）画像（カラー画像，とい意味）;
400-450 nm (Ocean Blue)
450-520 nm (Blue)
520-600 nm (Green)
630-690 nm (Red)
705-745 nm (Red Edge)
760-860 nm (NIR)
・分解能;
PAN画像： 0.5m
MS画像： 2m
・観測幅;
10km
ASNARO-1の観測画像 諸元
欧米の有名な商用人工衛星の観測画像との比較は以下となります．
これをみると，欧米の有名な商用衛星（WorldView,Peiades）と比べても，分解能（GSD)は大きな差はないことがわかります．ただ，観測幅（Swath Width）が半分ほであり，一度の観測運用での観測幅は狭くなります．観測対象が決まっているのであれば，それほど問題にはならないと思います．
Tellusで扱えるASNARO-1の観測画像は，上記のなかのMS画像になります．
次に， ASNARO-1の観測画像のTellusの総合開発環境（Jupyter Lab)での取得，そして解析方法として複数画像の組み合わせ（モザイク処理）の方法について紹介します．
ASNARO-1 観測画像の組み合わせ処理．
Tellusのアカウントや総合開発環境の取得およびその設定について，前回の記事をご参考ください．
・人工衛星画像のオープンフリーなプラットフォーム「Tellus」にてデータ解析．
リンク先をみながらコードを書くのは手間かと思いますので，画像取得などは重複しますがここでも紹介します．
Tellusの画像取得方法は，同じく以下のサイトを参考にしています．ありがとうございます．
参考：Tellusで超低高度衛星つばめの衛星画像タイルを取得して連結してみた
from IPython.display import HTML
HTML(r'<iframe width="960" height="480" src="https://tellus-map-test.herokuapp.com/" frameborder="0"></iframe>')
上記のサイトより，関心のある地域の緯度経度情報，およびタイル情報を取得します．
では，アプリケーションを実行します．
今回は，茨城県つくば市の画像を関心対象とします．
ASNARO-1のZoom範囲は10から18までとなります．18になると高分解の画像を取得することができますが，10km幅以上（目安）の広範囲の画像を取得し，その組わせをすると基本の総合開発環境では”メモリ不足”のエラーとなりました．そのため，今回はzoomを16としています．
（Zoomって何？と思われた方は，こちらのサイトをご参考ください．）
こちらで所得した関心域の緯度・経度の情報，および地理院タイルの座標情報を以下にコピペします．
AREA = {"min_lat":36.058,"min_lon":140.072,"max_lat":36.099,"max_lon":140.145}
lat= {"lat":[25712,25722]}
lon = {"lon":[58267,58280]}
zoom = {"zoom":16}
次に，この関心域のASNARO-1の観測画像の情報を取得します．以下のコードを実行します．
#ASNARO-1の光学観測画像の関心領域のメタ情報を取得する．
def getAsnaro1Scene(token,area):
url = 'https://gisapi.tellusxdp.com/api/v1/asnaro1/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()
#関心領域のASNARO-1の光学観測画像のメタ情報の取得．
res=getAsnaro1Scene(TOKEN,AREA)
print('number of scenes : {}'.format(len(res)))
for i in range(len(res)):
print(res[i]['acquisitionDate'])
こちらを実行すると，以下のアウトプットが出力されます．
number of scenes : 3
Wed, 11 Jul 2018 12:17:02 GMT
Sun, 30 Apr 2017 14:14:44 GMT
Tue, 10 Feb 2015 01:38:37 GMT
この出力は，関心域に重なる3つの画像があることを意味しています．重要なのは”これらの画像が関心域すべてを包含している（含んでいる）わけでない”ということです．
どういう意味かは，以下に説明します．
それでは，以下のコードを実行し，最新の観測画像の全体像（サムネイル）を見てみます．
thumbs = io.imread(res[0]['thumbs_url'])
io.imshow(thumbs)
これをみると，関心域のつくば市のつくば駅周辺のすべての情報がこの１枚の観測画像で包含されているように見えます．ただ，その関心域がこの左上の部分を含んでしまうと，情報がない白い部分も取得し，それを組み合わせ画像に使用してしまいます．そうなると，できあがった画像は白い空白部分を含むものとなり，意図しないものとなります．
また，右上部のピンクや白は雲になります．この雲もできるだけ利用したくありません．
前回の記事のレーダ衛星（ALOS-2：PALSAR-2）は雲があっても，夜間であっても画像を取得することができるため，上記の雲の懸念はないのが特徴です．ただ，そのかわりに画像の判読が難しいですが．．．
空白部分および雲画像を省く方法は，後ほど紹介します．
では，この画像のメタ情報を取得します．
res[0]
出力されたメタ情報は以下となります．
{'acquisitionDate': 'Wed, 11 Jul 2018 12:17:02 GMT',
'clat': 36.0569741131053,
'clon': 140.274741142563,
'cloudCover': 'N/A',
'entityId': '20181224062528808_AS1',
'max_lat': 36.0095766445904,
'max_lon': 140.417554331119,
'min_lat': 36.0078473769819,
'min_lon': 140.132107504947,
'path': 0,
'productId': '20181224062528808_AS1_D01_L1B',
'row': 0,
'thumbs_url': 'https://tile.tellusxdp.com/thums/ASNARO-1/20181224062528808_AS1.PNG'}
次に， 関心域について，空白部分が含まず，雲ができるだけないタイル画像を取得し組み合わ画像をつくります．一枚の画像を取得する方法は，前回の記事をご参考ください．（ご希望があれば，コメントいただければ追加します．）
#ASNARO-1の複数の光学観測画像より，空白画像がないもの，および雲が多くない画像を選択する．
def get_Asnaro1_scene_tile(zoom, xtile, ytile):
#zoom = 10-18
for i in range(len(res)):
scene_id = res[i]['entityId']
url = 'https://gisapi.{}/ASNARO-1/{}/{}/{}/{}.png'.format(
DOMAIN, scene_id, zoom, xtile, ytile)
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))
i = i + 1
#print(i)
else:
#print(i)
im = io.imread(BytesIO(r.content))
im_sum = im.sum()
if np.any(im == 0) or (im_sum > 50000000) : #空白画像もしくは雲が多い画像をスキップする．
i = i + 1
else:
break
return io.imread(BytesIO(r.content))
ここのポイントは，以下のコードになります．
if np.any(im == 0) or (im_sum > 50000000) : #空白画像もしくは雲が少ない画像をスキップする．
私は，光学画像の複数組み合わせ処理（モザイク処理）の知識はないため，この方法は試行錯誤した結果なのです．
空白部分があると，その部分は”0”を含んでいあるため，それを識別します．
また，対象画像に雲が含まれるとその値はダイナミックレンジの最大値である”255”を多く含んでいます．そのため，画像の数値をすべて足し算（積算）したとき，異常に大きいときは雲を多く含むと識別しました．（ここでは，50000000以上としました．）
取得した画像が上記のどちらかとなる場合は，それをスキップし別の画像を取得します．
以下のコードを実行すると，関心域のタイル画像を組みあわせた画像を取得することができます．
img_AS1 = getAsnaro1Img(zoom["zoom"],lat["lat"][0], lat["lat"][1], lon["lon"][0], lon["lon"][1])
望み通り，空白と雲を含まない組み合わせ画像を取得することができました．
ただ，見て分かる通り，右と左で明るさの異なる画像となっているのがわかります．
これは，光学観測画像は，観測するタイミングによって，雲や季節により写り方がかわるため，組み合わせるとこういった明暗が異なる画像となってしまいます．
それを認識して解析に用いることもできますが，その後の解析によっては統一化（標準化）したほうが便利なため，ここではこの画像の標準処理を試してみました．
pythonによる画像の標準処理はいくつか方法があり，例えば”画像の輝度を均一化する”という方法も試してみましたが，確かに均一にみえますが”白飛び”になってしまい，あまりよい画像を得られることができませんでした．
（色々試したおかげで，多くの知識を得ることができました．これがスキルを得る一番効果的な方法であり，この試行錯誤がエンジニアにとって一番楽しい時間と感じています．）
さて，どういう方法がよいのか試していたところ，以前の衛星画像から建物地図を作るときにもちいた画像の標準化処理が有効では，と思って使ってみました．
Deep Learning(深層学習)のための画像の標準化，リサイズ，分割，結合処理．
ここで用いた画像の標準化処理のコードは以下となります．
#画像の輝度の標準化処理．
def normalize(arr):
arr = arr.astype('float')
for i in range(3):
minval = arr[...,i].min()
maxval = arr[...,i].max()
if minval != maxval:
arr[...,i] -= minval
arr[...,i] *= (255.0/(maxval-minval))
return arr
これを複数画像を組み合わせするコードに追加します．
#ASNARO-1の光学観測画像の関心領域の画像の組み合わせ（モザイク処理．画像の標準化処理）
def getAsnaro1ImgC(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_Asnaro1_scene_tile( zoom, item_h, item_v)
img_h = normalize(img_h)
img_h = img_h.astype('int32')
else:
img_h1 = get_Asnaro1_scene_tile( zoom, item_h, item_v)
img_h1 = normalize(img_h1)
img_h1 = img_h1.astype('int32')
img_h=np.concatenate((img_h, img_h1), axis=1)
if (img is None):
img = img_h
else:
img = np.concatenate((img, img_h), axis=0)
io.imshow(img)
#print(res[2]['acquisitionDate'])
print(img.shape)
return img
上記を実行します．
img_AS1C = getAsnaro1ImgC(zoom["zoom"],lat["lat"][0], lat["lat"][1], lon["lon"][0], lon["lon"][1])
実行結果が以下となります．
標準処理を含まないものと比べると，画像内での明暗が少なくなっているのがわかります．
拡大して比較してみます．
複数タイル画像の組み合わせ画像（標準化処理なし）
複数タイル画像の組み合わせ画像（標準化処理あり）
拡大していみると，画像の明暗の差がかなり低下できているのがわかります．
これで，建物セグメンテーションのための前処理は終了です．
（追加更新） 衛星画像のアンシャープン処理．
3連休(9月14日から16日)で時間があったので，高解像度（Zoom=18）の画像の処理について試してみました．
衛星画像の利用方法のひとつとして，深層学習による駐車場の車のカウントがあげられます．
人工衛星画像とAI技術（㈱パスコ）
変わった解析方法として，レーダ衛星の画像を用いた駐車場の状況を推測する方法も提案されています．（マニアックだ．）
いつ空いてるの!? 無料衛星データでディズニーランドの混雑予想チャレンジ（前編）(宙畑)
まずは，王道の光学観測衛星の画像からの駐車場の車の数をかぞえるため，そのベースとなる画像の前処理について考えました．
対象とした駐車場は，建設中の新国立競技場の聖徳記念絵画館の駐車場になります．
上記に記載した方法にて対象地域の緯度・経度情報を取得します．
AREA = {"min_lat":35.674,"min_lon":139.711,"max_lat":35.682,"max_lon":139.722}
lat= {"lat":[103226,103233]}
lon = {"lon":[232806,232814]}
zoom = {"zoom":18}
上記の緯度・経度情報，およびその対象となる地理院タイルの座標をインプットし，TellusのASNARO-1の対象の画像を検索します．
#関心領域のASNARO-1の光学観測画像のメタ情報の取得．
res=getAsnaro1Scene(TOKEN,AREA)
print('number of scenes : {}'.format(len(res)))
for i in range(len(res)):
print(res[i]['acquisitionDate'])
検索結果が以下となります．
number of scenes : 22
Wed, 31 Oct 2018 12:31:16 GMT
Sat, 20 Oct 2018 12:52:04 GMT
Tue, 31 Jul 2018 12:42:49 GMT
Sat, 30 Jun 2018 12:17:46 GMT
Fri, 08 Jun 2018 12:25:21 GMT
Sat, 02 Jun 2018 12:36:46 GMT
Sat, 21 Apr 2018 12:38:02 GMT
Sun, 25 Mar 2018 12:27:00 GMT
Wed, 03 Jan 2018 15:30:59 GMT
Tue, 02 Jan 2018 02:16:36 GMT
Fri, 29 Dec 2017 15:24:58 GMT
Sun, 05 Nov 2017 13:55:39 GMT
Sat, 06 May 2017 15:29:27 GMT
Tue, 04 Apr 2017 15:27:01 GMT
Mon, 24 Oct 2016 14:10:19 GMT
Fri, 17 Jun 2016 15:39:32 GMT
Thu, 05 May 2016 14:08:24 GMT
Fri, 18 Mar 2016 15:31:58 GMT
Mon, 30 Nov 2015 14:31:57 GMT
Tue, 03 Nov 2015 14:32:11 GMT
Thu, 06 Aug 2015 03:12:35 GMT
Thu, 19 Feb 2015 12:40:48 GMT
東京都内ということで，多くの画像がストックされています．
最新の画像は，2018年10月31日に撮像された画像になります．
その画像を以下を実行し抽出します．また，取得画像の保存（ｐｎｇフォーマット）も実行します．
#ASNARO-1の光学観測画像の関心領域の画像の組み合わせ（モザイク処理）
def getAsnaro1Img(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_Asnaro1_scene_tile( zoom, item_h, item_v)
else:
img_h=np.concatenate((img_h, get_Asnaro1_scene_tile( 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[2]['acquisitionDate'])
print(img.shape)
return img
img_AS1 = getAsnaro1Img(zoom["zoom"],lat["lat"][0], lat["lat"][1], lon["lon"][0], lon["lon"][1])
Image.fromarray(img_AS1).save('asnaro1_tokyo_nationalstadium_18.png')
その結果が以下となります．
薄曇りの日だったため，全体的に暗く映っています．この画像の上記にある標準化処理を行います．
#画像の輝度の標準化処理．
def normalize(arr):
arr = arr.astype('float')
for i in range(3):
minval = arr[...,i].min()
maxval = arr[...,i].max()
if minval != maxval:
arr[...,i] -= minval
arr[...,i] *= (255.0/(maxval-minval))
return arr
img_AS1N = normalize(img_AS1)
img_AS1N = img_AS1N.astype('int32')
io.imshow(img_AS1N)
この出力が以下の画像となります．
先程よりも明るくコントラストがはっきりしています．
次に，この画像に映る対象物の輪郭をよりはっきりさせるために，画像処理をいくつか試しました．
画像のフィルタ処理は多くの手法があり，またそのパラメータの設定も重要になります．
ここでは，pythonの画像系モジュールであるPILに備わっている機能をまずは使ってみました．
Pillow (PIL) - 画像にフィルタを適用する方法
いくつか試した結果，一番良かったパンシャープン処理の方法を示します．
img_AS1N = Image.open("asnaro1_tokyo_nationalstudium_18N.png")
img_AS1NU = img_AS1N.filter(filter=ImageFilter.UnsharpMask())
img_AS1NU.save(("asnaro1_tokyo_nationalstudium_18NU.png"))
ここで，3つの画像を拡大し，聖徳記念絵画館を拡大してみます．
高倍率(zoom=18)の画像(元画像)であっても駐車場の車の識別はかなり難しいです．
その画像の標準化処理によってある程度見やすくなりましたが，それでも車の輪郭がボケています．
その画像にアンシャープン処理を行うことで，車の形がある程度識別できるようになりました．
今後，この処理を行った画像に，深層学習による画像処理を行い車の分布図を求めてみます．
Deep Learning で航空写真から自動車をカウントする
まとめと今後について．
以上の解析により，オープンフリーなプラットフォームである「Tellus」にて，高分解能な光学画像であるASNARO-1のデータ取得や，関心のある領域の複数画像の組み合わせ画像の取得方法，またその画像の標準処理ができるようになりました．
これでベースはできましたので，次はこの画像から深層学習処理により建物地図を作成してみたいと思います．
衛星画像のSegmentation（セグメンテーション）により建物地図を作成する．
ここでのタイル画像の組み合わせ処理や，画像の標準化処理は自己流になります．ここでの解析手法のコメントや，他によい方法の提案などのコメントをいただけるととても嬉しいです．
最後に，ご参考に全体のコードを示します．TOKENの*****部分はご自身のものに変更してください．
コメントがありましたら，お気軽にお知らせください．励みになります．
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
TOKEN = '****************************'
DOMAIN = 'tellusxdp.com'
#ASNARO-1の光学観測画像の関心領域のメタ情報を取得する．
def getAsnaro1Scene(token,area):
url = 'https://gisapi.tellusxdp.com/api/v1/asnaro1/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()
#ASNARO-1の複数の光学観測画像より，空白画像がないもの，および雲が多くない画像を選択する．
def get_Asnaro1_scene_tile(zoom, xtile, ytile):
#zoom = 10-18
for i in range(len(res)):
scene_id = res[i]['entityId']
url = 'https://gisapi.{}/ASNARO-1/{}/{}/{}/{}.png'.format(
DOMAIN, scene_id, zoom, xtile, ytile)
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))
i = i + 1
#print(i)
else:
#print(i)
im = io.imread(BytesIO(r.content))
im_sum = im.sum()
if np.any(im == 0) or (im_sum > 50000000) : #空白画像もしくは雲が多い画像をスキップする．
i = i + 1
else:
break
return io.imread(BytesIO(r.content))
#画像の輝度の標準化処理．
def normalize(arr):
arr = arr.astype('float')
for i in range(3):
minval = arr[...,i].min()
maxval = arr[...,i].max()
if minval != maxval:
arr[...,i] -= minval
arr[...,i] *= (255.0/(maxval-minval))
return arr
#ASNARO-1の光学観測画像の関心領域の画像の組み合わせ（モザイク処理）
def getAsnaro1Img(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_Asnaro1_scene_tile( zoom, item_h, item_v)
else:
img_h=np.concatenate((img_h, get_Asnaro1_scene_tile( 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[2]['acquisitionDate'])
print(img.shape)
return img
#ASNARO-1の光学観測画像の関心領域の画像の組み合わせ（モザイク処理．画像の標準化処理）
def getAsnaro1ImgC(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_Asnaro1_scene_tile( zoom, item_h, item_v)
img_h = normalize(img_h)
img_h = img_h.astype('int32')
else:
img_h1 = get_Asnaro1_scene_tile( zoom, item_h, item_v)
img_h1 = normalize(img_h1)
img_h1 = img_h1.astype('int32')
img_h=np.concatenate((img_h, img_h1), axis=1)
if (img is None):
img = img_h
else:
img = np.concatenate((img, img_h), axis=0)
io.imshow(img)
#print(res[2]['acquisitionDate'])
print(img.shape)
return img
#関心領域の国土地理院タイル座標の取得．
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":36.058,"min_lon":140.072,"max_lat":36.099,"max_lon":140.145}
lat= {"lat":[25712,25722]}
lon = {"lon":[58267,58280]}
zoom = {"zoom":16}
#関心領域のASNARO-1の光学観測画像のメタ情報の取得．
res=getAsnaro1Scene(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_AS1 = getAsnaro1Img(zoom["zoom"],lat["lat"][0], lat["lat"][1], lon["lon"][0], lon["lon"][1])
#取得した画像の保存
Image.fromarray(img_AS1).save('asnaro1_sample_tile_Tsukuba_16.png')
#複数画像の組み合わせ画像の取得（モザイク処理および標準化処理）
img_AS1C = getAsnaro1ImgC(zoom["zoom"],lat["lat"][0], lat["lat"][1], lon["lon"][0], lon["lon"][1])
#取得した画像の保存
io.imsave('asnaro1_sample_tile_tsukuba_16B.png',img_AS1C)
更新．
2019/9.15 誤記を修正しました．