【Python】シェープファイルを利用して衛星データを収集する【GIS】

前回の記事「【Python】シェープファイルデータを地図上にプロットする方法【GIS】」でシェープファイルをデータフレーム化し、マップ上にプロットする方法を紹介しました。

今回は作成したデータフレームから緯度・経度情報を抽出し、その地点の衛星データを収集する方法を紹介します。シェープファイルのデータフレーム化については以下の記事をご確認ください。
https://labo-code.com/python/shapefile/

動作検証済み環境

macOS Monterey(12.6), python3.10.7

Google Earth Engineを用いた衛星データ解析入門 ~タイムラプスを作成し地球の環境変動を観察~

Pythonを使った衛星データ解析の初心者向け技術書を販売中！

Google Earth Engineを使えば簡単に
自宅にいながら地球環境の変動を観察できます

技術書ページへ

シェープファイルから衛星データを収集してみる

シェープファイルから緯度経度を抽出する

前回の記事でも紹介しましたが、シェープファイルは図形情報と属性情報を持っています。この属性情報には座標系が含まれています。この座標系が緯度・経度を示す地理座標系の場合、そのまま衛星データの収集に利用することができます。

それでは、pythonでシェープファイルから衛星データを収集するプログラムをコーディングしてみます。以下に実装例を示します。

利用するデータ

前回の記事で作成したアカマツの植生データが入ったcsvファイルvegetation.csvを利用します。csvファイルは以下のディレクトリに格納します。

Desktop
|_ labcode
   |_ vegetation.csv

ソースコード

下記のコードをgis_satellite.pyという名前でDesktop/labcode/ディレクトリに保存します。

from functools import reduce

import ee
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

# GEEの認証
# ee.Authenticate()
ee.Initialize()

# 関心領域を切り取る関数を定義
def reduce_region_function(img):
  stat = img.reduceRegion(
    reducer=ee.Reducer.mean(),
    geometry=geo_point,
    scale=10
  )
  return ee.Feature(geo_point, stat).set({'millis': img.date().millis()})

# GEEデータを辞書型に変換する関数を定義
def fc_to_dict(fc):
  prop_names = fc.first().propertyNames()
  prop_lists = fc.reduceColumns(
    reducer=ee.Reducer.toList().repeat(prop_names.size()),
    selectors=prop_names).get('list')
  return ee.Dictionary.fromLists(prop_names, prop_lists)

# NDVI計算用関数を定義
def calc_ndvi(image):
  ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI')
  return image.addBands(ndvi)

# Timestamp作成用関数を定義
def add_date_info(df):
  df['Timestamp'] = pd.to_datetime(df['millis'], unit='ms')
  df['Year'] = pd.DatetimeIndex(df['Timestamp']).year
  df['Month'] = pd.DatetimeIndex(df['Timestamp']).month
  return df

# csvファイルの読み込み
gdf_tree = pd.read_csv('vegetation.csv')

# csvファイルをデータフレーム化
gdf_tree = gpd.GeoDataFrame(
  gdf_tree.loc[:, [c for c in gdf_tree.columns if c != "geometry"]],
  geometry=gpd.GeoSeries.from_wkt(gdf_tree["geometry"]),
  crs="epsg:4326" # 座標参照系の設定
)

# 緯度経度をx,y座標に変換
gdf_tree_geo = gdf_tree['geometry'].apply(lambda p:list(p.exterior.coords.xy))

# 辞書型を定義
ndvi_df = {}

# 衛星データのデータフレームを作成
for i in range(0,1639):
  x,y = gdf_tree_geo[i]
  cords = np.dstack((x,y)).tolist()
  double_list = reduce(lambda x,y: x+y, cords)
  single_list = reduce(lambda x,y: x+y, double_list)
  geo_point = ee.Geometry.Polygon(single_list)
  sentinel2 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED').filterBounds(geo_point).filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20))
  ndvi = sentinel2.map(calc_ndvi)
  ndvi_stat_fc = ee.FeatureCollection(ndvi.map(reduce_region_function)).filter(ee.Filter.notNull(ndvi.first().bandNames()))
  ndvi_dict = fc_to_dict(ndvi_stat_fc).getInfo()
  ndvi_df[i] = pd.DataFrame(ndvi_dict)

# 月ごとの中央値を計算
vg_df = pd.concat(ndvi_df)
vg_df['timestamp'] = pd.to_datetime(vg_df['millis'], unit='ms')
vg_df.set_index('timestamp', inplace=True)
vg_df.to_csv('ndvi_vg.csv')
vg_df = vg_df.resample('M').median()
vg_df = add_date_info(vg_df)

# カラムの並び替え
vg_df = vg_df.reset_index()

# 反射率のグラフを作成
vg_plot = vg_df.reindex(columns=['B1','B2','B3','B4','B5','B6','B7','B8','B8A','B9','B11','B12'])
vg_plot = vg_plot[0:1]
y = vg_plot.T.unstack().reset_index(level=0, drop=True)
x = np.array([443,490,560,665,705,740,783,842,865,945,1610,2190])
ax = plt.subplot()
ax.scatter(x,y)
ax.plot(x,y)
ax.set_xlabel("Wavelength [nm]")
ax.set_ylabel("BoA reflectance × 10,000")
ax.set_xlim(400, 2500)
plt.savefig('reflectance.png')

プログラムを実行する

ターミナルを開き、

$ cd Desktop/labcode/

と入力し、ディレクトリを移動します。あとは以下のコマンドを実行するだけです。（ $マークは無視してください）

$ python gis_satellite.py

実行結果

~/Desktop/labcode にndvi_vg.csvとreflectance.pngというファイルができていれば成功です！
ndvi_vg.csvには前回作成したcsvファイルの地点を測定した衛星データの値が格納されています。

また、作成されたpngファイルは以下のようになっています。

このコードでは2018年12月の衛星データを取得しています。アカマツは常緑の針葉樹林で、冬場も緑の葉を持っています。

今回作成した図は縦軸が反射率、横軸が波長率となっており、取得した衛星データの反射率をバンド毎にプロットしています。Sentinel-2では赤色光がB4、近赤外線がB8となっており、この図の左から4番目のプロットと8番目のプロットが対応しています。植物は赤色光を光合成に利用し、近赤外線を反射する特徴が見事に表れています。

反対に1610 nmのB11や2190 nmのB12といった短波長赤外の波長域は水の吸収が起きるため反射率の低下がみられます。

コードの解説

上に書いたソースコードで重要なところを解説をしていきます。

基本的には以前の記事（【Python】衛星リモートセンシングで時系列変化を捉える【NDVI】）のコードがベースになってるので、そちらも併せて参照してください。

from functools import reduce

import ee
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

GEE、geopandas、matlibplot、numpy、pandasのライブラリを使用します。ライブラリが存在しない場合はpip installでインストールしましょう。

# csvファイルの読み込み
gdf_tree = pd.read_csv('vegetation.csv')

# csvファイルをデータフレーム化
gdf_tree = gpd.GeoDataFrame(
  gdf_tree.loc[:, [c for c in gdf_tree.columns if c != "geometry"]],
  geometry=gpd.GeoSeries.from_wkt(gdf_tree["geometry"]),
  crs="epsg:4326" # 座標参照系の設定
)

# 緯度経度をx,y座標に変換
gdf_tree_geo = gdf_tree['geometry'].apply(lambda p:list(p.exterior.coords.xy))

csvファイルをデータフレーム化します。6行目と7行目でデータフレームのgeometryのデータ型をgeometryに変換しています。8行目ではgeometryの座標参照系を設定しています。epsg:4326はESPGコードの一つでEuropean Petroleum Survey Groupという団体が定めたGISの座標系に関連する要素を定義したものとなっています。4326はWGS84という測地系と緯度経度座標系が定義されています。

これらを設定することで12行目のexterior.coordsを用いたxy座標への変換が可能となります。

# 衛星データのデータフレームを作成
for i in range(0,1639):
  x,y = gdf_tree_geo[i]
  cords = np.dstack((x,y)).tolist()
  double_list = reduce(lambda x,y: x+y, cords)
  single_list = reduce(lambda x,y: x+y, double_list)
  geo_point = ee.Geometry.Polygon(single_list)
  sentinel2 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED').filterBounds(geo_point).filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20))
  ndvi = sentinel2.map(calc_ndvi)
  ndvi_stat_fc = ee.FeatureCollection(ndvi.map(reduce_region_function)).filter(ee.Filter.notNull(ndvi.first().bandNames()))
  ndvi_dict = fc_to_dict(ndvi_stat_fc).getInfo()
  ndvi_df[i] = pd.DataFrame(ndvi_dict)

3~7行目でxy座標に変換したgeometryをGEEで利用できるPolygonの形に変換しています。

8~12行目では、以前の記事と同じNDVIの算出から衛星データのデータフレーム化までを行っています。

# 反射率のグラフを作成
vg_plot = vg_df.reindex(columns=['B1','B2','B3','B4','B5','B6','B7','B8','B8A','B9','B11','B12'])
vg_plot = vg_plot[0:1]
y = vg_plot.T.unstack().reset_index(level=0, drop=True)
x = np.array([443,490,560,665,705,740,783,842,865,945,1610,2190])
ax = plt.subplot()
ax.scatter(x,y)
ax.plot(x,y)
ax.set_xlabel("Wavelength [nm]")
ax.set_ylabel("BoA reflectance × 10,000")
ax.set_xlim(400, 2500)
plt.savefig('reflectance.png')

1行目でデータフレームvg_dfから測定バンドのカラムだけを抽出し、以降の行の処理でそのデータを散布図にしています。今回は2行目でグラフ化する対象を0行目のデータに絞っています。グラフ化したいデータに応じて、この行を修正するといいでしょう。

最後にグラフをpngとして出力しています。