欧州宇宙機関が実施している地球観測ミッション「コペルニクス計画」の人工衛星Sentinel-2は、空間解像度10mの衛星画像データをフリーで提供している。
アカウント作成などに関しては、以下の記事が非常にわかりやすく参考になる。
ここでは追加の情報として、
- 行政区画などのShapeファイルの入手
- 地理座標系と投影座標系を変換して衛星画像を行政区画でクリッピング
- 正規化植生指数（陸生植物の活性度の指標）を近赤外域反射率、赤色光反射率から計算
する。
利用するPythonライブラリは以下のとおり。
# いつもの
import os
import numpy as np
# SentinelデータにアクセスするためのPython API
from sentinelsat import SentinelAPI
# 地理情報解析のためのPandasの拡張
import geopandas as gpd
# Leaflet.jsを利用した地図情報の可視化
import folium
# 地理情報システムで使われるGeoTiffを読み書きする
import rasterio
まずはCopernicus Open Access Hubで登録したユーザ名、パスワードを使ってAPIを用意。
user = 'USER'
password = 'PASSWORD'
api = SentinelAPI(user, password, 'https://scihub.copernicus.eu/dhus')
三島市形状データのダウンロード
描画したい行政区画（あるいは国立公園とか、なんらか興味のある区域）を決める。
日本国内の様々なエリアのShapeファイル（幾何学的形状や座標情報などを記述したデータ）は、国土交通省の国土数値情報のページでダウンロードできる。
ここでは静岡県三島市を描画することにする。（地元なので）
国土数値情報のページから、
国土数値情報ダウンロード
→ 政策区域
→ 行政地域
→ 行政区域（ポリゴン）
→ 東海：静岡県：令和二年（N03-200101_22_GML.zip）
とたどった先にzipのリンクがある。
一応、利用規約に基づく以下全般の原典表示。
出典：国土交通省国土数値情報ダウンロードサイト
「国土数値情報（シェープファイル形式データ）」（国土交通省）（ https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v2_4.html ）を加工して作成
ダウンロードしたデータを解凍すると、以下のようなファイルが入っている。
!tree ./data/N03-200101_22_GML/
./data/N03-200101_22_GML/
├── KS-META-N03-20_22_200101.xml
├── N03-20_22_200101.dbf
├── N03-20_22_200101.geojson
├── N03-20_22_200101.prj
├── N03-20_22_200101.shp
├── N03-20_22_200101.shx
└── N03-20_22_200101.xml
0 directories, 7 files
Shapeファイルは拡張子が.shpのやつ。これをGeoPandasで読み込む。
読み込まれたテーブルには静岡県全体の市区町村が含まれているので、三島市の項目を抜き出す。
shizuoka_shape = gpd.read_file('data/N03-200101_22_GML/N03-20_22_200101.shp',
encoding='shift-jis')
mishima_shape = shizuoka_shape[shizuoka_shape['N03_004'] == '三島市']
display(mishima_shape)
| N03_001 | N03_002 | N03_003 | N03_004 | N03_007 | geometry | |
|---|---|---|---|---|---|---|
| 137 | 静岡県 | None | None | 三島市 | 22206 | POLYGON ((138.92374 35.07665, 138.92378 35.076... |
"geometry"の部分が形状を示している。
このデータの座標参照系(CRS)はJGD2011（EPSGコードは6668）で記述されている。いまは使わないけど、あとで衛星画像に合わせてこの座標参照系を変換する。
mishima_shape.crs
<Geographic 2D CRS: EPSG:6668>
Name: JGD2011
Axis Info [ellipsoidal]:
- lon[east]: Longitude (Degree)
- lat[north]: Latitude (Degree)
Area of Use:
- undefined
Datum: Japanese Geodetic Datum 2011
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich
形状を可視化してみる。foliumで地図を呼び出して、三島市の形状を加えて描画する。中心の座標はどこでもいいけど適当にわかりやすい場所に。
latlon = [35.117655, 138.9382308]
m = folium.Map(latlon, zoom_start=11)
folium.GeoJson(mishima_shape).add_to(m)
m
衛星画像のダウンロード
さて、この形状をAPIに指定して該当する衛星画像をダウンロードするのだが、これをそのまま送信すると頂点の数が多すぎてAPIがクエリーを受け付けず検索に失敗する。
なので、三島市を包含する単純な四角形を計算して、その頂点データを送信することにする。
これはShapelyを使うと簡単。オブジェクトを包含する最小の四角形ポリゴンを返してくれる。
mishima_envelope = mishima_shape['geometry'].values[0].envelope
m = folium.Map(latlon, zoom_start=11)
folium.GeoJson(mishima_envelope).add_to(m)
m
この形状をAPIのクエリとして、三島市を含む衛星画像を検索する。
本当はもっと最近の画像を使いたかったんだけど、いろんな事情があってここでは2020年1月に撮影された画像から検索している。
いろんな事情：夏場は雲がかかっちゃうことが多くてあまり使える画像がない。あとそもそも三島市は画像タイルのちょうど「切れ目」にあたることが多く、単一の画像で三島市全域を含んでくれてる画像の選択肢があまりない。複数タイルを扱うのはめんどうなので、偶然全域を含んでくれてる時期で選んだ。
products = api.query(mishima_envelope,
date = ('20200101', '20200131'),
platformname = 'Sentinel-2',
processinglevel = 'Level-2A',
cloudcoverpercentage = (0,10)
)
len(products.keys())
5
結果は、各画像データセットのメタデータとして、データセットのユニークIDをキーとしたOrderedDictで返ってくる。今回の場合は5個の画像データセットについての情報が得られた。
for key in products.keys():
print(key)
print('\tBegin position: ', products[key]['beginposition'])
print('\tCloud percentage: ', products[key]['cloudcoverpercentage'])
cd6ea033-2588-4009-93d6-4b1a067ac937
Begin position: 2020-01-13 01:30:21.024000
Cloud percentage: 9.223427
037aa2a4-a3e1-402a-bfad-83371d7ae018
Begin position: 2020-01-13 01:30:21.024000
Cloud percentage: 0.301294
4342f23c-accc-425f-a6c3-5961b8f25022
Begin position: 2020-01-13 01:30:21.024000
Cloud percentage: 0.8914899999999999
2a1e5299-e3bf-40c2-b122-e314d06fedde
Begin position: 2020-01-13 01:30:21.024000
Cloud percentage: 2.562122
da95c8f3-9305-47be-9fe4-f6517b646939
Begin position: 2020-01-03 01:30:41.024000
Cloud percentage: 4.082906
このメタデータには多くの情報が含まれていて、いろいろな数値を見ながらどの画像を使うかを検討する。たとえば上記例では画像の撮影時期と、画像に占める雲の被覆率を表示している。
目的に応じて、雲ができるだけ少ない画像を使う、とかすると思うんだけど、今回の場合は前述の事情であまり選択肢がないので、"4342f23c-accc-425f-a6c3-5961b8f25022"のデータをダウンロードすることにする。
api.downloadにデータセットのIDを指定してダウンロードする。たいてい1ギガバイトを超えるサイズなので、出力先はある程度余裕のあるディスクにしておいた方がいい。
data_path = '/PATH/TO/OUTPUT/'
api.download('4342f23c-accc-425f-a6c3-5961b8f25022', directory_path=data_path)
ダウンロードしたデータを解凍すると、以下のようなファイルが含まれていることがわかる。
./S2A_MSIL2A_20200113T013021_N0213_R074_T54SUD_20200113T042045.SAFE/
├── AUX_DATA
├── DATASTRIP
│ └── DS_EPAE_20200113T042045_S20200113T013017
│ ├── MTD_DS.xml
│ └── QI_DATA
│ ├── FORMAT_CORRECTNESS.xml
│ ├── GENERAL_QUALITY.xml
│ ├── GEOMETRIC_QUALITY.xml
│ ├── RADIOMETRIC_QUALITY.xml
│ └── SENSOR_QUALITY.xml
├── GRANULE
│ └── L2A_T54SUD_A023809_20200113T013017
│ ├── AUX_DATA
│ │ └── AUX_ECMWFT
│ ├── IMG_DATA
│ │ ├── R10m
│ │ │ ├── T54SUD_20200113T013021_AOT_10m.jp2
│ │ │ ├── T54SUD_20200113T013021_B02_10m.jp2
│ │ │ ├── T54SUD_20200113T013021_B03_10m.jp2
│ │ │ ├── T54SUD_20200113T013021_B04_10m.jp2
│ │ │ ├── T54SUD_20200113T013021_B08_10m.jp2
│ │ │ ├── T54SUD_20200113T013021_TCI_10m.jp2
│ │ │ └── T54SUD_20200113T013021_WVP_10m.jp2
│ │ ├── R20m
│ │ │ ├── T54SUD_20200113T013021_AOT_20m.jp2
│ │ │ ├── T54SUD_20200113T013021_B02_20m.jp2
│ │ │ ├── T54SUD_20200113T013021_B03_20m.jp2
...
以下で使うのは、IMG_DATAディレクトリに含まれる解像度10m（R10m）の画像データ。
Sentinel-2の計測データは以下の表（ウィキペディアより引用）のように可視域、近赤外域、短波赤外域を含む13バンドのマルチスペクトルデータである。可視光スペクトルと近赤外の解像度10mは、そこまで高解像度とは言えないけど、無料で利用できるだけでありがたい。
| センチネル-2バンド | 中心波長(μm) | 解像度(m) |
|---|---|---|
| バンド1 沿岸のエアロゾル | 0.443 | 60 |
| バンド2 青 | 0.490 | 10 |
| バンド3 緑 | 0.560 | 10 |
| バンド4 赤 | 0.665 | 10 |
| バンド5 植生のレッドエッジ | 0.705 | 20 |
| バンド6 植生のレッドエッジ | 0.740 | 20 |
| バンド7 植生のレッドエッジ | 0.783 | 20 |
| バンド8 近赤外 | 0.842 | 10 |
| バンド8A 植生のレッドエッジ | 0.865 | 20 |
| バンド9 水蒸気 | 0.945 | 60 |
| バンド10 短波長赤外-巻雲 | 1.375 | 60 |
| バンド11 短波長赤外 | 1.610 | 20 |
| バンド12 短波長赤外 | 2.190 | 20 |
ダウンロードしたデータの"R10m"ディレクトリの中にそれぞれバンド２（青色、B02）、バンド３（緑色、B03）、バンド４（赤色、B04）、バンド８（近赤外、B08）の画像が含まれるのでこれらをrasterioで読み込んでみる。
R10_path = os.path.join(data_path, 'S2A_MSIL2A_20200113T013021_N0213_R074_T54SUD_20200113T042045.SAFE/GRANULE/L2A_T54SUD_A023809_20200113T013017/IMG_DATA/R10m')
b2 = rasterio.open(os.path.join(R10_path, 'T54SUD_20200113T013021_B02_10m.jp2'))
b3 = rasterio.open(os.path.join(R10_path, 'T54SUD_20200113T013021_B03_10m.jp2'))
b4 = rasterio.open(os.path.join(R10_path, 'T54SUD_20200113T013021_B04_10m.jp2'))
赤、緑、青の数値を書き込んでフルカラーの衛星画像を合成してみる。
RGB_File = os.path.join(data_path, 'satellite_image_rgb.tiff')
with rasterio.open(RGB_File,'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(b4.read(1),1) # band-1: Red
rgb.write(b3.read(1),2) # band-2: Green
rgb.write(b2.read(1),3) # band-3: Blue
rgb.close()
得られたファイルをQGISで見るとこんな画像になった。ぎりぎり三島が入ってる。
座標参照系変換と行政区画の切り抜き
この画像にたいして形状データを適用して、三島市のエリアだけ切り抜いてみる。
ここで、画像データ、形状データの座標参照系（coordinate reference system, CRS）が重要となる。CRSが一致していないと、形状の頂点を示す座標と、平面画像のどのピクセルがどの位置を示しているのかが対応しなくなってしまう。
ダウンロードした画像データ（のどれでもいいがここでは赤色）のCRSを調べてみると、この画像はEPSG:32654で記述されていることがわかる。
print(b4.crs)
EPSG:32654
一方で国土交通省のサイトからダウンロードした三島市の形状データは前述のようにJGD2011で記述されているので、これをEPSG:32654に変換する。
mishima_proj = mishima_shape.to_crs('EPSG:32654')
mishima_proj.crs
<Projected CRS: EPSG:32654>
Name: WGS 84 / UTM zone 54N
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: World - N hemisphere - 138°E to 144°E - by country
- bounds: (138.0, 0.0, 144.0, 84.0)
Coordinate Operation:
- name: UTM zone 54N
- method: Transverse Mercator
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
「EPSG:32654に投影された三島市の形状データ」を使って、さっき作った赤緑青合成画像を切り抜いてみる。あとでまた使うので関数にしておく。
import rasterio.mask
def clip_mishima(in_file, out_file):
with rasterio.open(in_file) as src:
out_image, out_transform = rasterio.mask.mask(src,
mishima_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(out_file, "w", **out_meta) as dest:
dest.write(out_image)
RGB_Masked_File = os.path.join(data_path, 'mishima_rgb.tiff')
clip_mishima(RGB_File, RGB_Masked_File)
切り抜かれた画像をQGISで開いてみる。きれいに三島市を切り抜いた衛星画像を生成できた。
NDVIの計算と可視化
NDVI（Normalized difference vegetation index, 正規化植生指数）と呼ばれる指標を計算すると、陸生植物の活性度を衛星画像から定量できる。
植物の葉っぱは近赤外域と可視光赤色域で反射率に大きな違いがあるらしい。
クロロフィルは700nm以下の可視光を強く吸収して光合成に利用するが、もっと長い波長の光はエネルギーも小さくて利用効率が悪く、また植物を加熱して損傷させてしまうので、近赤外光は葉の細胞構造で強く反射するように進化してきたらしい。
したがって、赤色と近赤外の反射率の違いを定量することで、その衛星画像に含まれる陸生植物の存在、活性度、光合成能力を見積もることができる。
Sentinel-2のデータは赤色（バンド４）と近赤外（バンド８）のスペクトルを含むので、それらを読み出してNDVIの計算式どおりそのまま計算すればオーケー。
b4 = rasterio.open(os.path.join(R10_path, 'T54SUD_20200113T013021_B04_10m.jp2'))
b8 = rasterio.open(os.path.join(R10_path, 'T54SUD_20200113T013021_B08_10m.jp2'))
red = b4.read()
nir = b8.read()
ndvi = (nir.astype(float)-red.astype(float))/(nir+red)
各ピクセルの値がNDVI値になったグレースケール画像を出力する。
out_meta = b4.meta
out_meta.update(driver='GTiff')
out_meta.update(dtype=rasterio.float32)
NDVI_File = os.path.join(data_path, 'satellite_image_ndvi.tiff')
with rasterio.open(NDVI_File, 'w', **out_meta) as dst:
dst.write(ndvi.astype(rasterio.float32))
このNDVI画像を同じように三島市の形状で切り抜く。
NDVI_Masked_File = os.path.join(data_path, 'mishima_ndvi.tiff')
clip_mishima(NDVI_File, NDVI_Masked_File)
QGISで見てみる。値が低いほど青色に、高いほど緑色になるように調整した。
高いNDVI値は陸生植物の存在、光合成活動を示唆する。前述の赤緑青画像と比較すると、やっぱり森っぽい場所とか、楽寿園と三嶋大社（左側市街の中の２つの緑エリア）で高い傾向にあることがわかる。
rasterioで読み込んだデータはNumpyのマトリックスになってるので、具体的な数値を計算するのも簡単。
ただし注意点として、さっきの手順で三島市を切り抜いた場合、切り抜かれた外側（三島市の外側）の値がゼロに設定された画像データになってしまう。NDVIは数値範囲([-1.0, 1.0])にゼロを含むので、そのままだと三島市の外側がNDVI値0のエリアになってしまう。平均NDVIとか計算するときにやっかいなので、切り抜きの際にnodata=numpy.nanパラメータを設定して、ゼロではなくNaNにしておいた方が扱いが楽。
with rasterio.open(NDVI_File) as src:
ndvi_mishima, _ = rasterio.mask.mask(src,
mishima_proj.geometry,
nodata=np.nan, # 対象エリアの外側をNaNに設定。
crop=True)
これで三島市内の平均NDVIが簡単に計算できるようになった。
np.nanmean(ndvi_mishima)
0.5356625