概要
これまで，欧州の地球観測の人工衛星（Sentinelシリーズ）の画像の取得方法について紹介しました．
無料で最新の衛星画像を入手する方法．
今回は，Sentinelが提供するAPIを用いて,jupyter lab などのPuythonのプラットフォームにSentinelの観測データを管理するCopernicusのサーバからの直接取得する方法，および取得したデータの画像処理，関心域の画像抽出，最後に抽出した画像のjpeg処理について紹介します．自己流ですので，よりスマートな方法がありましたらぜひコメントにてお教えください．
なお，この記事で紹介するコードは以下の記事を参考にさせていただきました．ありがとうございます．
Satellite Imagery Access and Analysis in python jupyter notebooks
@Midium
また，縁があって，わかりやすく有益な宇宙情報を発信されている宙畑にも記事として取り上げていただきました．応用例も含めて紹介していますので，そちらもご参照ください．
【コード付き】Sentinelの衛星データをAPI経由で取得してみた
全体のコードをgitにあげています．関心のある方はこちらよりダウンロードし，それぞれの環境にて試してみてください．
Sentinel2_API_access
RS_Python_Sentinel2_API.ipynb
開発環境について
私は，Tellusが提供している総合開発環境を用いて，Sentinelシリーズの画像取得，および画像処理を行っていました．ローカルの環境でも可能と思いますが，Tellusプラットフォームであれば，Tellusが提供している他の人工衛星の画像との組み合わせもできます．ご参考にしてください．
1.Sentinelシリーズ．
Sentinelシリーズは，欧州宇宙機関が開発・運用する地球観測シリーズで，すでに１から５まで打上げられ，運用されています．
今回は，その中でも多くの方に関心が高く馴染みのある，Sentinel-2の衛星画像の取得から解析までをとりあげます．
Sentinel-2 ©ESA 2000-2018.
Sentinel-2について，その観測計画や観測画像の詳細については，無料で最新の衛星画像を入手する方法．をご参考ください．
2.Sentinel衛星の画像取得
それでは，Sentinel衛星の画像を取得します．前回は，Sentinelの衛星データを配信しているCopernicus hubのウェブサイトより取得する方法を紹介しました．例えば，関心のある地域の時系列変化を取得したいなど多くの衛星画像を取得するときは，毎回衛星画像の検索，ダウンロードの作業が必要になります．
そこで，Sentinelが提供しているAPIを用いて，pythonのコードによりサーバから直接衛星データを取得し，表示処理をおこないます．Sentinel APIは以下のサイトをご参考ください．
Sentinel API python
2.1 関心域の情報取得
Tellusのプラットフォームにおける関心域の地理情報（緯度・経度）は，@t-matさんが提供されているTellusで超低高度衛星つばめの衛星画像タイルを取得して連結してみたのWebアプリケーションを用いて取得しました．ただ，こちらの適用範囲が日本域に限定されているため，今回は海外を含めた地理情報を取得するためのWebアプリケーションとして，米国キーン州立大学が提供する以下のサイトを利用します．
Polyline Tool@keen state colleage
#関心領域のポリゴン情報の取得．
from IPython.display import HTML
HTML(r'<iframe width="960" height="480" src="https://www.keene.edu/campus/maps/tool/" frameborder="0"></iframe>')
こちらを実行すると，以下の画面が表示されます．デフォルトが米国キーン州立大学になります．
左上のアイコン（＋ー）より画面を拡大し，関心域まで移動します．今回は東京近郊のSentinel-2の画像を取得します．
右クリックで関心域をセットし，ポリゴンを作ります．間違えた場合は”Reset”を押してやり直してください．
一筆書きでポリゴンができましたら，最後に”Colse Shape”をクリックしてください．これで始点と終点が同一点に補正されポリゴンが補正されます．
その後，ここで準備した緯度経度情報を選択しコピーします．
その後，コピーしたポリゴンの緯度・経度情報を以下のコードに貼り付けます．
その後の，関心域のポリゴンのGeojsonに変換し，foliumを用いて正しくセットできたか確認します．
これで，関心域の位置情報のポリゴンを準備することができました．次に，この関心域のSentinel-2の観測画像データをSentinel APIを用いてCopernicus hubのサーバーにリクエストしデータを取得します．
2.2 Sentinel APIを利用して画像取得
Copernicus hubにアクセスするには，アカウントが必要になります．アカウントの所得方法は，以前紹介しました以下の記事をご参考にしてください．
では，アカウントをセットし，Sentinel-2の衛星データを要求します．
#Copernicus hubのアカウントをセット．
user = '*******'
password = '********'
api = SentinelAPI(user, password, 'https://scihub.copernicus.eu/dhus')
#Sentinel-2の観測画像を要求
products = api.query(footprint_geojson,
date = ('20191201', '20200110'), #取得希望期間の入力
platformname = 'Sentinel-2',
processinglevel = 'Level-2A',
cloudcoverpercentage = (0,100)) #被雲率（0％〜100％）
footPrint_geojsonは先の準備した関心域のポリゴン情報になります．
その他に4つのパラメータを要求しています．それぞれの意味については，Sentinel APIの情報をご参考ください．
ここで取得するSentinel-2のプロダクトレベルは，L1-2Aを要求しています．
これは，衛星画像の位置合わせを終えたL-1Cレベルに，大気のゆらぎ（モヤ）を補正した画像になります．
みなさんの馴染みのあるGoogle mapの衛星画像と同じ，と考えてよいと思います．
Sentinel-2のプロダクトなどについては，以下のサイトをご参考ください．
Sentinel-2プロダクト
Sentinel-2 L1-2A
Leve1-2Aの観測画像は2019年1月以降の観測データがオンラインにアップされており，上記のコマンドにて取得することができます．
それ以前の画像については，個別にCopernicus hubに要求する必要があるらしく，以前紹介した記事を参考に個別に要求してください．（この場合は，自動取得はできないと思います．試していませんが...)
では，2019年12月1日から2020年1月10日までのSentinel-2のLevel1-2Aの衛星データを取得します．
len(products)
16
上記より，この期間の関心域の画像データは16個あることがわかりました．
Sentinel-2は決められた日時に自動的に観測されており，晴天でも荒天でも観測が実行されます．
そのため，16個の画像のなかにには雲しか映っていないものも含まれます．そこで，以下のコマンドを実行し，雲が一番少ない画像を選びます．
products_gdf = api.to_geodataframe(products)
products_gdf_sorted = products_gdf.sort_values(['cloudcoverpercentage'], ascending=[True])
products_gdf_sorted
結果が以下となります．
一部のみ載せていますが，16個の衛星画像のデータのそれぞれのメタ情報を見ることができます．
この時点で，一番上の衛星画像の情報が雲が一番少ない画像（被雲率が小さい）となっています．
この画像のプロダクト番号を取得し，Copernicus hubサーバより衛星画像を取得します．
uuid = products_gdf_sorted.iloc[0]["uuid"]
product_title = products_gdf_sorted.iloc[0]["title"]
api.download(uuid)
上記を実行すると，データのダウロードが開始します．
これで，Sentinel-2のデータを取得しました．次に，このデータより衛星画像を作成します．
2.3 取得データの画像処理（Geotiff)
衛星データの画像処理は，その専用モジュールであるRasterioを用いて行います．
Rasterioについては，以下のサイトをご参考にしてください．
Rasterio Quick Guide
以下のコマンドを実行し，取得したデータのzipファイルの解凍，衛星画像作成のためのB2,B3,B4バンドのデータを取得します．
#Zipファイルの解凍
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)
＃それぞれのファイルパスをセット．
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')
# Open Bands 4, 3 and 2 with Rasterio
b4 = rio.open(path_b4)
b3 = rio.open(path_b3)
b2 = rio.open(path_b2)
ここで，赤色であるB4の衛星画像を見てみます．
fig, ax = plt.subplots(1, figsize=(20, 20))
show(b4, ax=ax)
plt.show()
関東の画像であることはわかりましたが，どうも関心域の右側しか取得できていないような...
EO browserで関心域付近の衛星データのフットプリントを確認すると，関心域は両方の画像にまたいでいたことがわかりました．
Sentinel-2の画像は100km四方に区分され，それぞれが少し重なるようにセットされています．
これをみると，今回取得した画像の右側ではなく，左側の方がよいように思えます．
そのときは，先程の被雲率でソートした画像の２番目に雲が少ない画像に変更して，同じフローをして試して見てください．
そういった手間を省きたいときは，このEo Browserをみて，事前にそれぞれの観測画像のフットプリントが重ならない関心域を選ぶことで回避できると思います．（マニュアル操作になりますが…）
（２つ以上の観測画像のタイル処理をすればいいのですが，異なる観測日の画像を組み合わせるのは，絵が変わるので私はあまり好きではなくて）
次の画像とするときは，以下から実行してみてください．
uuid = products_gdf_sorted.iloc[1]["uuid"]
product_title = products_gdf_sorted.iloc[1]["title"]
api.download(uuid)
最後に，取得した衛星データの観測バンドを組み合わせて，観測画像のカラー処理をしデータを保存します．
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()
2.4 取得画像の関心域のみを抽出（Geotiff)
先に取得した画像は，Sentinel-2の観測画像の１シーン分（100km四方）のため，関心域以外も含まれています．
そこで，先に準備した関心域のポリゴンを用いて，上記の画像から関心域を切り出します．
nReserve_geo = gpd.read_file(str(object_name) +'.geojson')
#epsg情報の取得
epsg = b4.crs
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
先程確認した懸念のとおり，関心域と選択したSentnel-2の観測画像が一部しか重なっていなかったため，この部分だけが抽出されました．本来ならば，異なる衛星データを取得しなおし，同じ処理をするほうがよいのですが，その例は宙畑で記事になっていますので，ここではこのまま進めます．
2.5 関心域のJpeg画像の取得
次に，先程抽出した衛星の観測画像（geotiff)をjpeg画像に変換します．
jpeg画像に変換する理由は，データ容量が軽くなることと，Tellusの衛星画像の紹介にてもちいたkmz処理をすると，ローカルでGoogle Earthに重畳し閲覧することができるからです．
jpeg画像のkmz処理については，以下のサイトをご参考にしてください．
人工衛星画像の無料環境なプラットフォーム「Tellus」にてデータ解析． -Google Earth編-
では，以下を実行しjpegに変換します．jpegへの変換は gdalを用いました．
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)
作成したjpeg画像を以下を実行し閲覧します．
from PIL import Image
im = Image.open('Masked_' +str(object_name) +'.jpg')
im
あまりにも細い...
みなさんが実行されるときは，適切なものに変換してくださいね．
gitにあげたコードは，きちんと例になるものに関心域を変えておきます．
これで，SentinelのAPIを用いて，衛星データのCopernicus hubサーバからの取得，画像処理，抽出処理及びjpegへの変換の一連のフローを紹介しました．
3. おわりに
Sentinel APIの使用例を紹介しました．
ここで紹介した方法をベースに，関心域の時系列のGifアニメーションなどを作ることもできます．また，Sentinel-2ではなく，Sentinel-1のデータを取得することもできます．そのときは，APIのコマンドが少し異なりますので，Sentinel APIのサイトをご参考にしてください．
以上のアイデアはすでに実行しているため，評判がよければ記事化を考えますね．
参考資料
【コード付き】Sentinelの衛星データをAPI経由で取得してみた
無料で最新の衛星画像を入手する方法．
人工衛星画像の無料環境なプラットフォーム「Tellus」にてデータ解析． -Google Earth編-
Satellite Imagery Access and Analysis in python jupyter notebooks