解析ノートブック 2020/1/24

【コード付き】Sentinelの衛星データをAPI経由で取得してみた

Tellusにある衛星データとSentinelの衛星データを組み合わせた高度な解析を行いたい。そのファーストステップとして、まずはTellusの開発環境からSentinelの衛星データをAPI経由で取得してみました。

記事作成時から、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

Tellusで衛星データを触り始めた方の中には、もっと高頻度で更新される衛星データを触りたい！と思った方もいるのではないでしょうか？

そのような方に向けて、今回は欧州宇宙機関(ESA)が展開しているSentinel HubからSentinel衛星シリーズの1つであるSentinel-2のデータをAPIで取得する方法についてご紹介します。

Tellusの開発環境からSentinelの衛星データをAPI経由で取得できるようになることで、Tellusにある衛星データとSentinelの衛星データを組み合わせた高度な解析が可能になります。
※Sentinelの情報はTellus環境外からも取得できますが、Windowsなどを利用されている場合には、一部コマンドが動作しないことがあります。
※動作しない場合には、”動作しないライブラリ名”+”インストール”+”windows”等で検索してください。

1.Sentinel衛星シリーズとは

Sentinel衛星シリーズでは、ESAが展開しているCopernicus Programmeの一環として、地球の陸海空を観測しています。現行運用中の衛星は下記図の通り。

大気を観測する4号機、波高を観測する6号機も打ち上げに向けて準備が進められており、6種類の衛星で地球を観測していく予定です。

各衛星の観測結果は欧州宇宙機関などのプラットフォームから無償で公開されています。観測対象は欧州域だけでなく、全地球ということもあり、世界中のユーザが利用しています。

また、データを公開するのみならず、ユーザインタフェースを設けたプラットフォームも整備しており、以下、できることを簡単にご紹介します。

ひとつは、Tellus OSと同じようなグラフィックユーザインタフェース(GUI)です。下記のEO Browserでは各衛星の観測結果を日時・場所から絞り込むことで表示できます。例えば、面積を求めることが可能です。

Source : http://apps.sentinel-hub.com/eo-browser/

その他にも、以下キャプチャのように囲んだ範囲の植生の変化を求めることも可能です。秋ごろに植生指数の変化のピークが見られ、その後に下がっていることから、秋ごろに収穫する農作物を植えている地域なのだろうと推測できます。

最後に、本記事で実施する内容に関わりますが、アプリケーションプログラミングインタフェース(API)を叩くことでデータを取得できるようになっています。APIは開発環境に左右されることなく、自身の開発環境からでも叩くことができます。

では、いよいよ本記事の本丸へ。今回は2号機の衛星データをAPIから取得してみる方法をご紹介します。

SentinelのAPI詳細はこちらをご覧ください：https://scihub.copernicus.eu/userguide/

2.Sentinel Hubのユーザ登録方法

まず、APIを利用するためにはユーザ登録が必要です。ユーザ登録のためには下記サイトにアクセスします。

https://scihub.copernicus.eu/dhus/

アクセスするとこのような画面になります。

右上にある、右から3つ目の人型のアイコンをクリックします。

ポップアップ表示された画面の左下に表示されているSign upをクリックします。
そうすると下記のようにユーザ登録画面になります。必要事項記入の上、ユーザアカウントを作成します。

ユーザアカウント作成方法の詳細は特にこちらに記載されています：

https://scihub.copernicus.eu/userguide/SelfRegistratio

UsernameとPasswordはこのサイトにログインするときのみならず、APIを利用する際に使うので、メモしておきましょう。

3.環境構築

次に、APIを叩いて取得したデータを利用するために必要となるライブラリをインストールします。下記は必要十分なライブラリというわけではなく、あくまでこれから例示するプログラムで必要であったライブラリを示しています。ご自身の用途に合わせて、必要なライブラリをインストールしてください。

!pip install sentinelsat
!pip install rasterio
!apt install gdal-bin python-gdal python3-gdal 
!apt install python3-rtree 
!pip install git+git://github.com/geopandas/geopandas.git
!pip install descartes 
!pip install shapely
!pip install six
!pip install pyproj
!pip install descartes
!pip install geopandas

なお、冒頭にも記載しましたが、Sentinelのデータは公開データのため、当然Tellusの環境でなくても動作します。ローカル環境を用いてインストールする場合には、下記コマンドを実行するだけでは動作しないものもありますのでご注意ください。

例えば、
!pip install git+git://github.com/geopandas/geopandas.git
をローカル環境で実行するとエラーが出ます。”Git python インストール”等で検索すると、インストール方法が紹介されていますのでご覧ください。

4. APIでSentinel-2の画像を取得する

4.1. 関心領域の画像情報を取得する

では、Sentinel-2の画像を取得していきましょう。
まず、利用するライブラリをインポートします。データを地図で可視化するfoliumや、pandasの拡張機能であり衛星画像を処理する上で便利なgeopandas、衛星画像にマスク処理などをする際に便利なrasterioなどを利用します。

import os
import numpy as np

from sentinelsat import SentinelAPI, read_geojson, geojson_to_wkt 
import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

from shapely.geometry import MultiPolygon, Polygon
import rasterio as rio
from rasterio.plot import show
import rasterio.mask
import fiona
import folium

まず、データを取得したい領域の座標情報を取得します。今回はキーン州立大学が公開しているツールを利用し、画像を取得したい領域の座標を入手します。
“Close Shape”(下記図、赤枠)をクリックすることで、描画している図形を閉じることが可能です。例えば、右クリックで4点選択した後に“Close Shape”を押すと、四角形のポリゴンを入手できます。“Close Shape”を押さないと、後ほどの工程でエラーになることがあります。

from IPython.display import HTML
HTML(r’ src=”https://www.keene.edu/campus/maps/tool/” width=”960″ height=”480″ frameborder=”0″>
※実行する際には上記r’の次に、'<iframe’の中の文字列をコピー・挿入して実行してください

図形を描画したら下記でドラッグしている範囲の情報をコピーします。

上記でコピーした情報を下記にペーストしてください。これで、取得する範囲を定義できました。

AREA =  [
      [
        -220.291841,
        35.6593884
      ],
      [
        -220.2932143,
        35.4817801
      ],
      [
        -220.1380324,
        35.4817801
      ],
      [
        -220.1421523,
        35.6493456
      ],
      [
        -220.291841,
        35.6593884
      ]
    ]

定義した領域のポリゴン情報を生成し、ファイル名をつけます。

for i in range(len(AREA)):
    AREA[i][0] = AREA[i][0] +360
from geojson import Polygon
m=Polygon([AREA]) 

#ファイル名を定義. 好きな名称に設定してください. 
object_name = 'Tokyo_Bay'

import json
with open(str(object_name) +'.geojson', 'w') as f:
    json.dump(m, f)
footprint_geojson = geojson_to_wkt(read_geojson(str(object_name) +'.geojson'))

正しくポリゴンデータを定義できていない場合には、上記でエラーが出ます。エラーが出た場合には、コピペが正しくできているか、選択領域が広すぎないか、などを確認しましょう。

Sentinelからデータを取得するために、先ほど取得したSentinel Hubのユーザ情報を入力します。詳しくはこちらをご覧ください。

user = '自身のユーザ名を貼り付けてください' 
password = ' 自身のユーザ名を貼り付けてください' 
api = SentinelAPI(user, password, 'https://scihub.copernicus.eu/dhus')

定義した領域が合っているか、foliumの地図閲覧機能を用いて関心域を確認します。

m = folium.Map([(AREA[0][1]+AREA[len(AREA)-1][1])/2,(AREA[0][0]+AREA[len(AREA)-1][0])/2], zoom_start=10)

folium.GeoJson(str(object_name) +'.geojson').add_to(m)
m

欲しい領域であることが確認できたので、データを指定していきましょう。
今回は光学画像であるSentinel-2の画像を取得するため、既に指定している「場所」の情報以外である以下4つ、
・対象とする衛星
・期間
・データの処理レベル
・被雲率
を指定します。
Sentinel-2のプロダクトについてはSentinelの説明サイトをご参考にしてください。

products = api.query(footprint_geojson,
                     date = ('20191201', '20200110'), #取得希望期間の入力
                     platformname = 'Sentinel-2',
                     processinglevel = 'Level-2A',
                     cloudcoverpercentage = (0,100)) #被雲率（0％〜100％）

上記コマンドを実行した際に、例えば下記のようなエラーが表示される場合は、関心領域がSentinel-2の提供画像エリアよりも大きいか、適切な画像がない場合に多くみれます。関心領域を縮小するなどの対応にて、エラーが発生しない関心領域を設定してみてください。
SentinelAPIError: HTTP status 200 OK: Invalid query string. Check the parameters and format.

選択した範囲の取得希望期間、被雲率にあう観測履歴数を確認します。

len(products)

今回の例では、8枚画像があることが分かりました。

画像を取得する際に、雲がない方が良いですよね。そのため、要求期間で一番雲が少ない画像を選択します。取得できるファイルを一覧で表示しましょう。

products_gdf = api.to_geodataframe(products)
products_gdf_sorted = products_gdf.sort_values(['cloudcoverpercentage'], ascending=[True])
products_gdf_sorted

一覧で取得した画像を、被雲率が低い順に表示しましょう。

products_gdf_sorted.head()

ソートして先頭にある観測画像をダウンロードしましょう。

uuid = products_gdf_sorted.iloc[0]["uuid"]
product_title = products_gdf_sorted.iloc[0]["title"]
api.download(uuid)

取得したい画像情報を入手することができました。

なお、上記で取得できているURLにアクセスすることで、ローカル環境に直接データを取得することも可能です。ローカル環境で解析したい方はこちらからダウンロードしてみるのもよいでしょう。

4.2. 関心域の画像を抽出する

rasterioを用いて、取得したSentinel-2の観測画像を確認します。
今回は、Sentinel-2のBand4,3,2によるTrue Imageを作成しましょう。
波長の組み合わせによりどのような画像を作れるの？という方はこちらの記事をご覧ください。

file_name = str(product_title) +'.zip'
import zipfile
with zipfile.ZipFile(file_name) as zf:
 zf.extractall()

path = str(product_title) + '.SAFE/GRANULE'
files = os.listdir(path)

pathA = str(product_title) + '.SAFE/GRANULE/' + str(files[0])
files2 = os.listdir(pathA)

pathB = str(product_title) + '.SAFE/GRANULE/' + str(files[0]) +'/' + str(files2[1]) +'/R10m'
files3 = os.listdir(pathB)

画像を作るために、バンドを指定します。今回はバンド2をRに、バンド3をGに、バンド4をBに割り当てます。

path_b2 = str(product_title) + '.SAFE/GRANULE/' + str(files[0]) +'/' + str(files2[1]) +'/R10m/' +str(files3[0][0:23] +'B02_10m.jp2')
path_b3 = str(product_title) + '.SAFE/GRANULE/' + str(files[0]) +'/' + str(files2[1]) +'/R10m/' +str(files3[0][0:23] +'B03_10m.jp2')
path_b4 = str(product_title) + '.SAFE/GRANULE/' + str(files[0]) +'/' + str(files2[1]) +'/R10m/' +str(files3[0][0:23] +'B04_10m.jp2')

バンド4の画像を対象に、画像サイズや座標系を確認します。

b4 = rio.open(path_b4)
b3 = rio.open(path_b3)
b2 = rio.open(path_b2)

b4.count, b4.width, b4.height

b4.crs

今回の例の場合は、サイズは(1, 10980, 10980)、座標系はCRS.from_epsg(32654)と表示されました。

10980はピクセル数を表し、Sentinelは解像度が10 mGSD(つまり1ピクセルが10 m)なので、一辺が10980 pixel × 10 m/pixel≒100 kmの画像であることが分かります。

EPSGの326544というのは、地図投影する際によく利用される、WGS84のUTM座標系のことを表します。EPSGについて詳しく知りたい方はこちらをご参照ください。
では、画像を表示してみましょう。

fig, ax = plt.subplots(1, figsize=(20, 20))
show(b4, ax=ax)
plt.show()

バンド4の画像が表示されました。正しく表示できていそうなので、バンドを合成してTrue Colorの画像を表示していきましょう。

with rio.open(str(object_name) +'.tiff','w',driver='Gtiff', width=b4.width, height=b4.height, 
              count=3,crs=b4.crs,transform=b4.transform, dtype=b4.dtypes[0]) as rgb:
    rgb.write(b2.read(1),3) 
    rgb.write(b3.read(1),2) 
    rgb.write(b4.read(1),1) 
    rgb.close()

RGB_tokyo =rio.open(str(object_name) +'.tiff')
RGB_tokyo.crs

nReserve_geo = gpd.read_file(str(object_name) +'.geojson')
epsg = b4.crs

このまま画像をダウンロードしてもよいのですが、このファイルは見たい範囲以外も含んだSentinel-2の1シーンの画像となっています。ファイルサイズも大きいため、関心領域のみを抽出してダウンロードした方が効率的です。そのため、以下で抽出していきます。

nReserve_proj = nReserve_geo.to_crs({'init': epsg})

with rio.open(str(object_name) +'.tiff') as src:
    out_image, out_transform = rio.mask.mask(src, nReserve_proj.geometry,crop=True)
    out_meta = src.meta.copy()
    out_meta.update({"driver": "GTiff",
                 "height": out_image.shape[1],
                 "width": out_image.shape[2],
                 "transform": out_transform})
    
with rasterio.open('Masked_' +str(object_name) +'.tif', "w", **out_meta) as dest:
    dest.write(out_image)

msk = rio.open(r'Masked_' +str(object_name) +'.tif')
fig, ax = plt.subplots(1, figsize=(18, 18))
show(msk.read([1,2,3]))
plt.show

実行すると、下記のような画像が描画されると思います。画像自体は16ビットな中で、8ビットの画像として自動でヒストグラムを調整し、描画するため、このような見た目となります。

いちいちヒストグラム調整をして画像を表示するのは手間なので、ローカル環境でも表示しやすいjpg形式にファイルを変換してみましょう。

jpgに変換すると勝手にヒストグラムが調整される…？

from osgeo import gdal

scale = '-scale 0 250 0 30'
options_list = [
    '-ot Byte',
    '-of JPEG',
    scale
] 
options_string = " ".join(options_list)

gdal.Translate('Masked_' +str(object_name) +'.jpg',
               'Masked_' +str(object_name) +'.tif',
               options=options_string)
from PIL import Image
im = Image.open('Masked_' +str(object_name) +'.jpg')
im

希望する範囲の画像を取得できました。