植物 の分布 から田畑 の利用 状 況 を推測 する衛星 データの利 活用 方法
Tellusでは、
Tellusでは、
https://www.tellusxdp.com/ja/news/20220630_000502.html
Tellusの
なお、
衛星 データで分 かること・できること
光 の波長 って何 ? なぜ人工 衛星 は人間 の目 に見 えないものが見 えるのか
![](https://sorabatake.jp/wp-content/uploads/2022/07/image13-2.png)
これに
※
![](https://sorabatake.jp/wp-content/uploads/2022/07/image10-1300x449.jpg)
さらに、それぞれ
NDVI=(NIR-Red)/(NIR+Red)
![](https://sorabatake.jp/wp-content/uploads/2022/07/image6-3.png)
衛星 データプラットフォームTellus登録 と準備
TellusのAPIを
https://www.tellusxdp.com/market/login/
Tellusで
まずは
## 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
![](https://sorabatake.jp/wp-content/uploads/2022/07/image2-5-1300x669.png)
#各 個人 のAPIトークンを入力
token = 'XXXXXXXXXXXXXXXX'
衛星 データの検索 ・取得
Tellusの
https://www.tellusxdp.com/docs/travelers/
Tellusで
▼Tellusのデータ
データセットA
└ シーンa
└ ファイルa-i
└ ファイルa-ii
└ シーンb
└ ファイルb-i
└ ファイルb-i
・・・
「データセット」はある
「シーン」は
TellusのCOG形式 のAPIを使 って、ASNARO-2の画像 を重 ね合 わせて前後 の変化 抽出 を実施 してみた
https://www.tellusxdp.com/ja/traveler-service/guide/traveler-how-to/2_traveler.html
![](https://sorabatake.jp/wp-content/uploads/2022/07/image7-min-1-1300x708.png)
dataset_id='ea71ef6e-9569-49fc-be16-ba98d876fb73'
scenes = [
'202ce08d-ba4b-4ffe-8165-109fd3a8b917'
]
# データセット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を取得 する
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()
これらの
# データDL
download(scenes, dataset_id, token)
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_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
そこで、
「
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
"""
2.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_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)
そうすると
![](https://sorabatake.jp/wp-content/uploads/2022/07/image12-2.png)
※
# 青森 県 全域 の行政 区域 データの中 から藤崎 町 だけの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
![](https://sorabatake.jp/wp-content/uploads/2022/07/image5-7.png)
# マスク処理 済 み画像 の保存 先 のフォルダを指定
# ※自身 の環境 に合 わせて指定 先 は変更 してください
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画像 の表示
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)
そうすると
![](https://sorabatake.jp/wp-content/uploads/2022/07/image4-5.png)
耕作 が行 われていない可能 性 のある圃場 を抽出
ここまでで
これらの
しかし、
※なお、この
はじめに
なお、GeoPandasって?という
地理 空間 情報 を扱 うなら知 っておきたいPythonライブラリ、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
![](https://sorabatake.jp/wp-content/uploads/2022/07/image3-2.png)
こちらで
# 自身 の環境 に合 わせて指定 先 は変更 してください
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)
![](https://sorabatake.jp/wp-content/uploads/2022/07/image8-1.png)
ただし、このグラフだけではデータ
# もう少 しデータ数 を減 らしてみます
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)
![](https://sorabatake.jp/wp-content/uploads/2022/07/image9-3.png)
そこで、NDVIが
## 耕作 されていない可能 性 のあるデータを探 す
## 耕作 がされていない、つまり作付 けが行 われていないため夏場 でも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
このように
![](https://sorabatake.jp/wp-content/uploads/2022/07/image1-5.png)
このように
コードはGITHUBでも