衛星データとPythonを使って植物の活性状況を宇宙から可視化

#前提
Arcgisと呼ばれるGISツールを使いNDVI(植生指数）を計算するプログラム
植生指数とは植物の光の反射の特性を活かし、植物の量や活性状況を簡易的に示す指標
使用したエディタはJupyterで言語はPython

#準備
下記のURLを参考にArcgis Rest APIを使えるように登録
Client ID
Client Secret
の２つをのちに使用。

#コード

import pandas as pd
import requests
import urllib.request

# Pandasでの表示が省略されないように設定
pd.set_option("max_columns", 100)
pd.set_option('max_rows',1000)
pd.set_option('max_info_columns',100)

まずはライブラリをimportする

ARCGIS_CLIENT_ID = ''
ARCGIS_CLIENT_SECRET = ''

登録した際に取得したCLIENT_IDとCLIENT_SECRETを環境変数として登録

def get_arcgis_token():
    response = requests.get('https://www.arcgis.com/sharing/rest/oauth2/token/',
                            params = {'client_id' : ARCGIS_CLIENT_ID,
                                              'client_secret': ARCGIS_CLIENT_SECRET,
                                              'grant_type': 'client_credentials'})
    arcgis_token = response.json()['access_token']
    return arcgis_token

access_tokenを取得するための関数を定義
APIにアクセスするときに必要

field_center =[[120.095415,40.345695]]
satellite_date = '2019/05/20'
arcgis_token = get_arcgis_token()

・field_center　
NDVIを計測したい土地の座標
複数座標を渡すことも可能

・satellite_date　
データが欲しい時期

・arcgis_token　
APIにアクセスするためのtoken

#points
#複数地点の衛星波長データを取りたいときに使う
#ここではSentinelと呼ばれるヨーロッパの衛星を使う
def get_sentinel_sample(arcgis_token):
    response = requests.get('https://sentinel.arcgis.com/arcgis/rest/services/Sentinel2/ImageServer/getSamples',
                        params = urllib.parse.urlencode({'geometryType' : 'esriGeometryMultiPoint',
                           'geometry': { "points":field_center, 'spatialReference': { 'wkid': 4326}},
                           'mosaicRule': {'mosaicMethod': 'esriMosaicAttribute',
                                                  'where': 'category = 1 AND cloudcover <= 0.10',
                                                  'sortField': 'acquisitionDate',
                                                  'sortValue': satellite_date,
                                                  'ascending': 'true' 
                                                   },
                            'token': arcgis_token,
                            'returnFirstValueOnly': 'true',
                            'f': 'json' }))
    sent_sample = response.json()
    return sent_sample

・"points":field_center, 'spatialReference': { 'wkid': 4326}}
先ほど指定した座標
位置情報を表す際に様々な座標系があるのだがwkidでどの座標系を使用するのか定義する
4326は地理座標系でありGPSで得られる座標系のこと

・'category = 1 AND cloudcover <= 0.10'　
たまに背景画像である灰色の画像が返ってくるのでcategoryを指定し衛星画像しか返ってこないようにする
cloudcoverを指定し雲がか10%以下しかかっていない衛星画像を指定

・'sortValue': satellite_date,
先ほど指定した日付

・token': arcgis_token, 'returnFirstValueOnly': 'true',
先ほど取得したtoken
APIが非常に重いのでここでは初めの値だけ返すように指定

詳しいAPIの仕様については公式のドキュメントを参照(かなり分かりにくい)
https://developers.arcgis.com/rest/services-reference/get-samples.htm

sent_sample = get_sentinel_sample(arcgis_token)
sat_val = sent_sample['samples'][0][('value')]
df_sat_band = pd.DataFrame(columns=['band1','band2','band3','band4','band5','band6','band7','band8','band8a','band9','band10','band11','band12'])
band_val = pd.Series([int(strip_num.strip()) for strip_num in sat_val.split()],index=df_sat_band.columns)
#appendは戻り値で新しいインスタンスを返すため代入し直す
df_sat_band  = df_sat_band.append(band_val, ignore_index=True)
df_sat_band['NDVI'] = (df_sat_band['band8'] - df_sat_band['band4']) /(df_sat_band['band8'] + df_sat_band['band4'])

返ってきたデータを整形しpandasのデータフレーム型にしNDVIを計算しカラムに追加
Sentinelの場合バンド4が赤の波長域でバンド8が近赤外の波長域である
※NDVI計算式 NDVI=(IR-R)/(IR+R)

このようにNDVIを計算することができている

参考文献
https://developers.arcgis.com/rest/services-reference/get-samples.htm
https://qiita.com/boiledorange73/items/b98d3d1ef3abf7299aba