【時系列データ】衛星データの時系列分解を行う【STL分解】

import ee
import matplotlib.pyplot as plt
import pandas as pd
from statsmodels.tsa.seasonal import STL

# GEEの認証（必要に応じてコメントアウトを外す）
# ee.Authenticate()
ee.Initialize()

# 緯度経度、データを扱う衛星の設定
location_lonlat = [140.219017, 35.757631]
stat_region = ee.Geometry.Point(location_lonlat)
sentinel2 = (
    ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
    .filterDate('2019', '2024')
    .filterBounds(stat_region)
    .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
)

# 関心領域を切り取る関数を定義
def reduce_region_function(img):
    stat = img.reduceRegion(
        reducer=ee.Reducer.mean(),
        geometry=stat_region,
        scale=10
    )
    return ee.Feature(stat_region, 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)

# データフレーム成型用関数を定義
def add_date_info(df):
    # millis から Timestamp（日時型）を作成
    df['Timestamp'] = pd.to_datetime(df['millis'], unit='ms')
    
    # 年・月・日カラムを作成
    df['Year'] = df['Timestamp'].dt.year
    df['Month'] = df['Timestamp'].dt.month
    df['Day'] = df['Timestamp'].dt.day

    # YYYY-MM-DD 形式のカラムを追加
    df['Date'] = df['Timestamp'].dt.strftime('%Y-%m-%d')
    return df

# NDVIを計算
ndvi = sentinel2.map(calc_ndvi)
ndvi_stat_fc = (
    ee.FeatureCollection(ndvi.map(reduce_region_function))
    .filter(ee.Filter.notNull(["NDVI"]))
)

# 衛星データを辞書型に変換
ndvi_dict = fc_to_dict(ndvi_stat_fc).getInfo()
ndvi_df = pd.DataFrame(ndvi_dict)

# 日付情報を追加し、Date 列をインデックスに設定
ndvi_df = add_date_info(ndvi_df)
ndvi_df['Date'] = pd.to_datetime(ndvi_df['Date'])
ndvi_df.set_index('Date', inplace=True)

# 月末 (ME) を基準にリサンプリングし、平均を取る
ndvi_monthly = ndvi_df['NDVI'].resample('ME').mean()

# STLによる分解
result = STL(ndvi_monthly.dropna(), period=12, robust=True).fit()
observed = result.observed
trend = result.trend
seasonal = result.seasonal
resid = result.resid

# 可視化
fig, axes = plt.subplots(4, 1, figsize=(8, 8), sharex=True)

# 1) Observed
axes[0].plot(observed.index, observed.values, label='Observed')
axes[0].set_ylabel("NDVI")

# 2) Trend
axes[1].plot(trend.index, trend.values, color='C1', label='Trend')
axes[1].set_ylabel("Trend")

# 3) Seasonal
axes[2].plot(seasonal.index, seasonal.values, color='C2', label='Seasonal')
axes[2].set_ylabel("Season")

# 4) Residual
axes[3].plot(resid.index, resid.values, 'o', color='C3', label='Resid', markersize=3)
axes[3].axhline(y=0, color='black', linewidth=1)
axes[3].set_ylabel("Resid")

# Y 軸ラベル位置の調整
for ax in axes:
    ax.yaxis.set_label_coords(-0.1, 0.5)

plt.tight_layout()
plt.show()

# 緯度経度、データを扱う衛星の設定
location_lonlat = [140.219017, 35.757631]
stat_region = ee.Geometry.Point(location_lonlat)
sentinel2 = (
    ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
    .filterDate('2019', '2024')
    .filterBounds(stat_region)
    .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
)

# 関心領域を切り取る関数を定義
def reduce_region_function(img):
    stat = img.reduceRegion(
        reducer=ee.Reducer.mean(),
        geometry=stat_region,
        scale=10
    )
    return ee.Feature(stat_region, 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)

# データフレーム成型用関数を定義
def add_date_info(df):
    # millis から Timestamp（日時型）を作成
    df['Timestamp'] = pd.to_datetime(df['millis'], unit='ms')
    
    # 年・月・日カラムを作成
    df['Year'] = df['Timestamp'].dt.year
    df['Month'] = df['Timestamp'].dt.month
    df['Day'] = df['Timestamp'].dt.day

    # YYYY-MM-DD 形式のカラムを追加
    df['Date'] = df['Timestamp'].dt.strftime('%Y-%m-%d')
    return df

# NDVIを計算
ndvi = sentinel2.map(calc_ndvi)
ndvi_stat_fc = (
    ee.FeatureCollection(ndvi.map(reduce_region_function))
    .filter(ee.Filter.notNull(["NDVI"]))
)

# 衛星データを辞書型に変換
ndvi_dict = fc_to_dict(ndvi_stat_fc).getInfo()
ndvi_df = pd.DataFrame(ndvi_dict)

# 日付情報を追加し、Date 列をインデックスに設定
ndvi_df = add_date_info(ndvi_df)
ndvi_df['Date'] = pd.to_datetime(ndvi_df['Date'])
ndvi_df.set_index('Date', inplace=True)

# 月末 (ME) を基準にリサンプリングし、平均を取る
ndvi_monthly = ndvi_df['NDVI'].resample('ME').mean()

# STLによる分解
result = STL(ndvi_monthly.dropna(), period=12, robust=True).fit()
observed = result.observed
trend = result.trend
seasonal = result.seasonal
resid = result.resid