衛星データからは、植物の育成度合いを示す指標である植生指数が分かる、ということは「課題に応じて変幻自在？ 衛星データをブレンドして見えるモノ・コト #マンガでわかる衛星データ」で紹介しました。この植生指数は、下図に示すように、近赤外と赤の反射強度により求められます。

過去の記事ではNDVIの値から緑化度合いを算出していたのですが、今回は趣向を変えて、RGBの赤色(R)の反射強度の差を衛星データから見てみました。このRの反射強度が高い場合には「緑化度合いが強い」と仮定しています。夏ごろに実る農作物が多いと仮定すれば、耕作地面積は夏ごろの緑化度合いと相関があると言えます。つまり、夏ごろのRの反射強度の変化を見ることで、耕作地面積の変化を見ることができるのではないか、と考えました。
そして、その結果が正しいかどうか、地域経済分析システム(RESAS)が出している統計情報の変化を確認し、その関係性を確認したのが本記事になります。

(1)RESASデータとは

地域経済分析システム（RESAS）とは、経済産業省と内閣官房（まち・ひと・しごと創生本部事務局）が提供する、官民ビッグデータを集約し可視化するシステムです。TellusではRESASが整備する統計データをデータベースとして保持しており、開発環境からアクセスすることができます。

統計データは、農業系、漁業系、林業系、雇用系、観光系、財政系、人口系、まちづくり系、商取引系、医療・福祉系、産業系と、様々な分野で取り揃えがあります。詳しくはデータリファレンスをご覧ください。

※現在、本記事で取り上げている方法では、RESASのデータを取得することができなくなっております。ご利用の際はご注意ください。

(2)コード紹介

本記事では滋賀県を対象に、耕作地面積の変化を確認しました。
正確に比較しようと思うと田畑(圃場)のみ抽出し、その植生指数の変化を確認する必要があります。しかし、今回は簡単に、耕作地面積は緑化度合いの相対差と相関があり、「同時期であれば田畑以外の緑化状態は同じである」と仮定しました。そして、選択した衛星画像をFalse Colorとして赤色の割合を算出し、時系列に比較。加えて、上記で算出した値と、RESASから取得した耕作面積の合計値を比較し、その傾向が同じかどうかを確認してみました。

そのためのコードが以下になります。

from PIL import Image
from io import BytesIO
import matplotlib.pyplot as plt
import numpy as np
import psycopg2
import requests
%matplotlib inline
agriculture_land_total_25204_year = []
agriculture_land_total_25204_value = []

#RESASデータから耕作地合計面積を取得
conn = psycopg2.connect("postgresql://jovyan:password@postgres:5432/resas")
cur = conn.cursor()
cur.execute("SELECT year,city_name,agriculture_land_total FROM agriculture_land_city WHERE city_code=25204")
#取得したデータの表示
print('集計年 | 市町村名 | 耕作面積合計')
for row in cur:
    print(row)
    agriculture_land_total_25204_year.append(row[0])
    agriculture_land_total_25204_value.append(row[2])
avnir_year = []
avnir_value = []

#TellusAPIでavnir2のpngを取得
token = '自身のトークンを貼り付けてください'
scene_id1 = 'ALAV2A028882890' #2006/8/10
scene_id2 = 'ALAV2A247832890' #2010/9/19
tile = '/12/3596/1620' #滋賀県近江八幡市周辺
headers = { 'Authorization': 'Bearer %s' % token }

#scene_id1 タイル取得
#rspId取得
url = 'https://gisapi.tellusxdp.com/api/v1/av2ori/get_scene/' + scene_id1
response = requests.get( url, headers=headers )
data = response.json()
rspId1 = (data["rspId"])
#撮影年
avnir_year.append(int(data["acquisitionDate"].split()[3]))
#avnir取得
url = 'https://gisapi.tellusxdp.com/av2ori/'+ rspId1 + '/' + scene_id1 + tile + '.png'
payload = { 'preset':'false','opacity':0.85}
response1 = requests.get( url, params=payload, headers=headers )

#scene_id2 タイル取得
#rspId取得
url = 'https://gisapi.tellusxdp.com/api/v1/av2ori/get_scene/' + scene_id2
response = requests.get( url, headers=headers )
data = response.json()
print(data)
rspId2 = (data["rspId"])
#撮影年
avnir_year.append(int(data["acquisitionDate"].split()[3]))
#avnir取得
url = 'https://gisapi.tellusxdp.com/av2ori/'+ rspId2 + '/' + scene_id2 + tile + '.png'
payload = { 'preset':'false','opacity':0.85}
response2 = requests.get( url, params=payload, headers=headers )

im1 = Image.open(BytesIO(response1.content))
im2 = Image.open(BytesIO(response2.content))

#赤を強調・カウント
width,height = im1.size
cnt1 = 0
cnt2 = 0
for y in range(height):
    for x in range(width):
        r1,g1,b1,etc = im1.getpixel((x,y))
        r2,g2,b2,etc = im2.getpixel((x,y))
        if r1 > 100 & g1  100 & g2 < 100:
            cnt2=cnt2+1
avnir_value = [cnt1/(height*width), cnt2/(height*width)]
fig = plt.figure()
fig.add_subplot(1,2,1)
plt.imshow(np.array(im1))
fig.add_subplot(1,2,2)
plt.imshow(np.array(im2))

# 植生カウント表示
print('AVNIRから取得した植生カウント')
print(cnt1)
print(cnt2)
fig, ax1 = plt.subplots()
ax1.bar(agriculture_land_total_25204_year,agriculture_land_total_25204_value, color='slateblue', label="land total")
ax1.set_ylim(min(agriculture_land_total_25204_value)*0.9, max(agriculture_land_total_25204_value)*1.01)
ax2 = ax1.twinx()
ax2.plot(avnir_year,avnir_value,color='orangered', marker='o', label="red pixel rate")
ax2.set_ylim(min(avnir_value)*0.95, max(avnir_value)*1.05)
h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax1.legend(h1+h2, l1+l2, loc='lower right')
plt.show()

(3)整合性の確認結果

2章で示したコードを元にプログラムを実行すると、下図のような結果を得られます。確かに耕作地面積が減った年に、衛星画像から算出した緑化度合いも減っているのが分かります。このことから、NDVIの値を算出しなかったとしても、収穫量の変化傾向は捉えられそうなことが分かりますね。

(4)まとめ

今回は、耕作地面積を直接的に求めるのではなく、関係ありそうな値(衛星画像の赤の反射強度)から類推し、その類推はあっていそうか、ということを正解値(RESASのデータ)から確認する、ということを試みました。
RGBの衛星画像のRの反射強度とRESASのデータと比較したところ、確かに同じ変化傾向を示していそうだ、ということが示唆されました。

衛星の中には、近赤外(IR)の観測をしていないものもあります。そのような場合であったとしても、何かしらの変化傾向が見えそうだ、ということが分かったのではないでしょうか。

今回は2時期でしか比較していないのですが、本当に相関関係があるのか、ということを確認するために、複数時期で比較してみることや多地点で見てみる、他にはNDVIの値との相関を見てみる、なども試みてみても良いでしょう。