こんにちは

この記事はこちらのバックアップです

https://zenn.dev/fusic/articles/6807ccf192eda0

今回はSentinel-2の画像から筆ポリゴンを使って圃場ポリゴンの領域の画像を切り取る方法をメモしておきます

準備

準備するデータは

  • 筆ポリゴン
  • 衛星画像(Sentinel-2)

です。

筆ポリゴンについてはこちら

https://www.maff.go.jp/j/tokei/porigon/

ダウンロードの度にアンケートみたいなのに回答する必要があって面倒なので自分は全データダウンロードしました。

Sentinel-2の画像については適宜open access hubやAPIなどからダウンロードしてください。
今回はLevel-1Cのプロダクトを使用します

ちなみにsentinel-2のプロダクトとしては

  • Level-1C: 建物の倒れこみ等、幾何学的な歪みを補正したオルソ補正済みプロダクト
  • Level-2A: 1Cの補正に加え、大気の吸収・散乱の影響を軽減し、地表面の反射率に変換した大気補正済みプロダクト

自分は今回Notebookで処理を実行しています。

実装

ライブラリのインポート

import folium
import pyproj
import json
import rasterio
import geopandas as gpd
import matplotlib.pyplot as plt 

from pprint import pprint
from pathlib import Path
from rasterio.plot import show
from rasterio.mask import mask
from shapely.geometry import box
from skimage import io

# 佐賀の市町村コード https://www.j-lis.go.jp/spd/code-address/kyuusyuu/cms_16214187.html
json_path = './2022_41/2022_412015.json'

ちなみにファイル名は都道府県コードと紐づいててこちらを参照すると良い

https://makitani.net/shimauma/knowledge/japan-prefecture-code

筆ポリゴンのデータ確認

fude_polygon = gpd.read_file(json_path)
display(fude_polygon.head())
fude_polygon.plot(figsize=(12,12))

ついでに色々みてみる

# crs情報の表示
fude_polygon.crs
fude_polygon.shape

"""
EPSG:4326
(59928, 11)
"""

"""
田畑ごとの割合を確認
100 : 田んぼ
200 : 畑
"""
print(fude_polygon[fude_polygon['land_type'] == 100].shape[0] / fude_polygon.shape[0])
print(fude_polygon[fude_polygon['land_type'] == 200].shape[0] / fude_polygon.shape[0])

"""
0.71687691896943
0.28312308103057
"""

全部表示すると重いので一部だけ切り取って表示

# 筆ポリゴンの一部を可視化。全部は多いので一部だけ可視化
target_area_polygon = fude_polygon[
    (fude_polygon['point_lat'] > 33.16374436898851) & (fude_polygon['point_lat'] < 33.19595709521625) & \
    (fude_polygon['point_lng'] > 130.24408154375982)     & (fude_polygon['point_lng'] < 130.4628867046025)             
]

print(target_area_polygon.shape[0])

latlon = [33.16068647538113, 130.2977306972179]
m = folium.Map(latlon,zoom_start=15, control_scale=True)
folium.GeoJson(target_area_polygon).add_to(m)
m

NDVIの計算

# NDVIが未計算の時に実行する

sentinel_2_dir_list = [
    '../sentinel-2/S2A_MSIL2A_20230307T015651_N0509_R060_T52SFB_20230307T051354.SAFE',
    '../sentinel-2/S2B_MSIL2A_20230322T015659_N0509_R060_T52SFB_20230322T043851.SAFE',
    '../sentinel-2/S2B_MSIL2A_20230501T015659_N0509_R060_T52SFB_20230501T041707.SAFE',
    '../sentinel-2/S2A_MSIL2A_20230516T015651_N0509_R060_T52SFB_20230516T055754.SAFE',
    '../sentinel-2/S2A_MSIL2A_20230615T015701_N0509_R060_T52SFB_20230615T062501.SAFE'
]

def calc_ndvi_sentinel_2():
    for dir_path in sentinel_2_dir_list:
        # load img and calc ndvi
        b4_path = list(Path(dir_path).glob('**/*B04_10m.jp2'))[0]
        b8_path = list(Path(dir_path).glob('**/*B08_10m.jp2'))[0]
        b4 = rasterio.open(b4_path)
        b8 = rasterio.open(b8_path)
        red = b4.read()
        nir = b8.read()
        ndvi = (nir.astype(float)-red.astype(float))/(nir+red)

        # set meta
        out_meta = b4.meta
        out_meta.update(driver='GTiff')
        out_meta.update(dtype=rasterio.float32)

        # NDVI out path
        year = dir_path.split('_')[2][0:4]
        month_day = dir_path.split('_')[2][4:8]
        ndvi_out_path = f'./sentinel2_ndvi/{year}_{month_day}.tiff'

        with rasterio.open(ndvi_out_path, 'w', **out_meta) as dst:
            dst.write(ndvi.astype(rasterio.float32))
            
# calc_ndvi_sentinel_2()

ポリゴンを用いて切り取り、可視化

関数の用意

"""
ポリゴンの座標系をSentinel-2のデータと同じにする必要がある。
UTM座標系は32651〜32656まである。観測する地域に合わせて適宜調整する

https://tmizu23.hatenablog.com/entry/20091215/1260868350

EPSGの場合 326[ZoneID]となる
福岡、佐賀などはZone52なので32652
東京などだったらZone54なので32654

ちなみにEPSG 6688 ~ 6692も使えるけど、ZoneIDと一致しないので覚えづらいかも...
"""

def get_masked_image_for_sentinel_2(image_href, poly):
    data = rasterio.open(str(image_href))
    masked, mask_transform = mask(
		dataset=data, 
		shapes=poly.geometry.to_crs(crs={'init': 'epsg:32652'}), 
		crop=True
	)

    return masked[0,:,:]
    
def display_cropped_img(img_dir, polygon):
    for img_path in sorted(Path(img_dir).glob('**/*')):
        plt.figure(figsize=(5, 5))
        masked_img = get_masked_image_for_sentinel_2(str(img_path), polygon)
        plt.imshow(masked_img, cmap='bwr')
        plt.title(str(img_path).split('/')[1].split('.')[0])
## 佐賀周辺のポリゴンを手動で作成してNDVIの時系列を表示
bbox = box(minx=130.21971622772818, miny=33.071356339303456, maxx=130.50652945686997, maxy=33.30119764666647)
geo = gpd.GeoDataFrame({'geometry': bbox}, index=[0], crs='EPSG:4326')

display_cropped_img('./sentinel2_ndvi/', geo)

sentinel-2の地上分解脳は10mなので小さな圃場ではサイズが足りないため、比較的大きな圃場ポリゴンで実験してみる

# 佐賀周辺の筆ポリゴンでNDVIの時系列を表示
# サイズの大きいところだけ抽出する。面積に変換するためにepsg6689に変換する

fude_polygon_copy = fude_polygon.copy()
fude_polygon_copy.to_crs('epsg:32652',inplace=True)
fude_polygon_copy['area'] = fude_polygon_copy['geometry'].area

fude_polygon_copy_sorted_desc = fude_polygon_copy.sort_values(by='area', ascending=False)
fude_polygon_copy_sorted_desc
i = 24
sample = fude_polygon_copy_sorted_desc[i:i+1]
print(sample['area'])
display_cropped_img('./sentinel2_ndvi/', sample[["geometry"]])

圃場内でのNDVIの推移が直感的に見えますね
収穫作業が1日では終わらないのかな?などいろいろ考察ができます。

とまあ、こんな感じで筆ポリゴンを用いた衛星画像解析の基本ができました

あとはポリゴンの統計量(平均値とか)を計算して時系列データとして扱うとか、このあとの料理方法もいろいろありそうですね!