森林分野における衛星データ利用事例

この章で学習すること

本章では、茨城県日立市を対象に、衛星データからNDVIを算出し、森林の状態変化を可視化する方法を学びます。

  • numpy を用いた画像解析処理演算
  • rasterio を用いたカラー合成画像の作成とマスク処理
  • rasterio を用いたNDVIの計算
  • matplotlib を用いたグラフ作成

背景

持続可能な開発目標(SDGs)の目標14は「海の豊かさを守ろう」ですが、目標15では「陸の豊かさも守ろう」が定められています。陸域生態系の保護や土地劣化の阻止などと並び、森林は重要な要素であり、例えば以下のようなターゲットが含まれています(数字はターゲット番号)。

  • 森林を含む陸域生態系の保全・回復・持続可能な利用(15.1)
  • 持続可能な森林経営、森林減少阻止、劣化した森林の回復(15.2)
  • 持続可能な森林経営のための資金の調達と開発途上国への資源の動員(15.b)

適切な森林経営のためには現状の把握と監視が必要であり、大規模な森林を把握・監視するためには衛星が有効な手段です。

光学衛星で森林を監視する手段として、今回は植生指標に着目します。森林を含む全ての植物は近赤外線を反射し、可視光線のうち赤色の波長を吸収する特性があります。この近赤外と赤色の分光反射率を組み合わせることで植生指標を算出し、森林など植生の状態把握をすることが可能です。

この章では、一般的に用いられている植生指標としてNDVI(Normalized Difference Vegetation Index)を衛星データから計算し、それを時系列に可視化します。植生指数は、森林から他の土地被覆に変化した場合だけでなく季節によっても変化します。そのため、時系列で森林を解析するためには、まず、対象となる地域における季節と植生指標との関係を確認することがとても重要です。ここでは、森林の植生指標の季節変化を可視化する方法を学んでいきます。

準備

In [3]:
# google driveをマウント
# Tellusの開発環境やローカル環境を使用する想定であれば、このセルは削除可です。
from google.colab import drive 
drive = drive.mount('/content/drive')
Drive already mounted at /content/drive; to attempt to forcibly remount, call drive.mount("/content/drive", force_remount=True).

必要なライブラリのインストール

In [4]:
!pip install collection
!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 pandas
!pip install geopandas

print("done")
Collecting collection
  Downloading https://files.pythonhosted.org/packages/81/a7/12577601fd60036732cd5b078ed3aa9a2f888169ea70d84997b3607e63b4/collection-0.1.6.tar.gz
Building wheels for collected packages: collection
  Building wheel for collection (setup.py) ... done
  Created wheel for collection: filename=collection-0.1.6-cp36-none-any.whl size=5116 sha256=34611d31e518fea1ce0f1451e8c283698f24f3c076af569b89c1b3f280e58628
  Stored in directory: /root/.cache/pip/wheels/9e/f2/2b/a611b0dc83b770763e7962500ef158c09dc8161d3fce6e73de
Successfully built collection
Installing collected packages: collection
Successfully installed collection-0.1.6
Collecting rasterio
  Downloading https://files.pythonhosted.org/packages/c0/a8/63d45bb74c17c60e607b4beae77d68ad4c9ea6dff788534ce8c835d1d2f1/rasterio-1.2.0-cp36-cp36m-manylinux1_x86_64.whl (19.1MB)
     |████████████████████████████████| 19.1MB 50.0MB/s 
Requirement already satisfied: certifi in /usr/local/lib/python3.6/dist-packages (from rasterio) (2020.12.5)
Collecting cligj>=0.5
  Downloading https://files.pythonhosted.org/packages/42/1e/947eadf10d6804bf276eb8a038bd5307996dceaaa41cfd21b7a15ec62f5d/cligj-0.7.1-py3-none-any.whl
Collecting affine
  Downloading https://files.pythonhosted.org/packages/ac/a6/1a39a1ede71210e3ddaf623982b06ecfc5c5c03741ae659073159184cd3e/affine-2.3.0-py2.py3-none-any.whl
Collecting snuggs>=1.4.1
  Downloading https://files.pythonhosted.org/packages/cc/0e/d27d6e806d6c0d1a2cfdc5d1f088e42339a0a54a09c3343f7f81ec8947ea/snuggs-1.4.7-py3-none-any.whl
Requirement already satisfied: numpy in /usr/local/lib/python3.6/dist-packages (from rasterio) (1.19.5)
Requirement already satisfied: click<8,>=4.0 in /usr/local/lib/python3.6/dist-packages (from rasterio) (7.1.2)
Requirement already satisfied: attrs in /usr/local/lib/python3.6/dist-packages (from rasterio) (20.3.0)
Collecting click-plugins
  Downloading https://files.pythonhosted.org/packages/e9/da/824b92d9942f4e472702488857914bdd50f73021efea15b4cad9aca8ecef/click_plugins-1.1.1-py2.py3-none-any.whl
Requirement already satisfied: pyparsing>=2.1.6 in /usr/local/lib/python3.6/dist-packages (from snuggs>=1.4.1->rasterio) (2.4.7)
Installing collected packages: cligj, affine, snuggs, click-plugins, rasterio
Successfully installed affine-2.3.0 click-plugins-1.1.1 cligj-0.7.1 rasterio-1.2.0 snuggs-1.4.7
Reading package lists... Done
Building dependency tree       
Reading state information... Done
gdal-bin is already the newest version (2.2.3+dfsg-2).
python-gdal is already the newest version (2.2.3+dfsg-2).
python3-gdal is already the newest version (2.2.3+dfsg-2).
0 upgraded, 0 newly installed, 0 to remove and 13 not upgraded.
Reading package lists... Done
Building dependency tree       
Reading state information... Done
The following additional packages will be installed:
  libspatialindex-c4v5 libspatialindex-dev libspatialindex4v5
  python3-pkg-resources
Suggested packages:
  python3-setuptools
The following NEW packages will be installed:
  libspatialindex-c4v5 libspatialindex-dev libspatialindex4v5
  python3-pkg-resources python3-rtree
0 upgraded, 5 newly installed, 0 to remove and 13 not upgraded.
Need to get 671 kB of archives.
After this operation, 3,948 kB of additional disk space will be used.
Get:1 http://archive.ubuntu.com/ubuntu bionic/universe amd64 libspatialindex4v5 amd64 1.8.5-5 [219 kB]
Get:2 http://archive.ubuntu.com/ubuntu bionic/universe amd64 libspatialindex-c4v5 amd64 1.8.5-5 [51.7 kB]
Get:3 http://archive.ubuntu.com/ubuntu bionic/main amd64 python3-pkg-resources all 39.0.1-2 [98.8 kB]
Get:4 http://archive.ubuntu.com/ubuntu bionic/universe amd64 libspatialindex-dev amd64 1.8.5-5 [285 kB]
Get:5 http://archive.ubuntu.com/ubuntu bionic/universe amd64 python3-rtree all 0.8.3+ds-1 [16.9 kB]
Fetched 671 kB in 1s (655 kB/s)
Selecting previously unselected package libspatialindex4v5:amd64.
(Reading database ... 146374 files and directories currently installed.)
Preparing to unpack .../libspatialindex4v5_1.8.5-5_amd64.deb ...
Unpacking libspatialindex4v5:amd64 (1.8.5-5) ...
Selecting previously unselected package libspatialindex-c4v5:amd64.
Preparing to unpack .../libspatialindex-c4v5_1.8.5-5_amd64.deb ...
Unpacking libspatialindex-c4v5:amd64 (1.8.5-5) ...
Selecting previously unselected package python3-pkg-resources.
Preparing to unpack .../python3-pkg-resources_39.0.1-2_all.deb ...
Unpacking python3-pkg-resources (39.0.1-2) ...
Selecting previously unselected package libspatialindex-dev:amd64.
Preparing to unpack .../libspatialindex-dev_1.8.5-5_amd64.deb ...
Unpacking libspatialindex-dev:amd64 (1.8.5-5) ...
Selecting previously unselected package python3-rtree.
Preparing to unpack .../python3-rtree_0.8.3+ds-1_all.deb ...
Unpacking python3-rtree (0.8.3+ds-1) ...
Setting up libspatialindex4v5:amd64 (1.8.5-5) ...
Setting up python3-pkg-resources (39.0.1-2) ...
Setting up libspatialindex-c4v5:amd64 (1.8.5-5) ...
Setting up libspatialindex-dev:amd64 (1.8.5-5) ...
Setting up python3-rtree (0.8.3+ds-1) ...
Processing triggers for libc-bin (2.27-3ubuntu1.3) ...
/sbin/ldconfig.real: /usr/local/lib/python3.6/dist-packages/ideep4py/lib/libmkldnn.so.0 is not a symbolic link

Collecting git+git://github.com/geopandas/geopandas.git
  Cloning git://github.com/geopandas/geopandas.git to /tmp/pip-req-build-e8p4ne0u
  Running command git clone -q git://github.com/geopandas/geopandas.git /tmp/pip-req-build-e8p4ne0u
Requirement already satisfied: pandas>=0.24.0 in /usr/local/lib/python3.6/dist-packages (from geopandas==0.8.0+88.gb3db3be) (1.1.5)
Requirement already satisfied: shapely>=1.6 in /usr/local/lib/python3.6/dist-packages (from geopandas==0.8.0+88.gb3db3be) (1.7.1)
Collecting fiona>=1.8
  Downloading https://files.pythonhosted.org/packages/37/94/4910fd55246c1d963727b03885ead6ef1cd3748a465f7b0239ab25dfc9a3/Fiona-1.8.18-cp36-cp36m-manylinux1_x86_64.whl (14.8MB)
     |████████████████████████████████| 14.8MB 333kB/s 
Collecting pyproj>=2.2.0
  Downloading https://files.pythonhosted.org/packages/e4/ab/280e80a67cfc109d15428c0ec56391fc03a65857b7727cf4e6e6f99a4204/pyproj-3.0.0.post1-cp36-cp36m-manylinux2010_x86_64.whl (6.4MB)
     |████████████████████████████████| 6.5MB 46.2MB/s 
Requirement already satisfied: numpy>=1.15.4 in /usr/local/lib/python3.6/dist-packages (from pandas>=0.24.0->geopandas==0.8.0+88.gb3db3be) (1.19.5)
Requirement already satisfied: pytz>=2017.2 in /usr/local/lib/python3.6/dist-packages (from pandas>=0.24.0->geopandas==0.8.0+88.gb3db3be) (2018.9)
Requirement already satisfied: python-dateutil>=2.7.3 in /usr/local/lib/python3.6/dist-packages (from pandas>=0.24.0->geopandas==0.8.0+88.gb3db3be) (2.8.1)
Requirement already satisfied: cligj>=0.5 in /usr/local/lib/python3.6/dist-packages (from fiona>=1.8->geopandas==0.8.0+88.gb3db3be) (0.7.1)
Requirement already satisfied: certifi in /usr/local/lib/python3.6/dist-packages (from fiona>=1.8->geopandas==0.8.0+88.gb3db3be) (2020.12.5)
Collecting munch
  Downloading https://files.pythonhosted.org/packages/cc/ab/85d8da5c9a45e072301beb37ad7f833cd344e04c817d97e0cc75681d248f/munch-2.5.0-py2.py3-none-any.whl
Requirement already satisfied: attrs>=17 in /usr/local/lib/python3.6/dist-packages (from fiona>=1.8->geopandas==0.8.0+88.gb3db3be) (20.3.0)
Requirement already satisfied: six>=1.7 in /usr/local/lib/python3.6/dist-packages (from fiona>=1.8->geopandas==0.8.0+88.gb3db3be) (1.15.0)
Requirement already satisfied: click<8,>=4.0 in /usr/local/lib/python3.6/dist-packages (from fiona>=1.8->geopandas==0.8.0+88.gb3db3be) (7.1.2)
Requirement already satisfied: click-plugins>=1.0 in /usr/local/lib/python3.6/dist-packages (from fiona>=1.8->geopandas==0.8.0+88.gb3db3be) (1.1.1)
Building wheels for collected packages: geopandas
  Building wheel for geopandas (setup.py) ... done
  Created wheel for geopandas: filename=geopandas-0.8.0+88.gb3db3be-py2.py3-none-any.whl size=976378 sha256=20cdd7ab10c2365fc2fde1c0c04d4957caa90e363e8ad0984422fa39c088c63d
  Stored in directory: /tmp/pip-ephem-wheel-cache-wlbpbuwc/wheels/91/24/71/376c9c67192694168352afcccc2d264248f7e2cc6192997186
Successfully built geopandas
Installing collected packages: munch, fiona, pyproj, geopandas
Successfully installed fiona-1.8.18 geopandas-0.8.0+88.gb3db3be munch-2.5.0 pyproj-3.0.0.post1
Requirement already satisfied: descartes in /usr/local/lib/python3.6/dist-packages (1.1.0)
Requirement already satisfied: matplotlib in /usr/local/lib/python3.6/dist-packages (from descartes) (3.2.2)
Requirement already satisfied: numpy>=1.11 in /usr/local/lib/python3.6/dist-packages (from matplotlib->descartes) (1.19.5)
Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib->descartes) (2.4.7)
Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib->descartes) (2.8.1)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.6/dist-packages (from matplotlib->descartes) (0.10.0)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib->descartes) (1.3.1)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.6/dist-packages (from python-dateutil>=2.1->matplotlib->descartes) (1.15.0)
Requirement already satisfied: shapely in /usr/local/lib/python3.6/dist-packages (1.7.1)
Requirement already satisfied: six in /usr/local/lib/python3.6/dist-packages (1.15.0)
Requirement already satisfied: pyproj in /usr/local/lib/python3.6/dist-packages (3.0.0.post1)
Requirement already satisfied: certifi in /usr/local/lib/python3.6/dist-packages (from pyproj) (2020.12.5)
Requirement already satisfied: descartes in /usr/local/lib/python3.6/dist-packages (1.1.0)
Requirement already satisfied: matplotlib in /usr/local/lib/python3.6/dist-packages (from descartes) (3.2.2)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.6/dist-packages (from matplotlib->descartes) (0.10.0)
Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib->descartes) (2.8.1)
Requirement already satisfied: numpy>=1.11 in /usr/local/lib/python3.6/dist-packages (from matplotlib->descartes) (1.19.5)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib->descartes) (1.3.1)
Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib->descartes) (2.4.7)
Requirement already satisfied: six in /usr/local/lib/python3.6/dist-packages (from cycler>=0.10->matplotlib->descartes) (1.15.0)
Requirement already satisfied: pandas in /usr/local/lib/python3.6/dist-packages (1.1.5)
Requirement already satisfied: pytz>=2017.2 in /usr/local/lib/python3.6/dist-packages (from pandas) (2018.9)
Requirement already satisfied: numpy>=1.15.4 in /usr/local/lib/python3.6/dist-packages (from pandas) (1.19.5)
Requirement already satisfied: python-dateutil>=2.7.3 in /usr/local/lib/python3.6/dist-packages (from pandas) (2.8.1)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.6/dist-packages (from python-dateutil>=2.7.3->pandas) (1.15.0)
Requirement already satisfied: geopandas in /usr/local/lib/python3.6/dist-packages (0.8.0+88.gb3db3be)
Requirement already satisfied: pyproj>=2.2.0 in /usr/local/lib/python3.6/dist-packages (from geopandas) (3.0.0.post1)
Requirement already satisfied: shapely>=1.6 in /usr/local/lib/python3.6/dist-packages (from geopandas) (1.7.1)
Requirement already satisfied: fiona>=1.8 in /usr/local/lib/python3.6/dist-packages (from geopandas) (1.8.18)
Requirement already satisfied: pandas>=0.24.0 in /usr/local/lib/python3.6/dist-packages (from geopandas) (1.1.5)
Requirement already satisfied: certifi in /usr/local/lib/python3.6/dist-packages (from pyproj>=2.2.0->geopandas) (2020.12.5)
Requirement already satisfied: cligj>=0.5 in /usr/local/lib/python3.6/dist-packages (from fiona>=1.8->geopandas) (0.7.1)
Requirement already satisfied: six>=1.7 in /usr/local/lib/python3.6/dist-packages (from fiona>=1.8->geopandas) (1.15.0)
Requirement already satisfied: attrs>=17 in /usr/local/lib/python3.6/dist-packages (from fiona>=1.8->geopandas) (20.3.0)
Requirement already satisfied: click-plugins>=1.0 in /usr/local/lib/python3.6/dist-packages (from fiona>=1.8->geopandas) (1.1.1)
Requirement already satisfied: click<8,>=4.0 in /usr/local/lib/python3.6/dist-packages (from fiona>=1.8->geopandas) (7.1.2)
Requirement already satisfied: munch in /usr/local/lib/python3.6/dist-packages (from fiona>=1.8->geopandas) (2.5.0)
Requirement already satisfied: pytz>=2017.2 in /usr/local/lib/python3.6/dist-packages (from pandas>=0.24.0->geopandas) (2018.9)
Requirement already satisfied: numpy>=1.15.4 in /usr/local/lib/python3.6/dist-packages (from pandas>=0.24.0->geopandas) (1.19.5)
Requirement already satisfied: python-dateutil>=2.7.3 in /usr/local/lib/python3.6/dist-packages (from pandas>=0.24.0->geopandas) (2.8.1)
done

ライブラリのインポート

「基礎編」でご紹介したGDALライブラリは、地理空間データを扱う上で基礎となる重要なライブラリですが、GISの概念に慣れていないと理解しにくいパラメータがたくさん登場するため、衛星データを触り始めたばかりの頃は複雑に感じることもあるかと思います。

この教材では、地理空間データの処理がより便利になるPythonのライブラリをいくつかご紹介します。

  • numPy

    numpyは、数値計算を効率的に行うための拡張モジュールです。基本的な計算はpythonだけでも出来ますが、numpyを使うとその計算が高速になったり、より楽になったりします。多次元配列を必要とする処理において、特にその機能を発揮するライブラリです。

  • rasterio

    rasterioは、オンラインマッププロバイダーのMapboxが中心となって開発しているOSSのライブラリで、衛星データでは最も一般的なGeoTiffフォーマットのデータを扱う際に、シンプルな書き方ができるのが特徴のライブラリです。

  • matplotlib

    matplotlibは、Pythonでグラフや図を描画、作成するための包括的なライブラリです。折れ線グラフやヒストグラムといった2次元グラフだけなく3次元のグラフも描画することが可能です。

In [5]:
#必要ライブラリのインポート
import os
import numpy as np
import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
import rasterio as rio
import rasterio.mask
import fiona
import folium
import zipfile
import glob

from mpl_toolkits.axes_grid1 import make_axes_locatable
from shapely.geometry import MultiPolygon, Polygon
from rasterio import plot
from rasterio.plot import show
from rasterio.plot import plotting_extent
from rasterio.mask import mask
from PIL import Image
from osgeo import gdal

print("done")
done

衛星データの取得

この章では、ヨーロッパ宇宙機関(ESA)が運用するSentinel-2衛星が取得するデータを使用します。ESAのオープンデータポリシーに基づき、誰でも無償で観測データを利用することが可能です。

Sentinel-2について


Sentinel-2は陸域観測を主目的とした光学衛星であり、可視光、近赤外を含む12種類の異なる波長帯(バンド)で観測を行っています。

光学衛星は雲があると地表を観測できませんが、Sentinel-2は同じ設計の衛星2機を運用することで、5日に1回の観測を可能とし、晴天時のデータを入手できる可能性を高めています。

以下の2種類のプロダクトが提供されています。今回はLevel-1Cのプロダクトを使用します。

  • Level-1C: 建物の倒れこみ等、幾何学的な歪みを補正したオルソ補正済みプロダクト
  • Level-2A: 1Cの補正に加え、大気の吸収・散乱の影響を軽減し、地表面の反射率に変換した大気補正済みプロダクト

衛星の緒言についてはリモート・センシング技術センターのサイトを参照してください。

データの検索とダウンロード

データを取得するに際し、ESAのサイトCopernicus Open Access Hubでユーザ登録を行いましょう。ユーザー登録の詳細はこちらを参照してください。

Copernicus Open Access Hubの検索画面(図1)で、地図上で対象領域を指定した上で、以下の条件を入力して検索しましょう。今回は茨城県の日立市を対象としてデータを検索します。

  • 「Sensing period」(観測期間)に下記の期間を入力。

2017年 → 2017年1月1日から12月31日

2018年 → 2018年1月1日から12月31

  • 「Mission: Sentinel-2」のチェックボックスにチェック。
  • 「Product type」はプルダウンで「S2MSI1C」を選択。

Copernicus Open Access Hub.png

図1. Copernicus Open Access Hubの検索画面

検索したデータはサムネイル画像で雲の様子などが確認できるので、雲が無い良好なデータを見つけましょう。

今回は、2017年5月から2018年12月の7シーンのデータを使って年間の森林状態の変化を確認し、2017年および2018年の5月のデータを使って経年変化を確認します。

使用するデータ(全7シーン)の日付は以下の通りです。

  • 2017-05-08
  • 2017-10-27
  • 2017-12-26
  • 2018-04-20
  • 2018-05-20
  • 2018-10-02
  • 2018-12-19

Coperniqus Open Access Hubからデータをダウンロードすると、zip形式で圧縮されているので展開します。

In [6]:
# Tellusの開発環境やローカル環境を使用する想定であれば、適切なパスに書き換えてください。

a = glob.glob('/content/drive/MyDrive/s2/*.zip')

ziplis = []

for f in a:
  print(f[26:86])
  ziplis.append(f[26:86])

for file_title in ziplis:
  print("start unzip:"+ file_title)
  file_name = file_title +".zip"
  file_dir = "/content/drive/MyDrive/s2/"
  file_pass = os.path.join(file_dir + file_name)
  print(file_pass)
  with zipfile.ZipFile(file_pass) as zf:
    zf.extractall()
  print("done")
S2A_MSIL1C_20180420T011651_N0206_R031_T54SVF_20180420T030851
S2A_MSIL1C_20180520T011701_N0206_R031_T54SVF_20180520T030828
S2B_MSIL1C_20181002T011649_N0206_R031_T54SVF_20181002T030308
S2A_MSIL1C_20181219T013041_N0207_R074_T54SVF_20181219T031650
S2A_MSIL1C_20170508T012701_N0205_R074_T54SVF_20170508T013110
S2B_MSIL1C_20171226T012039_N0206_R031_T54SVF_20171226T024938
S2B_MSIL1C_20171027T011719_N0206_R031_T54SVF_20171027T030539
start unzip:S2A_MSIL1C_20180420T011651_N0206_R031_T54SVF_20180420T030851
/content/drive/MyDrive/s2/S2A_MSIL1C_20180420T011651_N0206_R031_T54SVF_20180420T030851.zip
done
start unzip:S2A_MSIL1C_20180520T011701_N0206_R031_T54SVF_20180520T030828
/content/drive/MyDrive/s2/S2A_MSIL1C_20180520T011701_N0206_R031_T54SVF_20180520T030828.zip
done
start unzip:S2B_MSIL1C_20181002T011649_N0206_R031_T54SVF_20181002T030308
/content/drive/MyDrive/s2/S2B_MSIL1C_20181002T011649_N0206_R031_T54SVF_20181002T030308.zip
done
start unzip:S2A_MSIL1C_20181219T013041_N0207_R074_T54SVF_20181219T031650
/content/drive/MyDrive/s2/S2A_MSIL1C_20181219T013041_N0207_R074_T54SVF_20181219T031650.zip
done
start unzip:S2A_MSIL1C_20170508T012701_N0205_R074_T54SVF_20170508T013110
/content/drive/MyDrive/s2/S2A_MSIL1C_20170508T012701_N0205_R074_T54SVF_20170508T013110.zip
done
start unzip:S2B_MSIL1C_20171226T012039_N0206_R031_T54SVF_20171226T024938
/content/drive/MyDrive/s2/S2B_MSIL1C_20171226T012039_N0206_R031_T54SVF_20171226T024938.zip
done
start unzip:S2B_MSIL1C_20171027T011719_N0206_R031_T54SVF_20171027T030539
/content/drive/MyDrive/s2/S2B_MSIL1C_20171027T011719_N0206_R031_T54SVF_20171027T030539.zip
done

解析結果を格納するための作業用ディレクトリを以下のように準備します。このディレクトリには、workディレクトリの下にTIF、NDVIおよびNDVIにマスク処理をしたデータを保存するためのディレクトリを作成します。

In [7]:
#作業用ディレクトリ(フォルダ)の作成

os.mkdir("work")

path = '/content/'
work_path = '/content/work'

RGB_dir = '/content/work/RGB_TIF'
os.mkdir(RGB_dir)

NDVI_dir = '/content/work/NDVI'
os.mkdir(NDVI_dir)

NDVI_mask_dir = '/content/work/NDVI_mask'
os.mkdir(NDVI_mask_dir)

コラム:API経由のLTAデータ取得


Sentinel-2のデータはAPI経由でも取得することが可能で、宙畑の記事で取得方法が紹介されています。

ただし、2019年8月以前に観測されたSentinel-2のLevel 2Aデータは全てLTA(Long Term Archives)に保管されており、それらデータをAPIで取得したい場合はESAにリクエストを送る必要があります。リクエストが承認されると、24時間以内にプロダクトがダウンロード可能となり、その時点から3日間利用可能となります。

2019年9月以降に観測されたデータであれば、常時取得が可能なのでリクエストを送らずにAPIで取得できます。

今回はAPIを使わずにデータを取得しましたが、APIでのデータ取得もぜひお試しください。

※2020年12月現在、複数のLTAデータを同時にリクエストすることができない状態となっています。リクエストは30分ごとに1回しか送信できないため、多くのシーンを取得したい場合は非常に時間がかかります。

カラー合成画像の作成

今回の解析では茨城県日立市の森林を対象としますので、それ以外の領域は不要です。対象領域のみを解析するため、国土地理院が公開しているポリゴンデータ(Shapefile)を使って不要な領域をマスク(隠す)します。

日立市のポリゴン取得

日立市のポリゴンデータは、国土地理院が公開している国土数値情報のうち行政区域データをダウンロードしましょう。ダウンロードは県単位になるため、茨城県から日立市のポリゴンを以下のように抜き出します。

抜き出したら、ポリゴンの座標系をSentinel-2のデータと同じUTM zone 54N (EPSG:32654) に変換します。座標系やEPSGについては、基礎分野の章や空間情報クラブのサイトを参照して下さい。

In [8]:
#ダウンロードしたデータを解凍して、ディレクトリを指定
shape_path = "/content/drive/MyDrive/N03-20200101_08_GML/"

# 茨城県のshpの読み込み
in_shape = gpd.read_file(os.path.join(shape_path + "N03-20_08_200101.shp"),encoding="shift-jis")

# 日立市でソート
shape_srt = in_shape[in_shape["N03_004"].isin(["日立市"])]
shape_srt = shape_srt.drop(columns=["N03_002","N03_003"])

#画像に合わせて投影変換
out_file = os.path.join(shape_path + "re_N03-20_08_200101.shp")
re_shape = shape_srt.to_crs({"init": "epsg:32654"})

#出力
re_shape.to_file(driver="ESRI Shapefile",filename=out_file)
print("done")

#GeoDataFrameを描画
f,ax = plt.subplots(1, figsize=(6,6))
ax = re_shape.plot(axes=ax)
plt.show
/usr/local/lib/python3.6/dist-packages/pyproj/crs/crs.py:53: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
  return _prepare_from_string(" ".join(pjargs))
done
/usr/local/lib/python3.6/dist-packages/geopandas/plotting.py:602: FutureWarning: 'axes' is deprecated, please use 'ax' instead (for consistency with pandas)
  FutureWarning,
Out[8]:
<function matplotlib.pyplot.show>

カラー合成とマスク処理

Sentinel-2は異なる12の波長帯(バンド)で観測していると説明しましたが、ダウンロードしたデータはバンド毎にjpeg2の形式で格納されています。

まずは、基礎分野で学習した「トゥルーカラー(True Color)」の組み合わせでカラー合成し、日立市以外の領域をマスクしてみましょう。

ここではrasterioで、バンド2(青)・バンド3(緑)・バンド4(赤)を組み合わせて合成し、ポリゴンデータを使って日立市以外をマスクします。

In [9]:
for file_title in ziplis:
  print("Start make RGB TIF image " +'<'+ file_title +'>')

  path_A = str(file_title) + '.SAFE/GRANULE/'
  f1 = os.listdir(path_A)
  path_B = str(file_title) + '.SAFE/GRANULE/' + str(f1[0])
  f2 = os.listdir(path_B)
  path_C = str(file_title) + '.SAFE/GRANULE/' + str(f1[0]) + '/IMG_DATA/'
  f3 = os.listdir(path_C)

  b2 = rio.open(str(file_title) + '.SAFE/GRANULE/' + str(f1[0]) +'/IMG_DATA/' +str(f3[0][0:23] +'B02.jp2'))
  b3 = rio.open(str(file_title) + '.SAFE/GRANULE/' + str(f1[0]) +'/IMG_DATA/' +str(f3[0][0:23] +'B03.jp2'))
  b4 = rio.open(str(file_title) + '.SAFE/GRANULE/' + str(f1[0]) +'/IMG_DATA/' +str(f3[0][0:23] +'B04.jp2'))
  
  RGB_path = os.path.join(RGB_dir,'sentinel-2_'+str(f3[0][7:15])+'_RGB.tif')
  

  RGB_colar = rio.open(RGB_path,'w',driver='Gtiff',
                       width=b4.width, height=b4.height,
                       count=3,
                       crs=b4.crs,
                       transform=b4.transform,
                       dtype=rasterio.uint16
                       )
  RGB_colar.write(b2.read(1),3) 
  RGB_colar.write(b3.read(1),2) 
  RGB_colar.write(b4.read(1),1) 
  RGB_colar.close()

  print("---masking---")

  with fiona.open(out_file, "r") as mask:
    masks = [feature["geometry"] for feature in mask]

  with rasterio.open(RGB_path) as src:
    out_image, out_transform = rasterio.mask.mask(src, masks, crop=True)
    out_meta = src.meta

  out_meta.update({"driver": "GTiff",
                 "height": out_image.shape[1],
                 "width": out_image.shape[2],
                 "transform": out_transform})

  with rasterio.open(RGB_path, "w", **out_meta) as dest:
    dest.write(out_image)
  
  #画像表示のため8bit形式で書き出し。
  scale = '-scale 0 255 0 15'
  options_list = ['-ot Byte','-of Gtiff',scale]
  options_string = " ".join(options_list)

  gdal.Translate(os.path.join(RGB_dir + "/" + 'sentinel-2_'+str(f3[0][7:15])+'.tif'),os.path.join(RGB_dir + "/" + 'sentinel-2_'+str(f3[0][7:15])+'_RGB.tif'),options = options_string)

  print("Done")
Start make RGB TIF image <S2A_MSIL1C_20180420T011651_N0206_R031_T54SVF_20180420T030851>
---masking---
Done
Start make RGB TIF image <S2A_MSIL1C_20180520T011701_N0206_R031_T54SVF_20180520T030828>
---masking---
Done
Start make RGB TIF image <S2B_MSIL1C_20181002T011649_N0206_R031_T54SVF_20181002T030308>
---masking---
Done
Start make RGB TIF image <S2A_MSIL1C_20181219T013041_N0207_R074_T54SVF_20181219T031650>
---masking---
Done
Start make RGB TIF image <S2A_MSIL1C_20170508T012701_N0205_R074_T54SVF_20170508T013110>
---masking---
Done
Start make RGB TIF image <S2B_MSIL1C_20171226T012039_N0206_R031_T54SVF_20171226T024938>
---masking---
Done
Start make RGB TIF image <S2B_MSIL1C_20171027T011719_N0206_R031_T54SVF_20171027T030539>
---masking---
Done
In [10]:
plt.figure(figsize=(18,9))

RGB_2018_spring = rio.open(RGB_dir+'/sentinel-2_20180420.tif')
RGB_2018_summer = rio.open(RGB_dir+'/sentinel-2_20180520.tif')
RGB_2018_autumn = rio.open(RGB_dir+'/sentinel-2_20181002.tif')
RGB_2018_winter = rio.open(RGB_dir+'/sentinel-2_20181219.tif')

ax1 = plt.subplot(1,4,1)
ax1.set_title("RGB_2018_spring")
ax2 = plt.subplot(1,4,2)
ax2.set_title("RGB_2018_summer")
ax3 = plt.subplot(1,4,3)
ax3.set_title("RGB_2018_autumn")
ax4 = plt.subplot(1,4,4)
ax4.set_title("RGB_2018_winter")

show(RGB_2018_spring.read([1,2,3]),ax=ax1)
show(RGB_2018_summer.read([1,2,3]),ax=ax2)
show(RGB_2018_autumn.read([1,2,3]),ax=ax3)
show(RGB_2018_winter.read([1,2,3]),ax=ax4)
Out[10]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f1e02fe62e8>

2018年の4シーン(4月、5月、10月、12月)のカラー合成画像です。見た目にはほとんど違いは分かりません。

森林の状態変化の可視化

森林地域のポリゴン取得

国土地理院の国土数値情報には、土地利用基本計画に基づき指定された森林地域のポリゴンも含まれています。このサイトから茨城県の森林地域のポリゴンをダウンロードし、日立市のポリゴンと同様に座標系をUTM zone 54N (EPSG:32654) に変換しましょう。

In [11]:
#森林マップポリゴン

#ダウンロードしたデータを解凍して、ディレクトリを指定
mori_shape_path = "/content/drive/MyDrive/A13-15_08_GML/"

# 森林マップポリゴンの読込(森林域のポリゴンは下2桁が07)
mori_shape = gpd.read_file(os.path.join(mori_shape_path + "a001080020160207.shp"),encoding="shift-jis")

#画像に合わせて投影変換
out_file_mori = os.path.join(mori_shape_path + "re_a001080020160207.shp")
re_mori_shape = mori_shape.to_crs({"init": "epsg:32654"})

#出力
re_mori_shape.to_file(driver="ESRI Shapefile",filename=out_file_mori)
print("done")

#GeoDataFrameを描画
f,ax = plt.subplots(1, figsize=(6,6))
ax = re_mori_shape.plot(axes=ax)
plt.show
/usr/local/lib/python3.6/dist-packages/pyproj/crs/crs.py:53: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
  return _prepare_from_string(" ".join(pjargs))
done
/usr/local/lib/python3.6/dist-packages/geopandas/plotting.py:602: FutureWarning: 'axes' is deprecated, please use 'ax' instead (for consistency with pandas)
  FutureWarning,
Out[11]:
<function matplotlib.pyplot.show>

茨城県の森林地域のポリゴンです。青い部分が森林になります。

正規化植生指数(NDVI)

植生の活性度を見るためには、正規化植生指数(Normalized Difference Vegetation Index: NDVI)という指数を使うことが有効です。NDVIは、近赤外線(NIR)と赤色(Red)の波長帯のデータを使って、以下の計算式で求めることができます。

$ NDVI=(NIR-Red)/(NIR+Red)$

</br>

Sentinel-2の場合、バンド4で赤色、バンド8で近赤外線を使っているので、以下の式でNDVIを算出します。

$NDVI=(Band 8-Band 4)/(Band 8+Band 4)$

</br>

なお、NDVIについては宙畑の記事でも紹介されているので参考にしてください。

日立市におけるNDVIを計算します。 Sentinel-2で雲の無いデータを取得できた4時期(4月、5月、10月、12月)全てでNDVIを計算します。

In [12]:
NDVI_dir = '/content/work/NDVI'

for file_title in ziplis:
  print("Start make NDVI " +'<'+ file_title +'>')

  path_A = str(file_title) + '.SAFE/GRANULE/'
  f1 = os.listdir(path_A)

  path_B = str(file_title) + '.SAFE/GRANULE/' + str(f1[0])
  f2 = os.listdir(path_B)

  path_C = str(file_title) + '.SAFE/GRANULE/' + str(f1[0]) + '/IMG_DATA/'
  f3 = os.listdir(path_C)

  b4 = rio.open(str(file_title) + '.SAFE/GRANULE/' + str(f1[0]) +'/IMG_DATA/' +str(f3[0][0:23] +'B04.jp2'))
  b8 = rio.open(str(file_title) + '.SAFE/GRANULE/' + str(f1[0]) +'/IMG_DATA/' +str(f3[0][0:23] +'B08.jp2'))
  
  red = b4.read()
  nir = b8.read()
  
  ndvi = np.where(((red != 0) & (nir != 0)),(nir.astype(float)-red.astype(float))/(nir.astype(float)+red.astype(float)), 0)

  out_meta = b4.meta
  out_meta.update(driver='GTiff')
  out_meta.update(dtype=rasterio.float32)

  NDVI_path = os.path.join(NDVI_dir,'sentinel-2_'+str(f3[0][7:15])+ '_ndvi.tif')

  with rio.open(NDVI_path, 'w', **out_meta) as dst:
    dst.write(ndvi.astype(rasterio.float32))

  print("---masking---")

  with fiona.open(out_file, "r") as mask:
    masks = [feature["geometry"] for feature in mask]

  with rio.open(NDVI_path) as src2:
    out_image_ndvi, out_transform = rasterio.mask.mask(src2, masks, crop=True)
    out_meta = src2.meta

  out_meta.update({"driver": "GTiff",
                 "height": out_image.shape[1],
                 "width": out_image.shape[2],
                 "transform": out_transform})

  with rasterio.open(NDVI_path, "w", **out_meta) as dst2:
    dst2.write(out_image_ndvi)

  print("Done")
Start make NDVI <S2A_MSIL1C_20180420T011651_N0206_R031_T54SVF_20180420T030851>
/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:21: RuntimeWarning: invalid value encountered in true_divide
---masking---
Done
Start make NDVI <S2A_MSIL1C_20180520T011701_N0206_R031_T54SVF_20180520T030828>
---masking---
Done
Start make NDVI <S2B_MSIL1C_20181002T011649_N0206_R031_T54SVF_20181002T030308>
---masking---
Done
Start make NDVI <S2A_MSIL1C_20181219T013041_N0207_R074_T54SVF_20181219T031650>
---masking---
Done
Start make NDVI <S2A_MSIL1C_20170508T012701_N0205_R074_T54SVF_20170508T013110>
---masking---
Done
Start make NDVI <S2B_MSIL1C_20171226T012039_N0206_R031_T54SVF_20171226T024938>
---masking---
Done
Start make NDVI <S2B_MSIL1C_20171027T011719_N0206_R031_T54SVF_20171027T030539>
---masking---
Done

次に画像を表示します。

In [13]:
fig,(ax1,ax2,ax3,ax4)=plt.subplots(1,4, figsize=(18,9))
plt.subplots_adjust(wspace=0.4)

NDVI_2018_spring = rio.open(NDVI_dir+'/sentinel-2_20180420_ndvi.tif')
NDVI_2018_summer = rio.open(NDVI_dir+'/sentinel-2_20180520_ndvi.tif')
NDVI_2018_autumn = rio.open(NDVI_dir+'/sentinel-2_20181002_ndvi.tif')
NDVI_2018_winter = rio.open(NDVI_dir+'/sentinel-2_20181219_ndvi.tif')

NDVI_2018_spring = NDVI_2018_spring.read(1)
NDVI_2018_summer = NDVI_2018_summer.read(1)
NDVI_2018_autumn = NDVI_2018_autumn.read(1)
NDVI_2018_winter = NDVI_2018_winter.read(1)

ndvi1=ax1.imshow(NDVI_2018_spring, cmap='RdYlGn')
ax1.set_title("NDVI_2018_spring")
ndvi1.set_clim(vmin=-1, vmax=1)
divider = make_axes_locatable(ax1)
ax_cb = divider.new_horizontal(size="5%", pad=0.1)
fig.add_axes(ax_cb)
plt.colorbar(ndvi1, cax=ax_cb)

ndvi2=ax2.imshow(NDVI_2018_summer, cmap='RdYlGn')
ax2.set_title("NDVI_2018_summer")
ndvi2.set_clim(vmin=-1, vmax=1)
divider = make_axes_locatable(ax2)
ax_cb = divider.new_horizontal(size="5%", pad=0.1)
fig.add_axes(ax_cb)
plt.colorbar(ndvi2, cax=ax_cb)

ndvi3=ax3.imshow(NDVI_2018_autumn, cmap='RdYlGn')
ax3.set_title("NDVI_2018_autumn")
ndvi3.set_clim(vmin=-1, vmax=1)
divider = make_axes_locatable(ax3)
ax_cb = divider.new_horizontal(size="5%", pad=0.1)
fig.add_axes(ax_cb)
plt.colorbar(ndvi3, cax=ax_cb)

ndvi4=ax4.imshow(NDVI_2018_winter, cmap='RdYlGn')
ax4.set_title("NDVI_2018_winter")
ndvi4.set_clim(vmin=-1, vmax=1)
divider = make_axes_locatable(ax4)
ax_cb = divider.new_horizontal(size="5%", pad=0.05)
fig.add_axes(ax_cb)
plt.colorbar(ndvi4, cax=ax_cb)
Out[13]:
<matplotlib.colorbar.Colorbar at 0x7f1e01fc2780>