農業分野における衛星データ利用事例

この章で学習すること

本章では、茨城県を対象として、過去約10年間のMODISデータを解析することで米の収量予測にチャレンジします。

  • pyMODIS を用いたMODIS衛星画像処理
  • gdal & numpy を用いた画像処理
  • rasterstat を用いた画像における統計情報の取得

背景

持続可能な開発目標(SDGs)の目標2「飢餓をゼロに」では、飢餓を終わらせるため、持続可能な農業を促進することが目標に定められています。直接農業に関わるターゲットは以下の通りです(数字はターゲット番号)。

  • 小規模食料生産者の農業生産性及び所得の倍増(2.3)
  • 持続可能な食料生産システムの確保とレジリエントな農業の実践(2.4)
  • 開発途上国の農業生産能力の向上(2.a)
  • 世界の農産物市場における貿易制限や歪みの是正と防止(2.b)

適切な農業経営や生産性向上のためには、作物の状況を正しく把握して判断することが重要です。既に先進的な農家はリモートセンシングの技術を導入しており、圃場単位ならドローン、広範囲に広がる場合は衛星が利用され始めています。

また、現在、農林水産省はICT技術等を活用した新たな農林水産統計調査の効率化を目指し、気象データと人工衛星から取得されるデータを用いて水稲収量の予測に取り組んでおり、現地調査に掛かる人員の効率化に取り組んでいます。詳細は農林水産省水稲の作柄に関する委員会の配布資料を参照してください。

この章では、2007年から2017年までのMODISのデータを解析し、茨城県全体を対象として水稲の収量を予測する式を作成し、2018年のデータで検証する方法を学んでいきます。

準備

In [ ]:
# google driveをマウント
# Tellusの開発環境やローカル環境を使用する想定であれば、このセルは削除可です。
from google.colab import drive 
drive = drive.mount('/content/drive')
Mounted at /content/drive

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

In [ ]:
!pip install pyMODIS
!pip install geopandas
!pip install pandas
!pip install matplotlib
!pip install numpy
!pip install scikit-learn
Collecting pyMODIS
  Downloading https://files.pythonhosted.org/packages/0b/e8/cb2d3deb1e94d4f141a3642abc030a47a4391eb5c466a4869554d10148d4/pyModis-2.1.0.tar.gz (3.0MB)
     |████████████████████████████████| 3.1MB 5.7MB/s 
Requirement already satisfied: numpy in /usr/local/lib/python3.6/dist-packages (from pyMODIS) (1.19.5)
Requirement already satisfied: future in /usr/local/lib/python3.6/dist-packages (from pyMODIS) (0.16.0)
Requirement already satisfied: requests in /usr/local/lib/python3.6/dist-packages (from pyMODIS) (2.23.0)
Requirement already satisfied: idna<3,>=2.5 in /usr/local/lib/python3.6/dist-packages (from requests->pyMODIS) (2.10)
Requirement already satisfied: urllib3!=1.25.0,!=1.25.1,<1.26,>=1.21.1 in /usr/local/lib/python3.6/dist-packages (from requests->pyMODIS) (1.24.3)
Requirement already satisfied: certifi>=2017.4.17 in /usr/local/lib/python3.6/dist-packages (from requests->pyMODIS) (2020.12.5)
Requirement already satisfied: chardet<4,>=3.0.2 in /usr/local/lib/python3.6/dist-packages (from requests->pyMODIS) (3.0.4)
Building wheels for collected packages: pyMODIS
  Building wheel for pyMODIS (setup.py) ... done
  Created wheel for pyMODIS: filename=pyModis-2.1.0-cp36-none-any.whl size=58346 sha256=a7789904c95172f0bbbd7cb29b3c68b08890512b42e574b928b5ae6284815fa1
  Stored in directory: /root/.cache/pip/wheels/a5/97/69/219a3bddbe564dc4c048f8ff8b90384e54ae9655c93d7cfeba
Successfully built pyMODIS
Installing collected packages: pyMODIS
Successfully installed pyMODIS-2.1.0
Collecting geopandas
  Downloading https://files.pythonhosted.org/packages/f7/a4/e66aafbefcbb717813bf3a355c8c4fc3ed04ea1dd7feb2920f2f4f868921/geopandas-0.8.1-py2.py3-none-any.whl (962kB)
     |████████████████████████████████| 972kB 5.5MB/s 
Requirement already satisfied: pandas>=0.23.0 in /usr/local/lib/python3.6/dist-packages (from geopandas) (1.1.5)
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 16.7MB/s 
Collecting fiona
  Downloading https://files.pythonhosted.org/packages/37/94/4910fd55246c1d963727b03885ead6ef1cd3748a465f7b0239ab25dfc9a3/Fiona-1.8.18-cp36-cp36m-manylinux1_x86_64.whl (14.8MB)
     |████████████████████████████████| 14.8MB 300kB/s 
Requirement already satisfied: shapely in /usr/local/lib/python3.6/dist-packages (from geopandas) (1.7.1)
Requirement already satisfied: pytz>=2017.2 in /usr/local/lib/python3.6/dist-packages (from pandas>=0.23.0->geopandas) (2018.9)
Requirement already satisfied: numpy>=1.15.4 in /usr/local/lib/python3.6/dist-packages (from pandas>=0.23.0->geopandas) (1.19.5)
Requirement already satisfied: python-dateutil>=2.7.3 in /usr/local/lib/python3.6/dist-packages (from pandas>=0.23.0->geopandas) (2.8.1)
Requirement already satisfied: certifi in /usr/local/lib/python3.6/dist-packages (from pyproj>=2.2.0->geopandas) (2020.12.5)
Collecting munch
  Downloading https://files.pythonhosted.org/packages/cc/ab/85d8da5c9a45e072301beb37ad7f833cd344e04c817d97e0cc75681d248f/munch-2.5.0-py2.py3-none-any.whl
Collecting click-plugins>=1.0
  Downloading https://files.pythonhosted.org/packages/e9/da/824b92d9942f4e472702488857914bdd50f73021efea15b4cad9aca8ecef/click_plugins-1.1.1-py2.py3-none-any.whl
Collecting cligj>=0.5
  Downloading https://files.pythonhosted.org/packages/42/1e/947eadf10d6804bf276eb8a038bd5307996dceaaa41cfd21b7a15ec62f5d/cligj-0.7.1-py3-none-any.whl
Requirement already satisfied: six>=1.7 in /usr/local/lib/python3.6/dist-packages (from fiona->geopandas) (1.15.0)
Requirement already satisfied: attrs>=17 in /usr/local/lib/python3.6/dist-packages (from fiona->geopandas) (20.3.0)
Requirement already satisfied: click<8,>=4.0 in /usr/local/lib/python3.6/dist-packages (from fiona->geopandas) (7.1.2)
Installing collected packages: pyproj, munch, click-plugins, cligj, fiona, geopandas
Successfully installed click-plugins-1.1.1 cligj-0.7.1 fiona-1.8.18 geopandas-0.8.1 munch-2.5.0 pyproj-3.0.0.post1
Requirement already satisfied: pandas in /usr/local/lib/python3.6/dist-packages (1.1.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: numpy>=1.15.4 in /usr/local/lib/python3.6/dist-packages (from pandas) (1.19.5)
Requirement already satisfied: pytz>=2017.2 in /usr/local/lib/python3.6/dist-packages (from pandas) (2018.9)
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: matplotlib in /usr/local/lib/python3.6/dist-packages (3.2.2)
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) (2.4.7)
Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib) (2.8.1)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib) (1.3.1)
Requirement already satisfied: numpy>=1.11 in /usr/local/lib/python3.6/dist-packages (from matplotlib) (1.19.5)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.6/dist-packages (from matplotlib) (0.10.0)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.6/dist-packages (from python-dateutil>=2.1->matplotlib) (1.15.0)
Requirement already satisfied: numpy in /usr/local/lib/python3.6/dist-packages (1.19.5)
Requirement already satisfied: scikit-learn in /usr/local/lib/python3.6/dist-packages (0.22.2.post1)
Requirement already satisfied: numpy>=1.11.0 in /usr/local/lib/python3.6/dist-packages (from scikit-learn) (1.19.5)
Requirement already satisfied: scipy>=0.17.0 in /usr/local/lib/python3.6/dist-packages (from scikit-learn) (1.4.1)
Requirement already satisfied: joblib>=0.11 in /usr/local/lib/python3.6/dist-packages (from scikit-learn) (1.0.0)

ライブラリのインポート

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

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

  • numPy

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

  • GeoPandas

    有名なデータ処理ライブラリpandasを、地理情報が取り扱いできるように拡張したライブラリです。pandasと同じ強力なデータ処理機能を持つだけでなく、他の地理空間データとの連携が可能です。

  • matplotlib

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

  • pyMOIDS

    pyMOIDSは、GNU GPL 2によって管理されているライブラリで、NASAが運用する衛星センサ"MODIS"のデータ処理専用のライブラリです。NASAのサーバから直接MODISデータをダウンロードして投影変換などの処理することが可能です。

  • scikit-learn

    scikit-learnは、オープンソースで公開されている機械学習用のライブラリです。誰でも無償で利用することが可能であり、分類/回帰/クラスタリングといったモデルを作成することが可能です。

In [ ]:
#必要ライブラリをインポート
import os
import glob
import numpy as np
import math
import matplotlib
import matplotlib.pyplot as plt
import shutil
import pandas as pd
import geopandas as gpd
import csv
import requests,urllib

from pymodis import downmodis
from pymodis.convertmodis_gdal import convertModisGDAL
from osgeo import gdal, gdalconst, gdal_array, osr
from sklearn.linear_model import LinearRegression
from sklearn import linear_model

print("done")
WxPython missing, no GUI enabled

衛星データの取得

この章では、アメリカ航空宇宙局(NASA)が運用するTerra衛星およびAqua衛星に搭載されている光学センサ「MODIS」が取得するデータを使用します。このデータは誰でも無償で利用することが可能です。

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

MODISについて

MODISは36種類の異なる波長帯(バンド)で地球を観測しており、プロダクトは以下の5種類に大別されます。

  • 放射輝度データ等(レベル1データ)
  • 大気プロダクト
  • 陸域プロダクト
  • 氷圏プロダクト
  • 海洋プロダクト

一般的にリモートセンシング技術を用いて植生を解析する際に用いられる代表的な指標としてNDVIとEVIが挙げられます。

今回は陸域プロダクトの1つであるMODIS Vegetation Index ProductよりEVIのデータを使用します。

コラム1:NDVIとEVI

光学衛星による植生の状況把握では、植生が強く反射する近赤外線および植生が吸収する可視域の赤色の波長がよく用いられます。この波長帯の組み合わせによって、植物の状態を表す正規化植生指数(Normalized Difference Vegetation Index: NDVI)を計算することができます。NDVIについては森林分野を参照して下さい。

NDVIの精度を向上させた植生の指標として、EVI (Enhanced Vegetation Index) も用いられています。

EVIは以下の式で求めることができます。

$EVI = \frac{G ×(IR - Red)}{(IR + C1 × Red)- (C2 × Blue + L)}$

ここで、Gはゲイン係数、C1とC2はエアロゾル補正係数、Lは背面効果補正係数を表しており、MODISのEVIプロダクトの場合、以下の数値が用いられています。

G = 2.5, C1 = 6, C2 = 7.5, L = 1

EVIは、近赤外線および赤色に加え、青色の波長帯を使うことで高バイオマス地域での感度を改善し、キャノピーのバックグラウンド信号の分離と大気の影響の低減を通じて植生モニタリングを改善して、植生信号を強化するように設計された「最適化された」植生指数です。

本章では、水稲の生育開始前と生育のピーク時(8月から9月)における差が明確に表れ、また、大気の影響に強く季節変化を捉えやすいと考えられるEVIを用いることとした。

通常、NDVIやEVIは-1から1までの値を取りますが、本章にて用いるMODISのプロダクトでは、符号付16bitのデータで提供されるているため、EVIの値は-2000から10000なっています。また、データがないエリアにおいては無効値として-3000が当てはめられています。

プロダクトの詳細はこちらを、EVIについてはこちらを参照してください。

API経由のデータ取得

APIを使ってMODISのデータをダウンロードするためには、NASAのEARTH DATAにアクセスして、ユーザー登録をしてください。

本章では、2007年から2017年までのMODISのデータを解析して収量を予測する式を作成し、2018年のデータで検証してみます。

そこで、NASAのLP DAAC (Land Processes Distributed Active Archive Center) よりAPI経由でデータをダウンロードし、対象期間のデータセットを作成します。

In [ ]:
#解析対象期間の指定とデータ保存用ディレクトリ(フォルダ)の作成

list = ["2007","2008","2009","2010","2011","2012","2013","2014","2015","2016","2017","2018"]

os.mkdir("average_EVI")

print("done")

本章では、Vegetation Indices 16-Day L3 Global 250mプロダクトを使用します。このプロダクトは、日別のEVIデータ16日分を合成することで雲の影響を減少させたものです。MODISのデータはほぼ毎日取得されますが、このプロダクトは16日ごとに提供されています。

また、プロダクト画像はシヌソイダル座標系で区切られ、各地域によってタイル番号が振られています。タイル番号についてはこちらを参照してください。

一般的に水稲収量を予測する際には植生指標の積算値が用いられますが、本章では水稲の植生指数が最も高くなる8月と9月で2時期のEVIを平均化して代表値とします。

また、8月と9月それぞれの代表プロダクトは画像を見て適切な組み合わせを選択しました。

ダウンロードしたデータにはそれぞれ異なるプロダクトがHDF形式で格納されているので、上記プロダクトを抜き出し、UTM zone 54Nの座標系を与え、GeoTIFF形式で保存します。

In [ ]:
for a in list:
  # 年ごとにデータ保存先のフォルダを作成
  os.mkdir(a)
  dest = "/content/" + a
  # タイル番号で取得するエリアを指定(茨城を取得するのでh29v05)
  tiles = "h29v05"
  # 代表プロダクトの組み合わせを下記のif文で指定
  if a == "2012":
    MODIS_date_After = "-09-20"
    MODIS_date_Before = "-08-25"
  elif a == "2013":
    MODIS_date_After = "-09-20"
    MODIS_date_Before = "-08-25"
  elif a == "2014":
    MODIS_date_After = "-09-20"
    MODIS_date_Before = "-08-25"
  else:
    MODIS_date_After = "-09-30"
    MODIS_date_Before = "-09-01"

  #データの検索期間を指定
  #検索期間のゴールを指定
  day = a + MODIS_date_After
  print(day)
  #検索期間のスタートを指定
  eday = a + MODIS_date_Before
  print(eday)
  #欲しいプロダクトを指定
  product = "MOD13Q1.006"
  #NASA Earthdataのアカウントユーザー名
  name = "ユーザー名を「””」の中に記入してください"
  #NASA Earthdataのアカウントパスワード
  path="アカウントパスワードを「””」の中に記入してください"
  #ダウンロード先のディレクトリを指定
  datapath = "/content/" + a
  #データのダウンロード
  modisDown = downmodis.downModis(destinationFolder=dest, product=product, user=name, password=path, tiles=tiles, today=day, enddate=eday)
  modisDown.connect()
  modisDown.downloadsAllDay()
  
  #使用するデータの選択
  hdf_list = glob.glob(os.path.join(datapath, 'MOD13Q1.*.hdf'))
  print(hdf_list)

  # EVIの抜出と投影変換
  for f in hdf_list:
    print(f)
    output_pref = "/content/" + a + "/" + str(f[14:])

    #HDFファイルに格納されているデータの数分のリスト。投影変換したいプロダクトを1,無視するものを0とする。
    #subsetのリストの長さはプロダクトによって異なります。本章の場合、EVI(250m)を使用するため左から2番目を"1"にする。
    subset = [0,1,0,0,0,0,0,0,0,0,0,0]

    convertsingle = convertModisGDAL(hdfname=f, prefix=output_pref, subset=subset, res=250, epsg=32654)
    convertsingle.run()

  #年ごとの平均値画像の作成
  EVI_list = glob.glob(datapath + "/" + "*EVI.tif")
  print(EVI_list) 


  src = gdal.Open(EVI_list[1], gdalconst.GA_ReadOnly)
  img = np.array(src.ReadAsArray()).astype(np.float32)

  pixel=img.shape[0]
  line=img.shape[1]


  base_array = np.zeros((pixel,line))

  for EVI in EVI_list:
    print("OPEN" + EVI)
    data = gdal.Open(EVI, gdalconst.GA_ReadOnly)
    imgs = np.array(data.ReadAsArray()).astype(np.float32)
    base_array = base_array + imgs
    base_array

  EVI_ave = base_array/len(EVI_list)

  array = EVI_ave

  xsize = src.RasterXSize
  ysize = src.RasterYSize
  geotransform = src.GetGeoTransform()

  print("make EVI average image")

  dst = os.path.join("/content/average_EVI/" + "EVI_average_" + a + ".tif")
  driver = gdal.GetDriverByName('GTiff')
  dst_raster = driver.Create(dst, xsize, ysize, 1, gdal.GDT_Int16)

  dst_raster.SetProjection(src.GetProjection())      #{出力変数}.SetProjection(座標系情報)
  dst_raster.SetGeoTransform(src.GetGeoTransform())  #{出力変数}.SetGeoTransform(座標に関する6つの数字)

  dst_band = dst_raster.GetRasterBand(1)
  dst_band.WriteArray(EVI_ave)

  dst_band.FlushCache()
  dst_raster = None
  print("done")
2007-09-30
2007-09-01
['/content/2007/MOD13Q1.A2007257.h29v05.006.2015167042547.hdf', '/content/2007/MOD13Q1.A2007273.h29v05.006.2015167045436.hdf']
/content/2007/MOD13Q1.A2007257.h29v05.006.2015167042547.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2007/MOD13Q1.A2007257.h29v05.006.2015167042547.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2007/MOD13Q1.A2007257.h29v05.006.2015167042547.hdf' reprojected
/content/2007/MOD13Q1.A2007273.h29v05.006.2015167045436.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2007/MOD13Q1.A2007273.h29v05.006.2015167045436.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2007/MOD13Q1.A2007273.h29v05.006.2015167045436.hdf' reprojected
['/content/2007/MOD13Q1.A2007257.h29v05.006.2015167042547.hdf_250m 16 days EVI.tif', '/content/2007/MOD13Q1.A2007273.h29v05.006.2015167045436.hdf_250m 16 days EVI.tif']
OPEN/content/2007/MOD13Q1.A2007257.h29v05.006.2015167042547.hdf_250m 16 days EVI.tif
OPEN/content/2007/MOD13Q1.A2007273.h29v05.006.2015167045436.hdf_250m 16 days EVI.tif
make EVI average image
done
2008-09-30
2008-09-01
['/content/2008/MOD13Q1.A2008273.h29v05.006.2015180221615.hdf', '/content/2008/MOD13Q1.A2008257.h29v05.006.2015179232557.hdf']
/content/2008/MOD13Q1.A2008273.h29v05.006.2015180221615.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2008/MOD13Q1.A2008273.h29v05.006.2015180221615.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2008/MOD13Q1.A2008273.h29v05.006.2015180221615.hdf' reprojected
/content/2008/MOD13Q1.A2008257.h29v05.006.2015179232557.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2008/MOD13Q1.A2008257.h29v05.006.2015179232557.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2008/MOD13Q1.A2008257.h29v05.006.2015179232557.hdf' reprojected
['/content/2008/MOD13Q1.A2008257.h29v05.006.2015179232557.hdf_250m 16 days EVI.tif', '/content/2008/MOD13Q1.A2008273.h29v05.006.2015180221615.hdf_250m 16 days EVI.tif']
OPEN/content/2008/MOD13Q1.A2008257.h29v05.006.2015179232557.hdf_250m 16 days EVI.tif
OPEN/content/2008/MOD13Q1.A2008273.h29v05.006.2015180221615.hdf_250m 16 days EVI.tif
make EVI average image
done
2009-09-30
2009-09-01
['/content/2009/MOD13Q1.A2009257.h29v05.006.2015194230815.hdf', '/content/2009/MOD13Q1.A2009273.h29v05.006.2015195131227.hdf']
/content/2009/MOD13Q1.A2009257.h29v05.006.2015194230815.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2009/MOD13Q1.A2009257.h29v05.006.2015194230815.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2009/MOD13Q1.A2009257.h29v05.006.2015194230815.hdf' reprojected
/content/2009/MOD13Q1.A2009273.h29v05.006.2015195131227.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2009/MOD13Q1.A2009273.h29v05.006.2015195131227.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2009/MOD13Q1.A2009273.h29v05.006.2015195131227.hdf' reprojected
['/content/2009/MOD13Q1.A2009273.h29v05.006.2015195131227.hdf_250m 16 days EVI.tif', '/content/2009/MOD13Q1.A2009257.h29v05.006.2015194230815.hdf_250m 16 days EVI.tif']
OPEN/content/2009/MOD13Q1.A2009273.h29v05.006.2015195131227.hdf_250m 16 days EVI.tif
OPEN/content/2009/MOD13Q1.A2009257.h29v05.006.2015194230815.hdf_250m 16 days EVI.tif
make EVI average image
done
2010-09-30
2010-09-01
['/content/2010/MOD13Q1.A2010273.h29v05.006.2015211141233.hdf', '/content/2010/MOD13Q1.A2010257.h29v05.006.2015211022501.hdf']
/content/2010/MOD13Q1.A2010273.h29v05.006.2015211141233.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2010/MOD13Q1.A2010273.h29v05.006.2015211141233.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2010/MOD13Q1.A2010273.h29v05.006.2015211141233.hdf' reprojected
/content/2010/MOD13Q1.A2010257.h29v05.006.2015211022501.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2010/MOD13Q1.A2010257.h29v05.006.2015211022501.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2010/MOD13Q1.A2010257.h29v05.006.2015211022501.hdf' reprojected
['/content/2010/MOD13Q1.A2010273.h29v05.006.2015211141233.hdf_250m 16 days EVI.tif', '/content/2010/MOD13Q1.A2010257.h29v05.006.2015211022501.hdf_250m 16 days EVI.tif']
OPEN/content/2010/MOD13Q1.A2010273.h29v05.006.2015211141233.hdf_250m 16 days EVI.tif
OPEN/content/2010/MOD13Q1.A2010257.h29v05.006.2015211022501.hdf_250m 16 days EVI.tif
make EVI average image
done
2011-09-30
2011-09-01
['/content/2011/MOD13Q1.A2011273.h29v05.006.2015222121912.hdf', '/content/2011/MOD13Q1.A2011257.h29v05.006.2015222094913.hdf']
/content/2011/MOD13Q1.A2011273.h29v05.006.2015222121912.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2011/MOD13Q1.A2011273.h29v05.006.2015222121912.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2011/MOD13Q1.A2011273.h29v05.006.2015222121912.hdf' reprojected
/content/2011/MOD13Q1.A2011257.h29v05.006.2015222094913.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2011/MOD13Q1.A2011257.h29v05.006.2015222094913.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2011/MOD13Q1.A2011257.h29v05.006.2015222094913.hdf' reprojected
['/content/2011/MOD13Q1.A2011273.h29v05.006.2015222121912.hdf_250m 16 days EVI.tif', '/content/2011/MOD13Q1.A2011257.h29v05.006.2015222094913.hdf_250m 16 days EVI.tif']
OPEN/content/2011/MOD13Q1.A2011273.h29v05.006.2015222121912.hdf_250m 16 days EVI.tif
OPEN/content/2011/MOD13Q1.A2011257.h29v05.006.2015222094913.hdf_250m 16 days EVI.tif
make EVI average image
done
2012-09-20
2012-08-25
['/content/2012/MOD13Q1.A2012257.h29v05.006.2015249161227.hdf', '/content/2012/MOD13Q1.A2012241.h29v05.006.2015249123812.hdf']
/content/2012/MOD13Q1.A2012257.h29v05.006.2015249161227.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2012/MOD13Q1.A2012257.h29v05.006.2015249161227.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2012/MOD13Q1.A2012257.h29v05.006.2015249161227.hdf' reprojected
/content/2012/MOD13Q1.A2012241.h29v05.006.2015249123812.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2012/MOD13Q1.A2012241.h29v05.006.2015249123812.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2012/MOD13Q1.A2012241.h29v05.006.2015249123812.hdf' reprojected
['/content/2012/MOD13Q1.A2012241.h29v05.006.2015249123812.hdf_250m 16 days EVI.tif', '/content/2012/MOD13Q1.A2012257.h29v05.006.2015249161227.hdf_250m 16 days EVI.tif']
OPEN/content/2012/MOD13Q1.A2012241.h29v05.006.2015249123812.hdf_250m 16 days EVI.tif
OPEN/content/2012/MOD13Q1.A2012257.h29v05.006.2015249161227.hdf_250m 16 days EVI.tif
make EVI average image
done
2013-09-20
2013-08-25
['/content/2013/MOD13Q1.A2013257.h29v05.006.2015268075618.hdf', '/content/2013/MOD13Q1.A2013241.h29v05.006.2015268074206.hdf']
/content/2013/MOD13Q1.A2013257.h29v05.006.2015268075618.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2013/MOD13Q1.A2013257.h29v05.006.2015268075618.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2013/MOD13Q1.A2013257.h29v05.006.2015268075618.hdf' reprojected
/content/2013/MOD13Q1.A2013241.h29v05.006.2015268074206.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2013/MOD13Q1.A2013241.h29v05.006.2015268074206.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2013/MOD13Q1.A2013241.h29v05.006.2015268074206.hdf' reprojected
['/content/2013/MOD13Q1.A2013241.h29v05.006.2015268074206.hdf_250m 16 days EVI.tif', '/content/2013/MOD13Q1.A2013257.h29v05.006.2015268075618.hdf_250m 16 days EVI.tif']
OPEN/content/2013/MOD13Q1.A2013241.h29v05.006.2015268074206.hdf_250m 16 days EVI.tif
OPEN/content/2013/MOD13Q1.A2013257.h29v05.006.2015268075618.hdf_250m 16 days EVI.tif
make EVI average image
done
2014-09-20
2014-08-25
['/content/2014/MOD13Q1.A2014241.h29v05.006.2015289183932.hdf', '/content/2014/MOD13Q1.A2014257.h29v05.006.2015291064259.hdf']
/content/2014/MOD13Q1.A2014241.h29v05.006.2015289183932.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2014/MOD13Q1.A2014241.h29v05.006.2015289183932.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2014/MOD13Q1.A2014241.h29v05.006.2015289183932.hdf' reprojected
/content/2014/MOD13Q1.A2014257.h29v05.006.2015291064259.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2014/MOD13Q1.A2014257.h29v05.006.2015291064259.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2014/MOD13Q1.A2014257.h29v05.006.2015291064259.hdf' reprojected
['/content/2014/MOD13Q1.A2014241.h29v05.006.2015289183932.hdf_250m 16 days EVI.tif', '/content/2014/MOD13Q1.A2014257.h29v05.006.2015291064259.hdf_250m 16 days EVI.tif']
OPEN/content/2014/MOD13Q1.A2014241.h29v05.006.2015289183932.hdf_250m 16 days EVI.tif
OPEN/content/2014/MOD13Q1.A2014257.h29v05.006.2015291064259.hdf_250m 16 days EVI.tif
make EVI average image
done
2015-09-30
2015-09-01
['/content/2015/MOD13Q1.A2015257.h29v05.006.2015306145336.hdf', '/content/2015/MOD13Q1.A2015273.h29v05.006.2015307053832.hdf']
/content/2015/MOD13Q1.A2015257.h29v05.006.2015306145336.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2015/MOD13Q1.A2015257.h29v05.006.2015306145336.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2015/MOD13Q1.A2015257.h29v05.006.2015306145336.hdf' reprojected
/content/2015/MOD13Q1.A2015273.h29v05.006.2015307053832.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2015/MOD13Q1.A2015273.h29v05.006.2015307053832.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2015/MOD13Q1.A2015273.h29v05.006.2015307053832.hdf' reprojected
['/content/2015/MOD13Q1.A2015257.h29v05.006.2015306145336.hdf_250m 16 days EVI.tif', '/content/2015/MOD13Q1.A2015273.h29v05.006.2015307053832.hdf_250m 16 days EVI.tif']
OPEN/content/2015/MOD13Q1.A2015257.h29v05.006.2015306145336.hdf_250m 16 days EVI.tif
OPEN/content/2015/MOD13Q1.A2015273.h29v05.006.2015307053832.hdf_250m 16 days EVI.tif
make EVI average image
done
2016-09-30
2016-09-01
['/content/2016/MOD13Q1.A2016257.h29v05.006.2016274144916.hdf', '/content/2016/MOD13Q1.A2016273.h29v05.006.2016292071240.hdf']
/content/2016/MOD13Q1.A2016257.h29v05.006.2016274144916.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2016/MOD13Q1.A2016257.h29v05.006.2016274144916.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2016/MOD13Q1.A2016257.h29v05.006.2016274144916.hdf' reprojected
/content/2016/MOD13Q1.A2016273.h29v05.006.2016292071240.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2016/MOD13Q1.A2016273.h29v05.006.2016292071240.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2016/MOD13Q1.A2016273.h29v05.006.2016292071240.hdf' reprojected
['/content/2016/MOD13Q1.A2016257.h29v05.006.2016274144916.hdf_250m 16 days EVI.tif', '/content/2016/MOD13Q1.A2016273.h29v05.006.2016292071240.hdf_250m 16 days EVI.tif']
OPEN/content/2016/MOD13Q1.A2016257.h29v05.006.2016274144916.hdf_250m 16 days EVI.tif
OPEN/content/2016/MOD13Q1.A2016273.h29v05.006.2016292071240.hdf_250m 16 days EVI.tif
make EVI average image
done
2017-09-30
2017-09-01
['/content/2017/MOD13Q1.A2017273.h29v05.006.2017290090926.hdf', '/content/2017/MOD13Q1.A2017257.h29v05.006.2017276133113.hdf']
/content/2017/MOD13Q1.A2017273.h29v05.006.2017290090926.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2017/MOD13Q1.A2017273.h29v05.006.2017290090926.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2017/MOD13Q1.A2017273.h29v05.006.2017290090926.hdf' reprojected
/content/2017/MOD13Q1.A2017257.h29v05.006.2017276133113.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2017/MOD13Q1.A2017257.h29v05.006.2017276133113.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2017/MOD13Q1.A2017257.h29v05.006.2017276133113.hdf' reprojected
['/content/2017/MOD13Q1.A2017273.h29v05.006.2017290090926.hdf_250m 16 days EVI.tif', '/content/2017/MOD13Q1.A2017257.h29v05.006.2017276133113.hdf_250m 16 days EVI.tif']
OPEN/content/2017/MOD13Q1.A2017273.h29v05.006.2017290090926.hdf_250m 16 days EVI.tif
OPEN/content/2017/MOD13Q1.A2017257.h29v05.006.2017276133113.hdf_250m 16 days EVI.tif
make EVI average image
done
2018-09-30
2018-09-01
['/content/2018/MOD13Q1.A2018273.h29v05.006.2018295112642.hdf', '/content/2018/MOD13Q1.A2018257.h29v05.006.2018282131938.hdf']
/content/2018/MOD13Q1.A2018273.h29v05.006.2018295112642.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2018/MOD13Q1.A2018273.h29v05.006.2018295112642.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2018/MOD13Q1.A2018273.h29v05.006.2018295112642.hdf' reprojected
/content/2018/MOD13Q1.A2018257.h29v05.006.2018282131938.hdf
Layer HDF4_EOS:EOS_GRID:"/content/2018/MOD13Q1.A2018257.h29v05.006.2018282131938.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days EVI reprojected
All layer for dataset '/content/2018/MOD13Q1.A2018257.h29v05.006.2018282131938.hdf' reprojected
['/content/2018/MOD13Q1.A2018273.h29v05.006.2018295112642.hdf_250m 16 days EVI.tif', '/content/2018/MOD13Q1.A2018257.h29v05.006.2018282131938.hdf_250m 16 days EVI.tif']
OPEN/content/2018/MOD13Q1.A2018273.h29v05.006.2018295112642.hdf_250m 16 days EVI.tif
OPEN/content/2018/MOD13Q1.A2018257.h29v05.006.2018282131938.hdf_250m 16 days EVI.tif
make EVI average image
done

サンプルとして2007年におけるEVI平均値画像を表示してみましょう。

In [ ]:
EVI_image = gdal.Open("/content/average_EVI/EVI_average_2007.tif")

EVI_array = EVI_image.ReadAsArray()

plt.figure(figsize=(20,10))
plt.imshow(EVI_array,vmin=-3000,vmax=10000,cmap='RdYlGn')

plt.title("EVI 2007")
plt.colorbar()
plt.show
Out[ ]:
<function matplotlib.pyplot.show>

上の画像ではEVIの範囲を-3000から10000で表示しています。 画像において緑色が強いほど植物の活性が高いことを示しており、赤が強いほど植生が少ないことを示しています。 また、水域が赤く染まっていますが、こちらは水域を無効値である-3000であることを示しています。

画像を見ると、東京都や神奈川県などの都市部においては赤色が強く、植生が少ないことが見て取れます。

水田域におけるEVIの集計

今回の解析では、茨城県を対象に水稲の収量を予測します。茨城県全域のEVIを解析すると、水田以外に畑や森林の植生も入ってしまうので、水田以外をマスク(隠す)した上で、水田域におけるEVIの平均値を計算します。

筆ポリゴンの取得

マスク処理には、農林水産省が提供している農地の区画情報である「筆ポリゴン」を使用します。筆ポリゴンは、都道府県あるいは市町村単位で提供されているGISデータ(Shapefile)です。

本来ならば、解析する年ごとに筆ポリゴンを入手する必要がありますが、過去のアーカイブは提供されていません。今回は、過去10年間の農地の変化は誤差の範囲として割り切り、提供されている最新版の筆ポリゴンを使って過去のEVIもマスクすることにします。

農林水産省ウェブサイトの筆ポリゴンダウンロードページから、茨城県の筆ポリゴンデータをダウンロードしましょう。zip形式で圧縮されたファイルを解凍すると、市町村単位のフォルダでファイルが格納されています。

筆ポリゴンの編集

市町村単位で提供されている筆ポリゴンを結合し、茨城県全体のデータを作成しましょう。

また、ポリゴン情報の中には「田」だけでなく「畑」の情報も格納されています。今回は水稲の収量を予測するので、畑のポリゴンを取り除く処理をします。

今回は、フリーのGISソフトであるQGISによる手順を紹介します。QGISはこちらからダウンロードすることができます。

  1. QGISを開き、各フォルダに格納されているShapefile (shp) を読み込む
  2. メニューバーにある「ベクタ」から「データ管理ツール」を開き、「ベクタレイヤのマージ」を起動する
  3. 「ベクタレイヤのマージ」ウィンドウが開かれるため、以下のように設定する
  • 「入力レイヤ」⇒ 各市町村のシェープフィルをすべて選択する
  • 「変換先の座標参照系(CRS)」⇒ EPSG:32654 WGS84/UTM zone 54N
  • 「出力レイヤ」⇒ 新規ファイル名を設定(「…」から保存先のディレクトリを選択可能)
  1. 「実行」をクリックすると茨城県 筆ポリゴンが出力される
  2. 上記4で作成した筆ポリゴンをQGISで開き、ポリゴンファイルを選択した状態でプラグインツールバーにある編集モード切替を選択する
  3. 上記4で作成したポリゴンファイルを選択し「属性テーブル」を開く
  4. 属性テーブル内にある「耕地の種類」から、「畑」の属性を持つポリゴンを選択し、削除する
  5. 属性テーブルを閉じ、プラグインツールバーにある「レイヤ編集内容の保存」を選択し、編集を終える

市町村単位の筆ポリゴンを結合し、水田を対象とした茨城県全体の筆ポリゴンを作成できました。Google Earthに重ねて確認してみましょう。

筆_all.jpeg

拡大して見てみましょう。水田の一筆一筆が黄色のポリゴンで囲われていることがわかります。

筆.jpeg

水田以外のマスク処理

作成した茨城県の水田の筆ポリゴンを利用し、再度QGISを使ってEVIのマスク画像を作成しましょう。手順は以下の通りです。

  1. QGISを開き、作成した茨城県の筆ポリゴンとEVIのGeoTIFF画像を表示する
  2. QGISのメニューバーにある「ラスタ」から「ベクタのラスタ化」を起動する
  3. 「ベクタのラスタ化」ウィンドウが開くので、以下のように設定して「実行」をクリックする
  • 「入力レイヤ」⇒ 茨城県の筆ポリゴン
  • 「固定値」⇒ 1.000
  • 「出力ラスタサイズの単位」⇒ 地理単位
  • 「水平方向の解像度」「鉛直方向の解像度」⇒ 250.000
  • 「出力領域」⇒ 「…」から「レイヤの領域を使う」を選択し、EVIの画像を選択
  • 「出力データ型」⇒ Byteを選択
  • 「追加のコマンドラインパラメータ」⇒ 「-at」と入力
  • 「出力ファイル」⇒ 新規ファイル名を設定(「…」から保存先のディレクトリを選択可能)

そうすると、下図のようなEVIの筆マスク画像が出力されます。

圃場マスク.jpeg

EVIの集計とグラフ表示

筆マスク画像を用いてEVIを集計しましょう。筆ポリゴン内で集計したEVIの値を平均値化し、その年のEVIの代表値とします。2007年から2018年まで、各年のEVI値を計算してみましょう。

In [ ]:
# mask_pathは適宜、データが格納されているディレクトリに変更してください。
# raster_pathは適宜、データを保存したいディレクトリに変更してください。


mask_path = "/content/drive/MyDrive/fude_raster_markK.tif"
raster_path = "/content/average_EVI"

mask = gdal.Open(mask_path)
mask_a = mask.ReadAsArray()

#保存用ディレクトリからEVIのGeoTiff画像を選択
EVI_year_datas = sorted(glob.glob(os.path.join(raster_path,"*.tif")))

#集計した値を入れるリストを作成
EVI_all = []

#EVIの平均値を算出
for data in EVI_year_datas:
  EVI = gdal.Open(data)
  EVI_a = EVI.ReadAsArray()
  evi = EVI_a * mask_a
  evi_2 = np.mean(evi[evi!=0])
  EVI_all.append(evi_2)

print(EVI_all)

#データフレーム化
df_EVI_all_year = pd.Series(EVI_all,index=["2007","2008","2009","2010","2011","2012","2013","2014","2015","2016","2017","2018"])
df_EVI_final = pd.DataFrame(df_EVI_all_year,columns=["EVI"])

df_EVI_final
[3300.8870090273213, 3288.4697196149937, 3273.2556744584604, 3100.4551523484984, 3271.328670213826, 3648.2165959227596, 3634.982025069249, 3618.293318188159, 3128.3178493852256, 3157.609067357513, 3564.8704489747115, 3293.9973894501904]
Out[ ]:
EVI
2007 3300.887009
2008 3288.469720
2009 3273.255674
2010 3100.455152
2011 3271.328670
2012 3648.216596
2013 3634.982025
2014 3618.293318
2015 3128.317849
2016 3157.609067
2017 3564.870449
2018 3293.997389

コラム2:ポリゴンデータによる空間統計

本章では、QGISで処理をしたマスク画像を用いて茨城県の農地におけるEVI値を推定しました。 "rasterstas"というライブラリを用いると、筆ポリゴンなどのポリゴンファイルを用いて推定が可能です。

In [ ]:
###sample###

"""
#ライブラリのインストール
!pip install rasterstats

from rasterstats import zonal_stats
import geopandas as gpd
import numpy as np

#raster_pathは作成したEVIのデータ、vector_pathは集計用のポリゴンファイルを指定してください。
raster_path = "/content/drive/MyDrive/aaa/"
vector_path = "/content/drive/MyDrive/fude_poly_cent_point/rice_cen_poly_32654.shp"

list_average_EVI = glob.glob(os.path.join(raster_path,"*.tif"))

list_1 = []
list_2 = []

for l in list_average_EVI:
  stats = zonal_stats(vector_path,l,stats="mean",all_touched=True,geojson_out=True)
  geostats = gpd.GeoDataFrame.from_features(stats)
  arr = geostats.values
  for n in range(geostats.shape[0]):
    g = (arr[n,3])
    h = f.append(g)
    ave = sum(list_1)/len(list_1)
  list_2.append(ave)
  print(list_2)
"""

rasterstasを用いる際、筆ポリゴンのような地物数の多いポリゴンデータで統計値を取得すると、非常に長い計算時間がかってしまうため、本章では処理の高速化のためマスク画像を用いました。

筆ポリゴンのように地物数の多いポリゴンデータを取り扱う際には、前処理としてバッファ(ポリゴン等から特定の距離の範囲のエリアを作成する機能)処理を施した後、重複したポリゴン同士を結合させることにより集約化したポリゴンを作成することで、処理の負荷を軽減させることができます。

2007年から2018年までのEVIの値を算出しました。年ごとの変化がわかりやすいようにグラフで表示してみましょう。

In [ ]:
fig= plt.figure(figsize=(12,8))
df_EVI_final.plot(title="EVI 2007 to 2018",color='g',grid=True,label="EVI",xlabel = "year",ylabel = "EVI")
plt.legend()
plt.show()
<Figure size 864x576 with 0 Axes>

日射量がEVIに及ぼす影響

EVIから収量を予測する前に、気象が水田のEVIに与える影響について見てみましょう。日本の水田は灌漑が整備されているので、降水量は作物の生育にあまり影響を及ぼしません。そこで、日射量とEVIを比較してみることにします。

AMeDASデータの取得

日本の気象データは、気象庁の地域気象観測システム「AMeDAS(Automated Meteorological Data Acquisition System)」が取得しています。詳細については気象庁のサイトを参照してください。日射量は「日照時間」として6分単位で観測されています。

気象庁の過去の気象データ・ダウンロードページ(下図参照)で、茨城県の2007年から2018年までの日照時間のデータをCSV形式でダウンロードします。

AMeDASキャプチャ.PNG

AMeDASは観測地点によって観測している対象が異なるため、日照時間を記録している地点だけを選択してください。

また、今回は8月と9月のEVIデータを取得したので、1ヶ月遡って7月から9月までの日照時間を取得し、その平均値を求め、その年の日射量の代表値とします。

コラム3:AMeDASデータの品質

ダウンロードしたCSVファイルには、過去の観測値に加え、選択されたオプションの組み合わせによってデータの品質情報、現象なし情報、均質番号などの値が付与されます。

  • 品質情報:データが正常かどうか、欠損が無いか等を分類した値(8:正常値、5:準正常値、4:資料不足値、2:疑問値、1:欠測、0:観測・統計項目でない)
  • 現象なし情報:観測点に常駐している観測員が現象の有無を記録した情報(1:現象なし、0:現象あり)
  • 均質番号:データの均質性を表したもの(観測点ごとに、観測時期の古いものから、1から順番に付与される)

本章では、品質情報における正常値のみを使用して解析を行っています。

また、茨城県以外の場所で試す場合は、AMeDAS観測地点と筆ポリゴン(田)の分布が空間的に偏りがないようにデータを選択してください。

AMeDASデータの詳細についてはこちらのページを参照してください。

欠損値があるデータの削除

取得したデータAMeDASデータには欠損値が含まれる場合があるので、CSVファイルを調べ、欠損値がある場合にはそのデータを削除します。

In [ ]:
#csvデータの読み込み
#csv_pathは適宜ダウンロードしたcsvファイルが保存されいているディレクトリに変更してください。

csv_path = "/content/drive/MyDrive/data.csv"
#csv_path = "/content/drive/data.csv"
df_ame = pd.read_csv(csv_path,encoding = "shift-jis")

#一行目には不要な空行があるので削除
df_ame_csv =df_ame.drop(0)

#観測地点の名称と列番号をCSVファイルを開いて確認し、リスト化
ibaraki_A_list = ["date","大子","笠間","水戸","古河","下妻","土浦","鹿嶋","日立","龍ヶ崎","筑波山","鉾田","常陸大宮","つくば(館野)","下館","北茨城"]
ibaraki_B_list = [0,1,4,7,11,14,17,20,23,26,29,32,35,38,42,45]

print("done")
done
In [ ]:
#CSVファイルから必要な要素を抜出して新しい表を作成
arr = []
for i in range(len(df_ame.count(1))-1):
    tmp = []
    for l in ibaraki_B_list:
        h = df_ame_csv.iloc[i,l]
        tmp.append(h)
    arr.append(tmp)

print("done")
done

取得したCSVファイルをデータフレーム形式で表示してみましょう。

In [ ]:
pd.set_option("display.max_rows",150)

df_ame_csv_2 = pd.DataFrame(arr)
df_ame_csv_2.set_index(0)
Out[ ]:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
0
NaN 大子 笠間 水戸 古河 下妻 土浦 鹿嶋 日立 龍ケ崎 筑波山 鉾田 常陸大宮 つくば(館野) 下館 北茨城
年月 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間)
NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
7-Jan 198.6 166.3 184.3 213.1 187.2 177.6 179.3 196.8 174.5 NaN 162.4 184.1 194.4 199.7 202.8
7-Feb 200.7 189.9 195.7 209.5 201.4 195.8 197.7 204.6 194.7 NaN 193.3 199.3 206.1 210.5 195.4
7-Mar 190.9 190.5 192.8 232.7 208.4 201.9 192.9 186.1 189.9 NaN 191 189.1 208.7 218.7 199.1
7-Apr 147 141.1 155.9 166.6 142.3 147.1 156 157.9 143.3 NaN 160.5 142.2 158.5 162.9 164.5
7-May 205.1 192.7 220.7 217.8 185.8 199.1 231.2 180 197 NaN 224 190.3 226.5 235.7 215.7
7-Jun 154.4 129.2 184.3 142.2 104.3 136.9 204.3 130.1 134 NaN 185.2 127.5 173.2 174.6 192.1
7-Jul 77.8 66.4 101.1 63.3 37.2 57.9 94.1 64.7 60.6 NaN 82 54.8 78.3 79.1 94.6
7-Aug 169.4 200.7 228.4 209.6 168.6 192.9 250.4 191.3 196.1 NaN 224.9 165.5 232.9 237.9 229.3
7-Sep 114.8 117.3 122 125 110.5 125.3 156.2 135.2 127.2 NaN 141.4 113.9 135.4 118.8 150.4
7-Oct 135.6 144.5 145 162.2 146.8 132.1 140.8 152.6 129.2 NaN 143.6 144.9 135.6 148.3 153.8
7-Nov 134.6 139.7 134.6 160.5 150.8 145.9 139.6 136.5 146.1 NaN 128.9 135 153.4 157.5 132.5
7-Dec 148.3 153.2 155.4 191 173.1 151.3 159.3 168.3 163.9 NaN 142.9 160.2 179.5 184.3 149.9
8-Jan 189.2 178.4 180.9 185.5 173.8 160.9 169.1 179 177.8 NaN 161.2 171.3 181.1 196.6 189.2
8-Feb 190.5 197.4 200.7 224.1 219.3 215.2 209 200.1 215 NaN 184.9 196.5 220.6 225.5 186.8
8-Mar 180.8 186.2 181.4 192.1 188.8 195.4 191.2 178 175.8 NaN 194.3 182.4 189.8 191.2 173.8
8-Apr 157.9 156.7 164.5 165.6 161.8 166.2 156.6 160.8 153.7 NaN 162.6 152.8 157.5 169.8 166.1
8-May 147.5 137.9 158.8 149.3 144.8 150.4 150.5 158.4 143.7 NaN 151.6 156.5 146.4 150.1 158.2
8-Jun 117 101.4 135.5 113.2 110.6 125.6 156.2 163.7 127.5 NaN 140.9 129.9 122.5 123.3 152.1
8-Jul 110.6 96.1 135.7 138.4 140.1 137 178.7 139.9 161.3 NaN 151.4 121.3 139.5 143 147
8-Aug 106.2 102.1 122.7 126.3 140.3 132.7 166.3 124.4 147.2 NaN 143.4 117.9 134.7 138.5 115.4
8-Sep 131 128.2 137.2 131 134.3 126.9 130.3 140.6 120.5 NaN 123.7 138.9 125.4 133.6 143.9
8-Oct 137.8 146.6 152.7 153 150 148.4 143.5 158.4 142.9 NaN 147.1 150.3 145.4 154.1 158.1
8-Nov 129 131.1 141.1 151 140.8 134.3 124.3 137.8 135.7 NaN 106.2 127.1 139.2 146.2 115.3
8-Dec 157.8 165.5 189.1 201.3 189.2 179.4 186.7 184.5 193.1 NaN 161 172.1 190.7 179 164.4
9-Jan 169.7 166.1 181.7 187 182.2 169.5 163.1 169.5 171 NaN 140.1 160.6 174.6 184.1 169.1
9-Feb 153 125.5 152.7 167.8 149.4 141.1 137.7 162.6 137.1 NaN 134 157 148.4 155.7 168.3
9-Mar 186.1 172.9 179.3 189.8 185.7 174.7 164.9 176.2 170.6 NaN 167.3 181.5 176.5 189.1 195.1
9-Apr 201.3 206.8 220.1 224.9 228.3 228.5 224.1 216.4 226.5 NaN 217.4 208.2 232.7 231 214.6
9-May 158.4 156.3 168.7 158.8 154.9 163.2 161 172 160.7 NaN 165 156.7 164.9 166.4 171.5
9-Jun 101.2 93.6 105.3 97.8 94.8 98.7 93.3 99.3 97.2 NaN 103.5 94.1 95 103.3 108.9
9-Jul 82.1 78.1 99.3 83 86 102.2 123.1 99.4 115.2 NaN 102.9 80.8 107.4 87.8 103.2
9-Aug 120.1 128 159 131.9 144.8 153.1 161.1 139.2 157.3 NaN 154.3 121 148.3 162.6 125.6
9-Sep 137.4 134.1 154.4 144.8 139.8 142 139.7 150.3 141.9 NaN 144.3 138.4 137.7 145.1 147.6
9-Oct 139.2 141.1 151.7 161.4 153.3 156.7 155.9 124.3 156.2 NaN 135.8 146.5 161.2 155.2 145.5
9-Nov 108.2 111.7 125.9 129.6 120.1 122.8 123 122.9 117.6 NaN 107.8 113.2 122.4 124.7 124
9-Dec 169.7 167.7 169.9 190.7 180.7 162.2 170.6 180.4 181.9 NaN 144.5 167 185.6 181.5 167.7
10-Jan 193.5 211.1 218.9 234.9 230.1 227.2 197.8 206.4 212.1 NaN 180.4 196.4 231.3 234.5 188.9
10-Feb 119.5 122.7 121.8 128.9 118.8 119.8 106.6 69.7 108.2 NaN 98.6 118.8 118.4 126.1 115.9
10-Mar 144.5 130.4 133.6 138.5 133.9 132.7 130.2 136.8 130.4 NaN 133 138.4 128.9 139.2 140.1
10-Apr 116.9 121.1 129.6 135.9 132.3 133 137.1 130.9 134.4 NaN 138.1 113.4 132.5 121.5 126.3
10-May 180.8 198.4 209 208.9 202.8 200.9 213.7 201.3 206.4 NaN 210.8 190.8 200.1 205.5 200.3
10-Jun 153.5 156.2 162.8 161.3 159.4 154.2 160.9 166.3 166.7 NaN 158.8 154.3 158.8 158.6 181.3
10-Jul 129.7 161.7 179.5 166.4 179.8 183.2 192.9 175.4 206.7 NaN 196.7 152 186.2 177.3 174.2
10-Aug 156.5 196.8 214 187 212.7 222.8 261.8 215.9 244.6 NaN 219.1 176.5 233.3 209.2 192.5
10-Sep 143.5 155.4 175.9 161 158.4 175.1 173.3 171.4 178.5 NaN 173.2 151.7 168.5 167.7 161.6
10-Oct 101.5 102.6 117.2 95.5 95.2 108.4 97 110.9 100.2 NaN 99.2 104.8 104.4 102 102.3
10-Nov 148 163.9 158.5 173 164.9 154.7 144.5 155.1 146.1 NaN 143.5 148 153 167.4 135.2
10-Dec 165 200.7 192.1 211.6 212.3 200.2 193.1 185.9 197.6 NaN 175 175.6 202.2 208.2 166.6
11-Jan 208.8 244.9 240.8 264.2 263.8 255.9 246.5 225.4 249.8 NaN 229 217.6 259.1 257.1 210.4
11-Feb 162.6 161.4 155.8 171.2 163 149.3 144.9 160.5 152.7 NaN 143.7 154.6 156.6 161.9 160.2
11-Mar 198.5 212.7 218.5 246 232.9 222.8 223.9 203.5 222.8 NaN 204.9 191.8 224.9 229.3 220.4
11-Apr 188.5 195.4 182.1 203.4 195 198.7 213.4 207.1 208.6 NaN 209.7 184.9 198.5 194.8 213.8
11-May 135.5 138.2 143.5 152.5 147.9 144.9 149.6 151.3 146.5 NaN 146.8 140.6 138 136.9 151.1
11-Jun 109.7 107.7 122.3 114.1 121.2 117.1 118.4 126.6 114.1 NaN 118.5 101.5 117.2 118.1 143.2
11-Jul 150.8 173.1 192.2 179 185.1 179.8 205.7 194.6 188.4 NaN 188.4 171.6 183.5 183.4 198.4
11-Aug 136.5 155.1 176.7 157.3 170.3 187.1 204.2 181.6 192.9 NaN 190 151.9 188.3 173 158.6
11-Sep 148.4 169.6 177.7 178.2 173.2 174.4 181.3 172.4 183.4 NaN 178.7 160.9 186.6 178.4 167
11-Oct 125.9 131.5 135.4 149.4 141.4 140.5 150.5 133.7 140.6 NaN 138.5 130.6 142.5 145.8 120.9
11-Nov 137.6 149.1 150.4 170.1 163.5 151.9 142.9 138 147.4 NaN 142 142.5 153.2 164.1 122.5
11-Dec 165.2 193.8 186.1 198.9 197.5 192.9 185.8 177.7 189.1 NaN 186.6 174.1 191.9 194.4 160.8
12-Jan 184.6 202.8 199.8 213.8 207.9 202.7 187.7 183.5 189.7 NaN 189.6 179.1 206.2 207.7 189.3
12-Feb 170.3 172.9 158.5 178.7 168.3 157.1 152.3 159 155.4 NaN 157.4 158.1 161.1 172.6 163.5
12-Mar 164.7 156.4 162 171.1 163.9 150.4 141.5 147.6 122.9 NaN 138.6 159.7 155.9 156.2 168.4
12-Apr 158.6 158.7 174.9 172.7 168.8 168.5 172.8 174.7 166.9 NaN 172.1 159.6 174 165.1 174.4
12-May 181.2 195.4 206.7 217.2 201.7 194.7 184.9 202.8 200.2 NaN 206.3 182.2 200.2 201.7 161.6
12-Jun 138 154.4 178.4 134.7 154.6 159.6 154.4 177.9 153.7 NaN 160.4 141.9 162.7 147.4 166.2
12-Jul 123 151.9 179.3 159.8 178.4 181.1 185.8 175.9 177.2 NaN 181.3 145.9 185.9 170.9 166.2
12-Aug 193.8 230.7 267.6 237.3 263.4 266.3 288 271.3 270.9 NaN 266.6 211.9 266.7 261.2 281.9
12-Sep 136.8 150.9 165.4 162.4 173.2 173.1 187.8 181.4 171.3 NaN 167.6 149.8 170.6 175.1 181.5
12-Oct 175.9 164.7 177.3 176.9 171.3 174.6 162.7 179.1 167.9 NaN 169.1 174.9 171.5 173.7 172.2
12-Nov 135.4 162.4 162.6 164.9 161.1 162.4 153.9 157.3 160.1 NaN 157.8 136.8 162.2 164 157.1
12-Dec 157.3 174.2 169.8 187.6 179.9 176.6 171.4 168.3 170.3 NaN 167.9 159.7 177.3 173.1 162.6
13-Jan 205.9 230.2 223.3 245.3 231.5 218.7 200 223.2 207.7 NaN 207.8 209.3 218.7 233.3 218.7
13-Feb 173.3 183.3 183 194.2 190.1 184 174.8 182.3 169.5 NaN 173.6 176.2 181.6 189.6 175.5
13-Mar 178 182 184.2 200.1 187.5 183.3 169.4 187.3 177.7 NaN 175.5 179.9 187.1 192.7 191.1
13-Apr 177.4 190.5 198.3 197.6 201.2 200.8 192.8 194.4 207 NaN 186.3 193 202.4 204 182.4
13-May 209.7 236.4 246.5 249.8 243.7 237 244.4 244.3 250.1 NaN 246.1 217.3 243.9 244.4 247.4
13-Jun 125.1 114.9 133.3 119.4 123.5 124.4 125.3 125.9 118.6 NaN 122.1 123.7 126.1 126.9 122.3
13-Jul 112.3 141.7 160.2 160.6 167.9 174.5 162.4 152 175.5 NaN 162.2 123.4 170.7 157.5 132.6
13-Aug 170.2 198.6 225.5 209.5 221 218.9 238.7 225.6 239.5 NaN 233.9 188.6 221.4 208.5 220.9
13-Sep 142.6 161.9 185 161.6 168.2 177.7 189.6 181.8 172.3 NaN 191 157 172.3 173.9 162.2
13-Oct 107.1 118 129.8 127.4 122.7 124 122.9 112.6 122.7 NaN 125.5 116.5 122.6 123.1 115.4
13-Nov 163 186.4 189.9 193.9 177.3 184.3 177.6 184.2 175.4 NaN 181.3 166.3 185.3 187 177.7
13-Dec 169.6 203.1 200 207 202.8 194.6 186.4 195.9 187.3 NaN 189.3 182.7 192.7 200.8 187.5
14-Jan 211.2 225.3 216.6 237.6 223.7 213.8 206.6 212.5 206.6 NaN 206.4 203.5 220.8 224.5 223.1
14-Feb 158 163.5 156.7 185.2 173.1 156 142.3 147.3 122 NaN 152.1 157.9 158 178.6 96.8
14-Mar 203.8 194.5 209.3 220.6 214.2 209.1 201.8 213.7 209.4 NaN 200.7 209.7 213.1 216.4 213.3
14-Apr 222.7 222.1 240.7 239.5 235.3 228.8 229.4 238.8 231.8 NaN 234.5 225.8 227.7 237.1 243.6
14-May 196.7 224.4 245.7 245.9 242.9 232.1 238.8 230.2 235.7 NaN 230.8 206 240.3 240.5 222.4
14-Jun 124.9 124.5 147.9 145.1 151.6 150.3 161.7 155.2 157.9 NaN 155.9 124.4 151.3 149.6 155.4
14-Jul 135.3 164.8 196 176.6 190.6 189.8 208.8 184.2 199.3 NaN 201.2 155.8 188.5 187.2 179.9
14-Aug 127.1 153.8 174.1 157.6 163.5 167.7 182.3 182.6 185.6 NaN 175.5 155.6 174.8 162.5 174
14-Sep 144.9 140.2 163.2 161 163.5 159.6 159.9 154.9 161.2 NaN 164 140.9 166 154.6 152.1
14-Oct 142.6 136 160 158.6 153.7 143.3 151.1 168.5 146.7 NaN 151.8 142.9 144.1 148.1 174.5
14-Nov 131.5 129.2 145.1 145 138 134.1 124.4 148.6 134.6 NaN 128 133.7 130.4 136.8 142.9
14-Dec 162.6 186.8 194.3 200.6 201.3 198.9 187 193.2 189.7 NaN 194.4 190.9 196.2 193.5 188.4
15-Jan 188.3 180.1 198 213 207.2 202.5 193.9 203.8 196.4 NaN 194.8 192.8 202.6 200.5 199.4
15-Feb 171.6 133.2 171.9 187.2 182 169.6 167.7 161.3 166.7 NaN 163.1 164.4 172.2 177.9 164.3
15-Mar 192.3 187 206.1 202.6 202 202.7 196.9 210.6 198.4 NaN 202.8 200.1 200.6 203.5 210.7
15-Apr 148.8 164.9 174 158 155.5 163 166.4 166.2 161.1 NaN 169.1 154 161.6 160.7 171.6
15-May 234.8 243.2 256.6 244.4 241.4 244.9 236.9 256.1 242.4 NaN 248.1 242 245.7 244.8 259.3
15-Jun 121.6 124.8 153.4 130.8 133.8 145.4 161 158 154.7 NaN 159.7 124.3 134.6 142.9 140.1
15-Jul 147.8 161.3 181 170.5 172.2 177.6 191.6 180.5 192.7 NaN 178.5 162.4 185.5 178.1 173.5
15-Aug 118.3 126.9 151.7 138.2 143.1 142.7 143.5 152.1 157.6 NaN 137.9 127 146.1 148 139.8
15-Sep 94.1 104.8 119.8 118 120.3 114.8 129.3 119.1 122.8 NaN 124 101.8 117.4 119.2 113.7
15-Oct 173.9 181 198.4 195.3 193.4 199.3 182 189.5 193.3 NaN 199.4 182.9 199 197.3 183.9
15-Nov 109.8 115.8 129.8 126 114.1 120.9 117.8 131.7 121.8 NaN 123.7 113.8 118.7 115.3 125.7
15-Dec 170.1 172.7 176.9 193.8 185.4 171.5 158.4 178.9 159 NaN 164.6 179.4 174.9 184.5 178.6
16-Jan 175.1 189.6 203.2 222.6 210.5 209.5 205.3 189.5 202.5 NaN 203.3 185.9 211.1 209.5 182.8
16-Feb 179.4 163.4 177.1 187.8 181.1 173.3 167.6 177.2 172.8 NaN 168 170.4 179.5 179.3 189
16-Mar 170.2 159.1 179.3 182.9 170.2 167 158.2 176.6 156.3 NaN 167.3 165.5 172.8 175 190.4
16-Apr 147.2 140.4 160.5 156.5 153.1 142.2 154.3 162.1 151.5 NaN 149.2 143.6 156.4 152 182.5
16-May 171.3 191.7 197 196.8 202.6 197.9 203.4 196.6 206.6 NaN 193.8 183.8 202.3 201.9 199
16-Jun 138.5 144.3 155.7 147.3 148 143.1 146.4 158.9 141.2 NaN 145.1 145.7 140.9 152.5 157
16-Jul 102.3 114.7 167.2 121.1 149.3 151 171 154.2 180.5 NaN 171.8 109.5 146.1 135.2 141.9
16-Aug 152.1 156.3 191.3 167.7 186.2 181.5 182.1 201.2 189.6 NaN 188.7 150.3 176.8 178.6 197
16-Sep 90.5 87.4 112.9 92.9 98.9 97 93 115.7 81.7 NaN 99.2 95.5 100.1 102.9 106.7
16-Oct 151.3 143.7 157.5 148.4 144.7 142.3 141.8 170.9 141 NaN 146.6 149.2 140.5 149.3 176.3
16-Nov 141 138.9 154.3 151.8 147.7 143 128.6 155.6 131.7 NaN 137.6 139.3 141.1 146.4 162.5
16-Dec 183.1 191.3 205.6 217.1 207.2 206.4 197 204.5 188.2 NaN 198.8 196.3 209.9 205.5 195.6
17-Jan 195.6 208.6 223.6 234.2 233.7 231.3 227.6 215.5 234.2 NaN 219.1 206 233.9 227.8 208.1
17-Feb 209.5 200.7 208 232.9 221.2 201.7 194.9 203.4 202.2 NaN 195.3 206.1 208.5 216.6 201.6
17-Mar 172.1 184.8 198.5 208.6 203 209.7 195.5 191.6 203.8 NaN 194.2 181.6 206.9 199.2 190.9
17-Apr 170.1 190.4 210.1 217.6 210.7 199.4 193 198.1 202.7 NaN 204.7 181.9 201 215.8 200.7
17-May 169.6 185.2 211.9 197.9 196.2 206.6 220.6 211 214.9 NaN 209.2 180.6 205.9 197.3 197.3
17-Jun 139.8 157.3 174.3 179.1 176.6 171.8 164.2 177.1 169.6 NaN 171.4 148.3 170.1 177.4 172.4
17-Jul 128.9 154.1 174 152.2 175.2 180.8 201.9 178.3 200.9 NaN 180 147.7 184.3 171 175.9
17-Aug 57.2 69.3 100 74.6 84.8 92.8 108.6 94.8 101.3 NaN 107.4 60.5 95 87.4 79.4
17-Sep 143.2 138.9 157.8 144.6 145.5 141.4 147.3 156.7 127.6 NaN 152.7 144.1 134.7 150.3 158.2
17-Oct 100 91.3 106.3 112.1 111.6 101.3 97.4 105.8 102.4 NaN 100.8 100.1 104.9 103.5 102.8
17-Nov 154.6 158.1 166.8 180 168.6 159.3 157.4 166.5 163 NaN 162.8 154.8 162.8 171.2 160.6
17-Dec 194.4 197.3 213.4 230.4 221.9 215.4 195.2 204.4 209.1 NaN 201.3 200 217.4 219.1 198.3
18-Jan 192.5 207.8 219.1 224.2 224.8 220.3 200.2 214.6 212.4 NaN 210.3 206.5 220.1 224 212.6
18-Feb 175.7 177.9 195.6 195 182.4 179.7 165.6 190.7 166.9 NaN 180.2 175.3 178.1 186.6 183.2
18-Mar 191.7 176.9 207.9 216 207.1 203 194.1 199 203.3 NaN 206.4 193.1 199.3 205.2 209
18-Apr 179.8 193.1 200.2 204.8 202.6 197.8 202.1 199.9 196.1 NaN 188.9 192.4 204.6 197.4 194.8
18-May 182.6 185.4 202.2 197.8 198.9 195 192.2 192.4 188.4 NaN 190.1 186.5 196.1 198.3 195.9
18-Jun 148.9 163.2 179.8 172 175.9 174.5 178.1 183.6 181.1 NaN 174.2 156.6 167.7 167.4 189.3
18-Jul 167.2 192.1 219.4 201.9 220.2 221.8 234 227.8 252.7 NaN 222.7 189.6 226.2 214.5 203.6
18-Aug 154.6 179.1 199.2 204.3 202.9 206.2 213 195.2 226.9 NaN 210.1 171.3 211.2 206.7 185.3
18-Sep 70.7 71.6 92.4 76.1 82.8 97.1 108.2 86.9 95.9 NaN 107.9 70.6 91.5 84.5 72.6
18-Oct 138.6 140.5 156.3 147.7 154.4 150.6 137.8 158.8 145.7 NaN 152.4 142.1 148.6 151.4 160.2
18-Nov 156.3 147.2 171.2 165 160 162.6 149.4 169.4 153.3 NaN 153.8 150.5 156.4 158.6 175
18-Dec 150.4 143.3 155.8 174.4 164.5 155 147.1 161.8 151.1 NaN 153.3 160.4 155.6 161.7 165.5

データフレームの一列目には観測年と月が"年-月"の形式で表記されており、各列には観測地点と、各年の月ごとの平均値が入っています。データフレーム中にある"NaN"は数値が入っていないことを表しています。この要因として、観測や統計を行っていないこと、正常にデータが取得されなかったことなどが挙げられます。

作成したデータフレームを見ると10列目(筑波山の列)と2行目がすべてNaNとなっています。

このままだと、正しい数値が算出できないため、10列目および2行目を削除します。

In [ ]:
#10列目の削除
df_ame_csv_3 =df_ame_csv_2.drop(columns=df_ame_csv_2.columns[10])

#2行目の削除
df_ame_csv_4 =df_ame_csv_3.drop(df_ame_csv_2.index[2])

print("done")
done

編集したAMeDASのデータを表示します。

In [ ]:
df_ame_csv_5=df_ame_csv_4.set_index(0)
df_ame_csv_5
Out[ ]:
1 2 3 4 5 6 7 8 9 11 12 13 14 15
0
NaN 大子 笠間 水戸 古河 下妻 土浦 鹿嶋 日立 龍ケ崎 鉾田 常陸大宮 つくば(館野) 下館 北茨城
年月 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間) 日照時間(時間)
7-Jan 198.6 166.3 184.3 213.1 187.2 177.6 179.3 196.8 174.5 162.4 184.1 194.4 199.7 202.8
7-Feb 200.7 189.9 195.7 209.5 201.4 195.8 197.7 204.6 194.7 193.3 199.3 206.1 210.5 195.4
7-Mar 190.9 190.5 192.8 232.7 208.4 201.9 192.9 186.1 189.9 191 189.1 208.7 218.7 199.1
7-Apr 147 141.1 155.9 166.6 142.3 147.1 156 157.9 143.3 160.5 142.2 158.5 162.9 164.5
7-May 205.1 192.7 220.7 217.8 185.8 199.1 231.2 180 197 224 190.3 226.5 235.7 215.7
7-Jun 154.4 129.2 184.3 142.2 104.3 136.9 204.3 130.1 134 185.2 127.5 173.2 174.6 192.1
7-Jul 77.8 66.4 101.1 63.3 37.2 57.9 94.1 64.7 60.6 82 54.8 78.3 79.1 94.6
7-Aug 169.4 200.7 228.4 209.6 168.6 192.9 250.4 191.3 196.1 224.9 165.5 232.9 237.9 229.3
7-Sep 114.8 117.3 122 125 110.5 125.3 156.2 135.2 127.2 141.4 113.9 135.4 118.8 150.4
7-Oct 135.6 144.5 145 162.2 146.8 132.1 140.8 152.6 129.2 143.6 144.9 135.6 148.3 153.8
7-Nov 134.6 139.7 134.6 160.5 150.8 145.9 139.6 136.5 146.1 128.9 135 153.4 157.5 132.5
7-Dec 148.3 153.2 155.4 191 173.1 151.3 159.3 168.3 163.9 142.9 160.2 179.5 184.3 149.9
8-Jan 189.2 178.4 180.9 185.5 173.8 160.9 169.1 179 177.8 161.2 171.3 181.1 196.6 189.2
8-Feb 190.5 197.4 200.7 224.1 219.3 215.2 209 200.1 215 184.9 196.5 220.6 225.5 186.8
8-Mar 180.8 186.2 181.4 192.1 188.8 195.4 191.2 178 175.8 194.3 182.4 189.8 191.2 173.8
8-Apr 157.9 156.7 164.5 165.6 161.8 166.2 156.6 160.8 153.7 162.6 152.8 157.5 169.8 166.1
8-May 147.5 137.9 158.8 149.3 144.8 150.4 150.5 158.4 143.7 151.6 156.5 146.4 150.1 158.2
8-Jun 117 101.4 135.5 113.2 110.6 125.6 156.2 163.7 127.5 140.9 129.9 122.5 123.3 152.1
8-Jul 110.6 96.1 135.7 138.4 140.1 137 178.7 139.9 161.3 151.4 121.3 139.5 143 147
8-Aug 106.2 102.1 122.7 126.3 140.3 132.7 166.3 124.4 147.2 143.4 117.9 134.7 138.5 115.4
8-Sep 131 128.2 137.2 131 134.3 126.9 130.3 140.6 120.5 123.7 138.9 125.4 133.6 143.9
8-Oct 137.8 146.6 152.7 153 150 148.4 143.5 158.4 142.9 147.1 150.3 145.4 154.1 158.1
8-Nov 129 131.1 141.1 151 140.8 134.3 124.3 137.8 135.7 106.2 127.1 139.2 146.2 115.3
8-Dec 157.8 165.5 189.1 201.3 189.2 179.4 186.7 184.5 193.1 161 172.1 190.7 179 164.4
9-Jan 169.7 166.1 181.7 187 182.2 169.5 163.1 169.5 171 140.1 160.6 174.6 184.1 169.1
9-Feb 153 125.5 152.7 167.8 149.4 141.1 137.7 162.6 137.1 134 157 148.4 155.7 168.3
9-Mar 186.1 172.9 179.3 189.8 185.7 174.7 164.9 176.2 170.6 167.3 181.5 176.5 189.1 195.1
9-Apr 201.3 206.8 220.1 224.9 228.3 228.5 224.1 216.4 226.5 217.4 208.2 232.7 231 214.6
9-May 158.4 156.3 168.7 158.8 154.9 163.2 161 172 160.7 165 156.7 164.9 166.4 171.5
9-Jun 101.2 93.6 105.3 97.8 94.8 98.7 93.3 99.3 97.2 103.5 94.1 95 103.3 108.9
9-Jul 82.1 78.1 99.3 83 86 102.2 123.1 99.4 115.2 102.9 80.8 107.4 87.8 103.2
9-Aug 120.1 128 159 131.9 144.8 153.1 161.1 139.2 157.3 154.3 121 148.3 162.6 125.6
9-Sep 137.4 134.1 154.4 144.8 139.8 142 139.7 150.3 141.9 144.3 138.4 137.7 145.1 147.6
9-Oct 139.2 141.1 151.7 161.4 153.3 156.7 155.9 124.3 156.2 135.8 146.5 161.2 155.2 145.5
9-Nov 108.2 111.7 125.9 129.6 120.1 122.8 123 122.9 117.6 107.8 113.2 122.4 124.7 124
9-Dec 169.7 167.7 169.9 190.7 180.7 162.2 170.6 180.4 181.9 144.5 167 185.6 181.5 167.7
10-Jan 193.5 211.1 218.9 234.9 230.1 227.2 197.8 206.4 212.1 180.4 196.4 231.3 234.5 188.9
10-Feb 119.5 122.7 121.8 128.9 118.8 119.8 106.6 69.7 108.2 98.6 118.8 118.4 126.1 115.9
10-Mar 144.5 130.4 133.6 138.5 133.9 132.7 130.2 136.8 130.4 133 138.4 128.9 139.2 140.1
10-Apr 116.9 121.1 129.6 135.9 132.3 133 137.1 130.9 134.4 138.1 113.4 132.5 121.5 126.3
10-May 180.8 198.4 209 208.9 202.8 200.9 213.7 201.3 206.4 210.8 190.8 200.1 205.5 200.3
10-Jun 153.5 156.2 162.8 161.3 159.4 154.2 160.9 166.3 166.7 158.8 154.3 158.8 158.6 181.3
10-Jul 129.7 161.7 179.5 166.4 179.8 183.2 192.9 175.4 206.7 196.7 152 186.2 177.3 174.2
10-Aug 156.5 196.8 214 187 212.7 222.8 261.8 215.9 244.6 219.1 176.5 233.3 209.2 192.5
10-Sep 143.5 155.4 175.9 161 158.4 175.1 173.3 171.4 178.5 173.2 151.7 168.5 167.7 161.6
10-Oct 101.5 102.6 117.2 95.5 95.2 108.4 97 110.9 100.2 99.2 104.8 104.4 102 102.3
10-Nov 148 163.9 158.5 173 164.9 154.7 144.5 155.1 146.1 143.5 148 153 167.4 135.2
10-Dec 165 200.7 192.1 211.6 212.3 200.2 193.1 185.9 197.6 175 175.6 202.2 208.2 166.6
11-Jan 208.8 244.9 240.8 264.2 263.8 255.9 246.5 225.4 249.8 229 217.6 259.1 257.1 210.4
11-Feb 162.6 161.4 155.8 171.2 163 149.3 144.9 160.5 152.7 143.7 154.6 156.6 161.9 160.2
11-Mar 198.5 212.7 218.5 246 232.9 222.8 223.9 203.5 222.8 204.9 191.8 224.9 229.3 220.4
11-Apr 188.5 195.4 182.1 203.4 195 198.7 213.4 207.1 208.6 209.7 184.9 198.5 194.8 213.8
11-May 135.5 138.2 143.5 152.5 147.9 144.9 149.6 151.3 146.5 146.8 140.6 138 136.9 151.1
11-Jun 109.7 107.7 122.3 114.1 121.2 117.1 118.4 126.6 114.1 118.5 101.5 117.2 118.1 143.2
11-Jul 150.8 173.1 192.2 179 185.1 179.8 205.7 194.6 188.4 188.4 171.6 183.5 183.4 198.4
11-Aug 136.5 155.1 176.7 157.3 170.3 187.1 204.2 181.6 192.9 190 151.9 188.3 173 158.6
11-Sep 148.4 169.6 177.7 178.2 173.2 174.4 181.3 172.4 183.4 178.7 160.9 186.6 178.4 167
11-Oct 125.9 131.5 135.4 149.4 141.4 140.5 150.5 133.7 140.6 138.5 130.6 142.5 145.8 120.9
11-Nov 137.6 149.1 150.4 170.1 163.5 151.9 142.9 138 147.4 142 142.5 153.2 164.1 122.5
11-Dec 165.2 193.8 186.1 198.9 197.5 192.9 185.8 177.7 189.1 186.6 174.1 191.9 194.4 160.8
12-Jan 184.6 202.8 199.8 213.8 207.9 202.7 187.7 183.5 189.7 189.6 179.1 206.2 207.7 189.3
12-Feb 170.3 172.9 158.5 178.7 168.3 157.1 152.3 159 155.4 157.4 158.1 161.1 172.6 163.5
12-Mar 164.7 156.4 162 171.1 163.9 150.4 141.5 147.6 122.9 138.6 159.7 155.9 156.2 168.4
12-Apr 158.6 158.7 174.9 172.7 168.8 168.5 172.8 174.7 166.9 172.1 159.6 174 165.1 174.4
12-May 181.2 195.4 206.7 217.2 201.7 194.7 184.9 202.8 200.2 206.3 182.2 200.2 201.7 161.6
12-Jun 138 154.4 178.4 134.7 154.6 159.6 154.4 177.9 153.7 160.4 141.9 162.7 147.4 166.2
12-Jul 123 151.9 179.3 159.8 178.4 181.1 185.8 175.9 177.2 181.3 145.9 185.9 170.9 166.2
12-Aug 193.8 230.7 267.6 237.3 263.4 266.3 288 271.3 270.9 266.6 211.9 266.7 261.2 281.9
12-Sep 136.8 150.9 165.4 162.4 173.2 173.1 187.8 181.4 171.3 167.6 149.8 170.6 175.1 181.5
12-Oct 175.9 164.7 177.3 176.9 171.3 174.6 162.7 179.1 167.9 169.1 174.9 171.5 173.7 172.2
12-Nov 135.4 162.4 162.6 164.9 161.1 162.4 153.9 157.3 160.1 157.8 136.8 162.2 164 157.1
12-Dec 157.3 174.2 169.8 187.6 179.9 176.6 171.4 168.3 170.3 167.9 159.7 177.3 173.1 162.6
13-Jan 205.9 230.2 223.3 245.3 231.5 218.7 200 223.2 207.7 207.8 209.3 218.7 233.3 218.7
13-Feb 173.3 183.3 183 194.2 190.1 184 174.8 182.3 169.5 173.6 176.2 181.6 189.6 175.5
13-Mar 178 182 184.2 200.1 187.5 183.3 169.4 187.3 177.7 175.5 179.9 187.1 192.7 191.1
13-Apr 177.4 190.5 198.3 197.6 201.2 200.8 192.8 194.4 207 186.3 193 202.4 204 182.4
13-May 209.7 236.4 246.5 249.8 243.7 237 244.4 244.3 250.1 246.1 217.3 243.9 244.4 247.4
13-Jun 125.1 114.9 133.3 119.4 123.5 124.4 125.3 125.9 118.6 122.1 123.7 126.1 126.9 122.3
13-Jul 112.3 141.7 160.2 160.6 167.9 174.5 162.4 152 175.5 162.2 123.4 170.7 157.5 132.6
13-Aug 170.2 198.6 225.5 209.5 221 218.9 238.7 225.6 239.5 233.9 188.6 221.4 208.5 220.9
13-Sep 142.6 161.9 185 161.6 168.2 177.7 189.6 181.8 172.3 191 157 172.3 173.9 162.2
13-Oct 107.1 118 129.8 127.4 122.7 124 122.9 112.6 122.7 125.5 116.5 122.6 123.1 115.4
13-Nov 163 186.4 189.9 193.9 177.3 184.3 177.6 184.2 175.4 181.3 166.3 185.3 187 177.7
13-Dec 169.6 203.1 200 207 202.8 194.6 186.4 195.9 187.3 189.3 182.7 192.7 200.8 187.5
14-Jan 211.2 225.3 216.6 237.6 223.7 213.8 206.6 212.5 206.6 206.4 203.5 220.8 224.5 223.1
14-Feb 158 163.5 156.7 185.2 173.1 156 142.3 147.3 122 152.1 157.9 158 178.6 96.8
14-Mar 203.8 194.5 209.3 220.6 214.2 209.1 201.8 213.7 209.4 200.7 209.7 213.1 216.4 213.3
14-Apr 222.7 222.1 240.7 239.5 235.3 228.8 229.4 238.8 231.8 234.5 225.8 227.7 237.1 243.6
14-May 196.7 224.4 245.7 245.9 242.9 232.1 238.8 230.2 235.7 230.8 206 240.3 240.5 222.4
14-Jun 124.9 124.5 147.9 145.1 151.6 150.3 161.7 155.2 157.9 155.9 124.4 151.3 149.6 155.4
14-Jul 135.3 164.8 196 176.6 190.6 189.8 208.8 184.2 199.3 201.2 155.8 188.5 187.2 179.9
14-Aug 127.1 153.8 174.1 157.6 163.5 167.7 182.3 182.6 185.6 175.5 155.6 174.8 162.5 174
14-Sep 144.9 140.2 163.2 161 163.5 159.6 159.9 154.9 161.2 164 140.9 166 154.6 152.1
14-Oct 142.6 136 160 158.6 153.7 143.3 151.1 168.5 146.7 151.8 142.9 144.1 148.1 174.5
14-Nov 131.5 129.2 145.1 145 138 134.1 124.4 148.6 134.6 128 133.7 130.4 136.8 142.9
14-Dec 162.6 186.8 194.3 200.6 201.3 198.9 187 193.2 189.7 194.4 190.9 196.2 193.5 188.4
15-Jan 188.3 180.1 198 213 207.2 202.5 193.9 203.8 196.4 194.8 192.8 202.6 200.5 199.4
15-Feb 171.6 133.2 171.9 187.2 182 169.6 167.7 161.3 166.7 163.1 164.4 172.2 177.9 164.3
15-Mar 192.3 187 206.1 202.6 202 202.7 196.9 210.6 198.4 202.8 200.1 200.6 203.5 210.7
15-Apr 148.8 164.9 174 158 155.5 163 166.4 166.2 161.1 169.1 154 161.6 160.7 171.6
15-May 234.8 243.2 256.6 244.4 241.4 244.9 236.9 256.1 242.4 248.1 242 245.7 244.8 259.3
15-Jun 121.6 124.8 153.4 130.8 133.8 145.4 161 158 154.7 159.7 124.3 134.6 142.9 140.1
15-Jul 147.8 161.3 181 170.5 172.2 177.6 191.6 180.5 192.7 178.5 162.4 185.5 178.1 173.5
15-Aug 118.3 126.9 151.7 138.2 143.1 142.7 143.5 152.1 157.6 137.9 127 146.1 148 139.8
15-Sep 94.1 104.8 119.8 118 120.3 114.8 129.3 119.1 122.8 124 101.8 117.4 119.2 113.7
15-Oct 173.9 181 198.4 195.3 193.4 199.3 182 189.5 193.3 199.4 182.9 199 197.3 183.9
15-Nov 109.8 115.8 129.8 126 114.1 120.9 117.8 131.7 121.8 123.7 113.8 118.7 115.3 125.7
15-Dec 170.1 172.7 176.9 193.8 185.4 171.5 158.4 178.9 159 164.6 179.4 174.9 184.5 178.6
16-Jan 175.1 189.6 203.2 222.6 210.5 209.5 205.3 189.5 202.5 203.3 185.9 211.1 209.5 182.8
16-Feb 179.4 163.4 177.1 187.8 181.1 173.3 167.6 177.2 172.8 168 170.4 179.5 179.3 189
16-Mar 170.2 159.1 179.3 182.9 170.2 167 158.2 176.6 156.3 167.3 165.5 172.8 175 190.4
16-Apr 147.2 140.4 160.5 156.5 153.1 142.2 154.3 162.1 151.5 149.2 143.6 156.4 152 182.5
16-May 171.3 191.7 197 196.8 202.6 197.9 203.4 196.6 206.6 193.8 183.8 202.3 201.9 199
16-Jun 138.5 144.3 155.7 147.3 148 143.1 146.4 158.9 141.2 145.1 145.7 140.9 152.5 157
16-Jul 102.3 114.7 167.2 121.1 149.3 151 171 154.2 180.5 171.8 109.5 146.1 135.2 141.9
16-Aug 152.1 156.3 191.3 167.7 186.2 181.5 182.1 201.2 189.6 188.7 150.3 176.8 178.6 197
16-Sep 90.5 87.4 112.9 92.9 98.9 97 93 115.7 81.7 99.2 95.5 100.1 102.9 106.7
16-Oct 151.3 143.7 157.5 148.4 144.7 142.3 141.8 170.9 141 146.6 149.2 140.5 149.3 176.3
16-Nov 141 138.9 154.3 151.8 147.7 143 128.6 155.6 131.7 137.6 139.3 141.1 146.4 162.5
16-Dec 183.1 191.3 205.6 217.1 207.2 206.4 197 204.5 188.2 198.8 196.3 209.9 205.5 195.6
17-Jan 195.6 208.6 223.6 234.2 233.7 231.3 227.6 215.5 234.2 219.1 206 233.9 227.8 208.1
17-Feb 209.5 200.7 208 232.9 221.2 201.7 194.9 203.4 202.2 195.3 206.1 208.5 216.6 201.6
17-Mar 172.1 184.8 198.5 208.6 203 209.7 195.5 191.6 203.8 194.2 181.6 206.9 199.2 190.9
17-Apr 170.1 190.4 210.1 217.6 210.7 199.4 193 198.1 202.7 204.7 181.9 201 215.8 200.7
17-May 169.6 185.2 211.9 197.9 196.2 206.6 220.6 211 214.9 209.2 180.6 205.9 197.3 197.3
17-Jun 139.8 157.3 174.3 179.1 176.6 171.8 164.2 177.1 169.6 171.4 148.3 170.1 177.4 172.4
17-Jul 128.9 154.1 174 152.2 175.2 180.8 201.9 178.3 200.9 180 147.7 184.3 171 175.9
17-Aug 57.2 69.3 100 74.6 84.8 92.8 108.6 94.8 101.3 107.4 60.5 95 87.4 79.4
17-Sep 143.2 138.9 157.8 144.6 145.5 141.4 147.3 156.7 127.6 152.7 144.1 134.7 150.3 158.2
17-Oct 100 91.3 106.3 112.1 111.6 101.3 97.4 105.8 102.4 100.8 100.1 104.9 103.5 102.8
17-Nov 154.6 158.1 166.8 180 168.6 159.3 157.4 166.5 163 162.8 154.8 162.8 171.2 160.6
17-Dec 194.4 197.3 213.4 230.4 221.9 215.4 195.2 204.4 209.1 201.3 200 217.4 219.1 198.3
18-Jan 192.5 207.8 219.1 224.2 224.8 220.3 200.2 214.6 212.4 210.3 206.5 220.1 224 212.6
18-Feb 175.7 177.9 195.6 195 182.4 179.7 165.6 190.7 166.9 180.2 175.3 178.1 186.6 183.2
18-Mar 191.7 176.9 207.9 216 207.1 203 194.1 199 203.3 206.4 193.1 199.3 205.2 209
18-Apr 179.8 193.1 200.2 204.8 202.6 197.8 202.1 199.9 196.1 188.9 192.4 204.6 197.4 194.8
18-May 182.6 185.4 202.2 197.8 198.9 195 192.2 192.4 188.4 190.1 186.5 196.1 198.3 195.9
18-Jun 148.9 163.2 179.8 172 175.9 174.5 178.1 183.6 181.1 174.2 156.6 167.7 167.4 189.3
18-Jul 167.2 192.1 219.4 201.9 220.2 221.8 234 227.8 252.7 222.7 189.6 226.2 214.5 203.6
18-Aug 154.6 179.1 199.2 204.3 202.9 206.2 213 195.2 226.9 210.1 171.3 211.2 206.7 185.3
18-Sep 70.7 71.6 92.4 76.1 82.8 97.1 108.2 86.9 95.9 107.9 70.6 91.5 84.5 72.6
18-Oct 138.6 140.5 156.3 147.7 154.4 150.6 137.8 158.8 145.7 152.4 142.1 148.6 151.4 160.2
18-Nov 156.3 147.2 171.2 165 160 162.6 149.4 169.4 153.3 153.8 150.5 156.4 158.6 175
18-Dec 150.4 143.3 155.8 174.4 164.5 155 147.1 161.8 151.1 153.3 160.4 155.6 161.7 165.5

作成した表から茨城県における7月から9月の平均日射量を算出します。

In [ ]:
amedas_sunshine_all = []

#年によって"range(7,19)"のかっこの中を変更
#集計する期間をインデックスで指定(-Jan,-Feb…,-Dec)
for y in range(7,19):
  df_ame_csv_6 = df_ame_csv_5[str(y)+"-Jul":str(y)+"-Sep"]
  ame_arr = df_ame_csv_6.values.astype(float)
  ame_arr_ave=np.average(ame_arr)
  amedas_sunshine_all.append(ame_arr_ave)

df_ame_all_year = pd.Series(amedas_sunshine_all,index=["2007","2008","2009","2010","2011","2012","2013","2014","2015","2016","2017","2018"])
df_ame_final = pd.DataFrame(df_ame_all_year,columns=["average_sunshine"])

df_ame_final
Out[ ]:
average_sunshine
2007 135.790476
2008 133.419048
2009 127.483333
2010 183.800000
2011 176.850000
2012 197.311905
2013 180.271429
2014 168.588095
2015 143.933333
2016 140.228571
2017 134.792857
2018 165.916667

日射量とEVIの比較

日射量の代表値とEVIの関係をグラフで表してみましょう。

In [ ]:
#ax1をAMeDASの日射量、ax2をEVIとして描画
fig, ax1 = plt.subplots(1,1,figsize=(10,8))
ax2 = ax1.twinx()
ax1.bar(df_ame_final.index,df_ame_final["average_sunshine"],color="orangered",label="AMeDAS")
ax2.plot(df_EVI_final["EVI"],linestyle='solid',color='g',marker='^',label='EVI')
ax1.set_ylim(100,250)
ax1.set_ylabel("AMeDAS")
ax1.set_xlabel("year")
ax2.set_ylim(2500,4000)
ax2.set_ylabel("EVI")
handler1, label1 = ax1.get_legend_handles_labels()
handler2, label2 = ax2.get_legend_handles_labels()
ax1.legend(handler1+handler2,label1+label2,borderaxespad=0)
ax1.grid(True)
fig.show()

このグラフから、日射量とEVIとの間に類似した傾向があることが見て取れます。日照時間が長ければ長いほど、水稲が順調に生育しているようです。

水稲収量の予測

ここまで、MODISデータから2007年から2018年の茨城県の水田のEVIを作成し、日射量による影響を見てみました。次に、2007年から2017年における水稲収量の実測値とEVIを用いて単回帰分析を行います。

また、予測式の検証として2018年のEVIを作成した予測式に当てはめ、同年の収量を予測し、推定値と実測値で比較してみることにします。

収量の実測値の取得

米の収量の実測値は、総務省統計局が運営するポータルサイト「e-Stat」から入手できます。このサイトは各府省が公表する統計データを管理、提供しており、APIでデータを入手することもできます。

米の収量は、kg/10aの単位で提供されてています。1a(アール)は100平方メートル(10m x 10m)です。

APIを利用するために、まずe-Statのユーザ登録のページで登録してください。

今回は、作物統計調査の水稲のデータセットから、茨城県を対象に2007年から2018年の収量情報を取得します。

また、2007年から2010年における県全体の収量がデータベース上で提供されていないことから、各市町村における収量を集計し、全体で平均値化することでその年の県全体の収量としました。

TellusのHow to useのページでも、API経由でe-Statから収量の情報を入手する方法を説明しているので参考にしてください。

In [ ]:
def get_json(base_url,params):
	params_str = urllib.parse.urlencode(params)
	url = base_url + params_str
	json = requests.get(url).json()
	return json
In [ ]:
#APIを登録した際に発行されるURL
appID = "49f9da1e75db39153649b295a88a24d99b3bb56c"
base_url = "http://api.e-stat.go.jp/rest/3.0/app/json/getStatsData?"
params = {"appId":appID,"lang":"J","statsDataId":"0003293480","metaGetFlg":"Y","cntGetFlg":"N","sectionHeaderFlg":"1"}

json = get_json(base_url, params)

必要なエリアのデータを取得します。 各都道府県および市町村にIDが降られているので、e-Statのデータベースより調べリスト化します。

In [ ]:
# 茨城県の市町村IDをリスト化
ID_list = ["08201","08202","08203","08204","08205","08207",
             "08208","08210","08211","08212","08214","08215",
             "08216","08217","08219","08220","08221","08222",
             "08223","08224","08225","08226","08227","08228",
             "08229","08230","08231","08232","08233","08234",
             "08235","08236","08302","08309","08310","08341",
             "08364","08442","08443","08447","08521","08542","08546","08564"]
In [ ]:
#データのダウンロード

df_yield_base = pd.DataFrame()

for ID in ID_list:
  # valuesを探す
  all_values = pd.DataFrame(json['GET_STATS_DATA']['STATISTICAL_DATA']['DATA_INF']['VALUE'])
  # @cat01が110は収量(kg/10a) @areaはリストにあるIDで指定し、茨城県の市町村ごとのデータをループ処理で取得する。
  values = all_values[(all_values["@cat01"] == "110") & (all_values["@area"] == ID)]
  # 取得した時間の末尾に000000がつくので削除
  values = values.replace("0{6}$", "", regex=True)
  values = values.rename(columns = {"@time": "年", "$": "10aあたりの収量","@unit": "単位"})
  values = values.astype({"10aあたりの収量":"int32"})
  df_yield_base = df_yield_base.append(values)

各市町村の収量を集計し、平均値化することでその年における茨城県の収量を算出します。

In [ ]:
year = ["2007","2008","2009","2010","2011","2012","2013","2014","2015","2016","2017","2018"]

df_yield = pd.DataFrame(columns=["年","収量合計[kg/10a]"])

#茨城県内の収量の平均値を算出
for y in year:
  ind = df_yield_base[df_yield_base["年"] == y]
  yield_ave = (ind["10aあたりの収量"].mean())
  df_yield_sub = pd.Series([y,yield_ave], index=df_yield.columns)
  df_yield = df_yield.append(df_yield_sub, ignore_index=True)

#単位を変換した年ごとの収量合計表を表示
df_yield = df_yield.set_index("年")

df_yield
Out[ ]:
収量合計[kg/10a]
2007 501.659091
2008 529.386364
2009 517.500000
2010 517.318182
2011 518.863636
2012 534.681818
2013 537.363636
2014 544.295455
2015 501.954545
2016 516.704545
2017 520.613636
2018 519.090909
In [ ]:
#データフレームをリスト化
yield_all_df = df_yield.values.tolist()

yield_all = []

for ap in range(len(year)):
  yield_all.append(yield_all_df[ap][0])

EVI値との比較

収量の実測値をEVIの値と並べてみましょう。

In [ ]:
#図表用のインデックス
name=["EVI","Yield[kg/10a]"]
index_2=["2007","2008","2009","2010","2011","2012","2013","2014","2015","2016","2017"]

#集計した値のリスト
EVI = EVI_all[0:11]
y_10 = yield_all[0:11]
In [ ]:
s1 = pd.Series(EVI,index=index_2[0:11],name=name[0])
s2 = pd.Series(y_10,index=index_2[0:11],name=name[1])
df_res = pd.DataFrame({name[0]:EVI,name[1]:y_10},index=index_2[0:11]).astype(float)
df_res
Out[ ]:
EVI Yield[kg/10a]
2007 3300.887009 501.659091
2008 3288.469720 529.386364
2009 3273.255674 517.500000
2010 3100.455152 517.318182
2011 3271.328670 518.863636
2012 3648.216596 534.681818
2013 3634.982025 537.363636
2014 3618.293318 544.295455
2015 3128.317849 501.954545
2016 3157.609067 516.704545
2017 3564.870449 520.613636
In [ ]:
#ax1をAMeDASの日射量、ax2をEVIとして描画
fig, ax1 = plt.subplots(1,1,figsize=(10,8))
ax2 = ax1.twinx()
ax1.bar(df_res.index,df_res["Yield[kg/10a]"],color="gold",label="Yield")
ax2.plot(df_res["EVI"],linestyle='solid',color='g',marker='^',label='EVI')
ax1.set_ylim(500,600)
ax1.set_ylabel("Yield")
ax1.set_xlabel("year")
ax2.set_ylim(2500,4000)
ax2.set_ylabel("EVI")
handler1, label1 = ax1.get_legend_handles_labels()
handler2, label2 = ax2.get_legend_handles_labels()
ax1.legend(handler1+handler2,label1+label2,borderaxespad=0)
ax1.grid(True)
fig.show()

このグラフを見ると、収量とEVIの傾向は概ね同じであることがわかります。ただし、2007年や2008年のように傾向が一致しない年もあり、考察が必要です。2008年は、8月中下旬が日照不足で登熟が抑制されたようですが、その後の天候回復で収量が回復したようです(水稲の作柄に関する委員会の資料参照)。

EVIと収量の単回帰分析

これまでに取得したEVIの値と収量の実測値を使って、単回帰分析を行います。

はじめに、EVIと収量との関係性を評価するために相関係数を算出します。

In [ ]:
x = df_res["EVI"]
y = df_res["Yield[kg/10a]"]

mod = LinearRegression()
df_x = pd.DataFrame(x)
df_y = pd.DataFrame(y)

#決定係数(R^2)を評価
mod_lin = mod.fit(df_x, df_y)
y_lin_fit = mod_lin.predict(df_x)
r2_lin = mod.score(df_x, df_y)

print('R = ' + str(round(math.sqrt(r2_lin), 4)))
R = 0.7478

計算した結果、R(相関係数)= 0.7 となりました。

相関係数は一般的に0.4から0.7の間であれば、正の相関が認められることからEVIと収量の間に相関関係が確認することができました。

EVIと収量の間に相関関係が確認できたので、収量を目的変数、EVIを説明変数として単回帰分析で収量予測式を作成します。

In [ ]:
# 予測モデルを作成

clf = linear_model.LinearRegression()

x2 = s1.values.reshape(-1,1)
y2 = s2.values.reshape(-1,1)

clf.fit(x2, y2) 
a = clf.coef_[0][0]
b = clf.intercept_[0]
c = clf.score(x2,y2)
plt.figure(figsize=(8, 8), dpi=80)
plt.scatter(x2, y2 ,c='#00FFFF', edgecolors="#37BAF2")

# 回帰直線の表示

plt.title('Linear regression')
plt.plot(x2, clf.predict(x2))
plt.xlabel('EVI')
plt.ylabel('Yield')
plt.grid()
plt.show()

# 各パラメータの表示 /// 表示の書き方を統一。型ををもっとシンプルな
print("回帰係数= ", clf.coef_[0][0])
print("切片= ", clf.intercept_[0])
print("決定係数= ", clf.score(x2, y2))
print("収量予測式= ", "y= " + str(a) + "x + " + str(b))
回帰係数=  0.04772084617820682
切片=  361.3913616586648
決定係数=  0.5592706325140075
収量予測式=  y= 0.04772084617820682x + 361.3913616586648

収量の予測と検証

作成した収量予測式を用いて2018年における収量を推定します。

In [ ]:
#2018年の結果を予測

#集計した値のリスト
list_2018 = [EVI_all[11] ,yield_all[11]]

yield_est = (a * list_2018[0]) + b

print("2018年の収量予測値[kg/10a]:",yield_est)
print("2018年の収量実測値[kg/10a]:",list_2018[1])

#差分計算

yield_diff = yield_est - list_2018[1]

print("誤差[kg/10a]:",yield_diff)
2018年の収量予測値[kg/10a]: 518.5837043920321
2018年の収量実測値[kg/10a]: 519.0909090909091
誤差[kg/10a]: -0.5072046988769898

2007年から2017年のEVIの値を用いて作成した収量予測式より2018年の収量を予測した結果、実測値に比べて約0.5 kg/10aの誤差で茨城県全体の収量を予測することができました。シンプルな予測式としてはまずまずの結果だと思います。

現在、農林水産省では複数の気象データを変数とした多変量解析による予測式の作成が取り組まれています。また、作成した予測式より上位5式を選択することで、各地域に合った予測式が作成されています。

本章ではEVIのみを使った簡単な単回帰分析で予測式を作りましたが、実際には他の変数や作物ごとの収量予測モデルが使われています。今回の方法を発展させ、複数の変数を組み合わせた予測にも挑戦してみてください。

まとめ

この章では、衛星データを用いて水稲収量を予測する方法について学習しました。

過去の衛星画像を時系列的に処理し推定モデルを作成するような場合、解析に用いるデータセットの品質が肝になります。解析を行う際には、画像や統計データの品質に注意し、予測が合わない場合の事由について考察することが重要です。

また、本章は解像度の低いMODISのデータを使った県単位の解析を行いました。狭い範囲を対象とする場合は別の衛星を使う必要があります。