Tellusでは、利便りべん性せい（データのCOG化か）と検索けんさく精度せいど（APIのフォーマット統一とういつ）の向上こうじょうを目指めざして、データ配信はいしん用ようAPIの更新こうしんを行おこないました。それに伴ともない、今回こんかいは新あたらしい配信はいしん用ようAPIを利用りようしてAVNIR-2のデータを取得しゅとくする方法ほうほうについてご紹介しょうかいします。
https://www.tellusxdp.com/ja/news/20220630_000502.html

Tellusの解析かいせき環境かんきょうを契約けいやくしていない方ほうであっても利用りよう可能かのうなAPIを用もちいて解析かいせきを行おこないますので、皆みなさまお手元てもとの環境かんきょうでお試ためしください。

なお、本ほん記事きじは技術評論社ぎじゅつひょうろんしゃから出版しゅっぱんされている月刊げっかん誌しSoftware Designに宙ちゅう畑はたけが寄稿きこうした連載れんさい「衛星えいせいデータプラットフォームTellusハンズオン」の内容ないようから、一部いちぶ修正しゅうせいして掲載けいさいするものです。

衛星えいせいデータで分わかること・できること

宙ちゅう畑はたけでは過去かこにも何なん度どか紹介しょうかいしていますが、改あらためて「波長はちょう」とは文字通もじどおり波なみの長ながさのことです。色いろや通信つうしんに使つかう電波でんぱや赤外線せきがいせん、紫外線しがいせんなどすべて電磁波でんじはであり、それぞれ「波長はちょう」によって異ことなる性質せいしつを持もちます。

詳くわしくはこちらをご覧らんください：

人ひとの目めはこの電磁波でんじはの中なかで可視かし光こうといわれる限かぎられた範囲はんいの波長はちょう帯たいしか認識にんしき（観測かんそく）できません。この可視かし光こうの波長はちょう帯たいを青あお、緑みどり、赤あかの色いろの組くみ合あわせで捉とらえています。

これに対たいして、人工じんこう衛星えいせいのセンサーでは可視かし光こう以外いがいの波長はちょう帯たいも観測かんそくできます。観測かんそくしている波長はちょうの中なかには、紫外線しがいせんや赤外線せきがいせん、電波でんぱをとらえるセンサーを搭載とうさいしているものもあるので、私わたしたちが目めにする地球ちきゅうとは異ことなった視点してんで地球ちきゅうの姿すがたを見みることができます。

波長はちょうによってなぜ見みえ方かたが違ちがうのでしょうか。それは対象たいしょうによって吸収きゅうしゅう（透過とうか）あるいは反射はんしゃする太陽光たいようこうの波長はちょうが異ことなるからです。

例たとえば、植物しょくぶつは他たの物質ぶっしつに比くらべて近赤外線きんせきがいせん※の帯域たいいきで強つよく反射はんしゃすることが分わかっています。つまり、近赤外線きんせきがいせんで観測かんそくした結果けっかを利用りようすると、植物しょくぶつの分布ぶんぷを調しらべやすいと言いえます。
※近赤外線きんせきがいせん（きんせきがいせん）：可視かし光こうに近ちかい波長はちょうの赤外線せきがいせん

さらに、それぞれ個別こべつの波長はちょうで調しらべるだけではなく、複数ふくすうの波長はちょうのデータを組くみ合あわせて計算けいさんすることで調しらべたいものが判別はんべつしやすくなる場合ばあいがあります。

今回こんかいは、よく利用りようされる波長はちょうの組くみ合あわせである正規せいき化か植生しょくせい指数しすう（NDVI）という指数しすうについて解析かいせき事例じれいを交まじえて紹介しょうかいします。

正規せいき化か植生しょくせい指数しすう＝NDVI(Normalized Difference Vegetation Index)とは、植物しょくぶつの活性かっせい度どを調しらべる指数しすうのことです。近赤外線きんせきがいせんの光ひかりの強つよさ（NIR）と赤色あかいろの光ひかりの強つよさ（Red）を以下いかの式しきに当あてはめることで求もとめることができます。

NDVI＝（NIR－Red）/（NIR＋Red）

例たとえば、以下いかのように過去かこ何なん年ねんか分ぶんのNDVIの値ねを時とき系列けいれつに見みていくことで、農作物のうさくもつの生育せいいく状じょう況きょうを調しらべることができます。

衛星えいせいデータプラットフォームTellus登録とうろくと準備じゅんび

TellusのAPIを利用りようするためには、Tellusのアカウント登録とうろくが必要ひつようです。まだアカウントをお持もちではない方ほうは、以下いかのページで必要ひつよう情報じょうほうを入力にゅうりょくし、アカウントを作成さくせいします。
https://www.tellusxdp.com/market/login/

Tellusで公開こうかいしているデータの利用りよう規則きそくはそれぞれのデータの提供ていきょう元もとのデータポリシーに基もとづきます。中なかには、Tellusが提供ていきょうする開発かいはつ環境かんきょう内ないでの利用りように限定げんていされ、ローカル環境かんきょうへのダウンロードが禁止きんしされているデータもあるので注意ちゅういが必要ひつようです。

今回こんかいは、Tellus環境かんきょう外がいからも利用りよう可能かのうなデータを使用しようしますので、Tellusアカウント登録とうろくをした上うえでAPIキーの発行はっこうさえすれば、お手元てもとのJupyter NotebookやGoogle Colabなどでもお試ためしいただけます。

本ほん記事きじではGoogle Colabを例れいに紹介しょうかいします。

まずは自分じぶんのGoogle Driveと接続せつぞくします。

## google driveと接続せつぞく
from google.colab import drive
drive.mount('/content/drive')

解析かいせきに必要ひつようなライブラリをインポートします。

!pip install rasterio
!pip install geopandas

import os
import requests
import json
import glob
import rasterio.mask
import fiona

import numpy as np
import rasterio as rio
import matplotlib.pyplot as plt
import geopandas as gpd
from rasterio.plot import show
from rasterio.plot import show_hist

API利用りようの際さいに必要ひつようなトークンは、Tellusにログイン後ご、右みぎ上じょうのユーザー名めいのプルダウンメニューから[開発かいはつ環境かんきょう]＞[APIアクセス]を選択せんたくし、[トークンの発行はっこう]から新あたらしいトークンを発行はっこうします。

#各かく個人こじんのAPIトークンを入力にゅうりょく
token = 'XXXXXXXXXXXXXXXX'

衛星えいせいデータの検索けんさく・取得しゅとく

Tellusの衛星えいせいデータAPIのリファレンスは以下いかで公開こうかいしています。
https://www.tellusxdp.com/docs/travelers/

Tellusで公開こうかいしている衛星えいせいデータは「データセット」「シーン」「ファイル」という構造こうぞうになっています。

▼Tellusのデータ構造こうぞう
データセットA
　　└ シーンa
　　　　└ ファイルa-i
　　　　└ ファイルa-ii
　　└ シーンb
　　　　└ ファイルb-i
　　　　└ ファイルb-i
・・・

「データセット」はある衛星えいせいに搭載とうさいされたセンサや、前ぜん処理しょり方式ほうしきなどが揃そろったデータのまとまりです。

「シーン」は特定とくていの場所ばしょをあるタイミングで撮影さつえいしたまとまりで、撮影さつえいの方法ほうほうなどにより、例たとえば赤あか青あお緑みどりの波長はちょうの要素ようそ毎ごとに複数ふくすうの「ファイル」が存在そんざいします。

詳くわしくはこちらの記事きじ等とうをご覧らんください。

衛星えいせいデータをAPIで取得しゅとくする前まえに、Tellus Traveler（ブラウザ）上じょうで必要ひつような衛星えいせいデータに関かんする情報じょうほうを取得しゅとくします。Travelerの使つかい方かたは以下いかをご覧らんください。
https://www.tellusxdp.com/ja/traveler-service/guide/traveler-how-to/2_traveler.html

今回こんかいは、2006年ねん10月がつ31日にちにAVNIR-2によって撮影さつえいされた青森あおもり県けんの画像がぞうを利用りようします。ご自身じしんで興味きょうみのある領域りょういきを検索けんさくし、データセットIDとシーンIDを取得しゅとくして解析かいせきにチャレンジしてみてください。

上記じょうきの衛星えいせい画像がぞうのデータセットIDとシーンIDは以下いかの通とおりです。

dataset_id='ea71ef6e-9569-49fc-be16-ba98d876fb73'

scenes = [
    '202ce08d-ba4b-4ffe-8165-109fd3a8b917'
]

調しらべたデータセットIDとシーンIDから、そのシーンに含ふくまれるファイル名めいを取得しゅとくします。Tellus Travelerを利用りようしてデータを取得しゅとくする関数かんすうは以下いかの通とおりです。

# データセットIDとシーンIDから、そのシーンに含ふくまれるファイル名めいを取得しゅとくする
def get_dataset_data_by_id_files(base_path, token, dataset_id, data_id):
    print(base_path, token, dataset_id, data_id)
    url = '{}/api/traveler/v1/datasets/{}/data/{}/files/'.format(
        base_path, dataset_id, data_id)
    headers = {
        'Authorization': 'Bearer ' + token,
        'content-type': 'application/json'
    }
    r = requests.get(url, headers=headers)
    assert r.status_code == 200
    return json.loads(r.content)

次つぎに、ファイルIDを指定していしてデータのダウンロードURLを取得しゅとくします。

# ファイルIDを指定していして、ダウンロードURLを取得しゅとくする
def get_dataset_data_by_id_files_by_id_download_url(base_path, token, dataset_id, data_id, file_id):
    url = '{}/api/traveler/v1/datasets/{}/data/{}/files/{}/download-url/'.format(
        base_path, dataset_id, data_id, file_id)
    headers = {
        'Authorization': 'Bearer ' + token,
        'content-type': 'application/json'
    }
    r = requests.post(url, headers=headers)
    assert r.status_code == 200
    return json.loads(r.content)

ダウンロードURLを用もちいて衛星えいせいデータをダウンロードする関数かんすうは以下いかの通とおりです。

# 衛星えいせいデータのダウンロードを行おこなう
def download(scenes, dataset_id, token, dist='./', base_path='https://www.tellusxdp.com'):
    for scene_id in scenes:       
        files = get_dataset_data_by_id_files(
            base_path, token, dataset_id, scene_id)
        print(files)
        rawdata = files['results']
        path = os.path.join(dist, scene_id)
        if len(rawdata) > 0:
            for file in rawdata:
                file_id = file['id']
                file_name = file['name']

                file_path = os.path.join(path, file_name)

                download_url = get_dataset_data_by_id_files_by_id_download_url(
                        base_path, token, dataset_id, scene_id, file_id)['download_url']

                r = requests.get(download_url, stream=True)

                if not os.path.exists(path):
                    os.makedirs(path)

                with open(file_path, 'wb') as f:
                    for chunk in r.iter_content(chunk_size=1024):
                        if chunk:

                            f.write(chunk)
                            f.flush()

これらの関数かんすうを実行じっこうすることで、衛星えいせいデータをご自身じしんのGoogle Driveにダウンロードできます。

# データDL
download(scenes, dataset_id, token)

NDVIの算出さんしゅつ

さて、様々さまざまなバンドのデータがダウンロードできたらそのデータを処理しょりしてNDVIを算出さんしゅつしていきましょう。

本来ほんらいNDVIは輝度きど値ちではなく、輝度きど値ちに大気たいき補正ほせい処理しょりを行おこない算出さんしゅつする地上ちじょう反射はんしゃ率りつを利用りようしますが、今回こんかいは処理しょりが簡素かんそな輝度きど値ちに変換へんかんして抽出ちゅうしゅつしています。

# NDVIの算出さんしゅつ
def mk_ndvi_rio(red_img, nir_img, red_gain, red_offset, nir_gain, nir_offset, out_ndvi):
 red = rio.open(red_img)
  nir = rio.open(nir_img)

  red_arr = red.read().astype(float)
  nir_arr = nir.read().astype(float)

  kwargs = {
      "width" : red.width,
      "height" : red.height,
      "count" : 1,
      "crs" : red.crs,
      "transform" : red.transform,
      "dtype" : rio.float32}
      
  ## DN2輝度きど変換へんかん
  red_arr = red_arr * red_gain + red_offset
  nir_arr = nir_arr * nir_gain + nir_offset
  
  ## NDVI算出さんしゅつ
  ndvi = np.where(((red_arr != 0) & (nir_arr != 0)),((nir_arr - red_arr) / (nir_arr + red_arr) ),0)

  write_ndvi = rio.open(out_ndvi,"w",driver="Gtiff",**kwargs)
  write_ndvi.write(ndvi)
  write_ndvi.close()

  return(out_ndvi)

NDVI値ちの画像がぞうを作成さくせいしていきます。

DATA_DIR = "/content/drive/MyDrive/work/study/Sorabatake/02_data"
DIR_LIST = os.listdir(DATA_DIR)
NDVI_DIR = "/content/drive/MyDrive/work/study/Sorabatake/03_ndvi"
#print(DIR_LIST) 

#データを取得しゅとく
for dir_name in DIR_LIST:

  dir_path = os.path.join(DATA_DIR, dir_name)
  #print(dir_path)

  subdir_list = os.listdir(dir_path)
  red_path = glob.glob(dir_path + "/IMG-03-*.tif")
  nir_path = glob.glob(dir_path + "/IMG-04-*.tif")
  hdr_path = glob.glob(dir_path + "/HDR-*.txt")
  hdr_path = str(hdr_path[0])
 
  # HDR情報じょうほうから、ゲイン、オフセット情報じょうほうを抽出ちゅうしゅつ
  hdr_f = open(hdr_path,"r")
  hdr_data = hdr_f.read()
  hdr_f.close

  # gain, offset
  red_gain = float(hdr_data[1753:1761])
  red_offset = float(hdr_data[1761:1769])
  nir_gain = float(hdr_data[1769:1777])
  nir_offset = float(hdr_data[1777:1785])
  print(red_gain, red_offset, nir_gain, nir_offset)

  red_path = red_path[0]
  nir_path = nir_path[0]
  print(red_path, nir_path)
  
  #NDVI画像がぞうの作成さくせい
  base_name = os.path.basename(red_path)
  ndvi_file = "NDVI-" + base_name[7:]
  ndvi_path = os.path.join(NDVI_DIR, ndvi_file)

作成さくせいしたNDVI画像がぞうを表示ひょうじします。

NDVI_20061031 = rio.open(NDVI_DIR + "/NDVI-ALAV2A040842780-OORIRFU-D069P3-20061031-001.tif")

NDVI_20061031 = NDVI_20061031.read(1)

fig = plt.figure(figsize=(18,18))


img1 = ax1.imshow(NDVI_20061031, cmap='RdYlGn', vmin=-1, vmax=1)
ax1.set_title("NDVI_20061031")
cbar1 = fig.colorbar(img1,ax=ax1, aspect=40,pad=0.08,shrink=0.8)

plt.show(fig)

以下いかのような画像がぞうが表示ひょうじされます。

植物しょくぶつが活発かっぱつなところを濃こい緑みどりで、植物しょくぶつのないところ（湖みずうみなど）が赤あかで表示ひょうじされています。一般いっぱん的てきなカラー画像がぞう（右みぎ）と比較ひかくしても、植物しょくぶつが活性かっせい化かしている部分ぶぶんがより分わかりやすくなっています。

NDVIを時とき系列けいれつにグラフ化かする

ここまでは、1枚まいの衛星えいせいデータについて分析ぶんせきをしてきましたが、同おなじ場所ばしょの衛星えいせいデータを時とき系列けいれつで見みていくことで、田畑たはたの利用りよう状じょう況きょうを推定すいていしてみます。例たとえば、夏場なつばにNDVIが低ひくいところがあれば、その田畑たはたは作付さくづけが行おこなわれていない、つまりは遊休ゆうきゅう農地のうちである可能かのう性せいがあると考かんがえることができます。

そこで、農林水産省のうりんすいさんしょうが公開こうかいしている筆ふでポリゴンと呼よばれる田畑たはたの区画くかくデータを使つかって、区画くかくごとのNDVI値ちの時とき系列けいれつ推移すいいを求もとめてみましょう。

「衛星えいせいデータの検索けんさく・取得しゅとく」の節ふしのコードを用もちいてデータをダウンロードします。ダウンロードしたデータからNDVIを算出さんしゅつしていきます。「衛星えいせいデータの検索けんさく・取得しゅとく」の続つづきとして以下いかを実行じっこうしてください。

def mk_ndvi_rio(red_img, nir_img, red_gain, red_offset, nir_gain, nir_offset, out_ndvi):
  """
    NAME:
        make_ndvi_rio
    
    ARGUMENT:
        red_img   : 赤あかバンド画が象ぞう
        nir_img   : 赤あか外がいパンド画像がぞう
        red_gain  : 赤あかバンドゲイン
        red_offset: 赤あかバンドオフセット
        nir_gain  : 赤あか外がいバンドゲイン
        nir_offset: 赤あか外がいバンドオフセット
        out_ndvi  : NDVI画像がぞう

    DISCRIPTION:
　      The original code was written by N.Furuta 

        ■DN値ちから輝度きど値ちへの変換へんかん式しき
          L = a×O+b
            L：輝度きど（W/m2/sr/μみゅーm）
            O：校正こうせい済ずみデジタル値ち(DN値ち)
            a：絶対ぜったい校正こうせい係数けいすう(ゲイン値ち)
            b：オフセット値ち
            ※プロダクトフォーマット説明せつめい書しょ AVNIR-2編へん、表ひょう3.3-8 アンシラリ2（ラジオメトリック校正こうせい）レコード、p3-30 参照さんしょう　
        
        ■NDVI値ちの計算けいさん
          (近きん赤あか外がいバンド - 赤あかバンド)/(近きん赤あか外がいバンド + 赤あかバンド)
          ※本来ほんらいは輝度きど値ちをさらに地上ちじょう反射はんしゃ率りつへ変換へんかんしたものを利用りよう
          ※本関ほんせき数すうでは簡単かんたんのため輝度きど値ちを用もちいて計算けいさん（簡易かんいNDVI）
        
    RETURN:
        out_ndvi : NDVI画像がぞう（パス）
  """
  
  red = rio.open(red_img)
  nir = rio.open(nir_img)

  red_arr = red.read().astype(float)
  nir_arr = nir.read().astype(float)

  kwargs = {
      "width" : red.width,
      "height" : red.height,
      "count" : 1,
      "crs" : red.crs,
      "transform" : red.transform,
      "dtype" : rio.float32}
      
  ## DN2輝度きど変換へんかん
  red_arr = red_arr * red_gain + red_offset
  nir_arr = nir_arr * nir_gain + nir_offset
  
  ## NDVI算出さんしゅつ
  ndvi = np.where(((red_arr != 0) & (nir_arr != 0)),((nir_arr - red_arr) / (nir_arr + red_arr) ),0)

  write_ndvi = rio.open(out_ndvi,"w",driver="Gtiff",**kwargs)
  write_ndvi.write(ndvi)
  write_ndvi.close()

  return(out_ndvi)

NDVI画像がぞうを作成さくせいします。

"""
２．NDVI画像がぞうの作成さくせい
  ## HDR情報じょうほうは公式こうしきドキュメントを参照さんしょう
  https://www.eorc.jaxa.jp/ALOS/alos-ori/doc/AVNIR-2_ORI_format_jp.pdf

"""
# 衛星えいせい画像がぞう（AVNIR-2）を保存ほぞんしているフォルダを指定してい
# 自身じしんの環境かんきょうに合あわせて指定してい先さきは変更へんこうしてください
DATA_DIR = "/content/drive/MyDrive/work/study/Sorabatake/02_data"
DIR_LIST = os.listdir(DATA_DIR)

# NDVI画像がぞうの保存ほぞん先さきのフォルダを指定してい
NDVI_DIR = "/content/drive/MyDrive/work/study/Sorabatake/03_ndvi"
#print(DIR_LIST) 

#データを取得しゅとく
for dir_name in DIR_LIST:

  dir_path = os.path.join(DATA_DIR, dir_name)
  #print(dir_path)

  subdir_list = os.listdir(dir_path)
  red_path = glob.glob(dir_path + "/IMG-03-*.tif")
  nir_path = glob.glob(dir_path + "/IMG-04-*.tif")
  hdr_path = glob.glob(dir_path + "/HDR-*.txt")
  hdr_path = str(hdr_path[0])
  #print(hdr_path)

  # HDR情報じょうほうから、ゲイン、オフセット情報じょうほうを抽出ちゅうしゅつ
  hdr_f = open(hdr_path,"r")
  hdr_data = hdr_f.read()
  hdr_f.close

  # gain, offset
  red_gain = float(hdr_data[1753:1761])
  red_offset = float(hdr_data[1761:1769])
  nir_gain = float(hdr_data[1769:1777])
  nir_offset = float(hdr_data[1777:1785])
  print(red_gain, red_offset, nir_gain, nir_offset)

  red_path = red_path[0]
  nir_path = nir_path[0]
  print(red_path, nir_path)
  
  #NDVI画像がぞうの作成さくせい
  base_name = os.path.basename(red_path)
  ndvi_file = "NDVI-" + base_name[7:]
  ndvi_path = os.path.join(NDVI_DIR, ndvi_file)

まず、目視もくしで変化へんかを観察かんさつしてみましょう。

以下いかのコードを実行じっこうすることで、複数ふくすう時期じきの青森あおもり県域けんいきのNDVI画像がぞうを表示ひょうじします。青森あおもり以外いがいのデータを取得しゅとくしている方ほうは以下いかのファイル名めいを修正しゅうせいしてください。

#NDVI画像がぞうの表示ひょうじ
NDVI_20061031 = rio.open(NDVI_DIR + "/NDVI-ALAV2A040842780-OORIRFU-D069P3-20061031-001.tif")
NDVI_20070618 = rio.open(NDVI_DIR + "/NDVI-ALAV2A074392780-OORIRFU-D069P3-20070618-002.tif")
NDVI_20080805 = rio.open(NDVI_DIR + "/NDVI-ALAV2A134782780-OORIRFU-D069P3-20080805-003.tif")
NDVI_20081105 = rio.open(NDVI_DIR + "/NDVI-ALAV2A148202780-OORIRFU-D069P3-20081105-000.tif")
NDVI_20090508 = rio.open(NDVI_DIR + "/NDVI-ALAV2A175042780-OORIRFU-D069P3-20090508-003.tif")
NDVI_20100626 = rio.open(NDVI_DIR + "/NDVI-ALAV2A235432780-OORIRFU-D069P3-20100626-002.tif")
NDVI_20100926 = rio.open(NDVI_DIR + "/NDVI-ALAV2A248852780-OORIRFU-D069P3-20100926-001.tif")

NDVI_20061031 = NDVI_20061031.read(1)
NDVI_20070618 = NDVI_20070618.read(1)
NDVI_20080805 = NDVI_20080805.read(1)
NDVI_20081105 = NDVI_20081105.read(1)
NDVI_20090508 = NDVI_20090508.read(1)
NDVI_20100626 = NDVI_20100626.read(1)
NDVI_20100926 = NDVI_20100926.read(1)

fig = plt.figure(figsize=(18,18))
ax1 = fig.add_subplot(331)
ax2 = fig.add_subplot(332)
ax3 = fig.add_subplot(333)
ax4 = fig.add_subplot(334)
ax5 = fig.add_subplot(335)
ax6 = fig.add_subplot(336)
ax7 = fig.add_subplot(337)

img1 = ax1.imshow(NDVI_20061031, cmap='RdYlGn', vmin=-1, vmax=1)
ax1.set_title("NDVI_20061031")
cbar1 = fig.colorbar(img1,ax=ax1, aspect=40,pad=0.08,shrink=0.8)
img2 = ax2.imshow(NDVI_20070618, cmap='RdYlGn', vmin=-1, vmax=1)
ax2.set_title("NDVI_20070618")
cbar2 = fig.colorbar(img2,ax=ax2, aspect=40,pad=0.08,shrink=0.8)
img3 = ax3.imshow(NDVI_20080805, cmap='RdYlGn', vmin=-1, vmax=1)
ax3.set_title("NDVI_20080805")
cbar3 = fig.colorbar(img3,ax=ax3, aspect=40,pad=0.08,shrink=0.8)
img4 = ax4.imshow(NDVI_20081105, cmap='RdYlGn', vmin=-1, vmax=1)
ax4.set_title("NDVI_20081105")
cbar4 = fig.colorbar(img4,ax=ax4, aspect=40,pad=0.08,shrink=0.8)
img5 = ax5.imshow(NDVI_20090508, cmap='RdYlGn', vmin=-1, vmax=1)
ax5.set_title("NDVI_20090508")
cbar5 = fig.colorbar(img5,ax=ax5, aspect=40,pad=0.08,shrink=0.8)
img6 = ax6.imshow(NDVI_20100626, cmap='RdYlGn', vmin=-1, vmax=1)
ax6.set_title("NDVI_20100626")
cbar6 = fig.colorbar(img6,ax=ax6, aspect=40,pad=0.08,shrink=0.8)
img7 = ax7.imshow(NDVI_20100926, cmap='RdYlGn', vmin=-1, vmax=1)
ax7.set_title("NDVI_20100926")
cbar7 = fig.colorbar(img7,ax=ax7, aspect=40,pad=0.08,shrink=0.8)

plt.show(fig)

そうすると以下いかのような結果けっかを得えられるかと思おもいます。

今回こんかいは例れいとして、青森あおもり県けん藤崎ふじさき町まちを対象たいしょうにNDVI解析かいせきを行おこなっていきます。そのために上記じょうきで算出さんしゅつしたNDVI画像がぞうを藤崎ふじさき町まちの区画くかくに沿そって切きり出だします。
※藤崎ふじさき町まちの区画くかくポリゴンは国土こくど地理ちり院いんが公開こうかいしている国土こくど数値すうち情報じょうほうのうち、行政ぎょうせい区域くいきデータを利用りようしています。

# 青森あおもり県けん全域ぜんいきの行政ぎょうせい区域くいきデータの中なかから藤崎ふじさき町まちだけのDVI値ちを切きり抜ぬく
# 藤崎ふじさき町まちのポリゴンを抽出ちゅうしゅつ

# 行政ぎょうせい区域くいきデータの保存ほぞん先さきのフォルダを指定してい
# ※自身じしんの環境かんきょうに合あわせて指定してい先さきは変更へんこうしてください
AOMORI_PATH = "/content/drive/MyDrive/work/study/Sorabatake/04_polygon/Aomori_N03-20210101_02_GML/"
# そのうちshape形式けいしきのデータを指定してい
AOMORI_SHP = "N03-21_02_210101.shp"
df_aomori = gpd.read_file(os.path.join(AOMORI_PATH, AOMORI_SHP), encoding="shift-jis")

Fujisaki_machi = df_aomori[df_aomori["N03_004"].isin(["藤崎ふじさき町まち"])]
Fujisaki_machi = Fujisaki_machi.drop(columns=["N03_002","N03_003"])

#投影とうえい変換へんかん
reprojest_Fujisaki_machi = Fujisaki_machi.to_crs({"init":"epsg:32654"})
reproject_out = os.path.join(AOMORI_PATH, "epsg_32654-"+ AOMORI_SHP)
#出力しゅつりょく
reprojest_Fujisaki_machi.to_file(driver="ESRI Shapefile",filename=reproject_out)

#GeoDataFrameを描画びょうが
f,ax = plt.subplots(1, figsize=(6,6))
ax = Fujisaki_machi.plot(axes=ax)
plt.show

同おなじ場所ばしょで解析かいせきを実行じっこうした場合ばあい、上記じょうきを実行じっこうすると以下いかのような結果けっかを得えられます。

次つぎに、NDVI画像がぞうを藤崎ふじさき町まちのポリゴンでマスク処理しょりし、藤崎ふじさき市し以外いがいの範囲はんいのNDVI画像がぞうを削除さくじょします。

# マスク処理しょり済ずみ画像がぞうの保存ほぞん先さきのフォルダを指定してい
# ※自身じしんの環境かんきょうに合あわせて指定してい先さきは変更へんこうしてください
MASK_DIR = "/content/drive/MyDrive/work/study/Sorabatake/03_ndvi/02_fujisaki_masked"
ndvi_list = os.listdir(NDVI_DIR)
ndvi_list = [f for f in ndvi_list if os.path.isfile(os.path.join(NDVI_DIR, f))]

for ndvi_file_name in ndvi_list:

  print("処理しょり開始かいし：" + ndvi_file_name)

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

  ndvi_data = os.path.join(NDVI_DIR, ndvi_file_name)
  print(ndvi_data)
  print(masks)
  with rasterio.open(ndvi_data) as src:
    out_image, out_transform = rasterio.mask.mask(src, masks, nodata=np.nan, crop=True)
    out_meta = src.meta

  out_meta.update({"driver": "GTiff",
                 "height": out_image.shape[1],
                 "width": out_image.shape[2],
                 "transform": out_transform})
  
  with rasterio.open(MASK_DIR + "/Masked-" + str(ndvi_file_name), "w", **out_meta) as dst:
    dst.write(out_image)

次つぎに、目視もくしで確認かくにんするために藤崎ふじさき町まちのNDVI画像がぞうを一覧いちらんで表示ひょうじします。

#藤崎ふじさき町まちで切きり出だしたNDVI画像がぞうの表示ひょうじ
masked_NDVI_20061031 = rio.open(MASK_DIR + "/Masked-NDVI-ALAV2A040842780-OORIRFU-D069P3-20061031-001.tif")
masked_NDVI_20070618 = rio.open(MASK_DIR + "/Masked-NDVI-ALAV2A074392780-OORIRFU-D069P3-20070618-002.tif")
masked_NDVI_20080805 = rio.open(MASK_DIR + "/Masked-NDVI-ALAV2A134782780-OORIRFU-D069P3-20080805-003.tif")
masked_NDVI_20081105 = rio.open(MASK_DIR + "/Masked-NDVI-ALAV2A148202780-OORIRFU-D069P3-20081105-000.tif")
masked_NDVI_20090508 = rio.open(MASK_DIR + "/Masked-NDVI-ALAV2A175042780-OORIRFU-D069P3-20090508-003.tif")
masked_NDVI_20100626 = rio.open(MASK_DIR + "/Masked-NDVI-ALAV2A235432780-OORIRFU-D069P3-20100626-002.tif")
masked_NDVI_20100926 = rio.open(MASK_DIR + "/Masked-NDVI-ALAV2A248852780-OORIRFU-D069P3-20100926-001.tif")

masked_NDVI_20061031 = masked_NDVI_20061031.read(1)
masked_NDVI_20070618 = masked_NDVI_20070618.read(1)
masked_NDVI_20080805 = masked_NDVI_20080805.read(1)
masked_NDVI_20081105 = masked_NDVI_20081105.read(1)
masked_NDVI_20090508 = masked_NDVI_20090508.read(1)
masked_NDVI_20100626 = masked_NDVI_20100626.read(1)
masked_NDVI_20100926 = masked_NDVI_20100926.read(1)

fig = plt.figure(figsize=(16,12))
ax1 = fig.add_subplot(331)
ax2 = fig.add_subplot(332)
ax3 = fig.add_subplot(333)
ax4 = fig.add_subplot(334)
ax5 = fig.add_subplot(335)
ax6 = fig.add_subplot(336)
ax7 = fig.add_subplot(337)

img1 = ax1.imshow(masked_NDVI_20061031, cmap='RdYlGn', vmin=-1, vmax=1)
ax1.set_title("NDVI_20061031")
cbar1 = fig.colorbar(img1,ax=ax1, aspect=40,pad=0.08,shrink=0.8)
img2 = ax2.imshow(masked_NDVI_20070618, cmap='RdYlGn', vmin=-1, vmax=1)
ax2.set_title("NDVI_20070618")
cbar2 = fig.colorbar(img2,ax=ax2, aspect=40,pad=0.08,shrink=0.8)
img3 = ax3.imshow(masked_NDVI_20080805, cmap='RdYlGn', vmin=-1, vmax=1)
ax3.set_title("NDVI_20080805")
cbar3 = fig.colorbar(img3,ax=ax3, aspect=40,pad=0.08,shrink=0.8)
img4 = ax4.imshow(masked_NDVI_20081105, cmap='RdYlGn', vmin=-1, vmax=1)
ax4.set_title("NDVI_20081105")
cbar4 = fig.colorbar(img4,ax=ax4, aspect=40,pad=0.08,shrink=0.8)
img5 = ax5.imshow(masked_NDVI_20090508, cmap='RdYlGn', vmin=-1, vmax=1)
ax5.set_title("NDVI_20090508")
cbar5 = fig.colorbar(img5,ax=ax5, aspect=40,pad=0.08,shrink=0.8)
img6 = ax6.imshow(masked_NDVI_20100626, cmap='RdYlGn', vmin=-1, vmax=1)
ax6.set_title("NDVI_20100626")
cbar6 = fig.colorbar(img6,ax=ax6, aspect=40,pad=0.08,shrink=0.8)
img7 = ax7.imshow(masked_NDVI_20100926, cmap='RdYlGn', vmin=-1, vmax=1)
ax7.set_title("NDVI_20100926")
cbar7 = fig.colorbar(img7,ax=ax7, aspect=40,pad=0.08,shrink=0.8)

plt.show(fig)

そうすると以下いかのような結果けっかを得えられます。

耕作こうさくが行おこなわれていない可能かのう性せいのある圃場ほじょうを抽出ちゅうしゅつ

ここまでで準備じゅんびが整ととのいました。いよいよ遊休ゆうきゅう農地のうち、つまりは耕作こうさくが行おこなわれていない可能かのう性せいのある圃場ほじょうを抽出ちゅうしゅつしていきましょう。

これらの判定はんてい条件じょうけんとして、ここでは稲いねの作付さくづけがされていない圃場ほじょう、つまり周囲しゅういと比較ひかくしてNDVI値ちの季き節ぶし変動へんどうが少すくない圃場ほじょうや夏場なつばのNDVI値ちが小ちいさい圃場ほじょうを判定はんていのための簡易かんい的てきな条件じょうけんとします。

しかし、闇雲やみくもにシーン全体ぜんたいから上記じょうきの条件じょうけんを当あてはめてしまっても、圃場ほじょう以外いがいのエリアも大量たいりょうに抽出ちゅうしゅつされてしまうことが考かんがえられます。そのため、まず初はじめにどこが圃場ほじょうなのかを判定はんていする必要ひつようがあります。その上うえで、全すべての時期じきのNDVIデータにおいてNDVIの変化へんかが少すくない圃場ほじょうを抽出ちゅうしゅつすることとします。

圃場ほじょうデータとして、ここでは農林水産省のうりんすいさんしょうが公開こうかいしている筆ふでポリゴンデータを利用りようします※。これは圃場ほじょうの一ひとつ一ひとつが単一たんいつのポリゴンとして格納かくのうされており、これを用もちいて、それぞれの圃場ほじょうごとのデータを抽出ちゅうしゅつすることが可能かのうとなります。

※なお、この筆ふでポリゴンが公開こうかいされているのは2020年ねん以降いこうのデータであるため、実際じっさいは衛星えいせいデータの撮像さつぞう日びである2006年ねん～2010年ねんと一致いっちしませんので、あくまで参考さんこうとしての利用りようになります。

はじめに筆ふでポリゴンデータをGeoPandasのデータフレームとして読よみ込こみ、衛星えいせい画像がぞうと重畳ちょうじょうするために投影とうえい変換へんかんを行おこないます。

なお、GeoPandasって？という方ほうは以下いかの記事きじをご覧らんください。

## 筆ふでポリゴンの読よみ込こみ
## 藤崎ふじさき町まちの筆ふでポリゴンを抽出ちゅうしゅつ

# 筆ふでポリゴンの保存ほぞん先さきのフォルダを指定してい
# ※自身じしんの環境かんきょうに合あわせて指定してい先さきは変更へんこうしてください
POLYGON_PATH = "/content/drive/MyDrive/work/study/Sorabatake/06_fude_polygon/02青森あおもり県けん（2021公開こうかい）/02361藤崎ふじさき町まち（2021公開こうかい）"
POLYGON_SHP = "02361藤崎ふじさき町まち（2021公開こうかい）__aa.shp"
df_fude_polygon = gpd.read_file(os.path.join(POLYGON_PATH, POLYGON_SHP), encoding="shift-jis")

#データの確認かくにん
df_fude_polygon

Fujisaki_PaddyField = df_fude_polygon[df_fude_polygon["耕地こうちの種類しゅるい"].isin(["田た"])]
Fujisaki_PaddyField = Fujisaki_PaddyField.drop(columns=["公開こうかい年度ねんど","調製ちょうせい年度ねんど"])

#投影とうえい変換へんかん
reprojest_Fujisaki_PaddyField = Fujisaki_PaddyField.to_crs({"init":"epsg:32654"})
reproject_PaddyFiled_out = os.path.join(POLYGON_PATH, "epsg_32654-"+ POLYGON_SHP)

#出力しゅつりょく
reprojest_Fujisaki_PaddyField.to_file(driver="ESRI Shapefile",filename=reproject_PaddyFiled_out, encoding="shift-jis")

#抽出ちゅうしゅつした藤崎ふじさき町まちの筆ふでポリゴンを描画びょうが
f,ax = plt.subplots(1, figsize=(6,6))
ax = reprojest_Fujisaki_PaddyField.plot(ax=ax)
plt.show

同おなじ場所ばしょのデータを取得しゅとくしている場合ばあいには、以下いかのような結果けっかを得えられます。広ひろい範囲はんいで田畑たはたが広ひろがっていることがわかりますね。

こちらで可視かし化かした各かくポリゴン内ないにおいて、衛星えいせい画像がぞうから算出さんしゅつしたNDVI値ちの平均へいきんを時とき系列けいれつに紐ひも付づけしていきます。

# 自身じしんの環境かんきょうに合あわせて指定してい先さきは変更へんこうしてください
MASK_DIR = "/content/drive/MyDrive/work/study/Sorabatake/03_ndvi/02_fujisaki_masked"
ndvi_list = os.listdir(MASK_DIR)
ndvi_list = [f for f in ndvi_list if os.path.isfile(os.path.join(MASK_DIR, f))]
print(ndvi_list)

## NDVI画像がぞうの読よみ込こみ
for ndvi_img in ndvi_list:
  ndvi_path = os.path.join(MASK_DIR, ndvi_img) 
  ndvi_title = ndvi_img[-16:-8]
  print(ndvi_path, ndvi_title)

  #ポリゴン毎ごとのndvi値ちの平均へいきんを格納かくのう
  #画像がぞう毎ごとにリスト初期しょき化か
  ndvi_mean = []

  print(str(ndvi_path) + " : 処理しょり開始かいし")

  with rasterio.open(ndvi_path) as ndvi_src:

    ## 筆ふでポリゴンの読よみ込こみ
    for index, item in reprojest_Fujisaki_PaddyField.iterrows():
      masks = item.geometry

      # 各かく筆ふでポリゴンに沿そってマスク処理しょり 
      out_ndvi_img, out_transform = rasterio.mask.mask(ndvi_src, [masks], nodata=np.nan, crop=False)
      # マスク内ないのNDVI値ちの平均へいきんを算出さんしゅつ -> リストへ
      ndvi_mean.append(np.nanmean(out_ndvi_img))
      #print(np.nanmean(out_ndvi_img))
      #print("------------------")
    
    # NDVI平均へいきん値ちのリストをdfに結合けつごう
    reprojest_Fujisaki_PaddyField[ndvi_title] = ndvi_mean

  print(str(ndvi_path) + "処理しょり終了しゅうりょう")

print("---結果けっかの確認かくにん---")
reprojest_Fujisaki_PaddyField

算出さんしゅつした結果けっかに対たいして、グラフ呼よび出だしを簡単かんたんにするためにカラムの順番じゅんばんを入いれ替かえます。

## NDVI値ちの時とき系列けいれつをソート
recolums = ["id","耕地こうちの種類しゅるい","geometry", "20061031", "20070618", 
           "20080805", "20081105", "20090508", "20100626", "20100926"]
reprojest_Fujisaki_PaddyField = reprojest_Fujisaki_PaddyField.reindex(columns=recolums)
## dfを別名べつめい保存ほぞん
reprojest_Fujisaki_PaddyField.to_file(driver='ESRI Shapefile', 
                                      filename = POLYGON_PATH + "/ndvi_clip.shp", 
                                      encoding="shift-jis")

さて、算出さんしゅつした結果けっかをグラフ化かしてみましょう。

## 追加ついかしたデータフレームからNDVI値ちを取得しゅとくし、グラフ表示ひょうじ
# データの読よみ込こみ
# ※自身じしんの環境かんきょうに合あわせて指定してい先さきは変更へんこうしてください
ndvi_clip_path = "/content/drive/MyDrive/work/study/Sorabatake/06_fude_polygon/02青森あおもり県けん（2021公開こうかい）/02361藤崎ふじさき町まち（2021公開こうかい）/ndvi_clip.shp"
ndvi_clip_df =  gpd.read_file(ndvi_clip_path, encoding="shift-jis")

ndvi_clip_df_t = ndvi_clip_df.drop(columns=["耕地こうちの種類しゅるい",	"geometry"])
ndvi_clip_df_t = ndvi_clip_df_t.T

list_tmp = ndvi_clip_df_t.loc["id"]
ndvi_clip_df_t = ndvi_clip_df_t.drop(index=["id"])
ndvi_clip_df_t.columns = list_tmp

ndvi_clip_df_t.plot.line( 
                         figsize = (10, 6),
                         title = "NDVI time series change",
                         ylim = ([-1.0, 1.0]),
                         xlabel = "Date",
                         ylabel = "NDVI_value",
                         alpha=0.5,
                         legend = False)

結果けっかが表示ひょうじされました。

ただし、このグラフだけではデータ量りょう（対象たいしょうポリゴン数すう：6538個こ）が多おおすぎるため、細こまかいNDVIの変動へんどうが良よくわからず、遊休ゆうきゅう農地のうちを見みつけ出だすことができません。そこで、もう少すこしデータ数すうを絞しぼって、傾向けいこうを見みてみましょう。

# もう少すこしデータ数すうを減へらしてみます
plot_item = []

# rangeの数字すうじを変かえることで表示ひょうじデータ数すうを変更へんこう可か
for i in range(30):
  plot_item.append(i)

ndvi_clip_df_t.plot.line(y = plot_item, 
                         figsize = (10, 6),
                         title = "NDVI time series change",
                         ylim = ([-1.0, 1.0]),
                         xlabel = "Date",
                         ylabel = "NDVI_value",
                         alpha=0.5,
                         legend = False)

変化へんかがわかりやすくなりました。

夏場なつばになると植物しょくぶつが青々あおあおとしていることからも想像そうぞうできると思おもうのですが、農作物のうさくもつは夏場なつばに向むけてNDVIが高たかくなる傾向けいこうがあります。そのような中なかで、夏場なつばであってもNDVI値ちが高たかくならない区画くかくがあることも分わかります。

そこで、NDVIが変化へんかしていない、つまりは耕作こうさくがされていない可能かのう性せいのあるポリゴン（圃場ほじょう）を絞しぼってみます。これらの圃場ほじょうの特徴とくちょうとして、夏場なつばでもNDVIの値ねが低ひくいと仮定かていし、夏場なつばである”20080805”のデータでNDVI値ちが0.0以下いかのポリゴンを抽出ちゅうしゅつしてみましょう。

## 耕作こうさくされていない可能かのう性せいのあるデータを探さがす
## 耕作こうさくがされていない、つまり作付さくづけが行おこなわれていないため夏場なつばでもNDVIの値ねが低ひくいと仮定かてい
#ndvi_clip_df=ndvi_clip_df.drop(columns=["耕地こうちの種類しゅるい",	"geometry"])
ndvi_low_id = ndvi_clip_df[ndvi_clip_df["20080805"]<= 0.0]
ndvi_low_id

このように夏場なつばであってもNDVIの値ねが低ひくい圃場ほじょうがいくつか抽出ちゅうしゅつされました。

実際じっさいに、NDVIが変化へんかしない圃場ほじょうを時とき系列けいれつに確認かくにんしてみると、以下いかの通とおりでした。赤あかく囲かこっている範囲はんいはたしかに目視もくしでみても変化へんかがなさそうな圃場ほじょうに見みえます。

実際じっさいにはこの場所ばしょが実際じっさいに耕作こうさくされていないのか、それ以外いがいの理由りゆうでこのように見みえているのかは定さだかではありませんが、少すくなくとも2006年ねんから2010年ねんまでの間あいだにおいてはNDVIが低ひくい、つまり植生しょくせいが比較的ひかくてき少すくない状態じょうたいであったことがデータから示唆しさされます。

このように衛星えいせい画像がぞう、特とくに赤あか外がい情報じょうほうを含ふくめた複数ふくすうの波長はちょうをうまく組くみ合あわせ、時とき系列けいれつに解析かいせきすることで地上ちじょうの様子ようすをモニタリングし、場合ばあいによっては見みたい対象たいしょう物ぶつ（今回こんかいの場合ばあいは作付さくづけを行おこなっていない田畑たはた）を抽出ちゅうしゅつすることも可能かのうとなります。

コードはGITHUBでも公開こうかいしています。https://github.com/sorabatake/softwaredesign_1