Landsat-8の光学画像を利用し、水稲の収量減少を捉える

Tag

はじめに

今回のチュートリアルではLandsat 8の光学画像を利用します。Landsat 8はUSGS(United States Geological Survey)が運用している衛星です。複数の観測波長帯を持ち、可視光、近赤外・熱赤外領域を含めると11種のバンド情報を取得することができます。詳しくは、APIリファレンスをご覧ください。

センサの波長帯(電磁波の領域)が変化すると、観測したものの見え方が異なってしまうのは何故でしょうか。全ての物質は、その性質(例えば表面の構造)によって、波長帯に固有の反射の強さを示します。簡単に言ってしまえば、特定の色を反射しているということです。この性質を利用し、波長帯の組み合わせを変えることにより、水や植物といった特定の物質だけを際立たせる画像を作成できます。詳しくは、「課題に応じて変幻自在? 衛星データをブレンドして見えるモノ・コト」をご参照ください。

観測波長帯の組み合わせとして代表的なものにNDVI(Normalized Difference Vegetation Index: 正規化植生指標)があります。上記の記事でも言及していますが、NDVIとは植物が強く反射する近赤外領域と、植物が吸収する赤色の領域の違いを計量した指標です。NDVIは-1から1の範囲をとり、値が大きいほど植物が繁茂しており、値が低くなると植生がなくなっていきます(0に近くなると裸地、砂地、岩地、そして人工物となり、負の値になると水などとなります)。

NDVIを利用すれば、画像から植生のみを抽出することが可能です。反対に、低い値を閾値とすれば、人工物を抽出することもできるでしょう。本チュートリアルでは、代かき後と考えられる、2014年4月と2019年4月の衛星画像における、耕作放棄地(または畑地)と水田を区別することで、水田の面積がどれほど減少したのかを捉えてみます。

NDVIは実際に水稲の栽培管理に利用されています。例えば、青森県の「青天の霹靂」では、この指標を利用し、出穂後の適切な刈り取り時期を推定し、栽培者に知らせるというシステムを運用しているようです。

私たちも開発環境を活用して問題に取り組んでみましょう。

APIを用いて波長帯合成画像を取得する

Landsat-8のメタデータを取得するために、APIを利用します。APIの仕様はASNAROと同様になっています。

import os, json, requests, math, glob, re
import numpy as np
from skimage import io
from io import BytesIO
import dateutil.parser
from datetime import datetime
from datetime import timezone
import matplotlib.pyplot as plt
%matplotlib inline
TOKEN = "ここにトークンを貼る"
def get_landsat_scene(min_lat, min_lon, max_lat, max_lon):
	"""Landsat8のメタデータを取得します"""
	url = "https://gisapi.tellusxdp.com/api/v1/landsat8/scene" \
    	+ "?min_lat={}&min_lon={}&max_lat={}&max_lon={}".format(min_lat, min_lon, max_lat, max_lon)
	headers = {
    	"Authorization": "Bearer " + TOKEN
	}
	r = requests.get(url, headers=headers)
	return r.json()

青森県周辺のシーンを取得してみましょう。緯度や経度の情報については、前回のチュートリアルと全く同じように、Tellus OS上に図形を描き、GeoJSONファイルをダウンロードして調べることができます。

# 弘前周辺
scenes = get_landsat_scene(40.49, 140.502, 40.744, 140.797)
# {\min_lat\:40.49,\min_lon\:140.502,\max_lat\:40.744,\max_lon\:140.797}
print(len(scenes))
# 47

指定した座標の近傍から47個のデータを取得することができました。これらのメタデータは撮影時間の情報を含んでいますので、任意の期間だけのデータを抽出することにします。今回は2014年1月1日から2019年6月1日までを指定します。

def filter_by_date(scenes, start_date=None, end_date=None):
    	"""衛星の撮影時間を基にしてメタデータを抽出する"""
    	if(start_date == None):
        	start_date = datetime(1900,1,1, tzinfo=timezone.utc) # tellusにあるlandsat8のデータは2013年以降から現在まで
    	if(end_date == None):
        	end_date = datetime.now(timezone.utc) # 現在の時間をデフォルトで取得する
    	return [scene for scene in scenes if start_date <= dateutil.parser.parse(scene["acquisitionDate"]) and end_date > dateutil.parser.parse(scene["acquisitionDate"])]
# 2014年1月1日から2019年6月1日まで
landsat_scenes = filter_by_date(scenes, datetime(2014, 1, 1, tzinfo=timezone.utc), datetime(2019, 6, 1, tzinfo=timezone.utc))
print(len(landsat_scenes))
# 43

47個の内、43個が上で指定した期間内に含まれていました。データが範囲内に収まっているかをサムネイル画像と共に、その撮影時間を示して確認してみましょう。

ext_date = [re.search(r'\d{2} [A-z]{3} \d{4}',ext_date['acquisitionDate']) for ext_date in landsat_scenes]
period = [i.group() for i in ext_date]
print(period)

devLandsat8_proj_20200220_22.png

これらの全てのサムネイル画像を一括で表示してみましょう。

def get_tile_images(scenes, r_num, c_num, figsize_x=12, figsize_y=12):
    	"""サムネイル画像を指定した引数に応じて行列状に表示"""
    	thumbs = [imgs['thumbs_url'] for imgs in scenes]
    	f, ax_list = plt.subplots(r_num, c_num, figsize=(figsize_x,figsize_y))
    	for row_num, ax_row in enumerate(ax_list):
        	for col_num, ax in enumerate(ax_row):
            	if len(thumbs) < r_num*c_num:
                	len_shortage = r_num*c_num - len(thumbs) # 行列の不足分を算出
                	count = row_num * c_num + col_num
                	if count < len(thumbs):
                    	match = re.search(r'\d{2} [A-z]{3} \d{4}',scenes[row_num * c_num + col_num]['acquisitionDate'])
                    	D_M_Y = match.group()
                    	ax.label_outer() # サブプロットのタイトルと、軸のラベルが被らないようにする
                    	ax.imshow(io.imread(thumbs[row_num * c_num + col_num]))
                    	ax.set_title(D_M_Y)
                	else:
                    	for i in range(len_shortage):
                        	blank = np.zeros([100,100,3],dtype=np.uint8)
                        	blank.fill(255)
                        	ax.label_outer()
                        	ax.imshow(blank)
            	else:
                	match = re.search(r'\d{2} [A-z]{3} \d{4}',thumbs[row_num * c_num + col_num]['acquisitionDate'])
                	D_M_Y = match.group()
                	ax.label_outer()
                	ax.imshow(io.imread(thumbs[row_num * c_num + col_num]))
                	ax.set_title(D_M_Y)
    	return plt.show()

# サムネイル画像を取得
get_tile_images(landsat_scenes,8,6,16,16)

devLandsat8_proj_20200220_19.png

Credit: Landsat-8 image courtesy of the U.S. Geological Survey

データが範囲内に収まっていることがわかりました。この中から20 Apr 2014と18 Apr 2019のデータを選択します。

続いて、波長帯を指定した画像を抽出してみましょう。rgbにどの波長帯を指定するかにより、画像の見映えは変化します。初期値ではTrueカラー(人間の目で見たものと同じ画像となる)となっています。

def get_landsat_blend(path,row,productId,zoom, xtile, ytile, opacity=1, r=4, g=3, b=2, rdepth=1, gdepth=1, bdepth=1, preset=None):
	url = "https://gisapi.tellusxdp.com/blend/landsat8/{}/{}/{}/{}/{}/{}.png?opacity={}&r={}&g={}&b={}&rdepth={}&gdepth={}&bdepth={}".format\
	(path, row, productId, zoom, xtile, ytile, opacity, r, g, b, rdepth, gdepth, bdepth)
	if preset is not None:
    	url += '&preset='+preset
	headers = {
    	    	"Authorization": "Bearer " + TOKEN
	}
	r = requests.get(url, headers=headers)
	return io.imread(BytesIO(r.content))

2014年4月と2019年4月のデータを取得します。

Apr2014_land = landsat_scenes[-3]
Apr2019_land = landsat_scenes[0]

2014年の中泊町を含むタイル画像を二つ取得し、それらを結合しましょう。

Apr2014_land_b1 = get_landsat_blend(Apr2014_land['path'],Apr2014_land['row'],Apr2014_land['productId'],12, 3645, 1535)
Apr2014_land_b2 = get_landsat_blend(Apr2014_land['path'],Apr2014_land['row'],Apr2014_land['productId'],12, 3645, 1536)
plt.figure(figsize = (6,12))
comb_Apr2014_land = np.vstack((Apr2014_land_b1,Apr2014_land_b2))
plt.imshow(comb_Apr2014_land) # make a stacked image

devLandsat8_proj_20200220_11.png

Credit: Landsat-8 image courtesy of the U.S. Geological Survey

画像からNDVIを計算する

ここで、NDVIを導く関数を定義しましょう。NDVI=(IR-Red)/(IR+Red)で計算が可能です。加えて、指定した波長帯の値だけを取得する関数を作成します。

def band_extract_land(scene, zoom, xtile, ytile, opacity=1,\
             	r=4, g=3, b=2, rdepth=1, gdepth=1, bdepth=1, preset=None):
	"""rに指定されたLandsat8のバンドの情報を抽出する"""
	url = "https://gisapi.tellusxdp.com/blend/landsat8/{}/{}/{}/{}/{}/{}.png?opacity={}&r={}&g={}&b={}&rdepth={}&gdepth={}&bdepth={}".format\
	(scene['path'], scene['row'], scene['productId'], zoom, xtile, ytile, opacity,\
 	r, g, b, rdepth, gdepth, bdepth)
	if preset is not None:
    	url += '&preset='+preset
	headers = {
    	"Authorization": "Bearer " + TOKEN
	}
	r = requests.get(url, headers=headers)
	ext_band = io.imread(BytesIO(r.content))[:,:,0]
	return ext_band

def calc_ndvi(red, nir):
	red_float = red.astype(np.float32)
	nir_float = nir.astype(np.float32)
	# calculate NDVI
	# (IR-R)/(IR+R)
	numerator = (nir_float - red_float)
	denominator = (nir_float + red_float)
	flat_numerator = numerator.flatten()
	flat_denominator = denominator.flatten()
	with np.errstate(divide='ignore',invalid='ignore'):
    	ndvi =  np.where(flat_denominator != 0, flat_numerator / flat_denominator, np.nan)
	return ndvi

2014年の4月のデータに対して計算をします。

# rにband4(visible red lightを指定)
r_Apr201401 = get_landsat_blend(Apr2014_land['path'],Apr2014_land['row'],\
                      	Apr2014_land['productId'],12,3645,1535,r=4,g=3,b=2)[:,:,0]
# rにband5(NIRを指定)
nir_Apr201401 = get_landsat_blend(Apr2014_land['path'],Apr2014_land['row'],\
                        	Apr2014_land['productId'],12,3645,1535,r=5,g=3,b=2)[:,:,0]
Apr2014_land_ndvi1 = calc_ndvi(r_Apr201401,nir_Apr201401)
print(Apr2014_land_ndvi1)
# rにband4(visible red lightを指定)
r_Apr201402 = get_landsat_blend(Apr2014_land['path'],Apr2014_land['row'],\
                      	Apr2014_land['productId'],12,3645,1536,r=4,g=3,b=2)[:,:,0]
# rにband5(NIRを指定)
nir_Apr201402 = get_landsat_blend(Apr2014_land['path'],Apr2014_land['row'],\
                        	Apr2014_land['productId'],12,3645,1536,r=5,g=3,b=2)[:,:,0]
Apr2014_land_ndvi2 = calc_ndvi(r_Apr201402,nir_Apr201402)
print(Apr2014_land_ndvi2)

devLandsat8_proj_20200220_13.png

計算された結果を二次元配列に直します。二つの配列を結合し、再び同じ場所の画像を描写してみましょう。

Apr2014_land_ndvi1_wb = Apr2014_land_ndvi1.reshape(256,256)
Apr2014_land_ndvi2_wb = Apr2014_land_ndvi2.reshape(256,256)
Apr2014_land_ndvi_wb_stack = np.vstack((Apr2014_land_ndvi1_wb,Apr2014_land_ndvi2_wb))
plt.figure(figsize = (6,12))
plt.imshow(Apr2014_land_ndvi_wb_stack)
plt.colorbar()

ndvi1.png

Credit: Landsat-8 image courtesy of the U.S. Geological Survey

同じ場所から、異なる見た目の画像を得ることができました。波長の組み合わせを操作することにより、私たちの目で捉えたものとは別のものを得ることができます。NDVI値では、水面が最も暗く、樹林帯が最も明るい色で表示されていることがわかります。

青森北部では5月の中旬から下旬にかけてが一般的な田植えの時期になります。この時期には代かきが終わっているはずです。そのNDVIは裸地と似ていると考えられます。既往の研究例*1や、NASAの説明*2に基づき、代かき後水田のNDVIを0.01から0.08と設定します。

*1"Mapping rice planting areas in southern China using the ...." https://www.sciencedirect.com/science/article/pii/S0895717710005248

*2 "Measuring Vegetation (NDVI & EVI) - NASA ...." 30 8月. 2000, https://earthobservatory.nasa.gov/features/MeasuringVegetation。アクセス日: 9 12月. 2019。

def threshold_ndvi(stack_img,min_val,max_val):
	"""閾値から当該のNDVI値を抽出。画像サイズは512*256のみ"""
	stack_img_ext = np.where((stack_img > min_val) & \
                        	(stack_img < max_val), \
                        	stack_img, np.nan)
	stack_img_reshape = stack_img_ext.reshape(512,256)
	return stack_img_reshape
Apr2014_land_ndvi_wb_reshape = threshold_ndvi(Apr2014_land_ndvi_wb_stack,0.01,0.08)
plt.figure(figsize = (6,12))
plt.imshow(Apr2014_land_ndvi_wb_reshape)

devLandsat8_proj_20200220_18.png

Credit: Landsat-8 image courtesy of the U.S. Geological Survey

目分量ですが、ある程度分離できたように見えます。続いて、設定した閾値のNDVI値が画像中にどの程度含まれているのかを計算しましょう。

with np.errstate(invalid='ignore'):
	Apr2014_land_ndvi_sum =  np.where((Apr2014_land_ndvi_wb_stack > 0.01) & (Apr2014_land_ndvi_wb_stack < 0.08), 1, 0)
sum(Apr2014_land_ndvi_sum.flatten())/len(Apr2014_land_ndvi_sum.flatten())
# 0.6172714233398438

おおよそ、60%が水田であると見積もることができました。同様のことを2019年の4月のデータに対しても行います。

# rにband4(visible red lightを指定)
r_Apr201901 = band_extract_land(Apr2019_land,12,3645,1535)
# rにband5(NIRを指定)
nir_Apr201901 = band_extract_land(Apr2019_land,12,3645,1535,r=5)
Apr2019_land_ndvi1 = calc_ndvi(r_Apr201901,nir_Apr201901)
print(Apr2019_land_ndvi1)
# rにband4(visible red lightを指定)。タイルはy方向に移動
r_Apr201902 = band_extract_land(Apr2019_land,12,3645,1536)
# rにband5(NIRを指定)
nir_Apr201902 = band_extract_land(Apr2019_land,12,3645,1536,r=5)
Apr2019_land_ndvi2 = calc_ndvi(r_Apr201902,nir_Apr201902)
print(Apr2019_land_ndvi2)
Apr2019_land_ndvi1 = Apr2019_land_ndvi1.reshape(256,256)
Apr2019_land_ndvi2 = Apr2019_land_ndvi2.reshape(256,256)
Apr2019_land_ndvi_wb_stack = np.vstack((Apr2019_land_ndvi1,Apr2019_land_ndvi2))
plt.figure(figsize = (6,12))
plt.imshow(Apr2019_land_ndvi_wb_stack)
plt.colorbar()
plt.show()

ndvi2.png

Credit: Landsat-8 image courtesy of the U.S. Geological Survey

二つの画像を目で比べてみても、両者に大きな違いはないように見えます。閾値を指定して、値を抽出した後の画像で比べてみましょう。

Apr2019_land_ndvi_wb_reshape = threshold_ndvi(Apr2019_land_ndvi_wb_stack,0.01,0.08)
plt.figure(figsize = (6,12))
plt.imshow(Apr2019_land_ndvi_wb_reshape)

devLandsat8_proj_20200220_17.png

Credit: Landsat-8 image courtesy of the U.S. Geological Survey

同じ閾値でデータを分類したところ、2019年の方が水田の抽出率が悪いように見えます。実際に数値で確かめてみましょう。

with np.errstate(invalid='ignore'):
	Apr2019_land_ndvi_sum =  np.where((Apr2019_land_ndvi_wb_stack > 0.01) & (Apr2019_land_ndvi_wb_stack < 0.08), 1, 0)
sum(Apr2019_land_ndvi_sum.flatten())/len(Apr2019_land_ndvi_sum.flatten())
# 0.5083389282226562

50%程度が水田だと求められました。

NDVI画像から指定した領域の水田を抽出する

これらの作業では画像内の全ての水田(と思しき部分)を抽出してしまっています。正確な値を求めるためには、中泊町の範囲で画像を切り取る必要があります。そのためには、行政区域のデータを国土数値情報から入手する必要があります。国土数値情報行政区域データ青森県N03-110331_02を上のリンクからダウンロードしてください。具体的なダウンロード方法については、別記事をご参照ください。

ダウンロードしたファイルは開発環境の~/work/dataに保存します。

basepath = os.path.expanduser('~')
filepath = basepath + '/work/data/N03-19_02_190101.geojson'
# GeoJSONファイルから、属性データの読み込み
with open(filepath) as f:
	df = json.load(f)
# 中泊町を抽出
nakadomari_features = [f for f in df['features'] if f['properties']['N03_004']=='中泊町']
print(len(nakadomari_features))
print(nakadomari_features[0])

入手したGeoJSONファイルは青森県の全ての市町村データを含んでいます。ここから中泊町だけを選択し、新しいリストに代入しました。実際に区分を表示してみます。

from descartes import PolygonPatch
fig = plt.figure(figsize = (6,12))
ax = fig.gca()
ax.add_patch(PolygonPatch(nakadomari_features[0]['geometry'], fc="#cccccc", ec="#cccccc", alpha=0.5, zorder=2 ))
ax.axis('scaled')
plt.show()

devLandsat8_proj_20200220_2.png

この画像を作成したNDVI画像と重ね合わせることができれば、正確な範囲で切り取ることができます。先ずは、以下の関数を定義しましょう。

タイル画像にはタイル座標の値が含まれていますが、緯度経度の情報が入っているわけではありません。行政区域と重ね合わせるためには、タイル座標を緯度経度に変換する必要があります。座標やタイルの変換における詳しい話はこちらをご覧ください。

from skimage.draw import polygon
def get_tile_bbox(zoom, topleft_x, topleft_y, size_x=1, size_y=1):
	"""
	タイル座標からバウンディングボックスを取得
	https://tools.ietf.org/html/rfc7946#section-5
	"""
	def num2deg(xtile, ytile, zoom):
    	# https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
    	n = 2.0 ** zoom
    	lon_deg = xtile / n * 360.0 - 180.0
    	lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
    	lat_deg = math.degrees(lat_rad)
    	return (lon_deg, lat_deg)
	right_top = num2deg(topleft_x + size_x , topleft_y, zoom)
	left_bottom = num2deg(topleft_x, topleft_y + size_y, zoom)
	return (left_bottom[0], left_bottom[1], right_top[0], right_top[1])

行政区域のGeoJSONファイルには、緯度経度を示すための値が含まれています。ポリゴンはその緯度経度の情報を持ち、上図のように区域を描画することができます。画像の重ね合わせをするために、このポリゴンの座標情報をピクセルの座標情報に置き換えましょう。

def get_polygon_image(points, bbox, img_size):
	"""
	bboxの領域内にpointsの形状を描画した画像を作成する
	"""
	def lonlat2pixel(bbox, size, lon, lat):
    	"""
    	経度緯度をピクセル座標に変換
    	"""
    	dist = ((bbox[2] - bbox[0])/size[0], (bbox[3] - bbox[1])/size[1])
    	pixel = (int((lon - bbox[0]) / dist[0]), int((bbox[3] - lat) / dist[1]))
    	return pixel
	pixels = []
	for p in points:
    	pixels.append(lonlat2pixel(bbox, (img_size['width'], img_size['height']), p[0], p[1]))
	pixels = np.array(pixels)
	poly = np.ones((img_size['height'], img_size['width']), dtype=np.uint8)
	rr, cc = polygon(pixels[:, 1], pixels[:, 0], poly.shape)
	poly[rr, cc] = 0
	return poly

最後に、ポリゴンの範囲内だけのピクセル座標を抜き出す関数を定義します。

def get_mask_image(features, bbox, px_size):
	"""
	図形からマスク画像を作成する
	Parameters
	----------
	features : array of GeoJSON feature object
    	マスクするfeature objectの配列
	bbox : array
    	描画領域のバウンディングボックス
	px_size : array
    	描画領域のピクセルサイズ
	Returns
	-------
	img_array : ndarray
    	マスク画像
	"""
	ims = []
	for i in range(len(features)):
    	points = features[i]['geometry']['coordinates'][0]
    	ims.append(get_polygon_image(points, bbox, px_size))
	mask = np.where((ims[0] > 0), 0, 1)
	for i in range(1,len(ims)):
    	mask += np.where((ims[i] > 0), 0, 1)
	return np.where((mask > 0), 0, 1).astype(np.uint8)

関数を実行して、中泊町のみを抽出する画像を作成してみましょう。

bbox = get_tile_bbox(12, 3645, 1535, 1, 2)
px_size = {'width': Apr2019_land_ndvi_wb_stack.shape[1], 'height': Apr2019_land_ndvi_wb_stack.shape[0]}
nakadomari_features_mask = get_mask_image(nakadomari_features, bbox, px_size)
plt.figure(figsize = (6,12))
io.imshow(nakadomari_features_mask)

devLandsat8_proj_20200220_10.png

中泊町以外は黄色く表示されているのがわかります。

この画像とAPIで取得した画像を利用して、マスキングされた画像を作成します。

Apr2014_mask_img = np.where(nakadomari_features_mask != 1, Apr2014_land_ndvi_wb_stack, 0)
plt.figure(figsize = (6,12))
plt.imshow(Apr2014_mask_img)
plt.colorbar()

devLandsat8_proj_20200220_14.png

Credit: Landsat-8 image courtesy of the U.S. Geological Survey

中泊町の範囲だけが切り取られました。この範囲内で、先ほどと同様の条件で水田を抽出してみます。

with np.errstate(invalid='ignore'):
	Apr2014_mask_img_sum =  np.where((Apr2014_mask_img > 0.01) & (Apr2014_mask_img < 0.08), 1, 0)
sum(Apr2014_mask_img_sum.flatten())/len(Apr2014_mask_img_sum.flatten())
# 0.374542236328125

約4割が水田であることが示されました。2019年のデータでも同様の処理を行います。

Apr2019_mask_img = np.where(nakadomari_features_mask != 1, Apr2019_land_ndvi_wb_stack, 0)
plt.figure(figsize = (6,12))
plt.imshow(Apr2019_mask_img)
plt.colorbar()

devLandsat8_proj_20200220_24.png

Credit: Landsat-8 image courtesy of the U.S. Geological Survey

with np.errstate(invalid='ignore'):
	Apr2019_mask_img_sum =  np.where((Apr2019_mask_img > 0.01) & (Apr2019_mask_img < 0.08), 1, 0)
sum(Apr2019_mask_img_sum.flatten())/len(Apr2019_mask_img_sum.flatten())
# 0.29027557373046875

約3割が水田であると示されました。さらに、二つの画像から差分をとって、具体的にどの部分が水田でないと判断されたのかを確認しましょう。

# 2014年と2019年の差分画像を作成する
diff_ndvi = Apr2014_mask_img_sum - Apr2019_mask_img_sum
# 2019年では水田でなくなった部分を抽出
diff_ndvi_ext = np.where(diff_ndvi==1,1,0)
plt.figure(figsize = (6,12))
plt.imshow(diff_ndvi_ext)

devLandsat8_proj_20200220_23.png

Credit: Landsat-8 image courtesy of the U.S. Geological Survey

黄色い部分が2019年では水田と判断されなかった箇所になります。2014年と2019年のデータを比較したところ、中泊町では約1割程度が水田ではなくなったと予想されます。

ここまでは、あくまで衛星画像からの水田の変化を捉えたものでした。実際に米の収量は衛星画像のように減少をしているのでしょうか。これを調べるために、政府の統計情報を取得してみましょう。

e-Statから統計情報を取得する

米の収量データは、開発環境からe-StatのAPIを利用し、取得できます。さらに、得られたデータを表やグラフで表示し、衛星データで得られた結果がどれほど妥当であったのかを考察してみましょう。

可視化ライブラリ(matplotlib)で日本語を使うため、最初に日本語フォントを開発環境で使えるように準備しましょう。

devLandsat8_proj_20200220_8.png

「Launcher」から「Other」の項目にある「Terminal」を選択します。

日本語フォントをダウンロードするには、Terminalから次のコマンドを実行してください。

cd ~/work; wget https://ipafont.ipa.go.jp/IPAfont/ipag00303.zip; unzip ipag00303.zip

※ipafontのダウンロードサイトのURLが変化すると、wgetは失敗します。その場合はipaフォントのサイトからipag00303.zipを直接ダウンロードしてください。ダウンロード後は/work以下に解凍したフォルダをアップロードしてください。

devLandsat8_proj_20200220_21.png

以上で日本語を使用する準備が完了しました。

次にe-Statを使う準備を始めます。e-StatではAPIを提供しており、利用規約に従って、データベースから様々な統計データを取得することができます。

ログインページをクリックし、ユーザー登録を済ませましょう。

devLandsat8_proj_20200220_15.png

ユーザー登録ページは、このページの下部にあります。

devLandsat8_proj_20200220_4.png

仮登録後、本登録をするためのメールが届きます。手順に従い空欄に記入をし、登録を完了してください。e-StatのAPIを利用するためには、アプリケーションIDを発行する必要があります。トップページからマイページをクリックしてください。API機能(アプリケーションID発行)を選択し、IDを発行します。

devLandsat8_proj_20200220_16.png

devLandsat8_proj_20200220_5.png

URLは”http://localhost/”と入力してください。これでe-StatのAPIを使用する準備が完了しました。実際にAPIを使用して作物統計調査のデータを取得しましょう。

import requests,urllib
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
def get_json(base_url,params):
	params_str = urllib.parse.urlencode(params)
	url = base_url + params_str
	json = requests.get(url).json()
	return json
def main():
	appID = "ここに自分のアプリケーションIDを貼る"
	base_url = "http://api.e-stat.go.jp/rest/2.1/app/json/getStatsData?"
	params = {
    	"appId":appID,
    	"lang":"J",
    	"statsDataId":"0003293480",# e-Statで調べる
    	"metaGetFlg":"Y",
    	"cntGetFlg":"N",
    	"sectionHeaderFlg":"1"
	}
	json = get_json(base_url, params)
	# valuesを探す
	all_values = pd.DataFrame(json['GET_STATS_DATA']['STATISTICAL_DATA']['DATA_INF']['VALUE'])
 	# @cat01が120は収量(t) @areaは02387(中泊町)を指定
	values = all_values[(all_values["@cat01"] == "120") & (all_values["@area"] == "02387")]
	# 取得した時間の末尾に000000がつくので削除
	values = values.replace("0{6}$", "", regex=True)
	# 変数名を変更
	values_df = values.rename(columns = {"@time": "年", "$": "収量","@unit": "単位"})
	values_df = values_df.sort_values("年") # 年を基準に並べ替え
	return values_df
# 以下を参照
# https://qiita.com/kushikota/items/eed1d2b26e710fdc60f8
# https://irukanobox.blogspot.com/2018/07/pythone-statapi.html

(*このサービスは、政府統計総合窓口(e-Stat)のAPI機能を使用していますが、サービスの内容は国によって保証されたものではありません)

関数を実行して、e-Statから中泊町の米収量データを取得しましょう。

nakadomari = main()
nakadomari.head()

devLandsat8_proj_20200220_7.png

データが取得できました。次はmatplotlibを用いてグラフを描写してみましょう。

from pylab import *
from matplotlib import font_manager # 日本語をmatplotlibで使う
fontprop = matplotlib.font_manager.FontProperties(fname="/home/jovyan/work/ipag00303/ipag.ttf", size = 12) # フォントファイルの場所を指定
import seaborn as sns
import csv
# 変数がobjectなのでintに変換
nakadomari.loc[:,"年"] = nakadomari["年"].astype(int)
nakadomari.loc[:,"収量"] = nakadomari["収量"].astype(int)
x_range = range(len(nakadomari))
x_labels = nakadomari["年"]
y_range = nakadomari["収量"]
plt.figure(figsize = (14,6))
plt.plot(x_range, y_range)
plt.title('米の収量変化2007年〜2018年(中泊町)', fontdict = {"fontproperties": fontprop})
plt.xticks(x_range, x_labels, rotation='vertical')
plt.ylabel('収量(t)', fontdict = {"fontproperties": fontprop})
plt.show()

devLandsat8_proj_20200220_6.png

中泊町では2013年から2015年にかけて急激に米の収量が減少していることがわかります。2014年から2018年では、15,100トンから12,800トンへと減収しており、約1.5割の減収があったことが確認できました。

市浦観測所の気象データを取得する

収穫量減少の原因は耕作面積減少の他にも、気象要因が考えられます。気象庁から青森県市浦観測所のデータを取得してみましょう。気象庁はAPIを配布していませんので、Pythonからスクレイピングでデータを入手します。スクレイピングはAPIと異なり、開発者がデータを取得するウェブコンテンツのDOMを知っておかなければなりません。また、高頻度のスクレイピングは、サーバーに負荷を与えることもあります。詳しくは、こちらの情報をご覧ください。

from bs4 import BeautifulSoup
import seaborn as sns
import csv
# 市浦(prec_no=31, block_no=1119)を取得。年、月の指定は任意
base_url = "http://www.data.jma.go.jp/obd/stats/etrn/view/daily_a1.php?prec_no=31&block_no=1119&year=%s&month=%s&day=&view="

(*出典:気象庁ホームページ)

# strをfloatへ
def str2float(str):
  try:
	return float(str)
  except:
	return 0.0
if __name__ == "__main__":
	var_list = [['年月日','降水量合計','時最大降水量','分最大降水量','平均気温','最高気温','最低気温',\
         	'平均風速','最大風速','最大風速風向','最大瞬間風速','最大瞬間風速風向','最多風向','日照時間','降雪','最深積雪']]
	# 2014年から2017年の範囲
	for year in range(2007,2019): # 期間の指定(年)
  	# print(year)
  	# 当該年の1月から12月
  	for month in range(1,13): # 月指定
    	# 年と月を当てはめる。
    	r = requests.get(base_url%(year, month))
    	r.encoding = r.apparent_encoding
    	# テーブルを取得
    	soup = BeautifulSoup(r.text)
    	rows = soup.findAll('tr',class_='mtx') # タグ指定、クラス指定
    	# 表の最初の1から3行目をスライス
    	rows = rows[3:]
    	# 全ての行を読み込む(日データ取得)
    	for row in rows:
        	row_vals = row.findAll('td')
        	rowDate = [] # 初期化
        	rowDate.append(str(year) + "/" + str(month) + "/" + str(row_vals[0].string))
        	for i in range(15): # 15 variables are required
            	rowDate.append(str2float(row_vals[i+1].string))
        	var_list.append(rowDate)
# https://qiita.com/Cyber_Hacnosuke/items/122cec35d299c4d01f10を参照

取得したデータをcsvにして開発環境に保存しましょう。

# csvとして保存
with open('Ichiura.csv', 'w') as file:
	writer = csv.writer(file, lineterminator='\n')
	writer.writerows(var_list)

保存したcsvからpandasデータフレームを作成します。

## csvからデータフレームを作成
ichiura = pd.read_csv('Ichiura.csv', sep=',')
ichiura.head()

devLandsat8_proj_20200220_12.png

変数の年月日をインデックスとします。続いて新しい変数、年、月、日をそれぞれ作成します。

ichiura['年月日'] = pd.to_datetime(ichiura['年月日'], format='%Y-%m-%d')
ichiura = ichiura.set_index('年月日')
# 年、月、日、それぞれを変数としてデータフレームに格納する
ichiura['Year'] = ichiura.index.year
ichiura['Month'] = ichiura.index.month
ichiura['Day'] = ichiura.index.day
# 変数名の変更(日本語から英語へ)
ichiura.rename(columns={'日照時間':'Sunshine'}, inplace=True)

devLandsat8_proj_20200220_3.png

月別の積算日照時間を計算します。

# 日照時間の総和を求める(年ごとに)
sum_sunlight = ichiura.groupby(['Year','Month'],as_index=False).agg({"Sunshine":"sum"})
sum_sunlight['Date'] = pd.to_datetime(sum_sunlight[['Year', 'Month']].assign(Day=1),format = '%Y%m%d') # 年月の変数を追加(プロット用)
sum_sunlight.head()

devLandsat8_proj_20200220_25.png

これでグラフを描く準備ができました。可視化モジュールを使い、実際に描画してみましょう。

# 英語の月を追加
sum_sunlight['Month_en']= sum_sunlight['Date'].apply(lambda x: x.strftime('%B'))
plt.figure(figsize=(18, 10))
sns.set(font_scale=1.2)
sum_sunlight_p = sum_sunlight.pivot("Month", "Year", "Sunshine")
ax = sns.heatmap(sum_sunlight_p, yticklabels=sum_sunlight['Month_en'].unique(), cmap = 'OrRd')
ax.set_title('2007年1月から2018年12月までの中泊町月別積算日射量変化(時)', fontdict = {"fontproperties": fontprop}, fontsize = 14)
plt.show()

devLandsat8_proj_20200220_20.png

ヒートマップで2007年1月から2018年12月までの積算日照時間の変化を示しました。北国である中泊町では、田植えの開始時期が5月中旬から下旬に集中しています。出穂は8月の1週目になりますので、6月から8月までの日照は稲の生育に重要と考えられます。2009年の7月を除けば、夏季の日照時間に大きな差異はないように見えます。ただし、2018年は濃い赤色がなく、薄い赤色が目立つため、全体的な日照時間が平年より少ないと考えられます。さらに降水量も調べてみましょう。

# 月別の降水量を計算
# 月別の降水量を計算
# as_index = Falseで、集計のベースを変数として残す
sum_prec = ichiura.groupby(['Year',"Month"],as_index=False).agg({"降水量合計":"sum"})
# sum_prec2 = pd.concat((sum_prec,month_df),axis=1)
sum_prec['Date'] = pd.to_datetime(sum_sunlight[['Year', 'Month']].assign(Day=1),format = '%Y%m%d')
sum_prec.head()
fig, ax= plt.subplots(figsize=(18, 6))
# x_range = range(len(sum_prec))
x_range = sum_prec['Date']
y_prec = sum_prec['降水量合計']
# 7月と8月をハイライト
# ax.bar(x, y, width=10)
bar_list = ax.bar(x_range, y_prec,color = 'green', align='center',width = 20)
ax.xaxis_date()
i = 0
while i < 12:
	bar_list[5+12*i].set_color('yellow')
	bar_list[6+12*i].set_color('yellow')
	bar_list[7+12*i].set_color('yellow')
	i += 1
plt.title('2007年から2018年の中泊町月別降水量合計(mm)', \
      	fontdict = {"fontproperties": fontprop},fontsize = 16)
# set ticks
ax.xaxis.set_major_locator(mdates.MonthLocator(interval = 3))
ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y/%m"))
plt.setp(ax.xaxis.get_majorticklabels(), rotation=90, fontsize=12)
plt.ylabel('降水量合計(mm)', fontdict = {"fontproperties": fontprop}, fontsize = 14)
plt.show()

devLandsat8_proj_20200220_1.png

気象データは多少の増減はあるものの、一年を通してのパターンは一貫しています。日照時間であれば、冬場に少なく、春から夏へとかけて増加しており、降水量であれば、冬と夏にピークがあります。図から考えると2018年は年を通して降水量が多く、日照時間が平年より少ないように見えます。しかし気象データのみで、収量の減少の原因を捉えることは難しそうです(厳密には、収量が気象データによってどのように変化するのかを統計モデルにより考えることもできます。機械学習で収量の予測モデルを作成することもできるでしょう)。

中泊町は水稲単一経営が多い地域でしたが*3、米の消費が下がる中で、野菜や花卉などを合わせて行う複合経営が見られるようになりました。このような流れの中で、一部の水田は畑地へと転換されてしまったのではないでしょうか。

複合経営は労力がかかり、高齢の農業従事者には経営が難しいという欠点があります。そのため、中泊町では若手人材の育成に積極的に取り組んでいます。その流れの中で、2013年には「ばろかだる会」という若手経営者で構成される組織が出来上がりました。奇しくも、米の収量減少時期と重なっています。

米の収量減少における、はっきりとした因果関係はさておき、衛星画像から米収量の低下を見積もることはできたかと思います。

しかし、光学画像では雲の多い撮影シーンを利用できないなどの制限もあります。雲に対する解決策としては、SAR画像を用いることが考えられます。また、さらに高い精度で水田を抽出するためには機械学習を用いる、現地データも利用することが考えられます。収量の予測モデルを立てるのであれば、複数の変数を様々なデータベースから取得して、妥当なモデルを勘案することができます。

今回のチュートリアルでは、衛星データを活用し、水稲収穫量を見積もることができました。画像から得られたデータを、実際の統計データと比較することで、見積もりの妥当性も併せて確かめることができました。しかし、その原因については考察をするには、さらなる解析が必要と考えられます。

いずれにしても、地上での変化を衛星データから簡易的に定量できるということは、非常に魅力的なことです。APIを活用し、Tellusの開発環境内で様々なデータ解析を行ってみてください。

*3【中泊町】 農山漁村の「地域経営」取組事例|青森県庁ウェブ ...." 24 5月. 2019, https://www.pref.aomori.lg.jp/sangyo/agri/tiikikeiei_torikumi_23nakadomari.html。アクセス日: 10 12月. 2019。

今回利用したスクリプトです。

import os, json, requests, math, glob, re
import numpy as np
from skimage import io
from io import BytesIO
import dateutil.parser
from datetime import datetime
from datetime import timezone
import matplotlib.pyplot as plt
%matplotlib inline
TOKEN = "ここにトークンを貼る"
def get_landsat_scene(min_lat, min_lon, max_lat, max_lon):
	"""Landsat8のメタデータを取得します"""
	url = "https://gisapi.tellusxdp.com/api/v1/landsat8/scene" \
    	+ "?min_lat={}&min_lon={}&max_lat={}&max_lon={}".format(min_lat, min_lon, max_lat, max_lon)
	headers = {
    	"Authorization": "Bearer " + TOKEN
	}
	r = requests.get(url, headers=headers)
	return r.json()
# 弘前周辺
scenes = get_landsat_scene(40.49, 140.502, 40.744, 140.797)
# {\min_lat\:40.49,\min_lon\:140.502,\max_lat\:40.744,\max_lon\:140.797}
print(len(scenes))
def filter_by_date(scenes, start_date=None, end_date=None):
    	"""衛星の撮影時間を基にしてメタデータを抽出する"""
    	if(start_date == None):
        	start_date = datetime(1900,1,1, tzinfo=timezone.utc) # tellusにあるlandsat8のデータは2013年以降から現在まで
    	if(end_date == None):
        	end_date = datetime.now(timezone.utc) # 現在の時間をデフォルトで取得する
    	return [scene for scene in scenes if start_date <= dateutil.parser.parse(scene["acquisitionDate"]) and end_date > dateutil.parser.parse(scene["acquisitionDate"])]
landsat_scenes = filter_by_date(scenes, datetime(2014, 1, 1, tzinfo=timezone.utc), datetime(2019, 6, 1, tzinfo=timezone.utc))
ext_date = [re.search(r'\d{2} [A-z]{3} \d{4}',ext_date['acquisitionDate']) for ext_date in landsat_scenes]
# サムネイル画像を取得
thumbs_landsat = [imgs['thumbs_url'] for imgs in landsat_scenes]
def get_tile_images(thumbs, r_num, c_num, figsize_x=12, figsize_y=12):
    	"""サムネイル画像を指定した引数に応じて行列状に表示"""
    	f, ax_list = plt.subplots(r_num, c_num, figsize=(figsize_x,figsize_y))
    	for row_num, ax_row in enumerate(ax_list):
        	for col_num, ax in enumerate(ax_row):
            	if len(thumbs) < r_num*c_num:
                	len_shortage = r_num*c_num - len(thumbs) # 行列の不足分を算出
                	count = row_num * c_num + col_num
                	if count < len(thumbs):
                    	match = re.search(r'\d{2} [A-z]{3} \d{4}',landsat_scenes[row_num * c_num + col_num]['acquisitionDate'])
                    	D_M_Y = match.group()
                    	ax.label_outer() # サブプロットのタイトルと、軸のラベルが被らないようにする
                    	ax.imshow(io.imread(thumbs[row_num * c_num + col_num]))
                    	ax.set_title(D_M_Y)
                	else:
                    	for i in range(len_shortage):
                        	blank = np.zeros([100,100,3],dtype=np.uint8)
                        	blank.fill(255)
                        	ax.label_outer()
                        	ax.imshow(blank)
            	else:
                	match = re.search(r'\d{2} [A-z]{3} \d{4}',landsat_scenes[row_num * c_num + col_num]['acquisitionDate'])
                	D_M_Y = match.group()
                	ax.label_outer()
                	ax.imshow(io.imread(thumbs[row_num * c_num + col_num]))
                	ax.set_title(D_M_Y)
    	return plt.show()
get_tile_images(thumbs_landsat,8,6,16,16)
def get_landsat_blend(path,row,productId,zoom, xtile, ytile, opacity=1, r=4, g=3, b=2, rdepth=1, gdepth=1, bdepth=1, preset=None):
	url = "https://gisapi.tellusxdp.com/blend/landsat8/{}/{}/{}/{}/{}/{}.png?opacity={}&r={}&g={}&b={}&rdepth={}&gdepth={}&bdepth={}".format\
	(path, row, productId, zoom, xtile, ytile, opacity, r, g, b, rdepth, gdepth, bdepth)
	if preset is not None:
    	url += '&preset='+preset
	headers = {
    	"Authorization": "Bearer " + TOKEN
	}
	r = requests.get(url, headers=headers)
	return io.imread(BytesIO(r.content))
# 2014年4月と2019年4月のデータを取得
Apr2014_land = landsat_scenes[-3]
Apr2019_land = landsat_scenes[0]
Apr2014_land_b1 = get_landsat_blend(Apr2014_land['path'],Apr2014_land['row'],Apr2014_land['productId'],12, 3645, 1535)
Apr2014_land_b2 = get_landsat_blend(Apr2014_land['path'],Apr2014_land['row'],Apr2014_land['productId'],12, 3645, 1536)
plt.figure(figsize = (6,12))
comb_Apr2014_land = np.vstack((Apr2014_land_b1,Apr2014_land_b2))
plt.imshow(comb_Apr2014_land) # make a stacked image
def calc_ndvi(red, nir):
	red_float = red.astype(np.float32)
	nir_float = nir.astype(np.float32)
	# calculate NDVI
	# (IR-R)/(IR+R)
	numerator = (nir_float - red_float)
	denominator = (nir_float + red_float)
	flat_numerator = numerator.flatten()
	flat_denominator = denominator.flatten()
	with np.errstate(divide='ignore',invalid='ignore'):
    	ndvi =  np.where(flat_denominator != 0, flat_numerator / flat_denominator, np.nan)
	return ndvi
# rにband4(visible red lightを指定)
r_Apr201401 = band_extract_land(Apr2014_land,12,3645,1535)
# rにband5(NIRを指定)
nir_Apr201401 = band_extract_land(Apr2014_land,12,3645,1535,r=5)
Apr2014_land_ndvi1 = calc_ndvi(r_Apr201401,nir_Apr201401)
print(Apr2014_land_ndvi1)
# rにband4(visible red lightを指定)。タイルはy方向に移動
r_Apr201402 = band_extract_land(Apr2014_land,12,3645,1536)
# rにband5(NIRを指定)
nir_Apr201402 = band_extract_land(Apr2014_land,12,3645,1536,r=5)
Apr2014_land_ndvi1_wb = Apr2014_land_ndvi1.reshape(256,256)
Apr2014_land_ndvi2_wb = Apr2014_land_ndvi2.reshape(256,256)
Apr2014_land_ndvi_wb_stack = np.vstack((Apr2014_land_ndvi1_wb,Apr2014_land_ndvi2_wb))
plt.figure(figsize = (6,12))
plt.imshow(Apr2014_land_ndvi_wb_stack)
plt.colorbar()
def threshold_ndvi(stack_img,min_val,max_val):
	"""閾値から当該のNDVI値を抽出。画像サイズは512*256のみ"""
	stack_img_ext = np.where((stack_img > min_val) & \
                        	(stack_img < max_val), \
                        	stack_img, np.nan)
	stack_img_reshape = stack_img_ext.reshape(512,256)
	return stack_img_reshape
Apr2014_land_ndvi_wb_reshape = threshold_ndvi(Apr2014_land_ndvi_wb_stack,0.01,0.08)
plt.figure(figsize = (6,12))
plt.imshow(Apr2014_land_ndvi_wb_reshape)
with np.errstate(invalid='ignore'):
	Apr2014_land_ndvi_sum =  np.where((Apr2014_land_ndvi_wb_stack > 0.01) & (Apr2014_land_ndvi_wb_stack < 0.08), 1, 0)
sum(Apr2014_land_ndvi_sum.flatten())/len(Apr2014_land_ndvi_sum.flatten())
# Around 62% is paddy fields
# rにband4(visible red lightを指定)
r_Apr201901 = band_extract_land(Apr2019_land,12,3645,1535)
# rにband5(NIRを指定)
nir_Apr201901 = band_extract_land(Apr2019_land,12,3645,1535,r=5)
Apr2019_land_ndvi1 = calc_ndvi(r_Apr201901,nir_Apr201901)
print(Apr2019_land_ndvi1)
# rにband4(visible red lightを指定)。タイルはy方向に移動
r_Apr201902 = band_extract_land(Apr2019_land,12,3645,1536)
# rにband5(NIRを指定)
nir_Apr201902 = band_extract_land(Apr2019_land,12,3645,1536,r=5)
Apr2019_land_ndvi2 = calc_ndvi(r_Apr201902,nir_Apr201902)
print(Apr2019_land_ndvi2)
Apr2019_land_ndvi1 = Apr2019_land_ndvi1.reshape(256,256)
Apr2019_land_ndvi2 = Apr2019_land_ndvi2.reshape(256,256)
Apr2019_land_ndvi_wb_stack = np.vstack((Apr2019_land_ndvi1,Apr2019_land_ndvi2))
plt.figure(figsize = (6,12))
plt.imshow(Apr2019_land_ndvi_wb_stack)
plt.colorbar()
Apr2019_land_ndvi_wb_reshape = threshold_ndvi(Apr2019_land_ndvi_wb_stack,0.01,0.08)
plt.figure(figsize = (6,12))
plt.imshow(Apr2019_land_ndvi_wb_reshape)
with np.errstate(invalid='ignore'):
	Apr2019_land_ndvi_sum =  np.where((Apr2019_land_ndvi_wb_stack > 0.01) & (Apr2019_land_ndvi_wb_stack < 0.08), 1, 0)
sum(Apr2019_land_ndvi_sum.flatten())/len(Apr2019_land_ndvi_sum.flatten())
# Around 51% is paddy fields
basepath = os.path.expanduser('~')
filepath = basepath + '/work/data/N03-19_02_190101.geojson'
# Call district category from downloaded the geoJSON file
with open(filepath) as f:
	df = json.load(f)
# extract 中泊町
nakadomari_features = [f for f in df['features'] if f['properties']['N03_004']=='中泊町']
print(len(nakadomari_features))
print(nakadomari_features[0])
from descartes import PolygonPatch
fig = plt.figure(figsize = (6,12))
ax = fig.gca()
ax.add_patch(PolygonPatch(nakadomari_features[0]['geometry'], fc="#cccccc", ec="#cccccc", alpha=0.5, zorder=2 ))
ax.axis('scaled')
plt.show()
from skimage.draw import polygon
def get_tile_bbox(zoom, topleft_x, topleft_y, size_x=1, size_y=1):
	"""
	タイル座標からバウンディングボックスを取得
	https://tools.ietf.org/html/rfc7946#section-5
	"""
	def num2deg(xtile, ytile, zoom):
    	# https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
    	n = 2.0 ** zoom
    	lon_deg = xtile / n * 360.0 - 180.0
    	lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
    	lat_deg = math.degrees(lat_rad)
    	return (lon_deg, lat_deg)
	right_top = num2deg(topleft_x + size_x , topleft_y, zoom)
	left_bottom = num2deg(topleft_x, topleft_y + size_y, zoom)
	return (left_bottom[0], left_bottom[1], right_top[0], right_top[1])
def get_polygon_image(points, bbox, img_size):
	"""
	bboxの領域内にpointsの形状を描画した画像を作成する
	"""
	def lonlat2pixel(bbox, size, lon, lat):
    	"""
    	経度緯度をピクセル座標に変換
    	"""
    	dist = ((bbox[2] - bbox[0])/size[0], (bbox[3] - bbox[1])/size[1])
    	pixel = (int((lon - bbox[0]) / dist[0]), int((bbox[3] - lat) / dist[1]))
    	return pixel
	pixels = []
	for p in points:
    	pixels.append(lonlat2pixel(bbox, (img_size['width'], img_size['height']), p[0], p[1]))
	pixels = np.array(pixels)
	poly = np.ones((img_size['height'], img_size['width']), dtype=np.uint8)
	rr, cc = polygon(pixels[:, 1], pixels[:, 0], poly.shape)
	poly[rr, cc] = 0
	return poly
def get_mask_image(features, bbox, px_size):
	"""
	図形からマスク画像を作成する
	Parameters
	----------
	features : array of GeoJSON feature object
    	マスクするfeature objectの配列
	bbox : array
    	描画領域のバウンディングボックス
	px_size : array
    	描画領域のピクセルサイズ
	Returns
	-------
	img_array : ndarray
    	マスク画像
	"""
	ims = []
	for i in range(len(features)):
    	points = features[i]['geometry']['coordinates'][0]
    	ims.append(get_polygon_image(points, bbox, px_size))
	mask = np.where((ims[0] > 0), 0, 1)
	for i in range(1,len(ims)):
    	mask += np.where((ims[i] > 0), 0, 1)
	return np.where((mask > 0), 0, 1).astype(np.uint8)
bbox = get_tile_bbox(12, 3645, 1535, 1, 2)
px_size = {'width': Apr2019_land_ndvi_wb_stack.shape[1], 'height': Apr2019_land_ndvi_wb_stack.shape[0]}
nakadomari_features_mask = get_mask_image(nakadomari_features, bbox, px_size)
plt.figure(figsize = (6,12))
io.imshow(nakadomari_features_mask)
Apr2014_mask_img = np.where(nakadomari_features_mask != 1, Apr2014_land_ndvi_wb_stack, 0)
plt.figure(figsize = (6,12))
plt.imshow(Apr2014_mask_img)
plt.colorbar()
with np.errstate(invalid='ignore'):
	Apr2014_mask_img_sum =  np.where((Apr2014_mask_img > 0.01) & (Apr2014_mask_img < 0.08), 1, 0)
sum(Apr2014_mask_img_sum.flatten())/len(Apr2014_mask_img_sum.flatten())
# Around 37% is paddy fields in Nakadomari-machi
Apr2019_mask_img = np.where(nakadomari_features_mask != 1, Apr2019_land_ndvi_wb_stack, 0)
plt.figure(figsize = (6,12))
plt.imshow(Apr2019_mask_img)
plt.colorbar()
with np.errstate(invalid='ignore'):
	Apr2019_mask_img_sum =  np.where((Apr2019_mask_img > 0.01) & (Apr2019_mask_img < 0.08), 1, 0)
sum(Apr2019_mask_img_sum.flatten())/len(Apr2019_mask_img_sum.flatten())
# Around 30% is paddy fields in Nakadomari-machi
# 2014年と2019年の差分画像を作成する
diff_ndvi = Apr2014_mask_img_sum - Apr2019_mask_img_sum
# 2019年では水田でなくなった部分を抽出
diff_ndvi_ext = np.where(diff_ndvi==1,1,0)
plt.figure(figsize = (6,12))
plt.imshow(diff_ndvi_ext)
import requests,urllib
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
def get_json(base_url,params):
	params_str = urllib.parse.urlencode(params)
	url = base_url + params_str
	json = requests.get(url).json()
	return json
def main():
	appID = "ここに自分のアプリケーションIDを貼る"
	base_url = "http://api.e-stat.go.jp/rest/2.1/app/json/getStatsData?"
	params = {
    	"appId":appID,
    	"lang":"J",
    	"statsDataId":"0003293480",# e-Statで調べる
    	"metaGetFlg":"Y",
    	"cntGetFlg":"N",
    	"sectionHeaderFlg":"1"
	}
	json = get_json(base_url, params)
	# valuesを探す
	all_values = pd.DataFrame(json['GET_STATS_DATA']['STATISTICAL_DATA']['DATA_INF']['VALUE'])
 	# @cat01が120は収量(t) @areaは02387(中泊町)を指定
	values = all_values[(all_values["@cat01"] == "120") & (all_values["@area"] == "02387")]
	# 取得した時間の末尾に000000がつくので削除
	values = values.replace("0{6}$", "", regex=True)
	# 変数名を変更
	values_df = values.rename(columns = {"@time": "年", "$": "収量","@unit": "単位"})
	values_df = values_df.sort_values("年") # 年を基準に並べ替え
	return values_df

関数を実行して、e-Statから中泊町の米収量データを取得しましょう。

nakadomari = main()
nkadomari.head()
nakadomari = main()
nakadomari.head()
from pylab import *
from matplotlib import font_manager # 日本語をmatplotlibで使う
fontprop = matplotlib.font_manager.FontProperties(fname="/home/jovyan/work/ipag00303/ipag.ttf", size = 12) # フォントファイルの場所を指定
import seaborn as sns
import csv
# 変数がobjectなのでintに変換
nakadomari.loc[:,"年"] = nakadomari["年"].astype(int)
nakadomari.loc[:,"収量"] = nakadomari["収量"].astype(int)
x_range = range(len(hoge))
x_labels = nakadomari["年"]
y_range = nakadomari["収量"]
plt.figure(figsize = (14,6))
plt.plot(x_range, y_range)
plt.title('米の収量変化2007年〜2018年(中泊町)', fontdict = {"fontproperties": fontprop})
plt.xticks(x_range, x_labels, rotation='vertical')
plt.ylabel('収量(t)', fontdict = {"fontproperties": fontprop})
plt.show()
# JMAから気象データを取得し、グラフで表示する
mport requests, os
import matplotlib.pyplot as plt
import pandas as pd
import matplotlib.dates as mdates
from bs4 import BeautifulSoup
from pylab import *
from matplotlib import font_manager # 日本語をmatplotlibで使う
fontprop = matplotlib.font_manager.FontProperties(fname="/home/jovyan/work/ipag00303/ipag.ttf", size = 12) # フォントファイルの場所を指定
import seaborn as sns
import csv
# 出典:気象庁ホームページ
# 市浦(prec_no=31, block_no=1119)を取得。年の月の指定は任意
base_url = "http://www.data.jma.go.jp/obd/stats/etrn/view/daily_a1.php?prec_no=31&block_no=1119&year=%s&month=%s&day=&view="
# strをfloatへ
def str2float(str):
  try:
	return float(str)
  except:
	return 0.0
if __name__ == "__main__":
	var_list = [['年月日','降水量合計','時最大降水量','分最大降水量','平均気温','最高気温','最低気温',\
         	'平均風速','最大風速','最大風速風向','最大瞬間風速','最大瞬間風速風向','最多風向','日照時間','降雪','最深積雪']]
	# 2014年から2017年の範囲
	for year in range(2007,2019): # 期間の指定(年)
  	# print(year)
  	# 当該年の1月から12月
  	for month in range(1,13): # 月指定
    	# 年と月を当てはめる。
    	r = requests.get(base_url%(year, month))
    	r.encoding = r.apparent_encoding
    	# テーブルを取得
    	soup = BeautifulSoup(r.text)
    	rows = soup.findAll('tr',class_='mtx') #タグ指定、クラス指定
    	# 表の最初の1から3行目をスライス
    	rows = rows[3:]
    	# 全ての行を読み込む(日データ取得)
    	for row in rows:
        	row_vals = row.findAll('td')
        	rowDate = [] #初期化
        	rowDate.append(str(year) + "/" + str(month) + "/" + str(row_vals[0].string))
        	for i in range(15): # 15 variables are required
            	rowDate.append(str2float(row_vals[i+1].string))
        	var_list.append(rowDate)
# csvとして保存
with open('Ichiura.csv', 'w') as file:
	writer = csv.writer(file, lineterminator='\n')
	writer.writerows(var_list)
## csvからデータフレームを作成
ichiura = pd.read_csv('Ichiura.csv', sep=',')
ichiura['年月日'] = pd.to_datetime(ichiura['年月日'], format='%Y-%m-%d')
ichiura = ichiura.set_index('年月日')
# Add columns with year and month name
# Japanese characters may cause error in following processes; therefore, they should be English
ichiura['Year'] = ichiura.index.year
ichiura['Month'] = ichiura.index.month
ichiura['Day'] = ichiura.index.day
ichiura.rename(columns={'日照時間':'Sunshine'}, inplace=True)
# 日照時間の総和を求める(年別)
sum_sunlight = ichiura.groupby(['Year','Month'],as_index=False).agg({"Sunshine":"sum"})
sum_sunlight['Date'] = pd.to_datetime(sum_sunlight[['Year', 'Month']].assign(Day=1),format = '%Y%m%d')
sum_sunlight.head()
sum_sunlight['Month_en']= sum_sunlight['Date'].apply(lambda x: x.strftime('%B'))
# heatmap
plt.figure(figsize=(18, 10))
sns.set(font_scale=1.2)
sum_sunlight_p = sum_sunlight.pivot("Month", "Year", "Sunshine")
ax = sns.heatmap(sum_sunlight_p, yticklabels=sum_sunlight['Month_en'].unique(), cmap = 'OrRd')
ax.set_title('2007年1月から2018年12月までの中泊町月別積算日射量変化(時)', fontdict = {"fontproperties": fontprop}, fontsize = 14)
plt.show()
# 月別の降水量を計算
sum_prec = ichiura.groupby(['Year',"Month"],as_index=False).agg({"降水量合計":"sum"})
sum_prec.head()
sum_prec['Date'] = pd.to_datetime(sum_sunlight[['Year', 'Month']].assign(Day=1),format = '%Y%m%d')
sum_prec.head()
# barplot
bar_list = ax.bar(x_range, y_prec,color = 'green', align='center',width = 20)
ax.xaxis_date()
i = 0
while i < 12:
	bar_list[5+12*i].set_color('yellow')
	bar_list[6+12*i].set_color('yellow')
	bar_list[7+12*i].set_color('yellow')
	i += 1
plt.title('2007年から2018年の中泊町月別降水量合計(mm)',       	fontdict = {"fontproperties": fontprop},fontsize = 16)
# set ticks
ax.xaxis.set_major_locator(mdates.MonthLocator(interval = 3))
ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y/%m"))
plt.setp(ax.xaxis.get_majorticklabels(), rotation=90, fontsize=12)
plt.ylabel('降水量合計(mm)', fontdict = {"fontproperties": fontprop}, fontsize = 14)
plt.show()