🐈

Sentinel-2の画像から筆ポリゴンを使って圃場ポリゴンの領域の画像を切り取る

石橋　龍

2023/07/21に公開

こんにちは
福岡に来て4年目になりますが夏の暑さはやっぱりこたえる北海道出身の石橋です。

今回は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'

まずは都道府県コードと紐づいててこちらを参照する

その後、市町村コードを調べる（多分 xxx県市町村コードとかで検索すると出てくる）
佐賀の場合はこんな感じで公開されてる

筆ポリゴンのデータ確認

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に変換することもついでにやっておく

# 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日では終わらないのかな？などいろいろ考察ができます。

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

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

Fusic 技術ブログPublication

さまざまな個性を受け入れて有機的につなぐ社内環境を整える。あらゆる事業機会の創出と実現を繰り返し、世の中に対する視点を絶えず増やして成長していく。あっと驚くような角度から発展できるポイントを見つけ、そこにいい感じにフィットする形でテクノロジーを組み込んで、世の中をちょっとずつ、時には大胆にアップデートしつづけていく。

準備

実装

NDVIの計算

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

Discussion